Unsupervised Training of Convex Regularizers using Maximum Likelihood Estimation
Authors: Hong Ye Tan
Ziruo Cai
Marcelo Pereyra
Subhadip Mukherjee
Junqi Tang
Carola-Bibiane Schönlieb
Venue: Accepted by TMLR ()

Abstract: Imaging is a canonical inverse problem, where the task of reconstructing a ground truth from a noisy measurement is typically ill-conditioned or ill-posed. Recent state-of-the-art approaches for imaging use deep learning, spearheaded by unrolled and end-to-end models and trained on various image datasets. However, such methods typically require the availability of ground truth data, which may be unavailable or expensive, leading to a fundamental barrier that can not be addressed by choice of architecture. Unsupervised learning presents a powerful alternative paradigm that bypasses this requirement by allowing to learn directly from noisy measurement data without the need for any ground truth. A principled statistical approach to unsupervised learning is to maximize the marginal likelihood of the model parameters with respect to the given noisy measurements. This paper proposes an unsupervised learning approach that leverages maximum marginal likelihood estimation and stochastic approximation computation in order to train a convex neural network-based image regularization term directly on noisy measurements, improving upon previous work in both model expressiveness and dataset size. Experiments demonstrate that the proposed method produces image priors that are comparable in performance to the analogous supervised models for various image corruption operators, maintaining significantly better generalization properties when compared to end-to-end methods. Moreover, we provide a detailed theoretical analysis of the convergence properties of our proposed algorithm.

Published in Transactions on Machine Learning Research (12/2024)
Unsupervised Training of Convex Regularizers using Maximum Likelihood Estimation
Hong Ye Tan
hyt35@cam.ac..u
Deportment of Applied Mathematics and Theoretical Physics
University of Cambridge, United Kingdom
Ziruo Cai
_ cwiziruo@@sjtu.cco edu.ccm
School of Mathematical Sciences Shanghai Jiao Tong University, China
Marcelo Pereyra
m.pereyra@huu.cc.om
School of Mathematical & Computer Sciences Heriot- Watt University, Edinburgh
Subhadip Mukherjee
mmukherjee@ccc..ithgg ac.in
Deportment of Electronics and Electrical Communication Engineering Indian Institute of Technology, Kharogpur, India
Junqi Tang
j.tang.20hhaam 20hham..ccmmm ac. uk
School of Mathematics University of Birmingham, United Kingdom
Cabila-Bibiane Schönlieb
cbs31 @cam.ac. uk
Deportment of Ayplied Mathematics ond Theoretical Physics University of Combridge, United Kingdom
Reviewed on OpenReview: https: openreview. net/ forum? id=375tU77h05
Abstract
Imaging canonical inverse problem, where the task of reconstructing ground truth from noisy measurement is typically ill-conditioned or ill-posed. Recent state-of-the- deep learning, spearheaded by unrolled and end-to-end
art approaches for imaging models and trained on various image datasets. However, such methods typically require the availability of ground truth data, which may be unavailable or expensive, leading to a fundamental barrier that can not be addressed by choice of architecture. Unsupervised learning presents a powerful alternative paradigm that bypasses this requirement by allowing to learn directly from noisy mcasurement data without the need for any ground truth. A principled statistical approach to unsupervised learning is to maximize the marginal likelihood of the model parameters with respect to the given noisy measurements. This paper proposes unsupervised learning approach that leverages maximum marginal likelihood estimation and stochastic approximation computation in order to train convex neural network-based image regularization term directly on noisy measurements, improving upon previous work in both model expressiveness and dataset size. Experiments demonstrate that the proposed method produces image priors that are comparable performance the analogous supervised models for various image corruption operators, maintaining significantly better generalization propertics when compared to end-to-end methods. Morcover, provide detailed theoretical analysis of the convergence properties of our proposed algorithm.

---

Published in Transactions on Machine Learning Research (12/2024)
1 Introduction
Image reconstruction often requires solving high-dinenniooal inverse problems, ranging from image denoising and deconvolution, to phase retrieval or computed tomography. These problems are usually formulated as the recovery of a signal x E X from a noisy measurement y = Ax + w E Y, with forward operator A : -+ y and some noise w. A common trait of these problems is ill- ill-posedness, where the solution of Ax = y can have wero or infinitely many solutions, necessitating the use of regularized reconst ion operators (Benning & Burger 2018 Engl et al. 1996). A classical method of solving this is using variational regularization, combining a data fidelity term with priors such as wavelet priors total variation priors (Chambolle & Pock 2016 Arridge et al.) 2019 Recently, various end-to-end approaches utilizing neural networks have pushed the state-of-the-art for imaging, achieving highly detailed reconstructions at the cost of well-crafted training data and computational resources. These include direct neural network approaches, as well as iterative plug-and-play priors (Zhang et al., 2021 Venkatalarishnan et al. 2013).
Variational models. One way to define the reconstruction of a measurement y as the solution of a minimization problem arg ming fa (Ax) + (xx (Arridge et al., 2019: Chambolle & Pock, 2016). Here, fy:Y + R is data fidelity function measuring how well the reconstruction fits the measurement y E V, and Rx : X + RRR is regularization functional that promotes regularity of the reconstructions. Examples include computed tomography, where A is the Radon transform and RA may enforce some type of wavelet sparsity (Antoniadis & Fan 2001 Donoho 2006). In the particular case of Gaussian denoising, A=1, I, the fidelity can be taken to be fz(z) = | y 7/2 and regularization Rx(x) = V/| commonly known as the total variation (TV) regularized solution. Gaussian noise allows for another probabilistic interpretation using Tweedic's identity, which provides a relationship between the the posterior mean estimator and a ground-truth density p(x)): it states that if pe is given by convolving p with Gaussian kernel of variance 6, then the minimum mean- -quared estimator (MMSE) denoiser Dr satisfies (D: - ())x) = eV log (x)) (Efron 2011 Laumont et al. 2022). That is, the denoising step is in the direction of the posterior score, preciscly given by gradient descent. on a convolved version of (x) + g(x). This begins to endow learnable variational
models with probabilistic interpretation.
Learned variational models are a middle ground between fully classical model-based reconstructions and fully learned end-to-end models (Lunz et al., 2018 Mukherjee et al. 2020 Kobler et al., 2020). Heuristically, the fidelity term (Ax) should depend on the observation process involving forward operator A and the noise model, while the regularization Rx ignores the observation process and considers only the data. One advantage of the decoupling is that the learned regularizer can be applied to different forward models without retraining. This is exploited in ald--and methods, which implicitly use generic denoisers as regularizers for a varicty of non-denoising tasks such as demosaicing and super-resolution (Zhang et al. 2021 Hurault et al. 2022a). Moreover, the variational structure allows for the application of classical regularization theory under regularity assumptions on RA (Burger & Osher 2004 Hauptmann et al. 2024 Schwab et al. 2019)
Learned regularization. Notable examples of learned regularization include adversarial regularization and
network Tikhonov (Mukherjee et al. 2021 Lunz et. al. 2018 Li et al.) 2020a Alberti et al.
2021 In the
supervised regime, L.e presence of clean and noisy training pairs, one method is to minimize some distance metric between the clean image and the reconstructed noisy image as given by an approximate solution of the variational problem (Diamond et al. 2017 Additional regularization is often added to enforce regularity on the learned regularizer, such as having small Lipschitz constant, which is linked to Wasserstein robustness (Kuhn et al.] 2019). In the regime of weak supervision, given unpaired clean and noisy data, adversarial regularization aims to penalize noisy images and promote clean images by maximizing and minimizing the regularizing terms respectively (unz et al. 2018 Mukherjee et al. 2020). Using the structure of a fidelity and regularizer allows for more explicit description of the final reconstruction, and can improve the performance of a learned denoiser when used in an unrolled scheme (Zhang el. al., 2021). We note that it is possible to indirectly learn R2 by instead learning mappings such as the proximal proxn, or gradient VR,, while still maintaining well-posodness of the variational formulation through convex analysis identities (Gribonval & Nikolova 2020 Morcau 1965). For example, under certain Lipschitz conditions, denoisers can be interpreted as proximal operators of weakly convex functions, allowing for another degree of regularity (Hurault et al.
2

---

Published in Transactions on Machine Learning Research (12/2024)
2022b
Shumaylov et al. 2023 This denoiser-prior connection has been recently extended to diffusion models Graikos et al. 2022].
Feng et al. 2023
Bayesian approaches. Approaches rooted in the Bayesian statistical paradigm usually aim to characterize
the posterior distribution (density) of x conditioned on observing the noisy measurement y (Stuart 2010
Calvetti & Somersalo, 2018). Using Bayes' rule, this density can be expressed as p(arly) = f(yl=)(()) Myy,, where 2(y[z) = Pw (v, Ax) ) is the likelihood arising from the noise distribution of w, p(x) is the prior density on 2, and Z(y) is a normalizing constant. The modelled distribution is thus characterized by choice of family for the prior p(z). Computational Bayesian inference has been active area of development for the past 40 years, including hicrarchical image models (Geman & Geman 1986; Carlin et al. 1992), diffusions (Roberts & Stramer 2001), and dynamical systems (Friston 2002; Wilkinson 2018), and we refer to Green et al. (0015); Mukherjee et al. (0023) for more comprehensive overviews. A widely-used modern method to describe the
distrriutiin p y) is to draw samples from it using Markov Chain Monte Carlo (MCMC) methods (Gilks 2007). One distinct advantage of the Bayesian approach is the option of uncertainty quantification rather than point estimation, given sufficient computational resources (Carioni
ct al. 1995 Robert et al.)
1999;
et al.
2023
Nagel & Sudret 2016 Durmus et al. 2018 Holden et al. 2022). This is used in various imaging
methods such as Wasserstein generative adversarial networks [Arjovsky et al. 2017), and Cycle-GANs ZZhu
ct al. 2017). , which utilize optimal transport to flow from noisy to clean images with a learned prior. Since
Bayesian approaches use the entire distribution, another possible advantage over other training methods is that paired training examples (such as clean-noisy pairs) are often not needed, requiring only access to samples from the marginal distributions and not the joint distribution. This can be useful in tasks such as unsupervised image registration, where paired data is labor-intensive (Balakrishnan et al. 2018 Chen et al., 2022b
The variational formulation can be interpreted as a special case of the Bayesian formulation. Given the negative prior log-density g(x)= = - log p(z) and negative log- likelihood fu(x) = log ((ylz), the maximum a- a-posteriori estimator (MAP) #MAP = arg max, p(xly) xx fx(( - g(x)) is equal to the minimizer min, fy(x) + g(x). In Bayes' rule, the likelihood can thus be interpreted as a fidelity term, and the prior as regularization. Pereyra (019) further explores this connection, demoustrating that for log-concave distributions, the MAP estimator is optimal with respect to Bregman divergence in terms of expectation, and is the dual problem of
MMSE problem. This connects variational regularization with the statistical problem of MAP estimation.
Unsupervised learning. Many end-to-end approaches as well as approaches described above require clean ground truths I and measurements y, which can be limiting in cases where ground truth data is unavailable or expensive, such as in medical imaging. We are instead concerned with the situation where we have a finite amount of noisy data, such as if a clean dataset corrupted. One such approach is the deep image prior (DIP), in which the neural network is fed a random input, and subsequently has its parameters optimized to recover noisy mcasurement with early stopping (Ulyanov et al. 2018). This has shown remarkable performance without any supervised data, using only the implicit regularization of the neural network architecture. Theoretical results for DIP include optimization dynamics analysis (Heckel & Soltanolkotabi 2020; Van Veen et al., 2018) and Bayesian interpretations using Gaussian process limits [Cheng et al. 2019). Pajot et al. (2018) propose using the discriminator of a GAN as a regularizer, using Stein's unbiased risk estimator to train using only corrupted measurements (Metzler et al. 2018 Stein 1981]. ).
Noise2Inverse considers the specific case where one noisy measurement can be considered as the average of multiple noisier measurements to artificially generate more data, such in limited-angle tomography (Hendriksen et al. 2020). We note that some unsupervised approaches such as Noise2Noise generate an unlimited number of noisy measurements, where averaging out the empirical risks with the noisy targets leads to an approximation of the empirical risks when training with clean targets (Lehtinen et al. 2018 but we do not consider this simpler setting.
Equivariant imaging (EI) is another unsupervised framework that learns reconstruction functions from ouly compressed measurements (Chen et al. 2021 Pereyra & Tachella| 2024 Scanvie et al. 2023 The framework utilizes the fact that the set of plausible signals in X' are invariant to certain group of transformations g. By "rotating" the range of the adjoint operator, constructing virtual operators corresponding to invariances of the forward operator, the learned operator is able to learn reconstructions that are approximately equivariant
3

---

Published in Transactions on Machine Learning Research (12/2024)
to the transformations. This greatly mitigates problems arising from nontrivial millspaces and removes artifacts that are not invariant under transformations in g.
We briefly remark on unsupervised versions of diffusson models, which is an area of intense research efforts. Diffusion models trained in a supervised manner have been successfully deployed in a plug-and- play manner with a data guidance term, delivering state-of-the-art reconstructions across a variety of challenging imaging problems (Chung et al., 2023; [Batzolis et al. 2021 Mardani al. 2024). However, from Bayesian computation viewpoint, these "guided" diffusion models are highly inaccurate (e.g., they significantly underestimate uncertainty quantification results, see (Pereyra & Tachella, 2024
Thong et al,, 2024)). ) Alternatively, for some problems it is also possible to employ Stein's unbiased risk estimator to train diffusion models directly from measurement data, without ground truth (Kawar et al. 2024 Aali et al. 2023) However, in problems that are ill-posed, issues arise with the null-space of the forward operator, leading to improper training (Chen et al.] 2022a). Other approaches consider transferring powerful pre-trained diffusion models from one setting to another, such as change of noise, forward operator, or datasct (Wang al.) 2023 Zhu al. 2023 As discussed previously, we are instead interested in the case where do not. have access to multiple copies of corrupted measurements nor pre-trained model. In this sense, our proposed method is thus able to train priors in strictly harder scenarios.
Statistical estimation of parameters. The problem of learning model parameters sits naturally within the statistical estimation framework and can be tackled by using powerful statistical inference techniques, such as maximum likelihood estimation. In many cases, such techniques do not require ground truth data and they are known to outperform common heuristics such as the L-curve criterion (Hansen & O'Leary 1993) and Morozov's discrepancy principle (Morozov] 1966) (e-g., see Vidal et al.] (2020) for detailed comparison of maximum likelihood estimation and alternative strategies in the context of regularization parameters).
In addition, Bayesian statistical approaches allow for statistical interpretations of reconstructions and their associated parameters. Alberti et al. (2021) provides a statistically motivated method of learning a Tikhonov preconditioner based on covariance estimators. Expectation maximization allows for finding locally optimal parameters when there are unobserved latent variables in an iterative manner (Dempster al.] 1977 Robert et al. 1999), which can be applied to tasks such as image segmentation and tomography (Carson et 2002 Levitan & Herman 1987). Vidal et al. (2020) proposes stochastic approximation method of estimating the optimal regularisation parameters for TV regularized denoising or hyperspect ral unmixing for a single image using maximum marginal likelihood estimation, which can then be used for post-hoc reconstruction. The method of Vidal et al. (2020) was recently extended by Mbakam et al. (2024) to include the estimation of unknown parameters in the forward model for semi-blinn image deconvolution problems, and to calibrate deep generative priors within a Bayesian image estimation framework Melidonis et al.] (2024). This paper builds on this develop a stochastic approximation scheme train a convex neural network regularization terms in fully unsupervised manner.
1.1 Contributions
In this work, we propose a principled unsupervised method of training convex regularizers with the following objectives.
1. (Unsupervised.) Training a regularizing prior using only noisy measurements. This is achieved using maximum marginal likelihood estimation. In particular, we do not assume the existence of multiple noisy copies of the data, such as in Noise2Noise schemes (Lehtinen et al. 2018
2. (Eapressive..) We increase the dimensionality of the regularizer from the order of 10 to 10 Existing MLE-based methods consider reconstruction for single images by tuning scalar parameters for hand- crafted priors. By using a more expressive parameterization for the regularizer, we can achieve better performance, and demonstrate that the proposed MCMC-based method is applicable in high dimensions.
3. (Efficient.) Unsupervised methods have strictly less available information than supervised methods. By comparing the proposed unsupervised method with the corresponding supervised method with the
4

---

Published in Transactions on Machine Learning Research (12/2024)
same regularizer architecture, we demonstrate that this loss of information is almost mitigated. We observe that the regularizers trained in the proposed unsupervised manner have strong regularixing effect, and have reasonable performance gaps compared to the same model under supervised training.
4. (Convergent training.) Theoretical ergodic convergence of the regularizers under the proposed method. This is done using a particular choice of regularizer architecture, which is used to verify the assumptions of convergence results in previous works (De Bortoli et al. 2021 2020
5. (Convergent estimation.) A convex architecture also implies nice properties of the statistical estimation. In the context of sampling, a convex regularizer leads to a log-concave density and provides existence of all posterior moments, which implies that the estimation of desired statistical quantities well-posed.
This work is organized follows. Section 2 casts the variational regularization formulation into the Bayesian framework, and motivates estimating regularization parameters using maximum marginal likelihood estimation. This problem is then decomposed and solved using a stochastic gradient scheme. The proposed method is then given in Algorithm [ with its mini-batched version Algorithm 2 for use in larger datasets. This extends the work of Vidal et al. (2020), from training regularization parameters on single images, to training convex ridge regularizers on entire datasets in an unsupervised manner. Section 3 presents convergence theory for the proposed method using the convex ridge regularizer architecture. In particular, Theorem 2 states that the proposed stochastic gradient scheme has ergodie convergence. Section evaluates the performance of our proposed method against various supervised and unsupervised baseline methods.
||
2 Unsupervised Learning by Maximum Marginal Likelihood Estimation
In this section, we first reformulate the variational regularization framework a maximum a-posteriori problem. From this Bayesian interpretation, we consider finding the regularizer parameters using maximum marginal likelihood estimation, which uses only noisy measurements. Motivated by applying gradient ascent on the marginal likelihood, we then introduce the decomposition of the marginal likelihood into posterior and prior expectations as given in Vidal al. (2020).
Given a fidelity term fy and regularization term go with trained parameters 0, the standard variational reconstruction for given measurement is given follows:
c arg min (f)()) + go()))
We can reformulate this maximum a-posteriori problem mentioned in the introduction. Setting fy(z) log ((ylx) as the negative log- liglliill and go(z) - p(z(() the variational reconstruction is equivalent to
TMAP = arg max log [[)]==)(=((0] arg max pogy(,, 0),
(1)
22
where we use Bayes' rule: p(zly,0) 0) x ((v)x)p(z((). To apply this relation for reconstruction, we still need to learn 0. In the unsupervised paradigm, samples from the joint distribution p(x,y) y) and even the prior distribution p(z) may be unavailable, which leaves us with learning an estimator 0* directly from noisy mcasurements sampled from p(y).
Adopting empirical Bayesian approach, the regularization parameter 0 can be accurately estimated from noisy measurements y by maximum marginal likelihood estimation (Vidal et al. 2020
0t arg max log p(yl0) BEE
(2)
Here, p(u(0) is the marginal likelihood given by marginalizing over data p(vl0) [()])))(z(00ddr and p(x|0) & (x)(-go(x)) is the density of the prior Gibbs measure. Standard assumptions include the negative logeliielioood (x) = - log Eyyl=) being convex and C³ with TypTiischitz gradient (in x), as well as the admissible parameter set e being compact.
One method of maximizing the marginal likelihood is by gradient ascent on 8 which requires knowledge of Vo log p(y)0). However, this is intractable and requires approximation. We now detail how to approximately
5

---

Published in Transactions on Machine Learning Research (12/2024)
compute this gradient using only noisy measurements y using double MCMC approach Vidal et al. 2020), extending the aforementioned work from single measurements to datasets. The stochastic gradient approximator can then be used inside a gradient-step optimizer, referred to in the literature as the stochastic optimization with unadjusted Langevin (SOUL) algorithm (De Bortoli et al.
2021
2020)
We first decompose the marginal likelihood p(yl0) as follows (Vidal et al. 2020 Denote the normalizing
constant of the prior Gibbs measure p=z|0) (ep(-go(x)) as
7(0) exp(( dz.
3
Then p(x[0) - exp(-g((=)/ZZ() the prior density. While the marginal likelihood + p(y is intractable, under cerlain regularity assumptions] detailed in Proposition | and Vidal et al. (2020) Prop. A1), , it can be replaced with noisy estimate and decomposed in the following manner using Fisher's identity (Vidal et al, 2020
Douc et al. 2014):
/ ploly,, V) logp(r, ul0) dr
Ve log (og((y)0) =
= =- V (ogo(c))((((1)) 0) dr Vo log(((0))-
(4)
Noticing that VoZ(0) = Vogo exp( go (2)) de, write
v log Z(0) = (v,Z(0)) = Vogs(z))(al))
(5)
Putting
4 and
together, we have the following integral expression for the marginal likelihood:
V 1og(())0) Eo ( Exy V [ogo
Using this decomposition, allowing perform noisy gradient ascent 0 using only noisy observations y (and access to the likelihood f). In this work, we extend the setup of Vidal et al. 2020 in the following manner:
able to use MCMC methods to estimate the prior and normalizing constant,
Experimentally, the prior distributions are loosened from being of the form p(x(0) x (xp(-OTg(x)) where go = TT g(x) is homogeneous, to the more general parameterized form p(x(0) o (x)))))(r)) for some convex go-
The dimensions of the parameter space are increased from de I 10 to the order of 105, the dimensionality of small convex regularizer.
2.1 Jointly Sampling Prior and Posterior
To compute the gradient of the marginal likelihood, we need to evaluate the expectations in 66 In general, this is analytically intractable and needs to be done numerically, though propertics such as positive homogeneity can make this easier (Vidal al. 2020 Sec. 3.3).
First assume that both the likelihood E(y]z) and regularizer go are C² (in x)). Observe that (6) is a difference of expectations of Vogo(=) over two probability distributions, namely over the prior p(z(0) and posterior p(rly, 0). Following the setup of Vidal et al. 22020), consider two Markov chains based on the unadjusted Langevin algorithm (ULA) to approximate these expectations (Dalalyan 2017 Durmus & Moulines 2017).
For step-size parameters 7,7" > 0, the chains are given as follows, where Zk, Zi are i.i.d. standard Gaussians with identity covariance matrix:
R : X = = V + (Xx) + 22
(7a)
Ry X Xx 7 2 (Xx) / Z
(7b)
It sufficient. that p(y) f(vlz) positive p.d.f.s, x()) is differentiable, and that any 0 int(@), there noog around @ such that [Ver(y, uniformly majarized by some ğ(x) that integrable w.r.t. f(y)=).

---

Published in Transactions on Machine Learning Research (12/2024)
The above Markov kernels R and Ryy target the posterior p(x 0) cxp(-ff(x)- - go(x)) and the prior p(x[0) xx exp(-go(x)) respectively. This allows us to use Monte Carlo integration in 100, and approximate the marginal posterior gradient Va logp()))) We note that it is possible to lift the smoothness requirement on or 90 by using a different Markov kernel such as the Morcau-Yoside unadjusted Langevin algorithm (MYULA, Durmus et a 2018) or the proximal unadjusted Langevin algorithm (PULA, (urmuss et al.) 2019). With the approximate marginal posterior gradient, we can apply stochastic gradient descent on 0 with step-sizes 6n yielding the stochastic approximation proximal gradient (SAPG ULA) Algorithm | presented in Vidal
et al. 2020
Algorithm SAPG ULA
Require: Initial (0o, Xi, X0), (ōn, MA)NNN, O MC step-size parameters 7, 7,2,, iterations N
1: for = 0 to N- --1 1 do
2: if n > 0 then
3:
Sct X X=-1
D. Postorior chain b. Prior chain
m., -1
Set X xn-1 mn-1
4:
5:
end if
6:
for k = 0 to mn - 1 do
Sample X ~ 8, (Xp:)) 8: Sample X ~ ,, (XP :
7:
b Update posterior chain
D Update prior chain
9: end for
10: Set 0n+1= = To (On + 5ut1 Eien (Vogo(X)) - Vogo (XT))
b SOUT. step
11:
end for
In Step 10, the estimated parameter 0₂ is updated with gradient ascent on p(y]0), where the expectations
6 are approximated with Monte Carlo integration with my samples, done in Steps to 9. A projection II onto the compact set θ is then imposed after each update of 0. The following result states that convergence is attained when the Markov chain updates and the updates for ge happen in an alternating fashion. Moreover, single sample for the Monte Carlo integration over the prior and posterior is sufficient for convergence. This is key to computational efficiency, by reducing the number of samples needed from the Markov chain.
Theorem 1 (De Bortoli et al. 2020 Theorem 6). Suppose that - 00) is convea w.r.t. 0, and that A Asumptions / to 3 hold. Under certain technical Lipschitz conditions and decaying step-sizes, single sample is sufficient, i.e,, mn lends to almost sure convergence of (0m)neN to some mazimizer 0_ arg max ply 0).
We restate another version of this theorem in more formal manner in Section 3 In particular, we verify that the required assumptions hold when using the convex ridge regularizer architecture for go. This is done in the constant step-size setting, which instead gives convergence in expectation.
2.2 Reflecting Markov Chains
While the Markov chains 7 sample from biased versions of the target posterior and prior distributions in the case of unrestricted domains, this does not translate to constraincd sampling. Indeod, the natural domain for images requires that pixel values be non-negative. Negative pixels may cause problems in particular for the case of Poisson imaging, where the likelihood is zero for negative measurements. Melidonis et al. (0222) consider modifying the chains to force the Markov chain samples to adhere to the non-negativity constraint by projecting or reflecting into the target domain. minor modification of Equation 7 leads to the following reflected Markov chains, where the absolute value of vector |o| to be taken componentwise:
Ry X = | - VV,S + (Xx) + 2
(8a)
Ry X 18x V V (X+) + /2//1
(8b)
7

---

Published in Transactions on Machine Learning Research (12/2024)
We can rewrite the target densities in the following manner. the positive orthant C' = R
Denote (x(x) the characteristic function of
0,
x x C,
lo(x)
44 4R R4 = fx x > 0, i= 1,1-1,,11. l....dd. x x C.
(9)
HOC,
Let Eo(x) - [g(x)) + go (z)). From the above reflections on the boundary of R4, the original target distribution p(x(0) and p(zly. 0) are truncated and modified into the following densities:
(=((6())+((())
p(z²0)- = E (0(((())+(4
p(aly. 0) = Juq e- c-
(10)
Jn e=)o
It can be shown that under suitable assumptions, the reflected SDE admits a unique invariant measure on
Rd and the invariant density is given by (10) (Melidonis ct al., 2022 Thm 3.4). The Markov chains (8a) and
(8b) are the discrete counterparts of the reflected SDE on R constructed in (Melidonis et al. 2022
We note that reflection onto the positive orthant is preferred over projection since work with densities. In particular, projection will assign positive probability to the boundaries of the target domain Rd, which has zero Lebesgue measure. Therefore, the law of the projected SDE will not admit a density with respect to the Lebesgue measure, invalidating the requirements of presented theory. We leave the analysis of the reflected Markov chains to future work.
2.3 Extensions to Datasets - Mini-Batching
SAPG-ULA, Algorithm 1 is the basic method for solving the MLE problem (2) in the case where ge needs only be trained on single datapoints, where the problem sufficiently well-posed. The Markov kernels can be chosen to be either standard ULA (7) or reflected UTA ( as applied in Vidal et al.) (2020); Melidonis et al. [2022) respectively. However, these implementations are limited to reconstructions for single images, well as only one or two scalar parameters for 6. In the case where is a neural network, the MLE problem for a single image may be severely overparameterized, leading to slow convergence. Inspired by traditional machine learning, we aim to train a neural network regularizing prior with more data, such as on a standard image dataset like STL-10.
The first computational difficulty comes from Step 10 of Algorithm | where Vege needed. Since computing gradients of go over all images in the datasct is computationally infeasible, mini-batching is necessary in practice, similarly to how stochastic methods such SGD or Adam are applied rather than full-batch gradient descent. In this section, we present a straightforward extension of the SAPG scheme to the batched case in Algorithm 2 by using multiple posterior Markov chains and changing the order in which the prior Markov chains are updated. This has the upside of only having to store one mini-batch worth of gradients at any given time, which alleviates memory issues by updating go in a stochastic manner.
Suppose that we are given fixed batched noisy measurements (K)) such as when a noisy image dataset is partitioned. The batched SAPG scheme consists of multiple posterior Markov chains {xmm)) where the posterior chains X are updated according to Ry.0- Only one prior chain X" is employed, updated according to Ry,0. One chain is sufficient since the prior Markov kernel R does not depend on the measured data Y,, saving memory compared to running multiple prior Markov chains.
/6=1
The proposed batched SAPG-UTA method summarized in Algorithm 2 This extends the SAPG ULA scheme from single images datascts. The Markov kernels Ry,0., Ry.00 need only have stationary distribution approximating the target posterior p(aly, 0) and prior p(x[0) respectively, which can be of the form (7), another kernel which can be chosen based on the regularity properties of the target boglliielhoood and prior.
Compared lo classical supervised training, we make few approximations in Algorithm 2 The first approximation comes from the lack of ground-truth data, which is replaced with the prior given by the regularizer ge- Another approximation comes from approximating the expectation over prior priorpposterior with only a single sample, motivated by Theorem [ There will also be bias induced by the choice of ULA for Markov chain, well as the joint updates of parameters 0 and Markov chains.
8

---

Published in Transactions on Machine Learning Research (12/2024)
Algorithm Batched SAPG ULA
Require: Number of batches B, batched measurements 18 Initial {6o. X0,6, X010=1> Ma)eN, O parameters 7, , iterations N,
for -0 0 to N W-1 do
if 0 then
3:
Set X m.
1,b for b
> Posterior chain for each batch
Set X X
Single prior chain
end if
for k = 0 to m, - do
7:
Sample X ~ R 00 (X for b to B do
> Update prior chain
8:
Sample k+1 ~ 0y (Xb Y)) 10: end for
9:
b Update all posterior chains
11: end for
12: 0n+1 To (0,
duti B
2bml Em (xgg((X)) (Xgo
Accumulated-grad SOUL step
mm
13:
end for
2.4
Connection with Adversarial Regularization
We additionally note a connection with adversarial regularization. Under integrability conditions and by Fubini's theorem, (6) shows that the problem maximizing the marginal likelihood p(y(0) in 2 is equivalent maximizing
sup Eejolg eg
(11)
9
The Bayesian MLE problem 11 can be thought of as maximizing the value the regularizer weighted over all possible æ (given by the prior expectation), while minimizing the regularizer for good reconstructions (given by the posterior expectation). This should be contrasted with adversarial regularization (Lunz et al., 2018), which has loss functions of the form
sup En [g)( Exwerue [ge(1)].
(12)
Comparing, the expectations over the noisy and ground-truth distributions of the adversarial loss 12 are replaced with versions that depend the parameter 0, in particular, the prior and posterior distributions in 11 The prior can be interpreted as the family of all noisy reconstructions, and the posterior reconstructions given measurements, which replaces ground-truth images. The Bayesian MLE problem can thus be formulated as adversarial regularization between the prior and posterior distributions.
The difference of expectations form of [1 and 12 are also related to the Wasserstein-1 metric, which is a notion of distance between two finite measures. This connection is noted for example Arjovsky et al. (2017). For two probability measures Rd, the definition of the Wasserstein-1 metrie Wi and Sertoooovic--Rubensee duality result are follows (Ambrosio et al. 2013 Villani 2009):
/ ||x yl dy(2,)) sup gET().))
[ [fodu
Wi(H,v) = inf
g dy
(13)
gel gEl-Lip
where I(H,v) )) is the set of all couplings. i.c. probability measures on Rd x Rdd with marginals and v respectively. Therefore, if the suprema in 11 and 12 were instead over go that are 1-Lipschitz, the problems are reduced to finding the so-called Kantorovich potentials that achieve the supremum in (13 between the prior and posterior distributions, or the noisy and true distributions respectively.
3 Convergence Analysis
In this section, we first introduce the architecture used for go- We then state a summary of the main theorems of De Bortoli et al. (2021), which demonstrates convergence in expectation to the minimizer, with

---

Published in Transactions on Machine Learning Research (12/2024)
the associated assumptions provided afterwards. The required assumptions will then be verified for the convex ridge regularizer, which together show convergence of SAPG for our formulation. We provide short intuitive explanations for each assumption and how the regularizer satisfies them, with detailed proofs left to Appendix
The regularizer architecture used for these experiments is the convex ridge regularizer (CRR) Goujon et al. 2022 This takes the form of a one-layer network:
90 : 1
(w: x).
(14)
where with some abuse of notation, 1/4 R + IR are convex profile functions (ridges) applied element wise, and w, E R are learnable weights. The parameters 0 consist of the weights w; and parameters of the ridges the The gradient of the regularizer can be expressed as
Vgo(x) = WwowWoo
(15)
where W = [[w1... WelT E RCdd is linear, and a R + RC is a component-wise activation function with components o; = vif : + R. Note that due to the convexity of thi the derivatives o; are increasing. By restricting the parameters such that o, are strictly increasing, the CRR can be made to be uniformly strongly convex.
For case of notation, first suppose have sufficient regularity such that Fisher's identity 11 holds, which will be justified in Proposition 1 and Assumption Let the negative log joint density log p(z, v/0) be given by
Ex(z) = fx) + (z))..
(16)
Note that (z))(z) the diffusion term of the Langevin step (7a) from posterior p(zly, 8), while (ag is the diffusion term of the Langevin step (7b) from prior p(x: We will use ergodicity of both chains to demonstrate (biased) convergence in expectation of 0.
In the following theorem, we consider the scenario where a fixed mn = is used in Algorithm I for the Monte Carlo integration, and a fixed step-size > 0) is used for the Markov chains. Theorem 2 summarizes Theorems 4-6 of De Bortoli et al. (2021), , which concern the convergence properties the SOUL algorithm. An extension to nonconvex - logp(y)() (0) is given when θ a sufficiently regular manifold De Bortoli et al. [2021] Appendix B), achieving (a.s.) convergence some point 6 with -Vs log p(y(()))*) lying on the normal
cone of θ at 8*
Theorem 2. [De Bortoli et al. 2021, Theorem 4, 5, 6). Assume that Assumptions to 6 hold and log p(y conver w.r.l. O. Assume further that the slep-sizes (5a)EE satisfy supner on < IT, and that 7 S min(1, 2mz)) where ma is the tail coercivity constant in Assumption
5. Then for any N",
k= òy logp(ul0m) sim E max log p(y Er=
I <O(/+ VT).
(17)
n-000 dEE
2k OF
Theorem 2 bounds the bias of the SAPG method terms of the step-size of the Markov chains. This covers the Gaussian dcconvolution experiments the next section, but not the Poisson denoising experiments which employ reflected ULA. However, we believe that the SAPG algorithm is also convergent. when the Markov kernel is given by reflected ULA, and refer proof of convergence to future work.
We note that almost sure convergence can also be obtained by using decrcasing step-sizes or increasing batch-size m., (De Bortoli el al. 2021 Thm. 1-3). For the sake of simplicity and computational efficiency, we consider only constant step-sizes on and 7,7,, and fix m, = for our experiments. For an empirical study into the bias of the SAPG method for low-dimensional (where ground-truth maximum likelihood estimators are known), we refer to (De Bortoli et al.) 2021: Vidal et al.] (020)). These works find that the bias is not significant. if the hyperparameter guidelines are followed, as do in our experiments.
10

---

Published in Transactions on Machine Learning Research (12/2024)
4 Experiments
In the following image reconstruction experiments, we will compare the proposed unsupervised MLE method Algorithm 2 with the L-gradient-step supervised training method as detailed in Goujon et al. (2022). We use the same choice of architecture as Goujon et al. (2022) extended from grayscale images to color images. In 14 and (1), W modelled using a convolution, and the d, are parameterized using linear splines. This results in the being degree-2 splines, namely piecewise quadratic with continuous derivatives. We consider having 32 splines for our parameterization, and W being convolutions with 3 (color) imput channels and 32 output channels. Additional implementation details including the choice of W can be found in Appendix B
The supervised gradient-step (GS) training method of Goujon et al. (002)) involves performing gradient steps to minimize the variational functional induced by the convex regularizer, which is then trained to minimize the & distance between the reconstructed image and the ground truth. We compare the supervised GS training to our proposed unsupervised training (SAPG), which does not require noiso-free images.
Equivariant imaging poses another interesting baseline to compare against for unsupervised learning. We use the setup in Scanvic et al. (2023). which recently extends the El framework using scale transforms, used in cases where the forward operator and noise models are equivariant under the standard rotations and reflections. The learned reconstruction of the EI model takes the form of an end-to-end neural network, using the state-of-tho-art SwinIR architecture for image roconstruction, using 11.5M parameters (Liang et al., 2021). We train using the given parameters for up to 200 epochs, but use early stopping at 25 epochs due to divergence issues during training.
For basclines that do not rely on machine learning, compare with the total variation prior g(x) = V Rudin et al., 1992), , as well as using deep image prior (DIP) for reconstruction (Tachella et al. 2021 Ulyanov et al. 2018). For TV regularization, the regularization parameter is chosen via a grid search to maximize PSNR with respect to the ground truth. We use TV to additionally compare with the single-image self-supervised approach of (Vidal et al.) 2020). which automatically finds the regularization parameter using MLF.. Note that the reported quality for TV will be higher than if (Vidal et al. 2020) is used, due to the comprehensive search. For DIP, we compare using a U-Net as well as a convolutional neural network, optimized using Adam and gradient descent as in Tachella et al. (2021). The reported PSNR is taken to be the highest with respect to the ground truth along the optimization trajectory.
To evaluate our MLE-trained regularizers, we consider the posterior mean and MAP estimators under the posterior chain. The posterior mean is evaluated after 2 x 10 and 1 x 103 iterations of the posterior Markov chains with 5 x 1 iterations to warm-start, labeled 'MMSE' in the figures and tables with superscripts indicating the number of posterior samples. To compute the MAP estimators of our trained regularizers, we minimize the modified negative log-posterior (x = fy(x) (x) + Ags(ar), where fy is the negative likeliod (fidelity) function corresponding to the corrupted data y, and go is the learned convex regularizer, and À >> 0 is regularization parameter. The regularization parameter À is chosen via a grid search À E {0.1, 0.2, .... 1.0] to maximize PSNR after training. To compute the MAP estimate, the negative log-posterin posterior Pe is minimized using the Adam oplimizer for up to 10a iterations due to being faster than gradient descent, with learning rate 10 and other parameters as default. We report the training and testing time in Appendix C
To explore the strength of the unsupervisod SAPG training method along with the supervised gradient-step method of Goujon et al. (2022), we additionally compared with a weakly-supervised training regime, by warm-starting the regularizer using adversarial training based on Lunz et al. (2018). This uses the formulation [2] with the noisy distribution being corrupted images, and is labelled as WS-SAPG in following tables. The weakly-trainee regularizer is first provided ground truths of a fixed 2% of the training images, and the CRR is trained to be maximized on the corresponding corrupted images and minimized on the ground truth. The loss for warm-starting can be expressed as follows, where the last term controls the Lipschitz constant of go-
M
min Cwarmstart (0)
1 M
go(atrue) g0 ) + mix(()Lip(s))) Lip(ge)).
(18)
After warm-starting, the adversarially learned CRR is passed as an initialization into the standard SAPG method, where only the corrupted training dataset is provided for use in Algorithm
11

---

Published in Transactions on Machine Learning Research (12/2024)
Table 1: Table of PSNR (dB) the convex ridge regularizer (CRR) trained for Gaussian deconvolution. The initial PSNR of the corrupted image 22.38dB and the imitial SSIM is 0.612. The MAP regularization parameter for the MLE is À 0.6, and values are averaged over 50 test images. We observe only small drop in performance of the unsupervised SAPG and weakly-supervised WS-SAPG compared to the supervised gradient-step (GS) method of around 0.5dB, and outperform total variation. Equivariant imaging (EI) performs very well, possibly due to a more powerful neural network parameterization.
CitR-based
End- condd-t-eddd/aae
GS
SAPG
WS-SAPG
DIP
MAP
MMSE
MMSB
MAP
MMSE34
MMSB
UNet-Adam
UNet-GD
CNN-Adom
CNN-GD
PSNR
25.74
25.19
25.09
25.25
25.05
25.06
27.11
21.75
25.50
21.80
24.67
24.30
SSTM
0.803
0.795
0.779
0.780
0.502
0.708
0.767
0.846
0.796
0.784
0.565
0.723
0.739
Parana
15,537
11,504,163
1,036,032
187,992
(a) Ground truth
(b) Corrupted
(c) TV
(d) DIP-UNet
(e) DIP-CNN
(f) EI
(8) GS
(h) SAPG MMSE (1) SAPG MAP
(i) Std dev
Figure 1: Visual comparison of various reconstructions of blurred test image with Gaussian noise. The unsupervised MMSE (h) and MAP (i) reconstructions of the proposed SAPG method contain visual artifacts compared to the gradient-step training method. The standard deviation (j) shows the uncertainty around edges. We observe that El has the sharpest looking image, which is due to explicit knowledge of the forward operator as blur kernel, as well having a much richer end-to-end parameterization.
4.1 Gaussian Deconvolution
The first set of experiments Gaussian deconvolution on the natural image dataset STL-10. The images consist of 96 96 color images, which are then corrupted with Gaussian blur with kernel size 5 and blur strength 1, followed by additive 5% Gaussian noise. The negative log-likelihood for noisy blurred image y is thus given by
fx(x)
20 AA
|A w112,
where A is the blur kernel. For SAPG, the step-sizes 7,76 for the likelihood and prior Markov chains Ry Ry,00 respectively are given by = = 1e-4.
Table D details the reconstruction quality in terms of PSNR of the proposed SAPG method as well as various supervised and unsupervised baselines. We observe that the proposed unsupervised SAPG method is able to perform closcly with the supervised t-gradient-step method, with gap of only 0.5dB. Morcover, we observe that the MAP estimate with the learned reconstruction priors are generally better than the posterior mean estimates. The SAPG-learned prior also is competitive with/better than other priors such TV DIP,
12

---

Published in Transactions on Machine Learning Research (12/2024)
Table 2: Table of PSNR (dB) of the convex ridge regularizer trained for Poisson denoising. We observe competitive performance of the learned CRR-based priors compared end-to-end reconstructions, and only minor drop in performance between the proposed unsupervised SAPG method and supervised GS method. The initial average PSNR of the corrupted images is 21.16dB and the initial SSIM is 0.641. The MAP regularization parameter for the MLE is À 0.6. Averaged over 50 test images.
CRR-based
Ead-to-end /model-based
SAPG
WS-SAPG
GS
EI
TV
MAP
M
MMSE166
MAP
MMSE
MMSE
UNet-Adam
UNet-GD
CNN-Adam
CNN-GD
PSNR
28.22 29
27.91
26.65
27.52
28.14
27.34
27.75
24.09
26.50
21.81 0.375
26.28
26.59
SSIM
0.858
0.806
0.813
0.810
0.878
0.836
0.816
0.861
0.722
0.818
0.783
0.799
(a) Ground truth
(b) Corrupted
(c) TV
(d) DIP-UNet
(e) DIP-CNN
(f)
(g) GS
(h) SAPG MMSE (i) SAPG MAP
(j) Std dev
Figure Visual comparison of reconstructions for Poisson denoising. The proposed unsupervised SAPG method, shown in subfigures (h) and (i) both show significant denoising, but with the presence of artifacts. EI has color artifacts while the supervised GS method has more textural artifacts. DIP has strong smoothing effect, but also induces strange visual arlifacts around the larget shown in subfigure fd
the latter of which still has a lack of theoretical properties. Adding weak supervision does not significantly change the performance here.
Figure shows visual comparisons of the various methods applied to a corrupted test image. El has the best perceptual quality out the tested images with the powerful end-to-end parameterization. The reconstruction of the proposed method is reasonable, but has artifacts compared to the supervised gradient step method, possibly due the ill-posedness of the problem. We are also able to plot the standard deviation of the Markov chain samples
in Figure clearly showing areas of interest. around the edges.
4.2 Poisson Denoising
Poisson denoising naturally arises in low-photon imaging, where the measurement takes integer pixel values. We can model the measurement using a Poisson random variable y > Pois(πæ) taken pixel-wise over an image space RR, where 7) parameter chosen such that the expected mean pixel intensity y fixed mean intensity value (MIV). The forward operator is set to be the identity. In this case, the Poisson negative log-likelihood for positive-integer-vlluss measurement. y Pois(1ιz) is given follows, where 44 is the
13

---

Published in Transactions on Machine Learning Research (12/2024)
Table Table of PSNR (dB) and drop in PSNR compared to Table | of the convex ridge regularizer trained for Gaussian deconvolution when applied to uniform deconvolution. We observe significantly better transfer properties of the learned CRR- based regularizers compared the end-to-end El method, with only a small drop in performance between the supervised and unsupervised versions. We also observe better transfer of the unsupervised SAPG training compared to supervised GS training. The initial PSNR of the corrupted image is 20B1 and the initial SSIM is 0.475. The MAP regularization parameter for the MLE is = 0.7. Averaged over 50 test images.
SAPG
WS-SAPG
GS
EI
TV
MAP
MMSE
MMSE
MAP
M
MMSE
PSNR
23.90
23.73
23.62
23.55
23.79
22.81
23.05
21.37
22.20
Drop
1.84
1.46
1.47
1.55
1.46
2.24
2.01
5.74
2.55
SSTM
0.739
0.720
0.691
0.697
0.732
0.645
0.659
0.637
0.676
Drop
0.064
0.075
0.088
0.083
0.07
0.121
0.108
0.209
0.11
characteristic function 7 in 19, y! denotes factorial, and division of images is to be taken pixel-wise:
a
fx(x) Σ [(nx); - y4 log(mas;) + log(u:!) + 4 (x),
(xfy(x) = n/1 - y/æ,
The Poisson log-likelihood does not have globally Lipschitz gradient, degenerating around the boundary of Rt when is near 0. This makes the Poisson problem more difficult as compared to the Gaussian case due to requiring both approximating the log-likelihood for zero-valued pixels, as well as needing small step-size. Therefore, we use the following modified negative log-likelihood, for some mollification parameler > MMelddonis et al. 2022),
(x)
Σ In(2; + b) - log(n(2;; b)) + log(u!)) 1=1 Vafy(x) = nll 1/(x/
The introduction of mollifies the likelihood near zeros and makes (xf Lipschitz on R7. We use the paramcter b MIV /100 suggested in Melidonis et al. (2022), where the mean intensity value (MIV) is the mean of y PPois(nz) The parameter 1) is chosen such that the mean intensity value is set to be AMrV 25 for light Poisson denoising. For SAPG, the step-sizes for the likelihood and prior Markov chains are given by 7 5e-6, 7 = le-5 respectively. For El, the best performance is obtained by modifying the training loss to use the Poisson unbiased risk estimator instead of Stein's 's unbiased risk estimator, and using rotation transforms instead of scale transforms (Chen et al. 2022ah.
Table compares the PSNR of the recovered images using the proposed method against the previously detailed baselines. We again observe a reasonable gap of 0.15dB between unsupervised method and the supervised method with the same architecture. Figure 2 however still demonstrates the presence of some visual artifacts of our unsupervised method when compared to its supervised t-gradient-step counterpart. Moreover, the CRR-based methods exhibit patchy artifacts that are similar to the ground-truth image, while El exhibits varying color artifacts. For this problem, we achieve much closer performance to equivariant imaging compared to Gaussian dcconvolution. This highlights the difficulty of Poisson denoising for end-to-end models, as observed previously due to the high Lipschitz constants and large data range, and an empirical benefit of variational regularization when dealing with high-variance noise.
14

---

Published in Transactions on Machine Learning Research (12/2024)
aa
Ground truth
(b) Corrupted
(c) ET
(d) TV
GS
(f) SAPG MMSE
(g) SAPG MAP
(h) Std dev
Figure 3: Visual comparison of reconstruction for uniform deconvolution when using models trained for Gaussian deconvolution. We observe more residual blur for El when compared to the model-basee priors in the bottom row.
4.3 Transfer to Different Forward Operator
In addition to the Gaussian deconvolution experiments, consider robustness of the learned models with respect to changes in forward operator. We include the change in forward operator to provide some empirical justification to the choice of learned regularizer opposed to end-to-end methods. We evaluate the trained models on a x 5 uniform blur kernel with the same additive 5% Gaussian noise. In particular, we consider supervised and unsupervised CRR, as well as the end-to-end El method. The models trained for Gaussian deconvolution Section 4.1 are applied directly to uniform doconvolution, allowing for a change of forward operator in the fidelity for the variational model-based CRR methods.
Table 3 shows the reconstruction results for transferring the learned models between different tasks. Compared to the learned regularizers and TV, we observe a significantly larger drop in performance from the El model, which is trained on a specific forward operator, which introduces more blur artifacts. This may be due to the scnsitivity of end-to-end networks to changes in distribution. On the other hand, the model-based gradient step reconstruction still produces reasonably good results, due to the underlying variational form. The unsupervised versions have similar performance, albeit with added visual artifacts. This may be due image prior components that are difficult learn from only noisy patches.
We note that within learning convex ridge regularizers, our SAPG method drops less PSNR than the supervised gradient-step method when changing forward operator, suggesting better generalization. While both methods utilize the forward operator during training, the gradient-step training unrolls through 10 optimization iterations, whereas the unsupervised method uses it as a consistency term. One possible explanation is that the Langevin diffusion the unsupervised method expands the available data to also cover some small change in forward operator, while the supervisod method sees only data that is provided from the original forward operator.
We additionally note that better generalization under change of forward operator can be achieved using unrolling schemes, which allow for the forward operator to be explicitly specified during the testing phaso MMonaa el al. 2021 Zhao el al. 2023). Such unrolling methods are typically trained in the supervised regime Li et al. 2020b). While we found that a direct application of the unsupervised end-to-end El loss in (Scanvic ct al. 2023) to an unrolled architecture led to unsatisfactory results in our deconvolution experiments, a more nuanced approach constructing more transferable architectures could be an interesting direction for future work.
15

---

Published in Transactions on Machine Learning Research (12/2024)
5 Conclusion
We proposed unsupervised method of training convex neural network-basod image prior based the stochastic approximate proximal gradient algorithm, extending previous work in both model expressiveness and dataset size. Prior works consider a family of total variation priors, here we consider a parameterization with significantly more paramcters. Experiments demonstrate that the proposed method produces priors that are near competitive when compared to the alternative supervised training method for various image corruption operators, maintaining significantly better generalization properties when compared to end-to-end methods. Moreover, we present convergence theory for the proposed Markov chains with the choice of convex ridge regularizer architecture.
The proposed method has two main limitations, namely runtime for training and inference, as well as the convexity assumption on the regularizer. The runtime is mainly due to slow convergence of the MCMC and additional stochasticity when training the CRR. A detailed comparison the train and test runtimes of the compared methods is given in Appendix [C] where our method takes an order maguitude longer. The convexity assumption of the regularizer allows for convergence of SAPG, but may be restrictive in terms of expressiveness.
Interesting future works could include convergence analysis in the discrete case or with accelerated versions of the training algorithm such as Nesterov accelcration or Adam instead of SGD, or alternative Markov kernels like SK-ROCK instead of ULA (Pereyra et al. 2020) Other possible directions are to add an informative prior on 0, which may change the sampling dynamics, or devising fully unsupervised methods for warm-starting training. We also believe that the strong convexity assumption of the convex ridge regularizer can be dropped, due the simple structure and coercive tails. promising alternative weak convexity, which has been shown to significantly outperform similar convex models (Shumaylov et al.) 2023). Future work could include theoretical analysis of the SAPG method in this setting using techniques from Durmus & Moulines (2017): De Bortoli et al. (2021). Another interesting work would be to extend this work to learning diffusion models from single copies of noisy measurements.
Acknowledgments
HYT was supported by GSK.ai and the Masason Foundation. MP was funded by the UK Research and Innovation (UKRT) Engineering and Physical Sciences Research Council (EPSRC) through grants EP/T007346/1, EP /LP/V06134/1 and EP/ /P0000767//1 CBS acknowledges support from the Philip Leverhulme Prize, the Royal Society Wolfson Fellowship, the EPSRC advanced career fellowship EP EP/ RP/VO29428/B EPSRC Grants EP /S/002644//L and EP/ T003553// 1,, EP EP/N014588/ 1, EP EPTT011961// the Wellcome Innovator Awards 215733/7/19/7 and 221633/Z/20/2 the European Union Horizon 2020 rescarch and innovation programme under the Marie Sklodowska-Cut Grant agreement No. 777826 NoMADS, the Cantab Capital Institute for the Mathematics of Information and the Alan Turing Institute.
References
Asad Aali, Marius Arvinte, Sidharth Kumar, and Jonathan I Tamir. Solving inverse problems with score-based generative priors learned from noisy data. In 2023 57th Asilomar Conference on Signals, Systems, and IEEE, 2023.
Compulers, pp.
Giovanni S Alberti, Ernesto De Vito, Matti Lassas, Luca Ratti, and Matteo Santacesaria. Learning the optimal Tikhonov regularizer for inverse problems. Advances in Neural Information Processing Systems, 24:25205-252116 2021.
Luigi Ambrosio, Alberto Bressan, Dirk Helbing, Axel Klar, Enrique Zuazua, Luigi Ambrosio, and Nicola Gigli. user's guide to optimal transport. Modelling and Optimisation of Flows on Networks: Cetraro, Italy 2009, Editors: Benedetto Piccoli, Michel Rascle, pp. 1-155, 2013.
Anestis Antoniadis and Jianqing Fan. Regularization of wavelet approximations. Journal of the Americon Stolistical Association, 96(455):933 199999999779919 2001.
16

---

Published in Transactions on Machine Learning Research (12/2024)
Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. In International conference on machine learning, pp. 214-223. PMLR, 2017.
Simon Arridge, Peter Maass, Ozan Öktem, and Dibiane Schönlieb. Solving inverse problems using data-driven models. Acta Numerica, 28:1-174, 2019.
Guha Balakrishnan, Amy Zhao, Mert R Sabuneu, John Guttag, and Adrian V Dalca. An unsupervised learning model for deformable medical image registration. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 9252-9260, 2018.
Georgios Batzolis, Jan Stanczuk, Marola-Bibians Schönlicb, and Christian Etmann. Conditional image generalion with score-based diffusion models. ar preprint arXiv:2111. 13666, 2021.
Martin Benning and Martin Burger. Modern regularization methods for inverse problems. Acta numerica, 27: 1-111, 2018.
Martin Burger and Stanley Osher. Convergence rates of convex variational regularization. Inverse problems, 20(5):1411, 2004.
Daniela Calvetti and Erkki Somersalo. Inverse problems: From regularization to Bayesian inference. Wiley Interdisciplinary Reviews: Computational Statástics, 10(3):c1427) 2018.
Marcello Carioni, Subhadip Mukherjee, Hong Ye Tan, and Junqi Tang. Unsupervised approaches based on optimal transport and convex analysis for inverse problems in imaging. arXiv preprint 0t::2311.0897. 2023.
Bradley P Carlin, Alan E Gelfand, and Adrian FM Smith. Hierarchical Bayesian analysis of changepoint problems. Journal of the royal statistical society: series C (applied statistics), 41(2):389 400,, 1992.
Chad Carson, Serge Belongie, Hayit Greenspan, and Jitendra Malik. Blobworld: Image segmentation using expectaton--maxiiniaatoon and its application image querying. IEFE Transactions on poltern analysis and machine intelligence, 24(8):1026 1038, 2002.
Antonin Chambolle and Thomas Pock. An introduction to contimuous optimization for imaging. Acta Numerica, 25:161-319, 2016.
Dongdong Chen, Julián Tachella, and Mike E Davies. Equivariant imaging: Learning beyond the range space. In Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 4379-4388, 2021.
Dongdong Chen, Julián Tachella, and Mike E Davies. Robust equivariant imaging: fully unsupervised framework for learning to image from noisy and partial measurements. In Proceedings of the IEEE/CVF Conference Computer Vision and Pottern Recognition, pp.
2022a.
Junyu Chen, Eric C Frey, Yufan He, William Segars, Ye Li, and Yong Du. Transmorph: Transformer for unsupervisod medical image registration. Medical image analysis, 82:102615, 2022b.
Zezhou Cheng, Matheus Gadelha, Subhransu Maji, and Daniel Sheldon. A bayesian perspective on the deep image prior. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 5443-5451, 2019.
Hyungjn Chung, Jeongsol Kim, Michael Thompson Mccann, Marc Louis Klasky, and Jong Chul Ye. Diffusion posterior sampling for general noisy inverse problems. In The Eleventh International Conference on Learning Representations, 2023.
Arnak Dalalyan. Theoretical guarantees for approximate sampling from smooth and log-concave densities. Journal of the Royal Statistical Society Series B: Statistical Methodology, 79(3):651 676, 2017.
Valentin De Bortoli, Alain Durmus, Marcelo Pereyra, and Ana Fernandez Vidal. Maximum likelihood estimation of regularization parameters in high-dimensional inverse problems: an empirical Bayesian approach. Part II: Theoretical analysis. STAM Journal on Imaging Sciences, 13(4):1990-21) 20200, 2020.
17

---

Published in Transactions on Machine Learning Research (12/2024)
Valentin De Bortoli, Alain Durmus, Marcelo Pereyra, and Ana F Vidal. Efficient stochastic optimisation by unadjusted Langevin Monte Carlo: Application to maximum marginal likelihood and empirical Bayesian estimation. Statistics and Computing, 31:1-18 18, 2021.
Arthur Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from incomplete data via the em algorithm. Journal of the royal statistical society: series B (methodological), 39(1):1-22, 1977.
Steven Diamond, Vincent Sitzmann, Felix Heide, and Gordon Wetzstein. Unrolled optimization with deep priors. Xiv preprint Xivv 705.08041, 2017.
David I. Donoho. Compressed sensing. JEEE Transactions on information theory, 52((4:1189-13300 2006.
Randal Douc, Fric Moulines, and David Stoffer. Nonlinear time series: Theory, methods and applications with R ezamples. CRC press, 2014.
Alain Durmus and Eric Moulines. Nonasymptotic convergence analysis for the unadjusted Langevin algorithm. The Annals of Applied Probability, 27(3):1551-1887, 2017.
Alain Durmus, Eric Moulines, and Marcelo Pereyra. Efficient Bayesian computation by proximal Markov chain Monte Carlo: when Langevin meets Moreau. STAM Journal on Imaging Sciences, 11(1):473-50) 4066 2018.
Alain Durmus, Szymon Majewski, and Blažej Miasojedow. Analysis of langevin monte carlo via convex oplimization. Journal of Machine Learning Research, 20(73):1 46, 2019.
Bradley Efron. Tweedic's formula and selection bias. Journal of the American Statistical Association, 106 (496):1602, 2011.
Heinz Werner Engl, Martin Hanke, and Andreas Neubauer. Regulorization of inverse problems, volume 375. Springer Science & Business Media, 1996.
Berthy T Feng, Jamie Smith, Michael Rubinstein, Huiwen Chang, Katherine L Bouman, and William T Freeman. Score-based diffusion models as principled priors for inverse imaging. In Proceedingss of the IEEE/CVF International Conference on Computer Vision, pp. 10520- 10520-10531, 2023.
Karl J Friston. Bayesian estimation of dynamical systems: an application to fmrr.. Neurolmage, 16(2):513-55) 530, 2002.
Donald Geman and Stuart Geman. Bayesian image analysis. In Disordered systems and biological organization, pp. 301 319. Springer, 1986.
Walter R Gilks, Sylvia Richardson, and David Spiegelhalter. Markow chain Monte Carlo in practice. CRC press, 1995.
Alexis Goujon, Sebastian Neumayer, Pakshal Bohra, Stanislas Ducotterd, and Michael Unser. A neural- network-based convex regularizer for image reconstruction. arXiv preprint arXim:2211.12400 2022.
Alexandros Graikos, Nikolay Malkin, Nebojsa Jojic, and Dimitris Samaras. Diffusion models as plug-and-play priors. Advances in Newral Information Processing Systems, 35:14715 14728, 2022.
Peter J Green, Krzysatof Eatuszyúski, Marcelo Pereyra, and Christian P Robert. Bayesian computation:
summary of the current state, and samples backwards and forwards. Statistics and Compuling, 25:835-866 2015.
Rémi Gribonval and Mila Nikolova. A characterization of proximity operators. Journal of Mathematical Imaging and Vision, 62(6-7):778--78) 2020.
Per Christian Hansen and Dianne Prost O'Leary. The use of the l-curve in the regularization of discrete ill-posed problems. STAM journal scientific computing, 14(6):1487-1503, 1993.
18

---

Published in Transactions on Machine Learning Research (12/2024)
Andreas Hauptmann, Subhadip Mukherjee, Carola-Bibiane Schönlieb, and Ferdia Sherry. Convergent regularization in inverse problems and linear plug-and-play denoisers. Foundations of Compulational Mathematics, pp. 1 34, 2024.
Reinhard Heckel and Mahdi Soltanolkotabi. Denoising and regularization via exploiting the structural bias of convolutional generators. In International Conference on Learning Representations, 2020.
Allard Adriaan Hendriksen, Daniël Maria Pelt, and K Joost Batenburg- Noisc2inverse: Self-supervised deep convolutional denoising for tomography. IEEE Transactions on Compulational Imaging, 6:1320 1335, 2020.
Matthew Holden, Marcolo Pereyra, and Konstantinos C Zygalakis. Bayesian imaging with data-driven priors encoded by neural networks. SIAM Journal on Imaging Sciences, :5(2):892 924, 2022.
Samuel Hurault, Arthur Leclaire, and Nicolas Papadakis. Gradient step denoiser for convergent plug-and-play. In International Conference on Learning Represenlations (2),) 2022a.
Samuel Hurault, Arthur Leclaire, and Nicolas Papadakis. Proximal denoiser for convergent plug-and-play optimization with nonconvex regularization. In International Conference on Machine Learning, pp- 9483 9505. PMLR, 2022b.
Bahjat Kawar, Noam Elata, Tomer Michaeli, and Michael Elad. Gsure-based diffusion model training with corrupted data. Transactions Machine Learning Research, 2024.
Erich Kobler, Alexander Effland, Karl Kunisch, and Thomas Pock. Total deep variation for linear inverse problems. In Proceedings of the IEEE/CVF Conference on computer vision and pottern recognition, pp. 7549-7558, 2020.
Daniel Kuhm, Peyman Mohajerin Esfahani, Viet Anh Nguyen, and Soroosh Shafierzadel-Abaddee Wasserstein distributionally robust optimization: Theory and applications in machine learning. In Operations research 8 management science in the age of analytics, pp. 130-166. Informs, 2019.
Rémi Laumont, Valentin De Bortoli, Andrés Almansa, Julie Delon, Alain Durmus, and Marcelo Pereyra. Baycsian imaging using plug play priors: when Langevin meets Tweedie. SIAM Journal on Imaging Sciences, 15(2):701 737, 2022.
Jaakko Lehtinen, Jacob Munkberg, Jon Hasselgren, Samuli Laine, Tero Karras, Milka Aittala, and Timo Aila. Noise2noise: Learning image restoration without clean data. In International Conference on Machine Learning, pp. 2965 2974. PMLR, 2018
Emanuel Levitan and Gabor T Herman. A maximum a posteriori probability expectation maximization algorithm for image reconstruction in emission tomography. IEEE transactions on medical imaging, 6(3): 185 192, 1987.
Housen Li, Johannes Schwab, Stephan Antholzer, and Markus Haltmeier. Nett: Solving inverse problems with deep neural networks. Inverse Problems, 36(6):065005, 2020a.
Yuelong Li, Mohammad Tofighi, Junyi Geng, Vishal Monga, and Yonina C Eldar. Efficient and interpretable deep blind image deblurring via algorithm unrolling. IEEE Transactions on Computational Imaging, 6: 666-681, 2020b.
Jingyun Liang, Jiezhang Cao, Guolei Sun, Kai Zhang, Lue Van Gool, and Radu Timofte. Swinir: Image restoration using swin transformer. In Proceedings of the IEEE/CVF international conference computer vision, pp. 1833-1844, 2021.
Sebastian Lunz, Ozan Öktem, and Carola-Bibiane Schönlieb. Adversarial regularizers in inverse problems. Advances in neurol information processing systems, 31, 2018.
Morteza Mardani, Jiaming Song, Jan Kautz, and Arash Vahdat. A variational perspective on solving inverse problems with diffusion models. In The Twelfth International Conference on Learning Representations, 2024.
19

---

Published in Transactions on Machine Learning Research (12/2024)
Charlesquin Kemajou Mbakam, Marcelo Pereyra, and Jean-François Giovannelli. Marginal likelihood estimation in semiblind image deconvolution: A stochastic approximation approach. SIAM Journal on Imaging Sciences, 17(2):1206 1254, 2024.
Savvas Mclidonis, Paul Dobson, Yoann Altmann, Marcelo Pereyra, and Konstantinos C Zygalakis. Efficient. Bayesian computation for low-photon imaging problems. arXiv preprint orXiv:2206.0533 05350, 2022.
Savvas Melidonis, Matthew Holden, Yoann Altmann, Marcelo Pereyra, and Konstantinos C Zygalakis. Empirical Bayesian imaging with large-scale push-forward generative priors. IEEE Signal Processing Letters, 2024.
Christopher Metzler, Ali Mousavi, Reinhard Heckel, and Richard G Baraniuk. Unsupervisod learning with Stein's unbiased risk estimator. orXiv preprint Xiv:v:1805. 10531, 2018.
Scan P. Meyn and R. L. Tweedic. Stability of Markovian processes i: criteria for discrete-time chains. Advances in Applied Probability, 24(3):542-55) 574, 1992. doi: 10.2307/14274753
Vishal Monga, Yuelong Li, and Yonina C Eldar. Algorithm unrolling: Interpretable, elficient. deep learning for signal and image processing. IEEE Signal Processing Magazine, 38(2):18 44, 2021.
Jean-Jacques Moreau. Proximité et dualité dans un espace hilbertien. Bulletin de la Société mathématique de France, 93:273-299, 1965.
Vladimir Alekseevich Morozov. Regularization of incorrectly posed problems and the choice of regularization parameter. USSR Computational Mathematics and Mathematical Physics, 6(1):242-251) 1966.
Subhadip Mukherjee, Sören Dittmer, Zakhar Shumaylov, Sebastian Lunz, Ozan Öktem, and Carola-Bibiane Bibiane Schönlieb. Learned convex regularizers for inverse problems. orXie preprint arXiv:2008.02833 028339 2020.
Subhadip Mukherjee, Marcello Carioni, Ozan Öktem, and Carola-Bibiane Schönlieb. End-to-end reconstruction moots data-driven regularization for inverse problems. Advances in Neural Information Processing Systems, 11:21413-2122,, 2021.
Subhadip Mukherjee, Andreas Hauptmann, Ozan Öktem, Marcelo Pereyra, and Carola-Bibiane Schönlieb. Learned reconstruction methods with convergence guarantees: A survey of concepts and applications. IEEE Signal Processing Magazine, 40(1):164-182, 2023.
Joseph Nagel and Bruno Sudret. A unified framework for multilevel uncertainty quantification Bayesian inverse problems. Probabilistic Engineering Mechanics, 43:68 84, 2016.
Yurii Nesterov. Introductory lectures on conver optimization: A basic course, volume 87. Springer Science & Business Media, 2013.
Arthur Pajot, Emmanuel De Bézenac, and Patrick Gallinari. Unsupervised adversarial image reconstruction. In International conference on learning representations, 2018.
Marcelo Pereyra. Revisiting maximuum---poserri estimation in log- log-concave models. SIAM Journal Imaging Sciences, 12(1):650-670, 2019.
Marcelo Pereyra and Julián Tachella. Equivariant bootstrapping for uncertainty quantification in imaging inverse problems. In International Conference on Artificial Intelligence and Statistics, pp. 4141-4149. PMLR, 2024.
Marcelo Pereyra, Luis Vargas Mieles, and Konstantinos C Zygalakis. Accelerating proximal Markov chain Monte Carlo by using an explicit stabilized method. SIAM Journal on Imaging Sciences, 13(2):905 935, 2020.
Christian P Robert, George Casella, and George Casella. Monte Carlo statistical methods, volume Springer, 1999.
20

---

Published in Transactions on Machine Learning Research (12/2024)
Christian P Robert et al. The Bayesian choice: from decision-theoretic foundations to computational implementation,, volume 2. Springer, 2007.
Gareth O Roberts and Osnat Stramer. On inference for partially observed nonlinear diffusion models using the Metropolis- Hastings algorithm. Biometrika, 88(3):603 621, 2001.
Leonid I Rudin, Stanley Osher, and Emad Fatemi. Nonlinear total variation based noise removal algorithms. Physica D: nonlinear phenomena, 60(1-4):259 268, 1992.
Jérémy Scanvic, Mike Davies, Patrice Abry, and Julián Tachella. Sule-supe learning for image super- resolution and deblurring. arXiv preprint arXiv:2312. 11232, 2023.
Johannes Schwab, Stephan Antholzer, and Markus Haltmeier. Deep null space learning for inverse problems: convergence analysis and rates. Inverse Problems, (52):0225008, 2019.
Zakhar Shumaylov, Jeremy Budd, Subhadip Mukherjee, and Carola-Bibiane Schönlieb. Provably convergent data-driven convex-nonconvex regularization. In NewrIPS 2023 Workshop on Deep Learning and Inverse Problems, 2023.
Charles M Stein. Estimation of the mean of a multivariate normal distribution. The Annals of Statistics, pp. 1135-1151, 1981.
Andrew M Stuart. Inverse problems: a Bayesian perspective. Acta numerica, 19:451-559, 2010.
Julián Tacholla, Junqi Tang, and Mike Davies. The neural tangent link between cmn denoisers and non-local filters. In Proceedings of the IEEE/CVE Conference on Computer Vision and Pollern Recognition, pp. (618-8627) 2021.
David YW Thong, Charlesquin Kemajou Mbakam, and Marcelo Pereyra. Do Bayesian imaging methods report trustworthy probabilities? arXiv preprint arXie:2405.cm 08179, 2024.
Dmitry Ulyanov, Andrea Vedaldi, and Victor Lempitsky. Deep image prior. In Proceedings of the IEEE conference on computer vision and patern recognition, pp. 9446 9454, 2018.
Dave Van Veen, Ajil Jalal, Mahdi Soltanolkotabi, Eric Price, Sriram Vishwanath, and Alexandros G Dimakis. Compressed sensing with deep image prior and learned regularization. arXiv preprint arXiv:1806.06 06438, 2018.
Singanallur V Venkatakrishnan, Charles A Bouman, and Brendt Wohlberg. Plug-and-play priors for model based reconstruction. In 2013 IEEE global conference on signal and information processing, pp. 945-948. IEEE, 2013.
Ana Fernandez Vidal, Valentin De Bortoli, Marcelo Pereyra, and Alain Durmus. Maximum likelihood estimation of regularization parameters in high-dimensional inverse problems: An empirical Bayesian approach Part I: Methodology and experiments. SIAM Journal on Imaging Sciences, 13(4):1945 1989, 2020.
Cédric Villani. Optimal transport: old and new, volume 338. Springer, 2009.
Yinhuai Wang, Jiwen Yu, and Jian Zhang, Zero-shot image restoration using denoising diffusion mull-space model. In The Eleventh International Conference on Learning Representations, 2023.
Darren Wilkinson. Stochastic modelling for systems biology. Chapman and Hall/CRC 2018.
Kai Zhang, Yawei Li, Wangmeng Zuo, Lei Zhang, Lue Van Gool, and Radu Timofte. Plug-and-play image restoration with deep denoiser prior. IEEE Transactions on Pattern Analysis and Machine Intelligence, 44 (10):6360-6376, 2021.
Yanan Zhao, Yuclong Li, Haichuan Zhang, Vishal Monga, and Yonina C Eldar. convergent neural network for non-blind image deblurring. In 2023 IEEE International Conference on Image Processing (ICIP), pp. 1505-1509 TEFE, 2023.
21

---

Published in Transactions on Machine Learning Research (12/2024)
Jun-Yan Zhu, Taesung Park, Phillip Isola, and Alexei A Efros. Unpaired image-to-image translation using cycle-consistent adversarial networks. In Proceedings of the TEFE international conference on computer vision, pp. 2223 2232, 2017.
Yuanzhi Zhu, Kai Zhang, Jingyun Liang, Jiezhang Cao, Bihan Wen, Radu Timofte, and Luc Van Gool. Denoising diffusion models for plug-and-play image restoration. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pottern Recognition, pp. 2219, 2023.
A Theorems and Proofs
A.1 Required Assumptions
We begin with stating some basic assumptions. These correspond to A1-3 in De Bortoli et al. (0011) )-
Assumption 1. (Compact )omain.) ② is conver compact set and θ C B (0, Re) with Re > 0.
Assumption (Lipschitz gradient.) There crists an open set U Rda and La Z such that O U, logp(ul.) C¹ (U, and for any 01, 02 6
(o V loggp(u S Lo 100 - 02l1.
Assumption 3. For any 0 E O there exist Ho, He Rd + Ro and two probability distributions Ro the Borel space (R B (R") satisfying for any A E θ
f Ho(x)dmo(z) + J Ho(x)ddo(x))
V Velogp(y10)
In addition, (0,x)= + Ho(x) and (0,x) Ha( are measurable.
Assumptions 1 and 2 are standard compactness and Lipschitz assumptions in the literature, and note that Assumption holds if Fisher's identity (6) holds. We additionally require some assumptions guarantee
ergodicity of the two Markov chains in Equation 7 The following assumptions correspond to L1-3 in
De Bortoli et al. (2021).
Assumption 4. (Uniformly Lipschitz gradients.) For any 0 E O, To admits probability density function
with respect to to the Lebesgue measure proportional to - exp (x)) and wo admils probabilily densily function with respect to to the Lebesgue measure proportional to x + exp (-))(z)). In addition, (0, x) + ge and (6, x) + Ex(x) are continuous, x + go(x) and r + Ea(x) are differentiable for all E O and there erists
Lg 2 0, LE 2 O such that for any I1, x2 E Rd
sup V (z)()))) 1100-VVaaae V0 s [, - x2ll
#C
sup
V Ex(1)) Ex LR xall
6E
and g/ : 0), (o : 0 E {} are bounded.
Assumption 5. (Twil condition.) There erists my > 0 and m2, c, Ri 2 0 such that for any 0 θ and eeb Rdd
(V,Eo(z), 2 mi (( + m2 VV Eo())] c.
Assumption 6. (Well-condittineed w.r.t. 0.) There exist constants MI, 2 0 and M2,e 2 such that for any B2 E O, Rd
V. - (x)) < M1 .h.e/li ||0 - 021 exp mi / 1 ||
VV.90, (x) - Vago, < M2 ||0, 021 exp
mi / + ||x11//8)
22

---

Published in Transactions on Machine Learning Research (12/2024)
Assumption is necessary condition for the SDEs corresponding to the ULA Algorithms 1 and 2 to have unique solutions. Assumptions 5 and G provide conditions for geometric ergodicity of ULA, which asserts geometric convergence to the invariant distribution in total variation (Meyn & Tweedie 1992]. Other sufficient
conditions for geometric ergodicity of ULA can be found in Durmus & Moulines (2017). We note that the exponential term in Assumption 6 is not necessary, and we show a lighter bound growing linearly |||]
similar to that used in (De Bortoli et al. 2020).
We provide here some intuition into how each condition is satisfied by the convex ridge regularizer. Detailed proofs are left to the appendix.
Proposition 1. Suppose lakes the form of a convext ridge regalarizer 14 where the conver profile functions ii are C¹, parameterized using piecewise quadratic splines, and the parameters I takes values in some compact O. Suppose further that fy(x) conver and also has Lipschitz gradient, and E(ylæ)) = exp(- has finite second moment. Then Assumptions to o all hold.
Proof sketch. Assumption optimization.
7
The compact O can be chosen a-priori to be sufficiently large to allow for
A Asumption Recalling the form of the CRR, go(x) x E (w, x), this is Lipschitz on compact sets by chain rule. The dominated convergence theorem together with Fisher's identity transfers this to the log marginal Vep(yl0). See Appendix B.1. which writes the derivative of the CRR in closed form.
A Assumption Follows from Fisher's identity (6. The required conditions for Fisher's identity follow (Vidal
et al. 2020 Prop. A.I.), verified by the CRR:
1. For any 0 E θ and & E Rd, (z, + p(x, vl0), y + p(ylz) are positive p.d.f.s
Follows since fy, 90 are finite and convex, and that go is constructed to grow quadratically for large x.
2. For any E 1R 0 E in)( 0 + p(v,=10) is differentiable.
Follows since CRRs are continuously differentiable w.r.t. for fixed x.
3. For any y E Rd and 0 E iO), there exists >> and measurable ğ such that for any Ā e (e) and Rd V log x 0) S Ø(x) with f ((=)p(vlr)de < +o0.
Since θ is compact, V/Vo ply, xz00ح| = [z)]] can be bounded by e||t:|| as in Appendix B.1 which is integrable under p(ylz) by assumption of having finite second moment.
We have ny = p()00), मe = p(rly, 0), and Ho = Ho = Vogo- This is also given in De Bortoli et al. 2020 Remark 2).
Assumption Since θ is compact and ge are piecewise quadratic, the ridge regularizers have uniformly Lipschitz gradient (Goujon et al. 2022 Prop. IV.1.). The condition for Eg holds as well since (x) has Lipschitz gradient. See Lemma
A ssumption For simplicity, we can enforce the convex ridge regularizer to strongly convex by a choice of spline parameters. In this case, Eg = fy + ge is also strongly convex, and thus for sufficiently large x, is bounded below by <1 ||0|]. Apply Young's inequality and taking sufficiently large c, R1 for the desired tail condition. In general, we need only strong monotonicity of the gradient for sufficiently large x w.r.t. the origin. See Lemma 2
Assumption 0 By continuity of the convex ridge regularizer, and that both go and Es grow at most quadratically in ||a]]] See Lemma 3
This rest of this section is divided into three subsections for clarity, addressing Assumptions 4-6. Recall the definition of the CRR and the form of its gradient:
go x +
sha x), Vgo(x)= = WTo(Wx))
(19)
23

---

Published in Transactions on Machine Learning Research (12/2024)
where di IR + R convex profile functions (ridges) applied clement-wise, and w; E R are learnable weights, W = [W1-. ... we]" E RCxd is linear, and : IR IRC is component-wie activation function with components o, = v/- R + R.
A.2 Assumption 4 for CRRs
Assumption 4. (Uniformly Lipschitz gradients.) For any 0 E θ To admits probability density function with respect to to the Lebesque measure proportional to x + exp (x)) and Tg admits a probability densily function with respect to to the Lebesgue measure proportional to x + exp (-)(())) In addition, (0, x.) + go and (0, x) + Ex(x) are continuous, I + go(x) and I + Es(ø) are differentiable for all E and there exists
Lg 2 0, L 2 0 such that for any 11. 22 E Rd
sup V (zi) (xg0(12) Lgllar x2
AE
sup V - s Le|la - x2
6E
and ([ 0 E θ}, (IV : 0 Q} are bounded.
Lemma 1. Conver ridge regularizers Equation 14 satisfy Assumption
Proof. From Equation and Appendix B since @ is compact and bounded, the slope of a is uniformly bounded and hence uniformly Lipschitz, say with constant La 0. Moreover, [|W]] is uniformly bounded (in operator norm). Given 21, 22 E Rd.
V= & W (TTTT (W2))) W WoWo)))) s [ L, [ |2] x2ll
Let Lg Lo supace (IWII*). Then Lg +oo and complete the inequalities with LE = Ly + Lg, where both V.9o and thus Eg are uniformly Lipschitz with respect to I, over E O. The rest follows from compactness of and the continuous parameterization of the CRRss
A.3 Assumption 5 for CRRs
A Assumption 5. (Tuil condition.) There erists my > O and m2, c, Rı 2 0 such that for any E θ and
x Rd
(V=Eo(z), 2 m + m [
Lemma 2. Conver ridge regularizers Equation
14
satisfy Assumption
Proof. Convex ridge regularizers are strongly convex, which can be enforced by choice of spline parameters, see Goujon et al. (2022) Appendix B Therefore, given 0 E O let Eg have strong convexity constant.
ma > m > 0 bounded below by m. Given x, E RRd Eg satisfies
(x VE(r)) - VEo(z') 2 mo || x'112.
(20)
From Lemma v (a(-) is uniformly Læ LipLichciitz over 0 E O. For all x, E Rd, we have the following inequalities (Nestervv 2013
Thm. 2.1.5):
(x x'. V...E. (x) V Eo('))) < Lg || - 2'112,
(21)
(e)(2) Eo(''') < Lg (x x', Ee(x)) - VEo(x'))
(22)
We compute as follows:
|V
parallelogram law
< 2Lg (I - -E,(2)))))((())))))))))))))))))) (
by 22
24

---

Published in Transactions on Machine Learning Research (12/2024)
Rearranging, we have
-x- Ba( - ( 2 2
LE
We combine this with 20 by halving both to obtain
( x' VE VEo('')) 2 mg
x'11
|VEo(>)]
[VE6(2')11
(23)
4Lg
2L
Now let x' 0. For any 2, we have the following by rearranging.
(r. VFo(x) 1 (O) 2
me
|x|
[V
|VE2(0))
2
4LE
2LE
(z, VEo()))
2
mg
||0||
[x)]
[O]]
+ (z. VE₂
2
4Lg
2Lg
me
|
1
[V
1
(24)
2
4LE
2LF
where the final inequality follows from Cauchy-Schwarz.
For ||||||| & 4 (s me, we have from 24 that
4
1
(x, VEo(x)) 2
[Eo(2))
|VE2(
ALE
mg
2LE
For | 2 (l]VEa(0))//mm, mo, we have
(x, VE+(T)) 2 | [=]]V 1
1
|V
m
4Lg
2L
1 |VE2(>))1
1
E (O)Ill)) +
[VE(
4Lg
2L
Thus we have
1 [V
(4+1 mg
(x, VEo())) 2
+
[O)
(25)
ALE
2Lg
We conclude by noting that my > m > lower-bounded and that the inequalities continue to hold the ball radius is taken to be R supe 4 [O)]l which is finite as is compact.
A.4 Assumption 6 for CRRs
Assumption 6. /Well-condilioned w.r.l. 0.) There exist constants M1,8 2 0 and M2 2 0 such that, for any 0, O2 E 0, E Rdd
V Eo, ( (x), < M1,e 1|0 021 exp mi + 1|/11/88
| >90, (x:) Vzgo, (x)) < M2 ||0 02l
exp m= + | /8)
Lemma Conver ridge regularizers Equation
14
satisfy Assumption
Proof. From Equation 15 since ④ compact and bounded, the slope of bounded and a is uniformly
Lipschitz, say with Lipschitz constant bounded by T.g.
25

---

Published in Transactions on Machine Learning Research (12/2024)
For any 81, 02 E O and R assume that W1 is part of 01 and W2 is a part of 02, then we have
|| Vago, (x) Vagu, (z) <
||W] 01 ( (W1x) w 01 + o[
W2 + || o2(2x)) - (W22)]
S Lo ( | W W2ll |||||| 2L, ||]]||] ||||||
La [W W2 W2 |x||
< M2 |10 - 02 |||||
where
M2 Lo sup (| Wall + [22 + |[ Well < +00. 0,00
(26)
Conclude the second inequality by using Assumption and noting 1 |||||| < C Cexp I ( for some sufficiently large C. The first incquality comes from fy having Lipschitz gradient.
]
B Convex Ridge Regularizer Implementation
In this section, detail the convex ridge regularizer architecture, and provide closed-form derivative with respect to its parameters in Appendix B.1 Recall from Section 4 that the convex ridge regularizer takes the form
R: Σ the (w. x).
(27)
where W = [W1 wcl E RCxd with learnable weights w, E Id, and v : R + are convex profile functions (ridges). In practice, we use matrices for w,, so that w. a represents a vector, and apply taa component-wise This can be justified by using duplicate vis and weight-tying.
Choice of W. As detailed in Section VI.A in Goujon et al. 2022 W is taken as the composition of two zero-padded convolutions with window size 7 x 7, with 8 and 32 output channels respectively. While the experiments in Goujon et al. (2022) are originally in black-and-white, we extend to color images by using 3 input channels instead of 1 for the first convolution.
Choice of ridges v.. As Goujon et al. (022)), we parametorize each component's derivative o; = v/. R + R as a linear spline, and subsequently define the convex profile function as an integral (),(1) = f (s( ds. To parameterize the linear splines o,, we define the values at equispaced knots, clamping during training to maintain strictly increasing slopes, ensuring strong convexity of the ridge vù; We use 21 equally distant knots (centered at 0), with distance A 0.01 between them. Values defined between knots are given by linear interpolation, and values past the endpoints given by linear extension. In particular, we have ox(t) c-.at as t too and (t) ct as t + -oo for some Ct 0.
By clamping the allowable knots, we get that for compact parameter set O, are uniformly Lipschitz, used in Lemma 1
B.1 Gradient of the CRR
We justify Assumption 2 and the third claim for Assumption 3 in Proposition I by computing the derivative of go with respect to its parameters 0.
Fix a CRR architecture, where the linear spline o, = v/f has (fixed) knots at t-k...., tk, taking (learnable) values c- ........ CK- We assume that the knots are equispaced fę = KA for h =K K,.., K, that co 0, relaxing these assumptions are simple and do not significantly change the proof. For strong convexily, we require that the sequence CK is strictly increasing. We can thus reparameterize in terms of spline differences k..... dx that are lower bounded by some constant m > 0 to enforce strong convexity:
dj
k 0
(1(14) = C =
k
d3
k < 0,
de 2 m > for k: - K,.,,,,, K.
(28)
0
k
26

---

Published in Transactions on Machine Learning Research (12/2024)
As use linear interpolation between the knots (and extrapolation outside), this can be rewritten sum of ReLUs:
Σ== (dj - dj- max (0, 1 (111)/
I E [t,, 4+11;; 1); h: 2
0i(1)=
Σ b d max t0, -
t ( tkl, S
k=0
Σ
و j=1 (dj d max t
K
t20
j=-K K (dg - dg+1) max - t)/
We can then compute (t( Sh o,(t). We note that since o; is sum of ReLUs (and vice versa in the negative direction), the integral (t) is sum of rectified quadratic units ReLU (t) f f1120:
Σ k (dj dj- ReL
t
(b(() -
(29)
ンj=-k (d, - ) R(T f)/A t
Given y Rd,, fix any ò e in) and r e Rd The parameters 0 consist of the weights w, and spline values C-K,b-..,K corresponding to spline o, o; viy/ In the following, we drop the index for notational convenience. We can compute derivatives with respect to spline parameters as follows. Suppose w x E [&r, (t41,,), where we define (x+1 = +00. We can differentiate with respect the spline parameters using 29
1 (ReLU?(wTT (1-1)- (e..U(( -1)]- A
(o/0a,)a(('T
(30)
Similar computations hold instcad w Tx O, and define f- -0. Note that these bounds also hold uniformly in small ball around 0, as w is continuous in w. Moreover, this grows linearly in w WT.. We can also differentiate with respect to the weight vector and get sum of ReLU units, which again grows linearly in æ: assuming w x is positive,
V Vwb(ww) 2)
Σ dg ReLi) tz-1)x
(31)
A similar expression holds if w T is negative. This expression is asymptotically quadratic in ||=||-
To verify Assumptions 3 and 6 in Proposition 1 we need only consider the derivatives (3 and (31 We have that ((Veg((>) is asymptotically bounded by a quadratic in ||ж]],, as required for Fisher's identity given finite second moment. Moreover, since the derivatives are sums and products of ReLU and squared ReLU units, Lipschitz continuity of Vg log ply(0) compact θ comes from the dominated convergence theorem.
C Time comparison
For TV, we report the time required to reconstruct 50 images using 5000 iterations of gradient descent for a fixed regularization parameter. The parameter is chosen via manual grid search to maximize PSNR.
For DIP, we compared 4 different training settings (2 architectures and initializations) and reported the best metrics for the Gaussian doconvolution and Poisson denoising experiments. We report the average time required to reconstruct using DIP by dividing the total time to reconstruct 50 images using all 4 methods by 4.
27

---

Published in Transactions on Machine Learning Research (12/2024)
Table 4: Approximate time required to train and test the compared methods Gaussian deconvolution and Poisson denoising tasks. Testing is done on 50 images, unrolled for 5000 iterations for energy-based methods such CRR-based and TV.
baseb
End-to-end/ seltto-enddmmeel-bbeee modeel-baed
SAPG
DIP
GS
E1
TV
MAP
MMSE (2 x 10")
MMSE (1 x 10")
Adam
GD
Train (Gaussian)
10 mins
days
hours
mins
days
Train (Poisson)
45 mins
week
day secs
15 mins
21 hours
Test
secs
mins
11 mins
52 mins
8 secs
28


Reviews:
[{"review_id": "XSe1ryIaTW", "reviewer": "Reviewer_trdP", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n\n- *Clear presentation*: the writing is clear; the motivation and literature is thoroughly grounded in the inverse problems literature.\n- *Theoretical results*: I did not check the proofs in detail, but their statements and conditions are clear.\n- *Aims*: expressive parameterizations seems to be key to reach state-of-the-art performance, hence the associated improvements to scalability are appreciated.\n\n**Weaknesses**\n\n*Major weaknesses*\n\nI see two related weaknesses due to which I don't think the current submission is \"supported by accurate, convincing and clear evidence\".\n1. *Diffusion models are notably absent.* \n    It is by now well-established that diffusion models can provide strong priors for inverse problems [1-7], but there is no discussion or comparison with this class of methods. First, this makes it hard to position the paper. Second, I believe that many state-of-the-art models for inverse problems are in fact based on diffusions, hence the lack of such baselines significantly weakens the performance claims of the proposed method. Most relevant would be unsupervised baselines such as [8] and [9], but to support the comparison to supervised methods in claim 1 there are also many supervised methods like [10].   \n\n2. *The experiments are only weakly aligned with claims.* \n\n    - Regarding claim 1 above, I might have misunderstood this, but I don't think the original work is included as a baseline, which would the obvious way to provide evidence for the claim.\n    - Regarding claim 2, there is no experimental evidence supporting the claims of improved scalability due to the batched training.\n    - Regarding claim 3, the authors acknowledge that batched SAPG-ULA is biased. Since this bias isn't analyzed theoretically, I would have liked to see an experimental evaluation of its size and effect.\n    - I think the experiment in section 4.3 is of some interest, but it has no bearing on the claims and, besides, would probably merit a different approach. \n\n*Minor weaknesses*\n\n- Only using PSNR is a bit meagre. Computer vision papers usually include also a perceptual metric like FID or LPIPS. In addition, I think it could make sense to compare the marginal likelihood as well in this case.\n- I think the paper would benefit from being shortened. Some suggestions:\n    - shorten each paragraph of the literature review, but add one on diffusion models\n    - move sections 2.4 and 3.1 to appendix\n\n*Questions:*\n1. Would it be possible to generalize the class of regularizers to Input convex neural networks [11]?\n2. Would it be possible to use variance reduction techniques like common random numbers?\n\n*References*\n\n[1] Feng et al., Score-Based Diffusion Models as Principled Priors for Inverse Imaging. ICCV, 2023.\n\n[2] Graikos et al., Diffusion Models as Plug-and-Play Priors. NeurIPS, 2022.\n\n[3] Cardoso et al., Monte Carlo guided Diffusion for Bayesian linear inverse problems. ICLR, 2024\n\n[4] Mardani et al., A Variational Perspective on Solving Inverse Problems with Diffusion Models. ICLR, 2024.\n\n[5] Dou and Song, Diffusion posterior sampling for linear inverse problem solving: A filtering perspective. ICLR, 2024.\n\n[6] Luo et al., \"Taming diffusion models for image restoration: A review.\" arXiv preprint 2409.10353, 2024.\n\n[7] Zhao et al. Conditional sampling within generative diffusion models. arXiv preprint 2409.09650, 2024.\n\n[8] Kawar et al., Denoising diffusion restoration models. NeurIPS, 2022.\n\n[9] Wang et al., Zero-shot image restoration using denoising diffusion null-space model. ICLR, 2023.\n\n[10] Luo et al., 2023. Image restoration with mean-reverting stochastic differential equations. ICML, 2023.\n\n[11] Amos et al., Input convex neural networks. ICML, 2017.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The papers extends previous work for unsupervised learning of convex regularizers for inverse problems. In my view. the main claims are as follows:\n1. *Performance*: it performs better than the closest prior work, and almost as well as supervised baselines.\n2. *Scalability*: it can be trained in a higher-dimensional setting than the closest prior work.\n3. *Theoretically convergent*: given a particular choice of parametrization, the authors verify that the convergence conditions in previous works apply.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- Relate the work to diffusion models\n- Add diffusion-based baselines to the experiments\n- Provide additional evidence for the claims, as suggested above.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "oWoksppY4N", "reviewer": "Reviewer_bTr2", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\nThe paper is well-written and clearly presents the problem statement, relevant background, and literature review. Also, the explanations of the proposed method and theoretical foundations are easy to follow. By incorporating a more general and expressive regularizer architecture, the proposed method extends traditional approaches like total variation, making it adaptive to different applications. The paper addresses various aspects by including a theoretical convergence proof for the proposed method, as well as proposing a scalable and efficient approach that handles large datasets using mini-batching and MCMC techniques.\n\nWeaknesses:\n The experiments are somewhat limited in terms of dataset variety and evaluation metrics. Only PSNR has been considered as an evaluation metric, while other metrics such as SSIM (structural similarity index measure) could provide a more comprehensive assessment of recovered image quality. Also, while DIP (Deep Image Prior) is used as a baseline, it is relatively old compared to recent state-of-the-art methods. Considering more recent methods, especially diffusion-based models, would provide a better comparison for the effectiveness of the proposed approach. For example, the paper misses including recent works such as: \n\n- Feng, Berthy T., et al. (2023): \"Score-based diffusion models as principled priors for inverse imaging.\"\n- Chung, Hyungjin, et al. (2023): \"Diffusion Posterior Sampling for General Noisy Inverse Problems.\"\n- Zhu, Yuanzhi, et al. (2023): \"Denoising diffusion models for plug-and-play image restoration.\"\n\nMoreover, the method builds upon prior work by Vidal et al., but the original work from Vidal has not been considered as a baseline. It is important to compare against this to show the improvements clearly. The experiments focus primarily on Gaussian deconvolution and Poisson denoising. Including more diverse inverse problems (such as MRI reconstruction or super-resolution) would be interesting to show the broader applicability of the proposed method. Regarding the convergence proof, there is a question of novelty in this aspect. Similar proofs for convergence exist in the literature for algorithms that involve stochastic gradient descent (SGD) or stochastic optimization methods that combine gradient updates with Markov Chain Monte Carlo (MCMC) techniques. The proof of convergence in this paper also builds on existing work, and the authors specifically cite De Bortoli et al. (2021) as the foundation for their proof, which provides convergence guarantees for Langevin algorithms that use gradient-based methods combined with MCMC for sampling from complex distributions.\nAdditionally, the method itself is not entirely new, as mini-batching for large datasets is a well-established approach. So, while the paper is well-written, the overall novelty appears to be limited, and further clarification on what sets this work apart from existing methods is required. Moreover, the experiments could be extended by incorporating more diverse datasets, evaluation metrics, and stronger baselines to show the method’s effectiveness better.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper addresses the inverse problem of reconstructing images from noisy measurements. By leveraging maximum likelihood estimation (MLE), the authors propose an unsupervised method to train convex neural network-based image regularizers using only noisy measurements without requiring ground-truth clean data. This way, they attempt to extend existing unsupervised methods with a more expressive model for the regularizer. In terms of efficiency and scalability, the paper introduces a scalable method that extends the stochastic approximation proximal gradient (SAPG) algorithm that handles large datasets via mini-batching and efficient Markov Chain Monte Carlo (MCMC) sampling to overcome computational challenges. The authors also provide theoretical convergence guarantees for their method. The paper provides proof of ergodic convergence for the proposed SAPG algorithm using a convex ridge regularizer architecture. Their experimental results show that the proposed method performs comparably to state-of-the-art supervised approaches across two image corruption tasks, Gaussian deconvolution and Poisson denoising. They also showed better generalization compared to end-to-end supervised methods.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Critical:\n\n- Incorporating additional evaluation metrics\n- Including recent state-of-the-art methods for comparison (diffusion-based models)\n- Using Vidal et al. as a baseline\n- Clarifying the novelty in convergence proof\n- Clarifying novelty in mini-batching approach\n\nStrengthening:\n\n- Expanding the diversity of datasets used in experiments\n- Expanding the diversity of inverse problems (e.g., MRI reconstruction, super-resolution)", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "I3ISiu8RJ9", "reviewer": "Reviewer_HnW1", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths:**\n- The work tackles a timely and important problem, namely addressing how to learn data-driven regularizers in situations when one only has access to noisy measurements.\n- The paper is well-written and easy to read, while also providing a nice background on learned variational regularization and its connection to bayesian methods.\n- Under the assumptions of parameterizing $g_{\\theta}$ as a CRR, the authors can obtain nice theoretical guarantees on the asymptotic convergence in expectation of their sampling procedure. The result is mostly an application of previous work, but some work is necessary to show that the CRR indeed satisfies the assumptions necessary for the Theorem.\n- The performance is not far behind comparable supervised methods, with nice unforeseen properties such as stronger generalization properties with respect to the forward operator.\n\n**Weaknesses:**\n- Through the use of two MCMC methods to optimize $p(y |\\theta)$, this method can be potentially very slow. The theoretical results are also asymptotic in nature, which do not necessarily give guarantees in terms of how long it would take to obtain a high-quality sample. Based on my reading, there is no discussion on how long it takes to train the regularizer, which I think is important.\n- Focusing on convexity, while certainly nice in terms of the induced posterior's theoretical properties, does come with insurmountable limitations. In particular, in ill-posed problems where the posterior can be multi-modal (e.g., when the forward model is nonlinear), such convex regularizers can never capture this structure since log-concave distributions are necessarily uni-modal.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work considers solving inverse problems with learned regularization in situations where one does not have access to samples of ground-truth images and only access to noisy measurements. The authors consider parameterizing a regularizer via a neural network-based convex ridge regularizer (CRR) $g_{\\theta}$ and use an empirical Bayesian approach to maximize $p(y | \\theta)$, where $y$ constitute our measurements. To overcome the intractability of computing $\\nabla_{\\theta} p(y|\\theta)$, the authors show that gradient ascent on $p(y|\\theta)$ amounts to computing expectations of the gradient of $g_{\\theta}$ over the prior and posterior. To compute empirical averages of such expectations, an additional MCMC loop is applied. Once trained, the regularizer can then be used to compute an MAP or MMSE estimate for downstream inverse problems. Theoretically, the authors show that, when using a CRR, the bias of their sampling approach is asymptotically bounded by the step size parameters. They experimentally verify their approach on Gaussian deconvolution and Poisson denoising, showing that they obtain solid performance when learning only on noisy measurement data.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "**Requested changes:**\nOverall, I am positive on the paper and think that this is a solid submission to TMLR. The only requested change I have would be to include further discussions on limitations (e.g., runtime of the training algorithm and focusing on convexity), as this was not thoroughly discussed in the main paper or in the appendix.\n\n**Comments/questions:**\n- One of the main questions I have is in regards to runtime and how long it takes to train such a regularizer. How does this compare to, e.g., the adversarial regularization framework of Lunz et al? A more in-depth discussion about this would be nice and improve the quality of the paper.\n- Under the assumptions on the CRR, you effectively have a strongly log-concave density. Can the authors comment on whether one can obtain stronger sampling guarantees (e.g., non-asymptotic ones) with this additional structure?\n- One should note the structure of the noisy images $x_i^{\\mathrm{noisy}}$ in equation (18). Are these simply backprojected images, as in Lunz et al, or are these clean images contaminated with noise?\n- It is noted that DIP has a lack of theoretical guarantees. While I agree in general it is true one cannot achieve convergence, it could be good to reference a few papers in the intro that have addressed this in simpler settings (using NTK-based techniques), e.g.,\n\n\"Denoising and Regularization via Exploiting the Structural Bias of Convolutional Generators\", R. Heckel & M. Soltanolkotabi, ICLR 2020\n\n\n**Typos:**\n- Should be a parentheses around $E_{\\theta} + \\iota_{+}$.\n- The notation for Borel sets $\\mathcal{B}(\\mathbb{R}^d)$ should be defined.\n- In assumption 5, does the notation refer to the complement of the ball $B(0,R_1)$? I would potentially choose a different notation for the constant $c > 0$.\n- Top of page 12: the weight matrix is not bolded.\n- There is an exclamation mark in equations describing the Poisson denoising problem.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper considers inverse imaging problems, where ground-truth data is unavailable, and only noisy measurements are provided. The authors propose an unsupervised learning approach, where they train a convex neural network-based image regularizer, using maximum marginal likelihood estimation. They extend the stochastic approximation proximal gradient (SAPG) algorithm of (Vidal et al., 2020), from single measurements to datasets, leveraging mini-batching and efficient Markov Chain Monte Carlo (MCMC) sampling to overcome computational challenges. 
Theoretical convergence results of the proposed method are provided when using the convex ridge regularizer (CRR).  Experiments demonstrate that their method achieves performance comparable to analogous supervised models for image reconstruction tasks such as Gaussian deconvolution and Poisson denoising, while offering better generalization than end-to-end models.

Below is a summary of the main strengths and weaknesses highlighted by reviewers and based on my own assessment. 

Strengths:
- The problem considered is important; namely learning data-driven regularizers from only noisy measurements 
- The paper is well written and clear. It provides a clear motivation, background, and explanation of the method and its theoretical results.
- Nice theoretical guarantees on the asymptotic convergence when using CRR. The result follows from prior work, but requires some work to verify that the assumptions necessary hold.
- Proposed method improves over existing MLE-based methods in both model expressiveness and dataset size
- Experimental results demonstrate competitive performance with supervised methods while showing better generalization.

Weaknesses:
- Proposed method is slow. No theoretical bound on the time required to obtain good sample. 
- Focus on convex regularizers is limiting. Though it does allow for nice theoretical properties. 

The official recommendations from the reviewers were mixed, with two reviewers recommending to reject and one to accept. The authors have addressed in their responses most of the reviewers' concerns. A key concern raised by two reviewers is the lack of comparison with diffusion-based methods. The authors added a discussion noting that existing diffusion models require either clean data or multiple noisy measurements, which is different from the setting here where only single copies of noisy data are available. Given this difference in setting, I agree with Reviewer HnW1 that a comparison with diffusion models, while interesting, is not necessary. The other concern raised by Reviewer bTr2 is about the novelty of the work, which is not a necessary criteria for acceptance at TMLR.  Hence, I am recommending to accept this paper with minor revision.

Requested revisions: 
- In the paragraph added on p.10 about empirical studies of bias of SAPG, add the conclusion from these studies, i.e., this comment from your response: "The papers find that the bias is not significant if the hyperparameter guidelines are followed, as we do in our work."
- Add the motivation for the experiments in Section 4.3 mentioned in the response to Reviewer trdP
- Add a discussion of the limitations of the proposed method (runtime, convexity).
Decision: Accept (poster)