summaries: This paper proposes a scalable method for conformal Bayesian prediction, by estimating the modified posterior predictive density using "add-one-in" importance sampling.  On several tabular datasets with sparse and hierarchical models, estimation is shown to be efficient for non-extreme miscoverage levels, and the produced credible intervals have good coverage.  **Post-rebuttal update:** Thank you for the response.  I am keeping my score unchanged: on the plus side, the method is sensible and easy to implement; on the other hand, there are some (understandable) limitations when scaling to large datasets or applied on approximate posterior samples; and some reviewers raised questions about novelty, which I'm not most suited to judge, as I'm not an expert in this field.  In aggregate, the rating of "marginally above the acceptance threshold" still describes my opinions well.

summaries: The authors present an approach to using a Bayesian model to provide predictive intervals with  frequentist coverage within the framework of conformal prediction.  This work generalizes existing Bayesian methods for conformal inference that apply only to conjugate models.  The authors provide a straightforward generalization of their approach to partially exchangeable data.  Empirical validation on several datasets, and comparison to a baseline approach is provided.

summaries: This article proposes a novel approach to constructing prediction intervals with approximately correct frequentist coverage, even on misspecified models.  The basic idea is to employ the conformal prediction technique using the Bayesian posterior predictive density as the conformity (goodness-of-fit) measure.  The main innovation of the article is to use an importance sampling approximation to the posterior predictive density at candidate values of outcome $y$.  This leads to a computationally fast and generally applicable algorithm for using samples from the standard Bayesian posterior to construct well-calibrated prediction intervals, even under misspecification of the Bayesian model.  An extension to hierachical models is also provided, and empirical results are presented demonstrating the performance of the method.

summaries: This paper introduces Bayesian and importance sampling techniques into a previously known conformal prediction framework.  I had troubles identifying what exactly is the novel contribution of this paper. It could use a brief section stating "our own contributions are A, B , C". In my understanding the conformal prediction framework was known before and the novel items are using the Bayesian posterior predictive density as a conformity score and additionally estimating the predictive distribution given new data from a weighted importance sampling where the weights are just the likelihoods of the prameter draws based on the new data point.  There are modelling extensions from completely exchangeable models to partial exchangeability, i.e., where we have groups and complete exchangeability is still valid in the groups, but not outside the groups. From this flat hierarchical structure there is a further extension to hierarchical models that may model deeper hierarchies and via hyperpriors may share information between groups.

summaries: The authors develop a fast conformal Bayesian computation approach.  The original conformal prediction method can be computationally intensive:  It often requires refitting the model for every augmented value y. The authors avoid this refitting step by using importance sampling. The authors also show the extension to conditionally exchangeable data, which is especially useful for hierarchical data. The resulted conformal interval has finite sample correct coverage and is immune to model misspecification or distribution assumptions.

abstractText: We develop scalable methods for producing conformal Bayesian predictive intervals with finite sample calibration guarantees. Bayesian posterior predictive distributions, p(y | x), characterize subjective beliefs on outcomes of interest, y, conditional on predictors, x. Bayesian prediction is well-calibrated when the model is true, but the predictive intervals may exhibit poor empirical coverage when the model is misspecified, under the so calledM-open perspective. In contrast, conformal inference provides finite sample frequentist guarantees on predictive confidence intervals without the requirement of model fidelity. Using ‘add-one-in’ importance sampling, we show that conformal Bayesian predictive intervals are efficiently obtained from re-weighted posterior samples of model parameters. Our approach contrasts with existing conformal methods that require expensive refitting of models or data-splitting to achieve computational efficiency. We demonstrate the utility on a range of examples including extensions to partially exchangeable settings such as hierarchical models.

authors:
R.F. Barber
E.J. Candes
A. Ramdas
Tibshirani
R. J

title: Predictive inference with the jackknife+

venue: Annals of Statistics,

title: The concept of exchangeability and its applications

venue: Far East Journal of Mathematical Sciences, 4:111–122.

authors:
J.M. Bernardo
A.F. Smith

title: Bayesian theory, volume 405

venue: John Wiley & Sons.

authors:
J. Bradbury
R. Frostig
P. Hawkins
M.J. Johnson
C. Leary
D. Maclaurin
S. WandermanMilne

title: JAX: composable transformations of Python+NumPy programs

authors:
E. Burnaev
V. Vovk

title: Efficiency of conformalized ridge regression

venue: Conference on Learning Theory, pages 605–622.

authors:
E.J. Candès
L. Lei
Z. Ren

title: Conformalized survival analysis

venue: arXiv preprint arXiv:2103.09763.

authors:
B. Carpenter
A. Gelman
M.D. Hoffman
D. Lee
B. Goodrich
M. Betancourt
M.A. Brubaker
J. Guo
P. Li
A. Riddell

title: Stan: a probabilistic programming language

venue: Grantee Submission, 76(1):1–32.

authors:
W. Chen
Chun
K.-J.
R.F. Barber

title: Discretized conformal prediction for efficient distribution-free inference

authors:
W. Chen
Z. Wang
W. Ha
R.F. Barber

title: Trimmed conformal prediction for highdimensional models

venue: arXiv preprint arXiv:1611.09933.

title: The well-calibrated Bayesian

venue: Journal of the American Statistical Association, 77(379):605–610.

authors:
R. Dunn
L. Wasserman
A. Ramdas

title: Distribution-free prediction sets with random effects

venue: arXiv preprint arXiv:1809.07441v2.

authors:
B. Efron
T. Hastie
I. Johnstone
R Tibshirani

title: Least angle regression

venue: Annals of statistics,

authors:
E. Fong
C. Holmes

title: On the marginal likelihood and cross-validation

authors:
Fraser
D. A

title: Is Bayes posterior just quick and dirty confidence

venue: Statistical Science,

authors:
A. Gammerman
V. Vovk
V. Vapnik

title: Learning by transduction

venue: Proceedings of the Fourteenth conference on Uncertainty in artificial intelligence, pages 148–155.

authors:
A. Gelman
J.B. Carlin
H.S. Stern
D.B. Dunson
A. Vehtari
D.B. Rubin

title: Bayesian data analysis

venue: CRC press.

authors:
A. Gelman
J. Hill

title: Data analysis using regression and multilevel/hierarchical models

venue: Cambridge university press.

authors:
T. Gneiting
A.E. Raftery

title: Strictly proper scoring rules, prediction, and estimation

venue: Journal of the American Statistical Association, 102(477):359–378.

title: Fast exact conformalization of the lasso using piecewise linear homotopy

authors:
J. Lei
M. G’Sell
A. Rinaldo
R.J. Tibshirani
L. Wasserman

title: Distribution-free predictive inference for regression

authors:
J. Lei
L. Wasserman

title: Distribution-free prediction bands for non-parametric regression

venue: Journal of the Royal Statistical Society: Series B: Statistical Methodology, pages 71–96.

title: Calibrated Bayes: a Bayes/frequentist roadmap

venue: The American Statistician, 60(3):213–223.

authors:
M. Magnusson
M. Andersen
J. Jonasson
A. Vehtari

title: Bayesian leave-one-out crossvalidation for large data

venue: International Conference on Machine Learning, pages 4244–4253. PMLR.

authors:
T. Melluish
C. Saunders
I. Nouretdinov
V. Vovk

title: Comparing the Bayes and typicalness frameworks

venue: European Conference on Machine Learning, pages 360–371. Springer.

title: Monte carlo theory, methods and examples

authors:
F. Pedregosa
G. Varoquaux
A. Gramfort
V. Michel
B. Thirion
O. Grisel
M. Blondel
P. Prettenhofer
R. Weiss
V. Dubourg
J. Vanderplas
A. Passos
D. Cournapeau
M. Brucher
M. Perrot
E. Duchesnay

title: Scikit-learn: Machine learning in Python

venue: Journal of Machine Learning Research, 12:2825–2830.

authors:
Y. Romano
E. Patterson
E.J. Candès

title: Conformalized quantile regression

venue: arXiv preprint arXiv:1905.03222.

authors:
J. Salvatier
T.V. Wiecki
C. Fonnesbeck

title: Probabilistic programming in Python using PyMC3

venue: PeerJ Computer Science, 2:e55.

authors:
S. Seabold
J. Perktold

title: Statsmodels: Econometric and statistical modeling with python

venue: Proceedings of the 9th Python in Science Conference, volume 57, page 61. Austin, TX.

authors:
G. Shafer
V. Vovk

title: A tutorial on conformal prediction

venue: Journal of Machine Learning Research, 9(Mar):371–421.

authors:
R. Tibshirani
R. Foygel

title: Conformal prediction under covariate shift

venue: Advances in neural information processing systems.

authors:
A. Vehtari
A. Gelman
J. Gabry

title: Practical Bayesian model evaluation using leave-oneout cross-validation and WAIC

venue: Statistics and computing, 27(5):1413–1432.

authors:
A. Vehtari
D. Simpson
A. Gelman
Y. Yao
J. Gabry

title: Pareto smoothed importance sampling

venue: arXiv preprint arXiv:1507.02646.

title: Cross-conformal predictors

venue: Annals of Mathematics and Artificial Intelligence, 74(1):9–28.

authors:
V. Vovk
A. Gammerman
G. Shafer

title: Algorithmic learning in a random world

venue: Springer Science & Business Media.

title: Frasian inference

venue: Statistical Science, pages 322–325.

authors:
W.H. Wolberg
O.L. Mangasarian

title: Multisurface method of pattern separation for medical diagnosis applied to breast cytology

venue: Proceedings of the national academy of sciences, 87(23):9193–9196.

authors:
G. Zeni
M. Fontana
S. Vantini

title: Conformal prediction: a unified review of theory and new challenges

venue: arXiv preprint arXiv:2005.07972.

heading: 1 Introduction

text:
We consider Bayesian prediction using training data Z1:n = {Yi, Xi}i=1:n for an outcome of interest Yi and covariates Xi ∈ Rd. Given a model likelihood fθ(y | x) and prior on parameters, π(θ) for θ ∈ Rp, the posterior predictive distribution for the response at a new Xn+1 = xn+1 takes on the form
p(y | xn+1, Z1:n) = ∫ fθ(y | xn+1)π(θ | Z1:n) dθ , (1)
where π(θ | Z1:n) is the Bayesian posterior. Asymptotically exact samples from the posterior can be obtained through Markov chain Monte Carlo (MCMC) and the above density can be computed through Monte Carlo (MC), or by direct sampling from an approximate model. Given a Bayesian predictive distribution, one can then construct the highest density 100× (1−α)% posterior predictive credible intervals, which are the shortest intervals to contain (1 − α) of the predictive probability. Alternatively, the central 100 × (1 − α)% credible interval can be computed using the α/2 and 1− α/2 quantiles. Posterior predictive distributions condition on the observed Z1:n and represent subjective and coherent beliefs. However, it is well known that model misspecification can lead Bayesian intervals to be poorly calibrated in the frequentist sense (Dawid, 1982; Fraser et al., 2011), that is the long run proportion of the observed data lying in the (1− α) Bayes predictive interval is not necessarily equal to (1− α). This has consequences for the robustness of such approaches and trust in using Bayesian models to aid decisions.
Alternatively, one can seek intervals around a point prediction from the model, ŷ = µ̂(x), that have the correct frequentist coverage of (1 − α) . This is precisely what is offered by the conformal prediction framework of Vovk et al. (2005), which allows the construction of prediction bands with finite sample validity without assumptions on the generative model beyond exchangeability of the
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
data. Formally, for Zi = {Yi, Xi}i=1:n, Zi iid∼ P and miscoverage level α, conformal inference allows us to construct a confidence set Cα(Xn+1) from Z1:n and Xn+1 such that
P(Yn+1 ∈ Cα(Xn+1)) ≥ 1− α. (2) Note that P is over Z1:n+1, that is we also treat Xn+1 as random; see Lei and Wasserman (2014) for a discussion on the difficulties of attaining coverage conditional on Xn+1. In this paper we develop computationally efficient conformal inference methods for Bayesian models including extensions to hierarchical settings. A general theme of our work is that, somewhat counter-intuitively, Bayesian models are well suited for the conformal method.
Conformal inference for calibrating Bayesian models was previously suggested in Melluish et al. (2001), Vovk et al. (2005), Wasserman (2011) and Burnaev and Vovk (2014), where it is referred to as “de-Bayesing”, “frequentizing” and “conformalizing”, but only in the context of conjugate models. Here, we present a scalable MC method for conformal Bayes, implementing full conformal Bayesian prediction using an ‘add-one-in’ importance sampling algorithm. The automated method can construct conformal predictive intervals from any Bayesian model given only samples of model parameter values from the posterior θ ∼ π(θ | Z1:n), up to MC error. Such samples are readily available in most Bayesian analyses from probabilistic programming languages such as Stan (Carpenter et al., 2017) and PyMC3 (Salvatier et al., 2016). We also extend conformal inference to partially exchangeable settings which utilize the important class of Bayesian hierarchical models, and note the connection to Mondrian conformal prediction (Vovk et al., 2005, Chapter 4.5). Previously, the extension of conformal prediction to random effects was introduced in Dunn et al. (2020) in a non-Bayesian setting, with a focus on prediction in new groups, as well as within-group predictions without covariates. We will see that the Bayesian hierarchical model allows for a natural sharing of information between groups for within-group predictions with covariates. We discuss the motivation behind using the Bayesian posterior predictive density as the conformity measure for both the Bayesian and the frequentist, and demonstrate the benefits in a number of examples.

heading: 1.1 Background

text: The conformal inference framework was first introduced by Gammerman et al. (1998), followed by the thorough book of Vovk et al. (2005). Full conformal prediction is computationally expensive, requiring the whole model to be retrained at each test covariate xn+1 and for each value in a reference grid of potential outcomes, e.g. y ∈ R for regression. This makes the task computationally infeasible beyond a few special cases where we can shortcut the evaluation along the outcome reference grid, e.g. ridge regression (Vovk et al., 2005; Burnaev and Vovk, 2014) and Lasso (Lei, 2019). Shrinking the search grid is possible, but still requires many refittings of the model (Chen et al., 2016). The split conformal prediction method (Lei et al., 2018) is a useful alternative method which only requires a single model fit, but increases variability by dividing the data into a training and test set that includes randomness in the choice of the split, and has a tendency for wider intervals. Methods based on cross-validation such as cross-conformal prediction (Vovk, 2015) and the jacknife+ (Barber et al., 2021) lie in between the split and full conformal method in terms of computation and average length of the intervals. A detailed comparison of various conformal methods are provided in Barber et al. (2021, Section 4). A review of recent advances in conformal prediction is given in Zeni et al. (2020), and interesting extensions have been developed by works such as Tibshirani and Foygel (2019); Romano et al. (2019); Candès et al. (2021).

heading: 2 Conformal Bayes

heading: 2.1 Full Conformal Prediction

text:
We begin by summarizing the full conformal prediction algorithm discussed in Vovk et al. (2005); Lei et al. (2018). Firstly, a conformity (goodness-of-fit) measure,
σi := σ(Z1:n+1;Zi),
takes as input a set of data points Z1:n+1, and computes how similar the data point Zi is for i = 1, . . . , n + 1. A typical conformity measure for regression would be the negative squared error arising from a point prediction −{yi − µ̂(xi)}2, where µ̂(x) is the point predictor fit to the augmented dataset Z1:n+1, assumed to be symmetric with respect to the permutation of the input
dataset. The key property of any conformity measure is that it is exchangeable in the first argument, i.e. the conformity measure for Zi is invariant to the permutation of Z1:n+1. Under the assumption that Z1:n+1 is exchangeable, we then have that σ1:n+1 is also exchangeable, and its rank is uniform among {1, . . . , n+ 1} (assuming continuous σ1:n+1). From this, we have that the rank of σn+1 is a valid p-value. If we now consider a plug-in value Yn+1 = y (where Xn+1 is known), we can define the rank of σn+1 among σ1:n+1 as
r(y) = 1
n+ 1 n+1∑ i=1 1 (σi ≤ σn+1) .
For miscoverage level α, the full conformal predictive set, Cα(Xn+1) = {y ∈ R : r(y) > α}, (3)
satisfies the desired frequentist coverage as in (2) under the exchangeability assumptions above. Intuitively, we are reporting the values of y which conform better than the fraction α of observed conformity scores in the augmented dataset. A formal proof can be found in Vovk et al. (2005, Chapter 8.7). For continuous σ1:n+1, we also have from Lei et al. (2018, Theorem 1) that the conformal predictive set does not significantly over-cover.
In practice, beyond a few exceptions, the function r(y) must be computed on a fine grid y ∈ Ygrid, for example of size 100, in which case the model must be retrained 100 times to the augmented dataset to compute σ1:n+1, with plug-in values for yn+1 on the grid. This is illustrated in Algorithm 1 below. We note here that the grid method only provides approximate coverage, as y may be selected even if it lies between two grid points that are not selected. This is formalized in Chen et al. (2018), but we do not discuss this further. In Appendix D.6, we provide an empirical comparison of the grid effects. This is also valid for binary classification where we now have a finite Ygrid = {0, 1}, and so the grid method for full conformal prediction is exact and feasible.
Observed data is Z1:n, Xn+1; Specify miscoverage level α for each y ∈ Ygrid do
Fit model to augmented dataset {Z1, . . . , Zn, {y,Xn+1}} Compute σ1:n and σn+1 Store the rank, r(y), of σn+1 among σ1:n+1
end Return the set Cα(Xn+1) = {y ∈ Ygrid : r(y) > α}.
Algorithm 1: Full Conformal Prediction

heading: 2.2 Conformal Bayes and Add-One-In Importance Sampling

text:
In a Bayesian model, a natural suggestion for the conformity score, as noted in Vovk et al. (2005); Wasserman (2011), is the posterior predictive density (1), that is
σ(Z1:n+1;Zi) = p(Yi | Xi, Z1:n+1) for i ∈ {1 . . . , n+ 1}, where
p(Yi | Xi, Z1:n+1) = ∫ fθ(Yi | Xi)π(θ | Z1:n+1) dθ. (4)
We highlight that the datum Zi = {Yi, Xi} also lies in the conditioning Z1:n+1. This is a valid conformity score, as we have π(θ | Z1:n+1) ∝ π(θ) ∏n+1 i=1 fθ(Yi | Xi), and so σ is indeed invariant to the permutation of Z1:n+1. We refer to this method as conformal Bayes (CB) - this is a special case of full conformal prediction, and we will see shortly that the exchangeability structure of Bayesian models is key to constructing conformity scores in the partial exchangeability scenario.
Beyond conjugate models, we are usually able to obtain (asymptotically exact) posterior samples θ(1:T ) ∼ π(θ | Z1:n), e.g through MCMC, where T is a large integer. Such samples are typically available as standard output from Bayesian model fitting. The posterior predictive can then be computed up to Monte Carlo error through
p̂(Yn+1 | Xn+1, Z1:n) = 1
T T∑ t=1 fθ(t)(Yn+1 | Xn+1).
The key insight is that refitting the Bayesian model with {Z1, . . . , Zn, {y,Xn+1}} is well approximated through importance sampling (IS), as only {y,Xn+1} changes between refits. This leads immediately to an IS approach to full conformal Bayes, where we just need to compute ‘add-one-in’ (AOI) predictive densities. Here AOI refers to the inclusion of {Yn+1, Xn+1} into the training set, named in relation to ‘leave-one-out’ (LOO) cross-validation. Specifically, for Yn+1 = y and θ(1:T ) ∼ π(θ | Z1:n), we can compute
p̂(Yi | Xi, Z1:n+1) = T∑ t=1 w̃(t)fθ(t)(Yi | Xi) (5)
where w̃(t) are our self-normalized importance weights of the form
w(t) = fθ(t)(y | Xn+1), w̃(t) = w(t)∑T t′=1 w (t′) · (6)
We see that the unnormalized importance weights have the intuitive form of the predictive likelihood at the reference point {y,Xn+1} given the model parameters θ(t). The use of AOI importance sampling has similarities to the computation of Bayesian leave-one-out (LOO) predictive densities for cross-validation (Vehtari et al., 2017), which is also used in accounting for model misspecification. An interesting aspect of AOI in comparison with LOO is that AOI predictive densities are less vulnerable to importance weight instability for the following reasons:
• In LOO, the target π(θ | Z−i) generally has thicker tails than the proposal π(θ | Z1:n), leading to importance weight instability. In contrast, AOI uses the posterior π(θ | Z1:n) as a proposal for the thinner-tailed π(θ | Z1:n+1). For LOO the importance weights are proportional to 1/fθ(y | x), in contrast to the typically bounded fθ(y | x) for AOI.
• For AOI, we are predicting Zi given Z1:n+1 which is always in-sample unlike in LOO where the datum is out-of-sample, so we can expect greater stability with AOI.
• The IS weight stability is governed by Yn+1 = y, which is not random as we select it for the grid, unlike in LOO cross-validation. For sufficiently large α, we will not need to compute the AOI predictive density for extreme values of y.
The above discussion is a natural analogy to the discussion provided in Vehtari et al. (2017, Section 2.1.1) for computing LOO predictives. Note that the LOO predictive density is also a valid conformity score; we provide more details and verify the above comparisons between LOO and AOI empirically in Appendix B. We also provide some IS weight diagnostics in Section 4.1.1 and find that they are stable. In difficult settings such as very high-dimensions, one can make use of the recommendations of Vehtari et al. (2015) for assessing and Pareto-smoothing the importance weights if necessary.

heading: 2.3 Computational Complexity

text:
Given the posterior samples, we must compute the likelihood for each θ(t) at Z1:n, as well at {y,Xn+1} for y ∈ Ygrid. The additional computation required for CB for each Xn+1 is thus T × (n + ngrid) likelihood evaluations, which is relatively cheap. This is then followed by the dot product of an (n + 1) × T matrix with a T vector for each y, which is O(nT ), so the overall complexity is O(ngridTn). The large matrices involved in computing the AOI predictives suggests we can take advantage of GPU computation, and machine learning packages such as JAX (Bradbury et al., 2018) are highly suitable for this application.
The values ngrid and T are constants we select, and we now briefly provide guidance on this. Selecting a grid is a general problem for full conformal prediction (beyond just CB), and this has been investigated in Chen et al. (2016, 2018). Here, we opt for a default choice of ngrid = 100 in Section 4.1 and we show in Appendix D.6 that this is sufficient for accurate coverage. For the selection of T , there are standard methods based on the effective sample size (ESS) which we compute in Section 4.1.1. One could then approximate the Monte Carlo standard error in the density estimates which is proportional to 1/ √ ESS, or estimate the IS variance directly. The ESS will depend both on MCMC mixing and IS stability, which may further depend on the dimensionality of the model. More details can be found in Owen (2013, Chapter 9.3)

heading: 2.4 Motivation

text:
Much has been written on the contrasting foundations and interpretation of Bayes versus frequentist measures of uncertainty (Little, 2006; Shafer and Vovk, 2008; Bernardo and Smith, 2009; Wasserman, 2011), and we provide a summary in Appendix A. Here we motivate CB predictive intervals from both a Bayesian and frequentist perspective.
The pragmatic Bayesian, aware of the potential for model misspecification in either the prior or likelihood, may be interested in conformal inference as a countermeasure. CB predictive intervals with guaranteed frequentist coverage can be provided as a supplement to the usual Bayesian predictive intervals. The difference between the Bayesian and conformal interval may also serve as an informal diagnostic for model evaluation (e.g. Gelman et al. (2013)). Posterior samples through MCMC or direct sampling are typically available, and so CB through automated AOI carries little overhead.
The frequentist may also wish to use a Bayesian model as a tool for constructing predictive confidence intervals. Firstly, the likelihood can take into account skewness, heteroscedasticity unlike the usual residual conformity score. Secondly, features such as sparsity, support, and regularization can be included through priors, while CB ensures correct coverage. Thirdly, CB provides a computationally feasible method to compute full conformal prediction intervals for a larger class of models. Finally, a subtle issue that arises in full conformal prediction is that we lose validity if hyperparameter selection is not symmetric with respect to Zn+1, e.g. if we estimate the Lasso penalty λ using only Z1:n before computing the full conformal intervals with said λ(Z1:n). For CB, a prior on hyperparameters induces weighting of the hyperparameter values by implicit cross-validation for each refit (Gneiting and Raftery, 2007; Fong and Holmes, 2020). We highlight here that this issue does not affect the split conformal method, e.g. Lei and Wasserman (2014, Section 5).

heading: 3 Partial Exchangeability and Hierarchical Models

text:
A setting of particular interest is for grouped data, which corresponds to a weakening of exchangeability often referred to as partial exchangeability (Bernardo and Smith, 2009, Chapter 4.6). Assume that we observe data from J groups, each of size nj , where again Zi,j = {Yi,j , Xi,j}. We write the full dataset as Z = {Zi,j : i = 1, . . . , nj , j = 1, . . . , J}. We may not expect the entire sequence Z to be exchangeable, instead only that data points are exchangeable within groups. Formally, this means that p(Z1:n1,1, . . . , Z1:nJ ,J) = p(Zγ1(1):γ1(n1),1, . . . , ZγJ (1):γJ (nJ ),J) (7) for any permutations γj of 1, . . . , nj , for j = 1, . . . , J . Alternatively, we can enforce the usual definition of exchangeability but only consider permutations of 1, . . . , n such that the groupings are preserved. A simple example of this partial exchangeability is if Zi,j
iid∼ Pj for i = 1, . . . , nj , j = 1, . . . , J , where Pj can now be distinct.
Partial exchangeability is useful in multilevel modelling, e.g. where Z1:nj ,j records exam results on students within school j, for schools j = 1, . . . , J . Students may be deemed exchangeable within schools, but not between schools. Further examples may be found in Gelman and Hill (2006).

heading: 3.1 Group Conformal Prediction

text:
Given a new Xnj+1,j belonging to group j for j ∈ {1, . . . , J}, we seek to construct a (1 − αj) confidence interval for Ynj+1,j . We define a within-group conformity score as
σi,j := σZ−j (Z1:nj+1,j ;Zi,j)
for i = 1, . . . , nj + 1. We write Z−j as the dataset without group j, and Z−j in the subscript indicates the dependence of the conformity score on this, which we motivate in the next subsection. For each Z−j , we require the score to be invariant with respect to the permutation of Z1:nj+1,j . For Znj+1,j = {y,Xnj+1,j}, the conformal predictive set is then defined
rj(y) = 1
nj + 1 nj+1∑ i=1 1 ( σi,j ≤ σnj+1,j ) , Cαj ( Xnj+1,j ) = {y ∈ R : rj(y) > αj) (8)
In other words, we rank the conformity scores σ1:nj+1,j within the group j, and compute the conformal interval as usual with Algorithm 1. The interval is valid from the following.
Proposition 1. Assume that {Z,Znj+1,j} is partially exchangeable as in (7), and the conformity measure σi,j for group j is invariant to the permutation of Z1:nj+1,j . We then have
P ( Ynj+1,j ∈ Cαj ( Xnj+1,j )) ≥ 1− αj
where Cαj ( Xnj+1,j ) is defined in (8), and P is over {Z,Znj+1,j}.
Proof. Conditional on Z−j , the observations Z1:nj+1,j are still exchangeable, and thus so are σ1:nj+1,j from the invariance of the conformity measure. The usual conformal guarantee then holds:
P ( Ynj+1,j ∈ Cαj ( Xnj+1,j ) | Z−j ) ≥ 1− αj .
Taking the expectation with respect to Z−j gives us the result.
It is interesting to note that the above group conformal predictor coincides with the attributeconditional Mondrian conformal predictor of Vovk et al. (2005, Chapter 4.5), with the group allocations as the taxonomy. Validity under the relaxed Mondrian-exchangeability of Vovk et al. (2005, Chapter 8.4) is key for us here.

heading: 3.2 Conformal Hierarchical Bayes

text:
Under this setting, a hierarchical Bayesian model can be defined of the form
[Yi,j | Xi,j , θj , τ ] iid∼ fθj ,τ (· | Xi,j) i = 1, . . . , nj , j = 1, . . . , J
[θj | φ] iid∼ π(· | φ) j = 1, . . . , J
φ ∼ π(φ), τ ∼ π(τ).
Here τ is a common parameter across groups (e.g. a common standard deviation for the residuals under homoscedastic errors). The desired partial exchangeability structure is clearly preserved in the Bayesian model (Bernardo, 1996). De Finetti representation theorems are also available for partially exchangeable sequences (when defined in a slightly different manner to the above), which motivate the specification of hierarchical Bayesian models (Bernardo and Smith, 2009, Chapter 4.6).
The posterior predictive is once again a natural choice for the conformity measure. Denoting Z̄y as the entire dataset augmented with Znj+1,j = {y,Xnj+1,j}, we have
σi,j = p(Yi,j | Xi,j , Z̄y) = ∫ fθj ,τ (Yi,j | Xi,j)π(θj , τ | Z̄y) dθj dτ (9)
for i = 1, . . . , nj + 1. The within-group permutation invariance follows as the likelihood is exchangeable within groups, and thus so is the posterior and resulting posterior predictive. Practically, this structure allows for independent coefficients θj for each group, but partial pooling through π(θ | φ) allows information to be shared between groups. A fully pooled model, whilst still valid, is usually too simple and predicts poorly, whereas a no-pooling conformity score ignores information sharing between groups. More details on hierarchical models can be found in Gelman et al. (2013, Chapter 5). We point out that we can select a separate coverage level αj for each group, which will be useful when group sizes nj vary - we provide a demonstration of this in Appendix D.4.2. Computation of σi,j is again straightforward, where MCMC now returns [θ (1:T ) 1:J , φ
(1:T ), τ (1:T )] ∼ π(θ1:J , φ, τ | Z). We can then estimate (9) using AOI importance sampling as in (5) and (6) using the marginal samples {θ(1:T )j , τ (1:T )} ∼ π(θj , τ | Z) and weights w(t) = fθ(t)j ,τ(t)(y | Xnj+1,j).
In the above, we consider predictive intervals within groups with covariates, extending the withingroup approach of Dunn et al. (2020). Predictive intervals for new groups are possible with the Bayesian model, but a conformal predictor would require additional stronger assumptions of exchangeability to ensure validity. The relevant assumptions and methods based on pooling and subsampling for new group predictions are discussed in Dunn et al. (2020), but would require rerunning MCMC in our case. We leave this for future work, noting that utilizing the Bayesian predictive density directly here seems nontrivial due to different group sizes.

heading: 4 Experiments

text: We run and time all examples on an Azure NC6 Virtual Machine, which has 6 Intel Xeon E5-2690 v3 vCPUs and a one-half Tesla K80 GPU card. We use PyMC3 (Salvatier et al., 2016) for MCMC and sklearn (Pedregosa et al., 2011) for the regular conformal predictor; both are run on the CPU. Computation of the CB and Bayes intervals is implemented in JAX (Bradbury et al., 2018), and run on the GPU. The code is available online1 and further examples are provided in Appendix D.

heading: 4.1 Sparse Regression

text:
We first demonstrate our method under a sparse linear regression model on the diabetes dataset (Efron et al., 2004) considered by Lei (2019). The dataset is available in sklearn, and consists of n = 442 subjects, where the response variable is a continuous diabetes progression and the d = 10 covariates consist of patient readings such as blood serum measurements. We standardize all covariates and the response to have mean 0 and standard deviation 1.
The Bayesian model we consider is
fθ(y | x) = N (y | θTx+ θ0, τ2) π(θj) = Laplace(0, b), π(θ0) ∝ 1, π(b) = Gamma(1, 1) π(τ) = N+(0, c)
(10)
for j = 1, . . . , d, and where b is the scale parameter andN+ is the half-normal distribution. Note that a hyperprior on b has removed the need for cross-validation that is required for Lasso. We consider two values of c for the hyperprior on τ , which correspond to a well-specified (c = 1) and poorlyspecified (c = 0.02) prior; in the latter case our posterior on τ will be heavily weighted towards a small value. This model is well-specified for the diabetes dataset (Jansen, 2013, Chapter 4.5) under a reasonable prior (c = 1). We compute the central (1 − α) credible interval from the Bayesian posterior predictive CDF estimated using Monte Carlo and the same grid as for CB.
To check coverage, we repeatedly divide into a training and test dataset for 50 repeats, with 30% of the dataset in the test split. We evaluate the conformal prediction set on a grid of size ngrid = 100 between [ymin − 2, ymax + 2], where ymin, ymax is computed from each training dataset. The average coverage, length and run-times (excluding MCMC) with standard errors are given in Table 1 for α = 0.2. MCMC induced an average overhead of 21.9s for a = 1 and 26.8s for c = 0.02 for the Bayes and CB interval, where we simulate T = 8000 posterior samples. The CB intervals are only slightly slower than the Bayes intervals, and still a small fraction of the time required for MCMC, and is thus an efficient post-processing step. For c = 1, the Bayesian intervals have coverage close to (1− α) with the smallest expected length, with CB slightly wider and more conservative. However, when the prior is misspecified with c = 0.02, the Bayes intervals severely undercover, whilst the CB coverage and length remain unchanged from the c = 1 case.
As baselines, we compare to the split and full conformal method using the non-Bayesian Lasso as the predictor, with the usual residual as the nonconformity score. For the split method, we fit Lasso with cross-validation on the subset of size ntrain/2 to obtain the Lasso penalty λ. For the full conformal method, we use the grid method for fair timing, as other estimators beyond Lasso would not have the shortcut of Lei (2019). As setting a default λ = 1 gives poor average lengths, we estimate λ = 0.004 on cross-validation on one of the training sets, and use this value over the 50 repeats. However, we must emphasize again that this is somewhat misleading, as discussed in Section 2.4. A fairer approach would involve fitting Lasso with CV for each of the 100 grid values and 133 test values, but this is infeasible as each fit requires around 80ms, resulting in a total run-time of 17 minutes. On the other hand, the AOI scheme of CB is equivalent to refitting b for each grid/test value. In terms of performance, the split method has wider intervals than CB/full, but performs well given the low computational costs. The full conformal method performs as well as CB, but is comparable in time as MCMC + CB, whilst not refitting λ. We note that the value of c does not affect the split/full method.
1https://github.com/edfong/conformal_bayes

heading: 4.1.1 Importance Weights

text: For the diabetes dataset, we look at the ESS of the self-normalized importance weights (6), which can be computed as ESS = 1/ ∑T t=1{w(t)}2 for each xn+1 and y. The ESS as a function of y for a single xn+1 is shown in Figure 1 for the two cases c = 1, 0.02, with the CB conformal bands given for α = 0.2, 0.5. We have scaled the ESS plots by ESSMCMC/T , where T = 8000 is the number of posterior samples and ESSMCMC is the minimum ESS out of all posterior parameters return by PyMC3. We observe the ESS is well behaved and stable across the range of y values. In both cases, the ESS for α = 0.2 is sufficiently large for a reliable estimate of the conformity scores. However, for c = 0.02, the ESS decays more quickly with y as the Bayes predictive intervals are too narrow, which the CB corrects for. Other values of xn+1 produce similar behaviour.

heading: 4.2 Sparse Classification

text:
In this section, we analyze the Wisconsin breast cancer (Wolberg and Mangasarian, 1990), again available in sklearn. The dataset is of size 569, where the binary response variable corresponds to a malignant or benign tumour. The 30 covariates consist of measurements of cell nuclei. Again, we standardize all covariates to have mean 0 and standard deviation 1.
We consider the logistic likelihood fθ(y = 1 | x) = [1 + exp {− (θTx+ θ0)}]−1, with the same priors for θ, θ0 as in (10). The Bayesian predictive set is the smallest set from {0}, {1}, {0, 1} that contains at least (1− α) of the posterior predictive probability. The conformal baselines are as above but with L1-penalized logistic regression, and for the full conformal method we have λ = 1. We again have 50 repeats with 70-30 train-test split, and set α = 0.2. The grid method is now exact, and the size of the CB intervals can take on the values {0, 1, 2}. The results are provided in Table 2, where MCMC required an average of 45.4s to produce T = 8000 samples. We see that even with reasonable priors, Bayes can over-cover substantially, which CB corrects in roughly the same amount of time as it takes to compute the usual Bayes interval. However, we point out that CB may produce empty prediction sets, whereas Bayes cannot, and we investigate this in Appendix D.3.2.

heading: 4.3 Hierarchical Model

text:
We now demonstrate Bayesian conformal inference using a hierarchical Bayesian model for multilevel data. We stick to the varying intercept and varying slope model (Gelman et al., 2013), that is for j = 1, . . . , J :
fθj ,τ (yi,j) = N (yi,j | θTjXi,j + θ0,j , τ2) π(θj) = N (φ, s2), π(θ0,j) = N (φ0, s20)
(11)
with hyperpriors N (0, 1) on the location parameters φ, φ0 and Exp(1) on the standard deviations s, s0, τ . We now apply this to a simulated example, and an application to the radon dataset of Gelman and Hill (2006) is given in Appendix D.4.2.
We consider two simulation scenarios, with J = 5 groups and nj = 10 elements per group:
1. Well-specified: We generate group slopes θj iid∼ N (0, 1) for j = 1, . . . , J . For each j, we
generate Xi,j ∼ N (0, 1) and Yi,j ∼ N (θjXi,j , 1).
2. Misspecified: We generate group slopes and variances θj iid∼ N (0, 1), τj iid∼ Exp(1) for j = 1, . . . , J . For each j, we generate Xi,j ∼ N (0, 1) and Yi,j ∼ N ( θjXi,j , τ 2 j ) .
The first scenario has homoscedastic noise between groups as assumed in the model (11) whereas the second scenario is heteroscedastic between groups. To evaluate coverage, we only draw θ1:J , τ1:J once (and not per repeat), giving us the values
θ1:J = [1.33,−0.77,−0.32,−0.99,−1.07], τ1:J = [1.24, 2.30, 0.76, 0.28, 1.11].
For each of the 50 repeats, we draw nj = 10 training and test data points from each group using the above θ1:J (and τ1:J for scenario 2), and report test coverage and lengths within each group. We use a grid of size 100 between [−10, 10]. The group-wise average lengths and coverage are given in Table 3 again with α = 0.2. Again run-times are given post-MCMC, where MCMC required an average of 90.1s and 78.4s for scenarios 1 and 2 respectively to generate T = 8000 samples. The Bayes interval is again the central (1− α) credible interval. The CB and Bayes methods have comparable run-times, likely due to the small n. As a reference, fitting a linear mixed-effects model in statsmodels (Seabold and Perktold, 2010) to the dataset takes around 200ms, so the full conformal method, which requires refitting for each of the 100 grid value and 50 test values, would take a total of 17 minutes. For scenario 1, both Bayes and CB provide close to (1− α) coverage, with the Bayes lengths being smaller. This is unsurprising, as the Bayesian model is well-specified. In scenario 2, the Bayes intervals noticeably over/under-cover depending on the value of τ1:J in relation to the Bayes posterior mean τ̄ ≈ 1.3. CB is robust to this, adapting its interval lengths accordingly (in particular for Groups 2 and 4) and providing within-group validity.

heading: 5 Discussion

text:
In this work, we have introduced the AOI importance sampling scheme for conformal Bayesian computation, which allow us to construct frequentist-valid predictive intervals from a baseline Bayesian model using the output of an MCMC sampler. This extends naturally to the partially exchangeable setting and hierarchical Bayesian models.
Under model misspecification, or theM-open scenario (Bernardo and Smith, 2009), CB can produce calibrated intervals from the Bayesian model. In the partially exchangeable case, CB can remain valid within groups. We find that even under reasonable priors, Bayesian predictives can over-cover, and CB can help reduce the length of intervals to get closer to nominal coverage. Diagnosing Bayesian miscalibration is in general non-trivial, but CB intervals may provide a diagnostic tool while also automatically correcting the miscalibration. When posterior samples of model parameters are available, AOI importance sampling is only a minor increase in computation, and interestingly is much faster than the split method (with the Bayesian model) which would require another run of MCMC. For the frequentist, CB intervals enjoy the tightness of the full conformal method, for a single expensive fit with MCMC followed by a cheap refitting process. We are also free to incorporate prior information, and use more complex likelihoods or priors, as well as automatically fitting hyperparameters.
There are however limitations to our approach, dictated by the realities of MCMC and IS. Firstly, the intervals are only approximate up to MC error, although we expect some robustness as the CB intervals are computed from the ranks and not the exact values. The stability of AOI importance sampling also depends on the posterior predictive being a good proposal, which may break down if the addition of the new datum {y,Xn+1} has very high effect on the posterior. Finally, if our Bayesian model is grossly misspecified and the predictive densities are too narrow, the IS weights may have low ESS and the AOI estimates may have high variance.
If only approximate posterior samples are available, e.g. through variational Bayes (VB), then an AOI scheme may still be feasible, where one includes an additional correction term in the IS weights for the VB approximation, e.g. in Magnusson et al. (2019). However, this remains to be investigated, as the IS weight stability is strongly reliant on the quality of the approximate posterior. Combining this with the Pareto-smoothed IS method of Vehtari et al. (2015) may lead to a feasible solution. In our experience, CB intervals tend to be a single connected interval, which may allow for computational shortcuts in adapting the search grid. It would also be interesting to pursue the theoretical connections between the Bayesian and CB intervals, in a similar light to Burnaev and Vovk (2014).

heading: Acknowledgements

text: We thank the anonymous reviewers for their feedback on our paper. Fong is funded by The Alan Turing Institute Doctoral Studentship, under the EPSRC grant EP/N510129/1. Holmes is supported by The Alan Turing Institute, the Health Data Research, U.K., the Li Ka Shing Foundation, the Medical Research Council, and the U.K. Engineering and Physical Sciences Research Council through the Bayes4Health programme grant EP/R018561/1.

title: Conformal Bayesian Computation


summaries: Summary:  - This submission deals with designing stable and smooth denoisers that could be directly used for regularization, or for plug-and-play (PnP) methods to solve inverse problems.  - This is an improvement upon existing denoisers such as RED and PnP which are not necessarily stable.  - The key idea is to craft neural networks that output scalar-valued potential function g whose gradient is the denoiser residual. To guarantee convergence and stability, g has to be smooth and h=1/2x^2 - g convex.  - Advantages: the denoisers lead to symmetric Jacobians, allowing for a MAP estimation; denoisers can be integrated into optimization schemes with backtracking step size, and thus no need for enforcing Lipschitz continuity. - To construct such a network, GraDnCNN stacks conv layers and nonlinear activations, and proposes DnICNNs with non-negative weights and convex activations. Variants of DnICNN are discussed that achieve better denoising/recovery performance for Gaussian deblurring, and image superresolution.  Contributions: - construct stable denoisers based on neural networks as gradients of smooth potential functions g - integrate the denoisers as a regularizer for solving inverse problems and test it out for gaussian image deblurring and superresolution which outperforms RED and traditional PnP methods

summaries: This paper considers the problem of designing an image denoiser that can be used as a prior within plug-and-play priors (PnP) or regularization by denosing (RED). Prior work has proved the convergence of PnP/RED for contractive, nonexpansive, and/or firmly nonexpansive denoisers. The design of denoisers that enable PnP/RED convergence is an active research area.  This paper focuses on denoisers that correspond to the gradients of some potential functions. By using such denoisers within RED, one can guarantee convergence via traditional optimization. The key idea in the paper is to use a deep neural network to parameterize a potential function (instead of the denoiser), then train that function for denoising by using its gradient. Additionally, when the weights of the network are nonlinear and activation functions are nondecreasing and convex, the denoiser corresponds to a convex function.  The key contribution of the paper is the formulation of a procedure for training a potential function that can be used to obtain a denoiser. The theoretical arguments presented in Lemma 2 and Theorem 1, are well known in optimization (see below for reference). Similarly, Algorithm 1 corresponds exactly to the RED-SD algorithm using a potential-driven denoiser. The empirical results compare three proposed variants of denoisers against PnP-PGD and RED-SD on problems of image deblurring and superresolution.

summaries: The paper proposes image denoising based CNNs as regularizers that can be constructed as the gradient of smooth potential functions $g$. The main promise is mentioned to be making an optimization algorithm has explicit regularization, when equipped with those denoisers. Unlike original regularization by denoising (RED) that the explicit regularizer exists if the denoisers satisfy several strict conditions, the proposed denoisers only rely on smoothness. Particularly, the denoisers consist of smooth activation functions to ensure the differential continuity and are trained in a residual denoising fashion using MSE loss. As a result, the proposed algorithm can be formulated as steepest descent variant of RED optimized by directly adopting the gradient of the objective function $F$. Standard back-tracking is used to ensure convergence. A convergence analysis is provided following the common optimization theory. Finally, the performance of the proposed method is evaluated over two image inverse problems such as gaussian deblur and super-resolution, with satisfactory results compared to existing methods based on PnP and RED.

summaries: Motivated by the RED algorithm for solving inverse problems, the authors propose a neural-net denoiser that is the gradient of an explicit potential function.  This guarantees that the Jacobian is symmetric definite, which then guarantees that the RED algorithm minimizes an explicit cost function, which in turn makes convergence straightforward to prove, thereby solving an important open problem about the convergence of RED and plug-and-play (PnP)-type algorithms.  The fact that there is an explicit cost also enables the use of adaptive stepsize selection within RED, with one particular scheme proposed by the authors.  Numerical comparisons to RED and PnP algorithms show improved PSNR and SSIM performance, as well as enhanced stability.

summaries: **Summary**: This paper addresses an important issue in the framework of regularization by denoising (RED), which is that the explicit RED regularizer $\frac{1}{2}x^T(x-D_\sigma(x))$ exists only if a set of strict conditions on the denoiser are met. Due to this problem, RED fails to formulate explicit regularizers for advanced denoisers (e.g. CNN) since they generally do not satisfy the conditions. Although some work has established some theoretical analysis for RED algorithms using deep denoisers, a linkage between the denoiser and the regularizer is generally missing.   In this paper, the authors proposed to formulate the regularizer as a *smooth potential functional* $g(x)$, which is parameterized by a CNN with smooth activation functions. To link the potential functional to a denoiser, the authors trained the CNN such that the gradient $\nabla g(x)$ outputs the noise residual for the AWGN. In this way, the proposed method can formulate the gradient-based RED update by directly taking the gradient of the objective function $F(x)$.  Since the objective function is available, the authors also proposed to use the standard back-tracking trick in their framework to guarantee convergence. Additionally, a convergence analysis is also introduced by following the common nonconvex optimization analysis given the handy $F(x)$.

abstractText: In recent years there has been increasing interest in leveraging denoisers for solving general inverse problems. Two leading frameworks are regularization-by-denoising (RED) and plug-and-play priors (PnP) which incorporate explicit likelihood functions with priors induced by denoising algorithms. RED and PnP have shown state-of-the-art performance in diverse imaging tasks when powerful denoisers are used, such as convolutional neural networks (CNNs). However, the study of their convergence remains an active line of research. Recent works derive the convergence of RED and PnP methods by treating CNN denoisers as approximations for maximum a posteriori (MAP) or minimum mean square error (MMSE) estimators. Yet, state-of-the-art denoisers cannot be interpreted as either MAP or MMSE estimators, since they typically do not exhibit symmetric Jacobians. Furthermore, obtaining stable inverse algorithms often requires controlling the Lipschitz constant of CNN denoisers during training. Precisely enforcing this constraint is impractical, hence, convergence cannot be completely guaranteed. In this work, we introduce image denoisers derived as the gradients of smooth scalar-valued deep neural networks, acting as potentials. This ensures two things: (1) the proposed denoisers display symmetric Jacobians, allowing for MAP and MMSE estimators interpretation; (2) the denoisers may be integrated into RED and PnP schemes with backtracking step size, removing the need for enforcing their Lipschitz constant. To show the latter, we develop a simple inversion method that utilizes the proposed denoisers. We theoretically establish its convergence to stationary points of an underlying objective function consisting of the learned potentials. We numerically validate our method through various imaging experiments, showing improved results compared to standard RED and PnP methods, and with additional provable stability.

authors:
M. Abadi
P. Barham
J. Chen
Z. Chen
A. Davis
J. Dean
M. Devin
S. Ghemawat
G. Irving
M Isard

title: Tensorflow: A system for large-scale machine learning

venue: In: Symposium on Operating Systems Design and Implementation

authors:
E. Agustsson
R. Timofte

title: Ntire 2017 challenge on single image super-resolution: Dataset and study

venue: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops

authors:
R. Ahmad
C.A. Bouman
G.T. Buzzard
S. Chan
S. Liu
E.T. Reehorst
P. Schniter

title: Plug-and-play methods for magnetic resonance imaging: Using denoisers for image recovery

venue: IEEE Signal Processing Magazine 37(1),

authors:
B. Amos
L. Xu
J.Z. Kolter

title: Input convex neural networks

venue: In: International Conference on Machine Learning

authors:
M. Bayati
A. Montanari

title: The dynamics of message passing on dense graphs, with applications to compressed sensing

venue: IEEE Transactions on Information Theory 57(2),

authors:
A. Beck
M. Teboulle

title: Gradient-based algorithms with applications to signal recovery. Convex optimization in signal processing and communications pp

authors:
R. Berthier
A. Montanari
P.M. Nguyen

title: State evolution for approximate message passing with nonseparable functions

venue: Information and Inference: A Journal of the IMA 9(1),

authors:
G.T. Buzzard
S.H. Chan
S. Sreehari
C.A. Bouman

title: Plug-and-play unplugged: Optimization-free reconstruction using consensus equilibrium

venue: SIAM Journal on Imaging Sciences

title: Performance analysis of plug-and-play ADMM: A graph signal processing perspective

venue: IEEE Transactions on Computational Imaging

authors:
S.H. Chan
X. Wang
O.A. Elgendy

title: Plug-and-play ADMM for image restoration: Fixed-point convergence and applications

venue: IEEE Transactions on Computational Imaging 3(1),

authors:
P. Chatterjee
P. Milanfar

title: Is denoising dead

venue: IEEE Transactions on Image Processing 19(4),

authors:
R. Cohen
M. Elad
P. Milanfar

title: Regularization by denoising via fixed-point projection (RED-PRO)

venue: arXiv preprint arXiv:2008.00226

authors:
P.L. Combettes
J.C. Pesquet

title: Fixed point strategies in data science

venue: IEEE Transactions on Signal Processing

authors:
A. Danielyan
V. Katkovnik
K. Egiazarian

title: Image deblurring by augmented Lagrangian with BM3D frame prior. In: Workshop on Information Theoretic Methods in Science and Engineering (WITMSE)

venue: Tampere, Finland. pp

authors:
D.L. Donoho
A. Maleki
A. Montanari

title: Message-passing algorithms for compressed sensing

venue: Proceedings of the National Academy of Sciences

title: Tweedie’s formula and selection bias

venue: Journal of the American Statistical Association 106(496),

authors:
K. Egiazarian
V. Katkovnik

title: Single image super-resolution via BM3D sparse coding

venue: IEEE European Signal Processing Conference (EUSIPCO). pp

authors:
A.K. Fletcher
P. Pandit
S. Rangan
S. Sarkar
P. Schniter

title: Plug-in estimation in high-dimensional linear inverse problems: A rigorous analysis

venue: Advances in Neural Information Processing Systems

title: Should penalized least squares regression be interpreted as maximum a posteriori estimation

venue: IEEE Transactions on Signal Processing 59(5),

authors:
R. Gribonval
M. Nikolova

title: On Bayesian estimation and proximity operators. Applied and Computational Harmonic Analysis

title: A characterization of proximity operators

venue: Journal of Mathematical Imaging and Vision 62(6),

authors:
U.S. Kamilov
H. Mansour
B. Wohlberg

title: A plug-and-play priors approach for solving nonlinear imaging inverse problems

venue: IEEE Signal Processing Letters 24(12),

authors:
B. Kawar
G. Vaksman
M. Elad

title: SNIPS: Solving noisy inverse problems stochastically

title: Stochastic image denoising by sampling from the posterior distribution

authors:
J. Koh
J. Lee
S. Yoon

title: Single-image deblurring with neural networks: A comparative survey

venue: Computer Vision and Image Understanding 203,

authors:
W.S. Lai
J.B. Huang
Z. Hu
N. Ahuja
M.H. Yang

title: A comparative study for single image blind deblurring

authors:
R. Laumont
V. De Bortoli
A. Almansa
J. Delon
A. Durmus
M. Pereyra

title: Bayesian imaging using plug & play priors: when Langevin meets Tweedie

authors:
H.A. Le Thi
T.P. Dinh

title: DC programming and DCA: thirty years of developments

venue: Mathematical Programming 169(1),

authors:
A. Levin
B. Nadler

title: Natural image denoising: Optimality and inherent bounds

venue: IEEE CVPR. pp

authors:
A. Levin
B. Nadler
F. Durand
W.T. Freeman

title: Patch complexity, finite pixel correlations and optimal denoising

venue: European Conference on Computer Vision

authors:
B. Lim
S. Son
H. Kim
S. Nah
K. Mu Lee

title: Enhanced deep residual networks for single image super-resolution

venue: Proceedings of the IEEE conference on computer vision and pattern recognition workshops

authors:
D. Martin
C. Fowlkes
D. Tal
J. Malik

title: A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics

venue: Proc. 8th Int’l Conf. Computer Vision

authors:
C.A. Metzler
A. Maleki
R.G. Baraniuk

title: From denoising to compressed sensing

venue: IEEE Transactions on Information Theory 62(9),

title: Symmetrizing smoothing filters

authors:
P. Nair
R.G. Gavaskar
K.N. Chaudhury

title: Fixed-point and objective convergence of plug-and-play algorithms

authors:
J.C. Pesquet
A. Repetti
M. Terris
Y. Wiaux

title: Learning maximally monotone operators for image recovery

venue: arXiv preprint arXiv:2012.13247

authors:
E.T. Reehorst
P. Schniter

title: Regularization by denoising: Clarifications and new interpretations

venue: IEEE Transactions on computational imaging 5(1),

authors:
J. Rick Chang
C.L. Li
B. Poczos
B. Vijaya Kumar
A.C. Sankaranarayanan

title: One network to solve them all: Solving linear inverse problems using deep projection models

authors:
Y. Romano
M. Elad
P. Milanfar

title: The little engine that could: Regularization by denoising (RED)

authors:
E. Ryu
J. Liu
S. Wang
X. Chen
Z. Wang
W. Yin

title: Plug-and-play methods provably converge with properly trained denoisers

authors:
S. Sreehari
S.V. Venkatakrishnan
B. Wohlberg
G.T. Buzzard
L.F. Drummy
J.P. Simmons
C.A. Bouman

title: Plug-and-play priors for bright field electron tomography and sparse interpolation

venue: IEEE Transactions on Computational Imaging 2(4),

authors:
Y. Sun
J. Liu

title: Kamilov, U.S.: Block coordinate regularization by denoising

venue: Advances in Neural Information Processing Systems. pp

authors:
Y. Sun
B. Wohlberg

title: Kamilov, U.S.: An online plug-and-play algorithm for regularized image reconstruction

authors:
A.M. Teodoro
J.M. Bioucas-Dias
M.A. Figueiredo

title: Scene-adapted plug-and-play algorithm with convergence guarantees

venue: IEEE 27th International Workshop on Machine Learning for Signal Processing (MLSP). pp

title: A convergent image fusion algorithm using sceneadapted Gaussian-mixture-based denoising

venue: IEEE Transactions on Image Processing

authors:
M. Terris
A. Repetti
J.C. Pesquet
Y. Wiaux

title: Building firmly nonexpansive convolutional neural networks

venue: IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)

authors:
T. Tirer
R. Giryes

title: Image restoration by iterative denoising and backward projections

venue: IEEE Transactions on Image Processing 28(3),

authors:
S.V. Venkatakrishnan
C.A. Bouman
B. Wohlberg

title: Plug-and-play priors for model based reconstruction

venue: IEEE Global Conference on Signal and Information Processing. pp

authors:
K. Wei
A. Aviles-Rivero
J. Liang
Y. Fu
C.B. Schönlieb
H. Huang

title: Tuning-free plug-and-play proximal algorithm for inverse imaging problems

authors:
Z. Wu
Y. Sun
J. Liu
U.S. Kamilov

title: Online regularization by denoising with applications to phase retrieval

venue: Proceedings of the IEEE International Conference on Computer Vision Workshops

authors:
J. Xing
U. Kamilov
W. Wu
Y. Wang
M. Zhang

title: Plug-and-play priors for reconstruction-based placental image registration. In: Smart Ultrasound Imaging and Perinatal, Preterm and Paediatric

venue: Image Analysis,

authors:
X. Xu
J. Liu
Y. Sun
B. Wohlberg

title: Kamilov, U.S.: Boosting the performance of Plug-and-Play Priors via denoiser scaling

authors:
X. Xu
Y. Sun
J. Liu
B. Wohlberg

title: Kamilov, U.S.: Provable convergence of plug-and-play priors with MMSE denoisers

venue: IEEE Signal Processing Letters 27,

authors:
K. Zhang
W. Zuo
Y. Chen
D. Meng
L. Zhang

title: Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising

authors:
K. Zhang
W. Zuo
S. Gu
L. Zhang

title: Learning deep CNN denoiser prior for image restoration

venue: Proceedings of the IEEE CVPR. pp

authors:
K. Zhang
W. Zuo
L. Zhang

title: Deep plug-and-play super-resolution for arbitrary blur kernels

heading: 1 Introduction

text:
General Recovery Inverse problems aim to recover an unknown signal from its observed corrupted or compressed version. This fundamental task appears in a variety of fields such as machine learning, signal processing, computer vision and medical imaging. When a model for the degradation process is available, one may derive a data fidelity term which relates the measurements to the unknown signal through a forward operator. However, this term alone typically leads to ill-posed optimization problems with unstable solutions, requiring the incorporation of additional knowledge of the desired
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
signal into the inversion process. This has given rise to diverse regularization methods which have been investigated extensively over the last few decades.
RED and PnP A very basic inverse problem is denoising, in which a clean signal is to be recovered from its measurements contaminated with additive noise. During the last period of time there has been an overwhelming advancement in the ability to perform denoising, creating techniques that approach optimal performance [11, 29, 30]. This striking progress in denoising has been the motivation for the plug-and-play priors (PnP) [48, 10, 41] methodology, which introduced the paradigm of using denoisers in order to solve general inverse problems. PnP incorporates implicit regularization into the inverse solutions by replacing proximal operators in common optimization algorithms with powerful denoisers. This simple modification has led to remarkable performance in diverse inverse tasks [41, 50, 38, 22, 56, 3, 51], as well as to the development of an arsenal of PnP-based techniques which includes proximal-type methods but also goes beyond them [33, 7, 18, 15, 5, 55, 14, 17]. However, the lack of a clearly formulated regularization term corresponding to the denoiser poses a serious challenge in ensuring the stability of PnP-based methods, and typically requires careful parameter-tuning [49].
To fill this void, the regularization by denoising (RED) [39] framework presents an explicit Laplacianbased regularization that relies on a denoiser. Under certain conditions on the denoiser [39, 37], the RED functional is convex and its gradient is given by the denoising residual, allowing for easy incorporation in standard optimization procedures. RED has demonstrated impressive recovery performance while providing global convergence guarantees. However, various common denoisers do not satisfy the conditions required by RED [37]. In such cases, there exists no objective function which justifies the RED method [37] and the corresponding convergence analysis does not hold. Thus RED, as well as PnP, shows state-of-the-art empirical results but without provable stability. A comprehensive overview of PnP and RED is given in [3].
Related Work The study of the stability of both RED and PnP methods has become an active line of research in recent years. The convergence of two PnP techniques was proven in [10] and [47] under the assumption of bounded denoisers. When minimum mean squared error (MMSE) denoisers are used, the PnP framework was shown to converge vis-a-vis an explicit objective function [53]. Similarly, [7, 18] proved the convergence of approximate message passing algorithms with MMSE denoisers for random froward operators. Multiple studies [41, 43, 44, 9, 45] derived the convergence of PnP by relying on nonexpansive denoisers. The latter has been the core assumption of [8], which explained the PnP framework, including its stability, through the notion of consensus equilibrium. A similar condition was used in [37], which adopted the score-matching framework to provide a new interpretation of the RED algorithms and their convergence. [42] presented a block-coordinate RED method and derived its fixed-point convergence for block-nonexpansive denoisers. The work of [36] introduced a convergent PnP-type method using a maximally monotone denoiser, created by training a nonexpansive neural network in a supervised fashion. Alternatively, [40] provided convergence guarantees for certain PnP variants that use trained denoisers whose residuals are nonexpansive. This type of denoiser was used in [27] to develop Bayesian inference with PnP priors, including the extension of the PnP framework to Monte Carlo sampling. The authors of [24, 23] employ denoisers within Langevin dynamics to sample from the posterior distribution of any linear inverse problem. A different strategy was studied in [35], which proposed a special inner product derived from the given linear denoiser, and thus proved the convergence of PnP based on the theory of averaged nonexpansive operators. Recent work [12] introduced RED via fixed-projection framework (RED-PRO), which utilizes the fixed-points of the denoisers as regularization. RED-PRO unifies RED and PnP and derives their convergence for general demicontractive denoisers, including nonexpansive ones.
Recent studies, mentioned above, are motivated by or rely on the assumption that the denoiser is an MMSE or maximum a posteriori (MAP) estimator. However, such estimators are in the form of gradients, hence, they must exhibit symmetric Jacobians [20, 19, 21] amongst other properties. Furthermore, Jacobian symmetry offers additional advantages as it often leads to stable algorithms and allows one to better understand the behavior of such non-linear estimators [34]. Unfortunately, many state-of-the-art denoisers do not fulfill this condition [37]. An additional drawback of current RED and PnP methods is that the Lipschitz constant of the denoiser often has to be controlled to create a nonexpansive operator, which in turn ensures convergence. This is typically done by learning a deep network and adding a regularization term to the training loss function, whose purpose is to bound the Lipschitz constant of the resultant denoiser. However, exactly enforcing this property of the denoiser is impractical and it is only approximated in practice.
Contributions In this work, we develop a provably convergent technique which exploits denoisers for solving general inverse problems. (1) We introduce image denoisers that are constructed as the gradients of smooth potential functions, represented by deep neural networks and trained in a supervised fashion. (2) Equipped with these learned potential networks, we formulate a general inverse framework where we utilize our denoising potentials as an explicit regularization. (3) We use the existence of an underlying objective function to incorporate the proposed denoisers into standard optimization schemes with backtracking step size. Specifically, we present a simple gradient descent algorithm where the step size is tuned automatically to provably guarantee convergence to stationary points of the objective function. (4) We numerically validate our method through imaging experiments, showing improved results compared to RED and PnP, with additional provable stability.

heading: 2 Background

text:
Inverse Problems: Formulation We consider the recovery of an unknown signal x∗ ∈ Rdx from its corrupted observation y ∈ Rdy . The inversion task is traditionally formulated as an optimization problem of the form
x̂ = argmin x∈Rdx f(x; y) + λg(x), (1)
where f(· ; y) is a fidelity term related to the observations y, g(·) is a regularization functional that incorporates prior information or known statistics of the true signal x∗ to avoid ill-posed problems, and λ > 0 balances between the terms. A common model for the measurements is y = Ax +w where A is a linear degradation (forward) operator and w is an i.i.d. white Gaussian noise, leading to the following log-likelihood function
f(x;y) = 1
2 ||Ax− y||22. (2)
The above model applies to a wide variety of inverse problems such as compressed sensing, super resolution, deblurring, inpainting and much more.
In the following, we consider a general functional f(·) where we suppress its dependency on the observations y for brevity. Moreover, both g(·) and f(·) are typically assumed to be convex where the former is generally non-differentiable, while the latter is a smooth function. By the definition of smooth functions, there exists a Lipschitz constant Lf > 0 such that the gradient∇f of f(·) satisfies
||∇f(x1)−∇f(x2)|| ≤ Lf ||x1 − x2||, ∀ x1, x2 ∈ Rdx . (3)
Proximal Methods The minimization problem (1) has been investigated extensively over the decades, leading to a series of diverse approaches. Many common optimization techniques for solving (1) rely on the proximal operator of g(·), defined as
Pλg(x) , argmin z∈Rdx
1 2 ||x− z||+ λg(z). (4)
The proximal operator, which is a MAP estimator, provides an efficient way to overcome the lack of differentiablity of g(·), especially when it admits a closed-form expression, for example in the case of the `1-norm. A simple representative of this family of algorithms is the proximal gradient descent (PGD) method which updates the solution iteratively as follows
xk+1 = Pλg ( xk − µk∇f(xk) )
(5)
where µk > 0 is the step size. PGD has shown great performance in various applications along with strong convergence guarantees.
The PnP Framework Notice that the proximal operator can be seen, by definition, as the solution to a denoising problem. This observation serves as the basis of the PnP framework which extends standard proximal methods by replacing the proximal operator with powerful denoisers. Following this rationale, (5) is modified to create the PnP-PGD method that performs the following update rule
xk+1 = D ( xk − µk∇f(xk) ) , (PnP-PGD)
where D : Rdx → Rdx is an off-the-shelf denoiser.
The RED Framework An alternative strategy for leveraging advanced denoisers is the RED framework which utilizes the denoiser to construct an explicit regularization g(x) = 12 〈x, x−D(x)〉. Under certain conditions on the denoiser [39, 37], g(·) is convex and its gradient is readily given by the denoising residual∇g(x) = x−D(x). This allows for the solution of (1) by standard algorithms, e.g., steepest (or gradient) descent
xk+1 = xk − µk
1 + λ
( ∇f(xk) + λ ( xk −D(xk) )) . (RED-SD)
Surprisingly, RED-based methods have been found to be effective also when the required conditions on the denoiser are not satisfied and the expression for the regularization term is invalid. In such cases, RED uses the denoiser as an implicit regularization, similar to PnP.
Challenges in RED and PnP Both RED and PnP provide extremely efficient means for exploiting powerful denoisers to solve inverse problems. Yet their analysis poses two closely related challenges:
1. Asymmetric Denoisers: Many studies on RED and PnP [53, 27, 43, 44, 9, 45, 52, 42] motivate or base their derivations on the assumption that the denoiser is an approximation of either MAP or MMSE estimator. It has been proven that when the Jacobian of such estimators is well-defined it must be symmetric [20, 19, 21]. However, RED and PnP typically reach state-of-the-art results by relying on denoisers which do not exhibit this property, as in the case of the popular DnCNN technique; as a result, they lose their interpretation as optimization solvers. In such cases, while there are alternative optimization-free explanations for RED and PnP [8, 37], it is challenging to properly analyze their behavior and characterize the nature of their solutions.
2. Convergence Issues: Continuing the above line of reasoning, often there exists no regularization that explains the denoiser, namely no functional g(·) for which∇g(x) = D(x) or∇g(x) = x−D(x) [37]. This lack of a potential function creates difficulties in studying the theoretical convergence of the RED and PnP algorithms. However, a common sufficient condition to guarantee convergence is for the denoiser D(·) to be (averaged) nonexpansive [13], implying that
||D(x1)−D(x2)|| ≤ L||x1 − x2||, ∀ x1, x2 ∈ Rdx , (6) for some Lipschitz constant 0 < L ≤ 1. However, most predefined denoisers do not satisfy this assumption. An alternative approach [40, 36, 46] is to learn a nonexpansive denoiser by adding a regularization term to the training loss which constrains its Lipschitz constant. Unfortunately, this remains an open challenge since in general there are no practical means to absolutely enforce the resultant denoiser to be nonexpansive.
Our work addresses both of these challenges. We first introduce a new class of denoisers which stem from an explicit potential functional and thus demonstrate symmetric Jacobian by their construction. We then derive an inverse technique that utilizes our denoisers as well as their potentials to ensure convergence. Three important remarks are to be made with regard to work presented below. First, we emphasize that the proposed denoisers are not intended to compete with state-of-the-art denoisers, but rather are developed to exhibit favorable properties that benefit other inversion tasks. Second, our proposed recovery method is not derived to compete with end-to-end learned solutions to specific inverse problems [31, 2, 25, 26]. Rather, our learned denoisers, along with our optimization framework, can be used for all inverse problems (of the form in Eq. (1)) without additional training. Lastly, while we present a present a specific inversion method, the potential-driven denoisers can benefit many RED and PnP algorithms, hence, our work complements these frameworks.

heading: 3 Recovery by Potential-Driven Denoisers

text: In this section, we present our denoisers driven by learned potential functions. Designed as gradient operators, the denoisers display symmetric Jacobians, allowing for a MAP or an MMSE estimator interpretation. More importantly, we use the denoising potentials as an explicit regularization for general inverse problems, leading to provably convergent optimization schemes which utilize our proposed denoisers. Unlike previous works, we do not require bounding or approximation of the Lipschitz constants of the denoisers.

heading: 3.1 Potential-Driven Denoisers

text:
MAP and MMSE Denoisers as Gradients Two fundamental forms of denoisers are MAP and MMSE estimators. The theoretical properties of MAP and MMSE denoisers are well-understood and
they serve as the basis of many works which study the convergence of RED and PnP methods. The following lemmas reveals that MAP and MMSE denosiers share a common structure:
Lemma 1 (Corollary 1 [21]). Assume the denoiser is a proximal operator D , Pφ of some scalarvalued function φ. Then, there exists a convex scalar-valued function ψ(·) such that D = ∇ψ.
Lemma 2 (Tweedie’s Formula [16]). Assume the denoiser is an MMSE estimator D of images contaminated with white Gaussian noise of variance σ2. Then, D , I + σ2∇ log p(·) where I is the identity operator and p(·) denotes a smooth probability density of images.
The above results imply that for any MAP or MMSE denoiser D there exists a scalar-valued function ρ such that the denoiser, or equivalently its residual, can be expressed as D = ∇ρ. This observation sets the main foundation of our gradient-driven denoisers, described next.
Construction of Gradient-Based Denoisers Now, we consider the task of training a deep network to perform denoising. While we focus on image denoising, the following derivations can be generalized to any signal of interest. Motivated by the above observation, we design a denoising deep network which is the gradient of scalar-valued function, acting as a potential function. Three possible approaches to construct gradient-driven denoisers are as follows:
1. Directly design a deep network which is the gradient of a certain potential function. This approach avoids building explicitly the potential network. In the Supplemental we describe certain conditions under which a given network is a gradient.
2. Compose a potential function using an arbitrary deep networkNθ whose input is an image. Simple examples of this approach are 〈x, Nθ(x)〉, ||Nθ(x)||2, ||x −Nθ(x)||2, etc. Note that in general the network Nθ does not equal to the gradient of the composite potential.
3. Model the potential function as scalar-valued deep network and use its gradient as the denoiser.
Here we focus on the last approach since it is the most general. We construct a potential function gθ(·) ∈ R defined as a deep neural network with parameters θ. As illustrated in Fig. 1, the potential network follows a typical structure of feed-forward network with skip-connections except from two main modifications:
• The output of the network is a single scalar. • The acitivations are continuously differentiable and smooth functions.
Mathematically, the potential network consists of M > 0 layers described recursively as follows
zm+1 = sm ( Wmzm + Vmx + bm ) , m = 0, ...,M − 1. (7)
Here, x ∈ Rdx is the input noisy image, zm ∈ Rdm are the hidden variables where z0 ≡ 0 and gθ(x) = zM is the output of the network. The parameters of the network are θ , {bm,Wm, Vm, }M−1m=0 where bm is a bias term and Wm,Vm are linear operators with W0 ≡ 0. The activation functions are denoted by sm : R→ R and applied element-wise. To allow the network to accept inputs of variable sizes, we choose the linear mappings to be standard convolutions. We ensure the output of the network is a scalar by setting the number of output channels of the last convolution to be one, and using global average pooling as the final activation function. Another important design choice is that all activations are continuously differentiable and smooth functions. This creates a potential network that is smooth as the composition of smooth functions, which is crucial for guaranteeing convergence of our inversion method, presented in Section 3.2.
Given the potential function gθ(·), we define our denoiser Dθ(·) as
Dθ(x) , x−∇gθ(x). (8)
Alternatively, we can write or Rθ = ∇gθ where Rθ is the denoising residual. This construction ensures that the resultant denoiser, termed gradient-driven denoising CNN (GraDnCNN), is indeed a gradient of an explicit function, and thus it possesses a symmetric Jacobian. This form of gradientdriven denoisers represents more faithfully MAP and MMSE estimators, hopefully leading to better
approximations. Finally, we perform supervised residual learning and train the denoiser (or equivalently the potential) by minimizing a chosen loss function L(·, ·) with respect to the parameters
θ∗ = argmin θ ∑ n L ( xn − x∗n,Rθ(xn) ) , (9)
where {x∗n, xn} are pairs of clean and noisy images respectively. Special Potentials: Convex Functions and Difference of Convex Functions An additional interesting path we follow is setting all {Wm} (blue blocks in Fig. 1) to be non-negative and all activation functions {sm} (red blocks in Fig. 1) to be convex and non-decreasing mappings. This creates a potential function gθ(·) which is an input convex neural network (ICNN) [4], leading to a denoiser called denoising ICNN (DnICNN). As we explain later in Section 3.2, the use of convex potentials leads to convex inverse problems, thus it facilitates the optimization process and allows to derive recovery methods with global convergence guarantees.
The down side of convex potentials is that the non-negative weights {Wm} may limit the expressivity of the network. Hence, we propose an alternative where the potential network is defined as the difference of two ICNNs gθ(x) , gθ1(x)− gθ2(x). (10) The resultant potential enjoys favorable properties as the difference of two convex functions [28] and it can approximate non-convex functions (unlike the convex potentials). Thus, we present a third denoiser, referred to as denoising by difference of ICNNs (DnDICNN), which acts as a compromise between GraDnCNN and DnICNN. In the following section, we study the properties of these three denoisers in the context of our recovery scheme, both in terms of stability and recovery performance.

heading: 3.2 Convergent Inversion Method

text:
Inversion by Descent Given a trained gradient-driven denoiser, we wish to exploit the denoising power within general recovery processes. The existence of an underlying learned potential gθ(·) allows us to formulate the recovery as the following problem:
x̂ = argmin x∈Rdx f(x) + λgθ(x). (11)
Here we assume f(·) is a continuously differentiable smooth function, typically the log-likelihood of the observations given by (2). The potential network acts as a regularization term or a smooth image prior which creates a well-posed problem. Furthermore, by our construction the potential gθ(·) is smooth, thus, the complete objective, denoted by F (·), is smooth with some unknown Lipschitz constant L > 0. This facilities the optimization process considerably and we can solve (11) via a simple gradient descent algorithm
xk+1 = xk − µk∇F (xk), (12) where we start with some initial solution x0. Notice the above iteration can be written as
xk+1 = xk − µk ( ∇f(xk) + λ ( xk −Dθ(xk) )) (13)
where Dθ(·) is our denoiser originating from gθ(·). Backtracking Step Size Note that the step size µk must be set properly to reach convergence. Typically, one only has access to the denoiser itself, hence the step size is tuned empirically which may lead either to divergence or to slow convergence due to an overly small step size. By contrast, we have complete access to the objective function, thus, we can properly set the step size to ensure convergence. To that end we rely on the Descent Lemma [6] in the Supplemental, which implies that
F (xk+1) ≤ F (xk)− µk 2 ||∇F (xk)||22, ∀µk : 0 ≤ µk ≤ 1 L . (14)
When ∇F (xk) 6= 0 we obtain that F (xk+1) < F (xk); otherwise, we reach a critical point. The challenge that remains is setting the step size according the Lipschitz constant L which is typically unknown. To address this, we set a backtracking step size defined as µk , 1/L̂ where L̂ is an approximation of the true constant L. We start with some initial guess L̂ = L0 and at each iteration update L̂ and µk accordingly until inequality (14) is satisfied, ensuring we perform a descent step.
Algorithm 1: Regularization by Potential-Driven Denoising Input: initial point x0 ∈ Rdx , L0 > 0, β > 1, λ > 0, f(·), gθ(·) and the related denoiser Dθ(·).
1 Set F (·) , f(·) + λgθ(·) and L̂ = L0. 2 for k=0,1,2,... do 3 µk = 1/L̂.
4 xk+1 = xk − µk ( ∇f(xk) + λ ( xk −Dθ(xk) )) .
5 while F (xk+1) > F (xk)− µk2 ||∇F (xk)|| 2 2 do 6 L̂ = β · L̂, µk = 1/L̂. 7 xk+1 = xk − µk ( ∇f(xk) + λ ( xk −Dθ(xk) )) .
The resulting inversion technique, outlined in Algorithm 1, can be viewed as RED-SD with backtracking step size which is enabled by the existence of an available prior function gθ(·). In addition, note that throughout the iterations it holds that L̂ ≤ max(L0, βL) where L is the true unknown Lipschitz constant. Also note that the total number of updates of L̂ is bounded by logβ(L/L0) ≈ O(log2 L). Convergence Guarantees We remark that other more sophisticated optimization techniques can be considered for solving (11), including acceleration schemes. However, rather than providing the most efficient recovery method, we focus here on simplicity to highlight the benefits of learned potential-driven denoisers. Moreover, there is a vast literature which provides convergence guarantees for Algorithm 1 under various conditions. In the following theorem we present a standard convergence analysis based on the smoothness of the objective function. The proof is given in the Supplemental. Theorem 1. Assume that f(·) and gθ(·) are continuously differentiable and smooth functions. Define F (·) , f(·)+λgθ(·) and assume it is bounded below by F ∗. Let {xk}k≥0 be the sequence generated by Algorithm 1. Then,
min m=0,1,...,k−1
||∇F (xm)||22 ≤ 2max(L0, βL)
k
( F (x0)− F ∗ ) .
The above theorem implies that ||∇F (xk)||22 → 0 as k →∞ and that all limit points of the sequence {xk}k≥0 are stationary points of the objective function F (·). Notice this result applies to all the three denoisers (or potentials) presented in Section 3.1. This emphasizes the importance of the smoothness of the potential functions which ensures the convergence of the recovery, regardless to the highly non-convex nature of the objective function.
Commonly, the function f(·) is a convex function, e.g. when it is the log-likelihood. In such scenarios, the use of convex potential networks leads to a special case where the recovery problem is convex and we obtain global convergence guarantees: Corollary 1. Assume that f(·) and gθ(·) are convex in addition to the conditions of Theorem 1. Then, Algorithm 1 converges to global minima of the objective function F (·).
The proof is straightforward since according to Theorem 1 the sequence generated by Algorithm 1 converges to stationary points of the objective function, and any stationary point of a convex function
is a global minimizer. Thus, in the context of inverse problems, our convex potential networks connect the well-known world of convex optimization and with the highly non-linear non-convex nature of deep neural networks.
Finally, for the completeness of our work, we study the case where the recovery problem includes additional constraints on the solution, leading to the following formulation
x̂ = argmin x∈Rdx f(x) + λgθ(x) + h(x), (15)
where the function h(·) is convex yet possibly non-smooth. We provide the full analysis consisting of algorithmic and theoretical results in the Supplemental.

heading: 4 Numerical Experiments

text:
Here we study the performance of Algorithm 1 using our three proposed potential-driven denoisers: GraDnCNN, DnICNN, DnDICNN. We compare ourselves to PnP-PGD and RED-SD, applied with the popular DnCNN denoiser [54], for the tasks of Gaussian deblurring and image super resolution.
Datasets For training the denoising networks for blind Gaussian denoising we use the public DIV2K dataset [2] that consists of a total of 900 high resolution images, 800 for training and 100 for validation. We perform data augmentation by randomly cropping 100 patches of size 50× 50 from each image of the DIV2K training set. We then contaminate each patch with an additive i.i.d. white Gaussian noise (WGN) with a noise level σn chosen uniformly at random from the range [1, 55], creating our 80,000 sample training set. For our validation set, we replicate each image of the DIV2K validation set 10 times, and we add WGN to each copy with random noise level in the specified range. As our test images for evaluating the inverse methods for deblurring and super resolution, we use two different datasets. The first set, provided by the authors of [39], consists of 11 common color images (butterfly, starfish, etc), and we refer to it as Set11. For a more through assessment, we utilize all 500 natural images from the public Berkeley segmentation dataset (BSD500) [32]. All pixel values of the images described above are in the range [0, 255].
Training Denoisers We experiment with four different denoisers. The first is a variant of the DnCNN-B architecture [54] with 7 layers, 128 feature maps, and 3 × 3 convolution kernels. For GraDnCNN we use the architecture shown in Fig. 1, consisting of M = 6 layers, ELU activation functions, 128 feature maps, and 3× 3 bias-free convolutions where the final two convolution layers have a single output channel, and the last activation is global average pooling to produce a scalar. A similar architecture is used for DnICNN, where we set certain weights to be non-negative. The latter architecture is then used to construct our DnDICNN as the difference of two DnICNNs.
Given the datasets detailed above, we train each of the networks using an Adam optimizer for 100 epochs with a constant learning rate of 10−3. For the training loss, we use a modified version of mean squared error (MSE) cost function: ∑ n 1 σ2n MSE ( xn − x∗n,Rθ(xn) )
where x∗n denotes a clean patch, xn is the corresponding noisy patch contaminated with noise of level σn, andRθ(xn) is the output of the network. Since samples with larger noise levels typically lead to larger training errors, we find that this modification of dividing the MSE by the noise level is required to prevent the denoisers from over-smoothing. All experiments are performed in Tensorflow [1] where each model is trained on a single NVIDIA Tesla 32GB V100 GPU.
Inverse Problems We focus here on two common inverse problems: Gaussian deblurring and single image super resolution. For the task of deblurring we perform the following process on the test datasets (Set11, BSD500) described above. Each image is convolved with a 2D Gaussian function 25× 25 with a standard deviation of 1.6. Then, we add to the result an i.i.d. WGN of level σ̂ = √ 2. For super resolution we follow the same procedure except we add an intermediate step, before noise contamination, in which we scale-down the image by a factor of 3 at each axis.
Recovery Given the degraded test images and the trained denoisers, we compare five different inverse methods: PnP-PGD and RED-SD, both with DnCNN, as well as Algorithm 1 employed with our three potential-driven denoisers. For Algorithm 1 we use the constrained update rule, described in the Supplemental, where the proximal operator is the projection of pixel values of the solution into the range [0, 255]. We apply the same projection on each intermediate solution of PnP-PGD
and RED-SD for improved stability. For properly setting the hyperparameters of the different tested methods we use Set11 to perform an extensive hyperparameter search, provided in the Supplemental.
We first study the convergence by running all the techniques for a fixed number of 400 iterations for both inverse tasks, where we measure convergence by ||xk+1 − xk||/µk. The results, averaged over each dataset separately, are given in Figures 2 and 3. The foremost observation is that in all the scenarios PnP-PGD diverges at the a certain point. We conjecture the reason of this instability of PnP-PGD lies in the fact that unlike the other methods PnP-PGD includes no hyperparameters for constraining the denoiser, thus, powerful denoisers may lead to divergence by analogy with a too large learning rate. This can be prevented by considering more sophisticated version of PnP or by performing early stopping, however, it will introduce additional hyperparameter which need to be properly tuned for each data instance, making the technique less robust. Moreover, notice that the convergence of Algorithm 1 does not depend on the values of λ and µ, regardless of the potential-driven denoiser in use. This is a desired property that most RED and PnP methods do not exhibit. We note that the convex potential leads to fastest convergence which is expected, since in this case the objective function is smooth and convex, facilitating the optimization process significantly.
Next, we quantitatively compare the inverse algorithms by computing the average peak signal to noise ratio (PSNR) and structural similarity index (SSIM) of their solutions with respect to the ground truth images. As seen in Table1, Algorithm 1 with GraDnCNN produces the best performance in both tasks and for both test sets. Algorithm 1 with DnICNN yields the worst results among our potential-driven denoisers, yet comparable to PnP-PGD and RED-SD with DnCNN. The performance of Algorithm 1 with DnDICNN is consistent with the fact that DnDICNN presents a compromise between GraDnCNN and DnICNN in terms of the expressivity of the network. To visually support our claims, we provide selected recovered images in Fig. 4, while additional results are given in the Supplemental. An important remark is that in contrast to PnP and RED methods, presented here and many others, the solutions of Algorithm 1 and its performance are completely independed of the choice of the step-size µ. Hence, Algorithm 1 is more robust.

heading: 5 Conclusion

text:
In this paper we developed a provably convergent method with explicit denoising regularization for solving inverse problems. To that end, we first introduced the notion of potential-driven denoisers, i.e. denoisers which are based on the gradient of learned scalar-valued networks. We then utilized our denoisers along with their potential functions to derive an inverse gradient descent method with backtracking step size, ensuring convergence. We demonstrate experimentally that the resulting methods provide improved performance compared to standard RED and PnP methods, while offering provable stability.
The developed framework complements RED and PnP as our potential-driven denoisers can be used within any RED or PnP techniques to guarantee convergence, regardless of the data being used. Furthermore, while our work focus on learning denoisers, the concept of potential-driven networks may be generalized to benefit other image-to-image operations.

title: It Has Potential: Gradient-Driven Denoisers for Convergent Solutions to Inverse Problems