Mario González
email: mgonzalez@unorte.edu.uy

Andrés Almansa
email: andres.almansa@parisdescartes.fr

Pauline Tan
email: pauline.tan@sorbonne-universite.fr

CenUR RN, † Dmel

Solving Inverse Problems by Joint Posterior Maximization with Autoencoding Prior *

Keywords: Image Restoration, Inverse Problems, Bi-convex Optimization, Bayesian Statistics, Generative Models, Variational Auto-encoders AMS subject classifications. 68U10, 65K10, 65D18, 68T05, 90C26, 90C25, 90C30, 1

In this work we address the problem of solving ill-posed inverse problems in imaging where the prior is a variational autoencoder (VAE). Specifically we consider the decoupled case where the prior is trained once and can be reused for many different log-concave degradation models without retraining. Whereas previous MAP-based approaches to this problem lead to highly non-convex optimization algorithms, our approach computes the joint (space-latent) MAP that naturally leads to alternate optimization algorithms and to the use of a stochastic encoder to accelerate computations. The resulting technique (JPMAP) performs Joint Posterior Maximization using an Autoencoding Prior. We show theoretical and experimental evidence that the proposed objective function is quite close to bi-convex. Indeed it satisfies a weak bi-convexity property which is sufficient to guarantee that our optimization scheme converges to a stationary point. We also highlight the importance of correctly training the VAE using a denoising criterion, in order to ensure that the encoder generalizes well to out-of-distribution images, without affecting the quality of the generative model. This simple modification is key to providing robustness to the whole procedure. Finally we show how our joint MAP methodology relates to more common MAP approaches, and we propose a continuation scheme that makes use of our JPMAP algorithm to provide more robust MAP estimates. Experimental results also show the higher quality of the solutions obtained by our JPMAP approach with respect to other non-convex MAP approaches which more often get stuck in spurious local optima.

is log-concave with respect to x. To illustrate this, let us consider the case where the negative log-conditional is quadratic with respect to x (1.1) F (x, y) = 1 2σ 2 ∥Ax -y∥ 2 . This boils down to a linear degradation model that takes into account degradations such as, white Gaussian noise, blur, and missing pixels. When the degradation operator A is noninvertible or ill-conditioned, or when the noise level σ is high, obtaining a good estimate of x requires prior knowledge on the image, given by p X (x) ∝ e -λU (x) . Variational and Bayesian methods in imaging are extensively used to derive MMSE or MAP estimators, based on (a) explicit hand-crafted priors like Tikhonov regularization [START_REF] Tikhonov | On the regularization of ill-posed problems[END_REF], total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF][START_REF] Louchet | Posterior expectation of the total variation model: Properties and experiments[END_REF][START_REF] Pereyra | Proximal Markov chain Monte Carlo algorithms[END_REF] and its higher order [START_REF] Bredies | Total generalized variation[END_REF] and non-local [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF] variants, sparsity in a transformed domain [START_REF] Donoho | De-noising by soft-thresholding[END_REF], or in redundant representations like wavelet frames or patch dictionaries [START_REF] Elad | Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing[END_REF], or (b) learning-based priors like patch-based Gaussian mixture models [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF][START_REF] Yu | Solving inverse problems with piecewise linear estimators: From gaussian mixture models to structured sparsity[END_REF][START_REF] Afonso | Scene-Adapted Plug-and-Play Algorithm with Guaranteed Convergence: Applications to Data Fusion in Imaging[END_REF].

Neural network regression. Since deep neural networks (NN) showed their superiority in image classification tasks [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] researchers started to look for ways to use this tool to solve inverse problems too. The most straightforward attempts employed neural networks as regressors to learn a risk minimizing mapping y → x from many examples (x i , y i ) either agnostically [START_REF] Dong | Learning a deep convolutional network for image super-resolution[END_REF][START_REF] Zhang | Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising[END_REF][START_REF] Zhang | Ffdnet: Toward a fast and flexible solution for cnn-based image denoising[END_REF][START_REF] Gharbi | Deep joint demosaicking and denoising[END_REF][START_REF] Schwartz | Deepisp: Toward learning an end-toend image processing pipeline[END_REF][START_REF] Gao | Dynamic scene deblurring with parameter selective sharing and nested skip connections[END_REF] or including the degradation model in the network architecture via unrolled optimization techniques [START_REF] Gregor | Learning fast approximations of sparse coding[END_REF][START_REF] Chen | Trainable Nonlinear Reaction Diffusion: A Flexible Framework for Fast and Effective Image Restoration[END_REF][START_REF] Diamond | Unrolled optimization with deep priors[END_REF][START_REF] Gilton | Neumann networks for inverse problems in imaging[END_REF].

Implicitly decoupled priors. The main drawback of neural networks regression is that they require to retrain the neural network each time a single parameter of the degradation model changes. To avoid the need for retraining, another family of approaches seek to decouple the NN-based learned image prior from the degradation model. A popular approach within this methodology are Plug & Play (or PnP) methods. Instead of directly learning the logprior -log p X (x) = U (x) + C, these methods seek to learn an approximation of its gradient ∇U [START_REF] Siavash | Deep Mean-Shift Priors for Image Restoration[END_REF][START_REF] Siavash | Image Restoration using Autoencoding Priors[END_REF] or proximal operator prox U [START_REF] Singanallur | Plug-and-Play priors for model based reconstruction[END_REF][START_REF] Meinhardt | Learning proximal operators: Using denoising networks for regularizing inverse imaging problems[END_REF][START_REF] Zhang | Learning Deep CNN Denoiser Prior for Image Restoration[END_REF][START_REF] Chan | Plug-and-play admm for image restoration: Fixed-point convergence and applications[END_REF][START_REF] Ulugbek S Kamilov | A plug-and-play priors approach for solving nonlinear imaging inverse problems[END_REF][START_REF] Ryu | Plug-and-play methods provably converge with properly trained denoisers[END_REF], by replacing it by a denoising NN. Then, these approximations are used in an iterative optimization algorithm to find the corresponding MAP estimator in Equation (1.2) or more generally some sort of consensus equilibrium among the data fitting term and the priors [START_REF] Gregery | Plug-andplay unplugged: Optimization-free reconstruction using consensus equilibrium[END_REF].

Taking an apparently different approach Romano et al. introduced the regularization by denoising (RED) algorithm [START_REF] Romano | The little engine that could: Regularization by denoising (red)[END_REF] which uses a denoiser D σ to construct an explicit regularizer U (x) = 1 2 x T (x -D σ (x)). Under certain conditions (see below) its gradient ∇U = I -D σ can be conveniently computed in terms of the denoiser, leading to a gradient descent scheme for the associated MAP estimator, which is very easy to implement.

Explicitly decoupled generative priors. In another series of works pioneered by Bora et al. [START_REF] Bora | Compressed sensing using generative models[END_REF] and followed by [START_REF] Shah | Solving linear inverse problems using gan priors: An algorithm with provable guarantees[END_REF][START_REF] Raj | GAN-Based Projector for Faster Recovery With Convergence Guarantees in Linear Inverse Problems[END_REF][START_REF] Menon | PULSE: Self-Supervised Photo Upsampling via Latent Space Exploration of Generative Models[END_REF][START_REF] Huang | A Provably Convergent Scheme for Compressive Sensing under Random Generative Priors[END_REF][START_REF] Hand | Global Guarantees for Enforcing Deep Generative Priors by Empirical Risk[END_REF] the Plug & Play prior is explicitly provided by a generative model, most often a generative adversarial network G that maps a latent variable z ∼ N (0, I) to an image x = G(z) with the desired distribution p X as represented by the learning dataset. More precisely these methods solve an optimization problem on the latent variable z and the reconstructed image is provided by xmap = G( ẑmap ). As we show in the following sub-section and in appendix B.1, this corresponds (when α = 1) to the Maximum A Posteriori (MAP) estimator with respect to the z variable. In this work we adopt this framework with some extensions that help avoid getting trapped in spurious critical points of the non-convex objective function.

Empirical success of Plug & Play and RED. Plug & Play and RED approaches became very popular because they allow to repurpose very powerful state of the art denoisers as regularizers of a large family of inverse problems in a quite straightforward manner. They have been successfully applied to many different problems in imaging and they have thus empirically proven their superiority (in terms of achievable reconstruction quality with respect to more classical regularization techniques), and opened the way for the solution of more difficult inverse problems in imaging.

Theoretical questions. The success of Plug & Play and RED approaches largely outpaced our understanding of why and when these techniques lead to algorithms that provably converge to well-posed statistical estimators with well known properties. This is not surprising because obtaining convergence guarantees for non-convex optimization problems under realistic conditions is quite challenging.

A notable exception where strong convergence results have been obtained is the particular case of compressed sensing, where the rows of the degradation operator (or sensing matrix) A are independent realizations of a zero-mean Gaussian distribution. For this problem, Hand et al. [START_REF] Hand | Global Guarantees for Enforcing Deep Generative Priors by Empirical Risk[END_REF][START_REF] Huang | A Provably Convergent Scheme for Compressive Sensing under Random Generative Priors[END_REF] show that the optimization objective (1.3) has almost no spurious stationary points when the generator is assumed to be a random ReLU network with Gaussian weights. As a consequence, a minor modification to the gradient descent algorithm in [START_REF] Bora | Compressed sensing using generative models[END_REF] converges with high probability to the global optimum.

In this paper we are interested in more general inverse problems, where the sensing matrix A is not necessarily random but deterministic and highly structured most often dictated by our modeling of the acquisition device. In this more general setting the hypotheses of the CS results are not necessarily satisfied, and the kind of convergence guarantees that could be established for PnP algorithms with non-convex priors are much weaker (typically only convergence to a stationary point or fixed point is provided, not necessarily a global optimum), and most works concentrate in the implicit case, where the prior is not explicitly provided by a generative model, but implicit in a denoising algorithm.

In such a case the actual prior is unknown, the existence of a density whose gradient or proximal operator is well approximated by a neural denoiser is most often not guaranteed [START_REF] Edward | Regularization by denoising: Clarifications and new interpretations[END_REF], and the convergence of the algorithm is not guaranteed either unless the denoiser is retrained with specific constraints like idempotence [START_REF] Gupta | Cnn-based projected gradient descent for consistent ct image reconstruction[END_REF][START_REF] Shah | Solving linear inverse problems using gan priors: An algorithm with provable guarantees[END_REF], contractive residual [START_REF] Ryu | Plug-and-play methods provably converge with properly trained denoisers[END_REF] or exact, invertible, smooth MMSE denoisers [START_REF] Xu | Provable Convergence of Plug-and-Play Priors with MMSE denoisers[END_REF].

The effect of such training constraints on the quality of the denoisers and the associated priors is yet to be explored in detail. But even when these constraints are satisfied, convergence conditions can be quite restrictive, either (a) requiring the data-fitting term F to be strongly convex [START_REF] Ryu | Plug-and-play methods provably converge with properly trained denoisers[END_REF] (thus excluding many important problems in computational imaging where A is not full rank like interpolation, super-resolution, deconvolution with a non-invertible kernel or compressive sensing), and/or (b) constraining the regularization parameter λ outside of its useful range [START_REF] Ryu | Plug-and-play methods provably converge with properly trained denoisers[END_REF][START_REF] Xu | Provable Convergence of Plug-and-Play Priors with MMSE denoisers[END_REF]. 1Similarly, an early analysis of the RED approach [START_REF] Edward | Regularization by denoising: Clarifications and new interpretations[END_REF] provides a convergence proof, but only under quite restrictive conditions (locally homogeneous denoisers with symmetric Jacobian) which exclude most state of the art denoisers like DnCNN, BM3D, NLMeans. A more recent analysis of a stochastic variant of the RED algorithm [START_REF] Laumont | On Maximum-a-Posteriori estimation with Plug & Play priors and stochastic gradient descent[END_REF] (called PnP-SGD) significantly expands the family of denoisers that provide convergence guarantees, including DnCNN and the doubly-stochastic variant of NLM [START_REF] Sreehari | Plug-and-Play Priors for Bright Field Electron Tomography and Sparse Interpolation[END_REF]. These guarantees come, however, at the expense of a very small descent step which leads to a very computationally expensive algorithm with slow convergence. In addition, the experiments with PnP-SGD show that this algorithm is extremely sensitive to the initial condition, and it can be easily get stuck on spurious local minima if not initialized very carefully.

Focus of this work. Very recent works focused on developing MAP estimation algorithms with convergence guarantees under more realistic conditions. The convergence analysis of the RED framework, and its RED-PRO variant was further refined under a demicontractive condition for the denoiser [START_REF] Cohen | Regularization by Denoising via Fixed-Point Projection (RED-PRO)[END_REF]. This condition is, however, difficult to verify according to Pesquet et al. [START_REF] Pesquet | Learning Maximally Monotone Operators for Image Recovery[END_REF] who provides an alternative convergence analysis framework based on firmly non-expansive denoisers for which explicit training procedures exist [START_REF] Terris | Building firmly nonexpansive convolutional neural networks[END_REF]. In this work we explore alternative new ways to bring theory and practice closer together, by proposing novel Plug & Play algorithms to compute the MAP estimator of an inverse problem with a neural regularizer. Unlike previous approaches which were based on implicit priors, or on GAN-based explicit priors, our approach is based on an explicit generative prior that has been trained as a Variational AutoEncoder (VAE). As we shall see later, the additional VAE structure provides: (i) powerful mechanisms to avoid getting stuck in spurious local minima of the associated non-convex functional, and (ii) convergence guarantees under much less restrictive conditions on the inverse problem F and regularization parameter λ.

The next Section 1.1 reviews previous work on similar approaches to compute a MAP estimator from a generative prior that was trained either as a VAE or a GAN. Section 1.2 briefly introduces our approach and how it relates to previous work. The section finishes with an overview of the rest of the paper.

1.1. Maximum a Posteriori meets Generative Models. Our approach focuses on PnP algorithms where the prior is provided by a generative model. For instance one could use a generative adversarial network (GAN) to learn a generative model for X = G(Z) with Z ∼ N (0, I) a latent variable. The generative model induces a prior on X via the pushforward measure p X = G♯p Z , which following [48, section 5] can be developed as

p X (x) = p Z G -1 (x) det S(G -1 (x)) δ M (x)
where S = ∂G ∂z T ∂G ∂z is the squared Jacobian and the manifold M = {x : ∃z, x = G(z)} represents the image of the generator G. With such a prior p X , the x-optimization (1.2) required to obtain xmap becomes intractable (in general), for various reasons:

• the computation of det S,

• the inversion of G, and • the hard constraint x ∈ M. These operations are all memory and/or computationally intensive, except when they are partially addressed by the use of a normalizing flow like in [START_REF] Helminger | Blind Image Restoration with Flow Based Priors[END_REF][START_REF] Whang | Compressed Sensing with Invertible Generative Models and Dependent Noise[END_REF].

Current attempts to use such a generative model as a prior, like the one proposed by Bora et al. [START_REF] Bora | Compressed sensing using generative models[END_REF] for GANs, circumvent these difficulties by performing an optimization on z (in the latent domain) instead of x. Instead of solving Equation (1.2), they solve

ẑmap = arg max z p Y |X (y | G(z)) p Z (z) = arg min z F (G(z), y) + 1 2 ∥z∥ 2 , (1.4)
by assuming a standard Gaussian prior. This problem is much more tractable, and the corresponding x-estimate is obtained as

(1.5) xmap-z = G( ẑmap ).
As we show in appendix B.1, this new estimator does not necessarily coincide with xmap but it does correspond to the MAP-estimator of x after the change of variable x = G(z),

namely xmap-z = G arg max z p Z|Y (z | y) .
Since G is non-linear, this problem (or its equivalent formulation (1.4)) is highly non-convex and difficult to solve with global optimality guarantees. Nevertheless, in the particular case where A is a random Gaussian matrix (compressed sensing case) or when F is strongly convex, recent work shows that the global optimum can be reached with linear convergence rates by a small modification of a gradient descent algorithm [START_REF] Huang | A Provably Convergent Scheme for Compressive Sensing under Random Generative Priors[END_REF][START_REF] Hand | Global Guarantees for Enforcing Deep Generative Priors by Empirical Risk[END_REF], or by an ADMM algorithm with non-linear constraints [START_REF] Latorre | Fast and Provable ADMM for Learning with Generative Priors[END_REF][START_REF] Benning | Preconditioned ADMM with nonlinear operator constraint[END_REF][START_REF] Valkonen | A primal-dual hybrid gradient method for nonlinear operators with applications to MRI[END_REF]. To the best of our knowledge, these results do not extend, however, to the more general case we are interested in here, where A is deterministic and rank-deficient, and F is consequently not strongly convex. In this more general setting, convergence guarantees for this optimization problem remain extremely difficult to establish, as confirmed by experimental results presented in Section 3.

A common technique to solve difficult optimization problems like the one in Equation (1.4) is to use (Half Quadratic) splitting methods (1.6) xβ = arg min

x min z F (x, y) + β 2 ∥x -G(z)∥ 2 + 1 2 ∥z∥ 2 J 1,β (x,z)
combined with a continuation scheme, namely:

(1.7) xmap-z = lim β→∞ xβ .
The convergence of the continuation scheme in the last line is a standard result in Γ-convergence (see [START_REF] Dal | An Introduction to Γ-Convergence[END_REF] and appendix C). The corresponding splitting algorithm is presented in Algorithm 1.1.

Algorithm

β := β k 3:
for n := 0 to maxiter do 4:

z n+1 := arg min z J 1,β (x n , z) // Nonconvex 5:
x n+1 := arg min x J 1,β (x, z n+1 ) // Quadratic 6:

end for 7:

x 0 := x n+1 8: end for 9: return x n+1 However, unlike most cases of HQS which include a linear constraint between the two variables, this splitting algorithm still contains (line 4) a difficult non-convex optimization problem2 .

1.2. Proposed method: Joint MAP x,z . In this work we propose to address this challenge by substituting the difficult non-convex sub-problem by a local quadratic approximation provided by the encoder of a variational autoencoder.

Indeed, as we show in Section 2, a variational autoencoder allows to interpret the splitting Equation (1.6) as the negative logarithm of the joint posterior density p X,Z|Y (x, z | y). Therefore, solving Equation (1.6) amounts to compute a joint map x,z estimator that we denote by xβ mapx,z . Moreover if the same joint conditional density p X,Z|Y (x, z | y) is decomposed in a different manner, it leads to an approximate expression that makes use of the encoder, and is quadratic in z. If this approximation is good enough then the maximization of the joint log-posterior becomes a bi-concave optimization problem or approximately so. And in that case, an extension of standard bi-convex optimization results [START_REF] Gorski | Biconvex sets and optimization with biconvex functions: a survey and extensions[END_REF] shows that the algorithm converges to a stationary point.

We also highlight the importance of correctly training the VAE in such a way that the encoder generalizes well to noisy values of x outside of the support of p X (x). This can be achieved by training the VAE to reconstruct their clean inputs with noise injected at the input level, as proposed by Im et al. [START_REF] Daniel | Denoising criterion for variational auto-encoding framework[END_REF]. We observe that this modified training does not degrade the quality of the generative model, but makes our quasi-bi-convex optimization procedure much more robust.

Finally we show that a continuation scheme allows to obtain the map z estimator as the limit of a series of joint map x,z optimizations. This continuation scheme, in addition to the quasi-bi-convex optimization, and the initialisation heuristic provided by the denoising encoder leads to a much more robust non-convex optimization scheme which more often converges to the right critical point than a straightforward gradient descent of the map z model.

The remainder of this paper is organized as follows. In Section 2 we derive a model for the joint conditional posterior distribution of space and latent variables x and z, given the observation y. This model makes use of a generative model, more precisely a VAE with Gaussian decoder. We then propose an alternate optimization scheme to maximize the joint posterior model, and state convergence guarantees. Section 3 presents first a set of experiments that illustrates the convergence properties of the optimization scheme. We then test our approach on classical image inverse problems, and compare its performance with state-of-theart methods. Concluding remarks are presented in Section 4.

2. From Variational Autoencoders to Joint Posterior Maximization. Recently, some generative models based on neural networks have shown their capability to approximate the complex image distribution in a data-driven fashion. In particular, Variational Autoencoders (VAE) [START_REF] Diederik | Auto-Encoding Variational Bayes[END_REF] combine variational inference to approximate unknown posterior distributions of latent variable models with the ability of neural networks to learn such approximations.

Consider a graphical model z → x in which we assume that a latent variable z is responsible of the observed image x. For example, in an image of a handwritten digit we can imagine which digit is represented in the image, width, angle (and so on) as latent variables. We choose a generative model

p θ (x, z) = p θ (x|z)p Z (z)
where p Z (z) is some simple distribution (which we can easily sample from) and p θ (x|z) is the approximation of the probability distribution of x given z parameterized by a neural network (with weights θ) known as stochastic decoder.

The intractability of p θ (x) = p θ (x|z)p Z (z) dz is related to the posterior distribution p θ (z|x) by (2.1) p θ (z|x) = p θ (x|z)p Z (z) p θ (x) .

The variational inference approach consists in approximating this posterior with another model q ϕ (z|x) which, in our case, is another neural network with parameters ϕ, called a stochastic encoder.

Following [START_REF] Diederik | Auto-Encoding Variational Bayes[END_REF], we consider the Evidence Lower BOund (ELBO) as

(2.2) L θ,ϕ (x) := log p θ (x) -KL(q ϕ (z|x) || p θ (z|x)) ≤ log p θ (x)
where KL is the Kullback-Leibler divergence. Thus, given a dataset D = {x 1 , . . . , x N } of image samples, maximizing the averaged ELBO on D means maximizing log p θ (D) which is the maximum likelihood estimator of weights θ and minimizing KL(q ϕ (z|x) || p θ (z|x)) which enforces the approximated posterior q ϕ (z|x) to be similar to the true posterior p θ (z|x).

It can be shown [START_REF] Diederik | Auto-Encoding Variational Bayes[END_REF] that the ELBO can be rewritten as

(2.3) L θ,ϕ (x) = E q ϕ (z|x) [log p θ (x|z)] -KL(q ϕ (z|x) || p Z (z)).
The first term in (2.3) is a reconstruction loss similar to the one of plain autoencoders: it enforces that the code z ∼ q ϕ (•|x) generated by the encoder q ϕ can be used by the decoder p θ to reconstruct the original input x. The second term is a regularization term that enforces the distribution q ϕ (z|x) of the latent code z (given x) to be close to the prior distribution p Z (z). It is common to choose an isotropic Gaussian as the prior distribution of the latent code:

p Z (z) = N (z | 0, I) ∝ e -∥z∥ 2 /2
and a Gaussian encoder q ϕ (z|x) = N (z | µ ϕ (x), Σ ϕ (x)), so that the KL divergence in (2.3) is straightforward to compute. For the decoder p θ (x|z) a Gaussian decoder is the most common choice and as we will see we benefit from that.

2.1. Learning approximations vs. encoder approximations. In this work we construct an image prior using a Variational Autoencoder (VAE). Like any machine learning tool VAEs make different kinds of approximations. Let's distinguish two types of approximations that shall be important in the sequel: Learning approximation: The ideal prior p * X can only be approximated by our VAE due to its architectural constraints, finite complexity, truncated optimization algorithms, finite amount of data and possible biases in the data. Due to all these approximations, after learning we have only access to an approximate prior p X ≈ p * X . VAEs give access to this approximate prior p X via a generative model: taking samples of a latent variable Z with known distribution N (0, I) in R l (with l ≪ d), and feeding these samples through a learned decoder network, we obtain samples of X ∼ p X . The approximate prior itself

(2.4) p X (x) = p θ (x|z) p Z (z) dz
is intractable because it requires computing an integral over all possible latent codes z. However the approximate joint distribution is readily accessible

p X,Z (x, z) = p θ (x|z) p Z (z) thanks to p X|Z (x | z) = p θ (x|z)
which is provided by the decoder network. Encoder approximation: In the previous item we considered the VAE as a generative model without making use of the encoder network. The encoder network

pZ|X (z | x) := q ϕ (z|x) ≈ p Z|X (z | x)
is introduced as an approximate way to solve the intractability of p Z|X (z | x) = p θ (z|x) (which is related to the intractability of p θ (x) as observed in equation (2.1)).

Using the encoder network we can provide an alternative approximation for the joint distribution pX,Z (x, z)

:= q ϕ (z|x) p X (x) ≈ p X,Z (x, z)
which shall be useful in the sequel.

Put another way, the ideal joint distribution p * X,Z is inaccessible, but can be approximated in two different ways:

The first expression denoted p X,Z (x, z) only uses the decoder and is only affected by the learning approximation

p * X,Z (x, z) ≈ p X,Z (x, z) := p θ (x|z) p Z (z) .
The second expression denoted pX,Z (x, z) uses both encoder and decoder and is affected both by the learning approximation and by the encoder approximation

p X,Z (x, z) ≈ pX,Z (x, z) := q ϕ (z|x) p X (x)
In the following subsection we shall forget about the ideal prior p * X and joint distribution p * X,Z which are both inaccessible. Instead we accept p X (with its learning approximations) as our prior model which shall guide all our estimations. The approximation symbol shall be reserved to expressions that are affected by the encoder approximation in addition to the learning approximation.

Variational Autoencoders as Image Priors.

To obtain the Maximum a Posteriori estimator (MAP), we could plug in the approximate prior p X in equation (1.2), but this leads to a numerically difficult problem to solve due to the intractability of p X . Instead, we propose to maximize the joint posterior p X,Z|Y (x, z | y) over (x, z) which is equivalent to minimizing (2.5)

J 1 (x, z) := -log p X,Z|Y (x, z | y) = -log p Y |X,Z (y | x, z) p θ (x | z)p Z (z) = F (x, y) + H θ (x, z) + 1 2 ∥z∥ 2 .
Note that the first term is quadratic in x (assuming (1.1)), the third term is quadratic in z and all the difficulty lies in the coupling term H θ (x, z) = -log p θ (x | z). For Gaussian decoders [START_REF] Diederik | Auto-Encoding Variational Bayes[END_REF], the latter can be written as

H θ (x, z) = 1 2 d log(2π) + log det Σ θ (z) + ∥Σ -1/2 θ (z)(x -µ θ (z))∥ 2 . (2.6)
which is also convex in x. Hence, minimization with respect to x takes the convenient closed form:

arg min x J 1 (x, z) = A T A + σ 2 Σ -1 θ (z) -1 × A T y + σ 2 Σ -1 θ (z)µ θ (z) .
(2.7)

Unfortunately the coupling term H and hence J 1 is a priori non-convex in z. As a consequence the z-minimization problem

(2.8) arg min z J 1 (x, z)
is a priori more difficult. However, for Gaussian encoders, VAEs provide an approximate expression for this coupling term which is quadratic in z. Indeed, given the equivalence

p θ (x | z) p Z (z) = p X,Z (x, z) = p Z|X (z | x) p X (x) ≈ q ϕ (z | x) p X (x) (2.9)
we have that (2.10)

H θ (x, z) + 1 2 ∥z∥ 2 ≈ K ϕ (x, z) -log p X (x) .
where

K ϕ (x, z) = -log q ϕ (z | x).
Therefore, this new coupling term becomes

K ϕ (x, z) = -log N (z; µ ϕ (x), Σ ϕ (x)) = 1 2 l log(2π) + log det Σ ϕ (x) + ∥Σ -1/2 ϕ (x)(z -µ ϕ (x))∥ 2 ,
which is quadratic in z. This provides an approximate expression for the energy (2.5) that we want to minimize, namely

(2.11) J 2 (x, z) := F (x, y) + K ϕ (x, z) -log p X (x) ≈ J 1 (x, z).
This approximate functional is quadratic in z, and minimization with respect to this variable yields (2.12) arg min

z J 2 (x, z) = µ ϕ (x).
In the case of linear VAEs,

p Z (z) = N (z; 0, I) (2.13) p X|Z (x | z) = N (x; V θ z + v θ , Σ θ ). (2.14)
It is easily shown that the posterior is also Gaussian [START_REF] Lucas | Don't blame the ELBO! A linear VAE perspective on posterior collapse[END_REF], namely

(2.15) p Z|X (z | x) = N (z; M V T θ (x -v θ ), Σ θ M ) where M = (V T θ V θ + Σ θ ) -1 .
Hence, the linear encoder q ϕ (z | x) that minimizes the ELBO is that of equation (2.15) so the approximation (2.9) is exact and then J 1 = J 2 .

Alternate Joint Posterior Maximization.

The previous observations suggest to adopt an alternate scheme to minimize -log p X,Z|Y (x, z | y) in order to solve the inverse problem. We begin our presentation by a simple version of the proposed algorithm, which aims at managing the case where the approximation of J 1 by J 2 is exact (at least in the sense given in Assumption 1 below); then we propose an adaptation for the more realistic non-exact case and we explore its convergence properties.

When J 1 = J 2 is a bi-convex function, Algorithm 2.1 is known as Alternate Convex Search. Its behavior has been studied in [START_REF] Gorski | Biconvex sets and optimization with biconvex functions: a survey and extensions[END_REF][START_REF] Aguerrebere | A Bayesian Hyperprior Approach for Joint Image Denoising and Interpolation, With an Application to HDR Imaging[END_REF]. Here we shall consider the following (strong) assumption, which includes the strictly bi-convex case (J 1 = J 2 ):

Assumption 1. For any x, if z * is a global minimizer of J 2 (x, •), then z * is a global mini- mizer of J 1 (x, •).
The proposed alternate minimization takes the simple and fast form depicted in Algorithm 2.1, which can be shown to converge to a stationary point of J 1 under Assumptions 1 and 2, as stated in Proposition 2.1 below. Note that the minimization in step 2 of Algorithm 2.1 does not require the knowledge of the unknown term -log p X (x) in Equation (2.11) since it does not depend on z. z n+1 := arg min z J 2 (x n , z) // Quadratic approx 3:

x n+1 := arg min x J 1 (x, z n+1 ) // Quadratic 4: end for 5: return x n+1 , z n+1

The convergence analysis of the proposed schemes requires some general assumptions on the functions J 1 and J 2 : Assumption 2. J 1 (•, z) is convex and admits a unique minimizer for any z. Moreover, J 1 is coercive and continuously differentiable.

The unicity of the minimizers of the partial function J 1 (•, z) can be dropped. In this case, the proof of the convergence of Algorithm 2.1 has to be slightly adapted.

The convergence property of Algorithm 2.1 will be investigated in a wider framework below (Proposition 2.1). Note that all the properties required in Assumption 2 are satisfied if we use a differentiable activation function like the Exponential Linear Unit (ELU) [START_REF] Clevert | Fast and Accurate Deep Network Learning by Exponential Linear Units (ELUs)[END_REF] with α = 1, instead of the more common ReLU activation function. More details can be found in Appendix A.

Approximate Alternate Joint Posterior Maximization.

When the autoencoder approximation in (2.11) is not exact (Assumption 1), the energy we want to minimize in Algorithm 2.1, namely J 1 may not decrease. To ensure the decay, some additional steps can be added. Noting that the approximation provided by J 2 provides a fast and accurate heuristic to initialize the minimization of J 1 , an alternative scheme is proposed in Algorithm 2.2.

In Algorithm 2.2, gd is a gradient descent scheme such that for any starting point z 0 , the output

z + satisfies ∂J 1 ∂z (x, z + ) = 0 and J 1 (x, z + ) ≤ J 1 (x, z 0 )
Hence, one can consider for instance a gradient descent scheme which finds a local minimizer of J 1 (x, •) starting from z 0 . Our experiments with Algorithm 2.2 (Section 3.5) show that during the first few iterations (where the approximation provided by J 2 is good enough) z 1 and z 2 reach convergence faster than z 3 . After a critical number of iterations the opposite is true (the initialization provided by the previous iteration is better than the J 2 approximation, and z 3 converges faster).

Algorithm 2.2 Joint posterior maximization -approximate case

Require: Measurements y, Autoencoder parameters θ, ϕ, Initial conditions x 0 , z 0 Ensure: x, ẑ = arg max x,z p X,Z|Y (x, z | y)

1: for n := 0 to maxiter do 2:

z 1 := arg min z J 2 (x n , z) // Equation (2.12)

3:

z 2 := gd z J 1 (x n , z), starting from z = z 1 4:
z 3 := gd z J 1 (x n , z), starting from z = z n 5:

for i := 1 to 3 do 6:

x i := arg min x J 1 (x, z i ) // Equation (2.7)

7:
end for 8:

i * := arg min i∈{1,2,3} J 1 (x i , z i ) 9:
(x n+1 , z n+1 ) := (x i * , z i * ) 10: end for 11: return x n+1 , z n+1

These observations suggest that a faster execution, with the same convergence properties, can be achieved by the variant in Algorithm 2.3, which avoids the costly computation of z 2 and z 3 when unnecessary. Hence, in practice, we will use Algorithm 2.3 rather than Algorithm 2.2. However, Algorithm 2.2 provides a useful tool for diagnostics. Indeed, the comparison of the evaluation of J 1 (x i , z i ) for i = 1, 2, 3 performed in step 8 permits to assess the evolution of the approximation of J 1 by J 2 .

Algorithm 2.3 is still quite fast when J 2 provides a sufficiently good approximation, since in that case the algorithm chooses i * = 1, and avoids any call to the iterative gradient descent algorithm. Even if we cannot give a precise definition of what sufficiently good means, the sample comparison of K ϕ and H θ as functions of z, displayed in Figure 3(a), shows that the approximation is fair enough in the sense that it preserves the global structure of J 1 . The same behavior was observed for a large number of random tests.

Note that Algorithm 2.1 is a particular instance of Algorithm 2.3 in the case where Assumption 1 holds, and n 1 = n 2 = 0 and if grad descent gives a global minimizer of the considered function (in this case, the computation of z 1 , z 2 , are skipped and only z 3 is computed). Proposition 2.1 (Convergence of Algorithm 2.3). Let {(x n , z n )} be a sequence generated by Algorithm 2.3. Under Assumption 2 we have that:

1. The sequence {J 1 (x n , z n )} converges monotonically when n → ∞.

2. The sequence {(x n , z n )} has at least one accumulation point.

3. All accumulation points of {(x n , z n )} are stationary points of J 1 and they all have the same function value.

Proof. Since we are interested in the behaviour for n → ∞, we assume n > n 2 in Algorithm 2.3.

1. Since n > n 2 the algorithm chooses i * = 3 and z n+1 = z 3 . According to the definition of grad descent, one has

J 1 (x n , z n+1 ) ≤ J 1 (x n , z n ) Algorithm 2.

Joint posterior maximization -approximate case (faster version)

Require: Measurements y, Autoencoder parameters θ, ϕ, Initial condition x 0 , iterations n 1 ≤ n 2 ≤ n max Ensure: x, ẑ = arg max x,z p X,Z|Y (x, z | y)

1: for n := 0 to n max do 2:

done := FALSE 3:

if n < n 1 then 4:

z 1 := arg min z J 2 (x n , z) // Equation (2.12)

5:

x 1 := arg min x J 1 (x, z 1 ) // Equation (2.7) if not done and n < n 2 then 12:

6: if J 1 (x 1 , z 1 ) < J 1 (x n , z n ) then 7: i * := 1 // J 2 is
z 1 := arg min z J 2 (x n , z) 13:
z 2 := gd z J 1 (x n , z), starting from z = z 1 14:

x 2 := arg min x J 1 (x, z 2 ) // Equation (2.7) if not done then 21:

15: if J 1 (x 2 , z 2 ) < J 1 (x n , z n ) then 16: i * := 2 // J 2 init
z 3 := gd z J 1 (x n , z), starting from z = z n 22:

x 3 := arg min x J 1 (x, z 3 ) // Equation (2.7)

23:

i * := 3

24:

end if

25:

(x n+1 , z n+1 ) := (x i * , z i * ) 26: end for 27: return x n+1 , z n+1 and by optimality one has

J 1 (x n+1 , z n+1 ) ≤ J 1 (x n , z n+1 ).
Hence, since J 1 is coercive (thus, lowerbounded), Statement 1 is straightforward.

2. Thanks to the coercivity of J 1 , the sequences {(x n , z n )} and {(x n , z n+1 )} are bounded, thus admit an accumulation point.

3. Using Fermat's rule and the definition of grad descent, one has

∂J 1 ∂z (x n , z n+1 ) = 0 and ∂J 1 ∂x (x n+1 , z n+1 ) = 0.
Let (x * , z * ) be an accumulation point of {(x n , z n )}. By double extraction, one can find two subsequences such that (x n j +1 , z n j +1 ) → (x * , z * ) and (x n j , z n j +1 ) → ( x * , z * ) By continuity of ∇J 1 , one gets that ∂J 1 ∂z ( x * , z * ) = 0 and

∂J 1 ∂x (x * , z * ) = 0
In particular, the convexity of J 1 (•, z * ) and Assumption 2 ensure that x * is a global minimizer of J 1 (•, z * ). Besides, the inequalities proved in Point 1 above show that

J 1 (x * , z * ) = J 1 ( x * , z * ) = lim n→∞ J 1 (x n , z n )
that is, x * is also a global minimizer of J 1 (•, z * ). Since J 1 (•, z * ) has a unique minimizer, one has x * = x * , and

∂J 1 ∂z (x * , z * ) = 0 namely (x * , z *
) is a stationary point of J 1 . Note that we have also proved that x n j and x n j +1 have same limit.

Remark 2.2. Note that if n 1 = n 2 = ∞ we cannot assume that i * = 3. In that case statements 1 and 2 are still valid but the third statement is not. The reason is that for i * ∈ {1, 2} we cannot guarantee the chain of inequalities

J 1 (x n+1 , z n+1 ) ≤ J 1 (x n , z n+1 ) ≤ J 1 (x n , z n ) but only J 1 (x n+1 , z n+1 ) ≤ J 1 (x n , z n ).
This is consistent with the design of the algorithm where iterations n < n 2 serve as an heuristic to guide the algorithm to a sensible critical point. However, convergence to a critical point is only guaranteed by the final iterations n > n 2 .

2.5. MAP-z as the limit case for β → ∞. If one wishes to compute the map-z estimator instead of the joint map-x-z from the previous section, one has two options:

1. Use your favorite gradient descent algorithm to solve equation (1.4). 2. Use Algorithm 2.3 to solve a series of joint map-x-z problems with increasing values of β → ∞ as suggested in Algorithm 1.1. In the experimental section we show that the second approach most often leads to a better optimum.

In practice, in order to provide a stopping criterion for Algorithm 1.1 and to make a sensible choice of β-values we reformulate Algorithm 1.1 as a constrained optimization problem arg min

x,z : ∥G(z)-x∥ 2 ≤ε F (x, y) + 1 2 ∥z∥ 2 .
The corresponding Lagrangian form is

(2.16) max β min x,z F (x, y) + 1 2 ∥z∥ 2 + β ∥G(z) -x∥ 2 -ε +
and we use the exponential multiplier method [START_REF] Tseng | On the convergence of the exponential multiplier method for convex programming[END_REF] to guide the search for the optimal value of β (see Algorithm 2.4)

Algorithm 2.4 map-z as the limit of joint map-x-z. Require: Measurements y, Tolerance ε, Rate ρ > 0, Initial β 0 , Initial x 0 , Iterations 0 ≤ n 1 ≤ n 2 ≤ n max Ensure: arg min x, z : ∥G(z)-x∥ 2 ≤ε F (x, y) + 1 2 ∥z∥ 2 . 1: β := β 0 2: x 0 , z 0 := Algorithm 2.3 starting from x = x 0 with β, n 1 , n 2 , n max . 3: converged := FALSE 4: k := 0 5: while not converged do 6:

x k+1 , z k+1 := Algorithm 2.3 starting from x = x k with β and n 1 = n 2 = 0 7: 3.1. Baseline algorithms. To validate our approach, we perform comparisons on several inverse problems with the following algorithms:

C = ∥G(z k+1 ) -x k+1 ∥ 2 -ε 8: β := β exp(ρC)
• CSGM (Bora et al. [START_REF] Bora | Compressed sensing using generative models[END_REF]) directly computes the zmap estimator as defined in Equation (1.4) using gradient descent. We run CSGM using the decoder of a VAE as generator G starting at random z 0 . In addition, as Bora et al. note that random restarts are important for good performance, we also compute the best result (as measured by (1.4)) among m = 10 different random initializations z 0 and refer to this variant as mCSGM.

• PULSE [START_REF] Menon | PULSE: Self-Supervised Photo Upsampling via Latent Space Exploration of Generative Models[END_REF] is very similar to CSGM but restricts the search of the latent code z to the sphere of radius √ l, arguing that it concentrates most of the probability mass of a Gaussian distribution N (0, I) on a high-dimensional space R l .

• PGD-GAN [START_REF] Shah | Solving linear inverse problems using gan priors: An algorithm with provable guarantees[END_REF] performs a projected gradient descent of F (x, y) wrt x:

w k = x k -ηA T (Ax k -y) x k+1 = G(arg min z ∥w k -G(z)∥). (3.1)
• In addition, we implement the splitting method of Algorithm 1.1 which is a simple continuation scheme for the zmap estimator of Equation (1.4). For a fair comparison we run all algorithms on the same prior, i.e. the same generator network G = µ θ where µ θ is the decoder mean from the VAE model that we trained for JPMAP.

Inverse problems.

Here, we briefly describe the inverse problems y = Ax + η, η ∼ N (0, σ 2 I) to be considered for validating our approach:

• Denoising: A = I and σ large.

• Compressed Sensing: the sensing matrix A ∈ R q×d has Gaussian random entries A ij ∼ N (0, 1/q), where q ≪ d is the number of measurements.

• Interpolation: A is a diagonal matrix with random binary entries, so masking a percentage p of the image pixels. • (Non-blind) Deblurring: Ax = h * x where h is a known convolution kernel.

• Super-resolution: A is a downsampling/decimation operator of scaling factor s.

AutoEncoder and dataset.

In order to test our joint prior maximization model we first train a Variational Autoencoder like in [START_REF] Diederik | Auto-Encoding Variational Bayes[END_REF] on the training data of MNIST handwritten digits [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF].

The stochastic encoder takes as input an image x of 28 × 28 = 784 pixels and produces as an output the mean and (diagonal) covariance matrix of the Gaussian distribution q ϕ (z|x), where the latent variable z has dimension 8. The architecture of the encoder is composed of 3 fully connected layers with ELU activations (to preserve continuous differentiability). The sizes of the layers are as follows: 784 → 500 → 500 → (8 + 8). Note that the output is of size 8 + 8 in order to encode the mean and diagonal covariance matrix, both of size 8.

The stochastic decoder takes as an input the latent variable z and outputs the mean and covariance matrix of the Gaussian distribution p θ (x|z). Following [START_REF] Dai | Diagnosing and Enhancing VAE Models[END_REF] we chose here an isotropic covariance Σ θ (z) = γ 2 I where γ > 0 is trained, but independent of z. This choice simplifies the minimization problem (2.8), because the term det Σ θ (z) (being constant) has no effect on the z-minimization. The architecture of the decoder is also composed of 3 fully connected layers with ELU activations (to preserve continuous differentiability). The sizes of the layers are as follows: 8 → 500 → 500 → 784. Note that the covariance matrix is constant, so it does not augment the size of the output layer which is still 784 = 28 × 28 pixels.

We also trained a VAE on CelebA [START_REF] Liu | Deep learning face attributes in the wild[END_REF] images cropped to 64 × 64 × 3, with latent dimension ranging from 64 to 512. We choose a DCGAN-like [START_REF] Radford | Unsupervised representation learning with deep convolutional generative adversarial networks[END_REF] CNN architecture as encoder and a symmetrical one as decoder with ELU activations, batch normalization and isotropic covariance as before. For more details, see the code3 .

We train these architectures using PyTorch [START_REF] Paszke | Automatic differentiation in pytorch[END_REF] with batch size 128 and Adam algorithm for 200 epochs with learning rate 0.0001 and rest of the parameters as default.

3.4. Need to train the VAE with a denoising criterion. It should be noted that when training our Variational Autoencoder we should be more careful than usual. Indeed in the most widespread applications of VAEs they are only used as a generative model or as a way to interpolate between images that are close to M, i.e. the image of the generator µ θ . For such applications it is sufficient to train the encoder µ ϕ , Σ ϕ on a training set that is restricted to M.

In our case however, we need the encoder to provide sensible values even when its input x is quite far away from M: the encoder has to actually fulfill two functions at the same time:

1. (Approximately) project x to its closest point in M, and 2. compute the encoding of this projected value (which should be the same as the encoding of the original x. Evaluating the effectiveness of JPMAP vs CGSM as a function of σDVAE (same setup of Figure 1). Without a denoising criterion σDVAE = 0 the JPMAP algorithm may provide wrong guesses z 1 when applying the encoder in step 2 of Algorithm 2.2. For σDVAE > 0 however, the alternating minimization algorithm can benefit from the robust initialization heuristics provided by the encoder, and it consistently converges to a better local optimum than the simple gradient descent in CSGM.

Traditional VAE training procedures do not ensure that the encoder generalizes well to x ̸ ∈ M. In order to ensure this generalization ability we adopt the training procedure of the DVAE (Denoising VAE) proposed by Im et al. [START_REF] Daniel | Denoising criterion for variational auto-encoding framework[END_REF], which consists in adding various realizations of zero-mean Gaussian noise of variance σ 2 DVAE to the samples x presented to the encoder, while still requiring the decoder to match the noiseless value, i.e. we optimize the parameters in such a way that

(3.2) µ θ (µ ϕ ( x)) ≈ x
where x = x + σ DVAE ε and ε ∼ N (0, I) for all x in the training set and for many realizations of ε.

More specifically, if we take a corruption model p( x | x) like above, it can be shown [START_REF] Daniel | Denoising criterion for variational auto-encoding framework[END_REF]] that is an alternative ELBO of (2.3). In practice, using Monte Carlo for estimating the expectation E p( x|x) in (3.3), we only need to add noise to x before passing it to the encoder q ϕ during training, as mentioned in (3.2).

(3.3) Lθ,ϕ (x) = E p( x|x) E q ϕ (z| x) [log p θ (x|z)] -KL(q ϕ (z| x) || p Z (z)) (
Our experiments with this denoising criterion confirm the observation by Im et al. [START_REF] Daniel | Denoising criterion for variational auto-encoding framework[END_REF] that it does not degrade the quality of the generative model, as long as σ DVAE is not too large (see Figure 1). As a side benefit, however, we obtain a more robust encoder that generalizes well for values of x that are not in M but within a neighbourhood of size ≈ σ DVAE around M. This side benefit, which was not the original intention of the DVAE training algorithm in [START_REF] Daniel | Denoising criterion for variational auto-encoding framework[END_REF] is nevertheless crucial for the success of our algorithm as demonstrated in Figure 2. The same figure shows that as long as σ DVAE ≥ 5 its value does not significantly affect the performance. In the sequel we use σ DVAE = 15.

3.5. Effectiveness of the encoder as a fast approximate minimizer. Proposition 2.1 shows that the proposed alternate minimization scheme in Algorithm 2.3 converges to a stationary point of J 1 . And so does the gradient descent scheme in [START_REF] Bora | Compressed sensing using generative models[END_REF]. Since both algorithms have to deal with non-convex energies, they both risk converging to spurious local minima. Also both algorithms solve essentially the same model when the variance γ of the coupling term tends to zero.

If our algorithm shows better performance (see next subsection), it is mainly because it relies on a previously trained VAE in two fundamental ways: (i) to avoid getting trapped in spurious local minima and (ii) to accelerate performance during the initial iterations (n < n min ). These two features are only possible if the autoencoder approximation is good enough and if the encoder is able to provide good initializations for the non-convex zoptimization subproblem in line 13 of Algorithm 2.3.

Figures 3 and4 illustrate these two properties of our VAE. We do so by selecting a random x 0 from MNIST test set and by computing z * (z 0 ) := gd z J 1 (x 0 , z) with different initial values z 0 . These experiments were performed using the ADAM minimization algorithm with learning rate equal to 0.01. Figure 4(a) shows that z * (z 0 ) reaches the global optimum for most (but not all) initializations z 0 . Indeed from 200 random initializations z 0 ∼ N (0, I), 195 reach the same global minimum, whereas 5 get stuck at a higher energy value. However these 5 initial values have energy values J 1 (x 0 , z 0 ) ≫ J 1 (x 0 , z 1 ) far larger than those of the encoder initialization z 1 = µ ϕ (x 0 ), and are thus never chosen by Algorithm 2.3. The encoder initialization z 1 on the other hand provides much faster convergence to the global optimum.

In addition, this experiment shows that we cannot assume z-convexity: The presence of plateaux in the trajectories of many random initializations as well as the fact that a few initializations do not lead to the global minimum indicates that J 1 may not be everywhere convex with respect to z. However, in contrast to classical works on alternate convex search, our approach adopts weaker assumptions and does not require convexity on z to prove convergence in Proposition 2.1.

In Figure 4(b) we display the distances of each trajectory to the global optimum z * (taken as the median over all initializations z 0 of the final iterates z * (z 0 )); note that this optimum is always reached, which suggests that z → J 1 (x 0 , z) has a unique global minimizer in this case. Finally, Figure 3 shows that the encoder approximation is quite good both in the latent space (Figure 3(a)) and in image space (Figures 3(b) and 3(c)). It also shows that the true posterior p θ (z|x) is pretty close to log-concave near the maximum of q ϕ (z|x).

3.6. Image restoration experiments. Choice of x 0 : In the previous section, our validation experiments used a random x 0 from the data set as initialization. When dealing with an image restoration problem, Algorithms 2.2 and 2.3 require an initial value of x 0 to be chosen. In all experiments we choose this initial value as A T y.

Choice of n 1 and n 2 : After a few runs of Algorithm 2.2 we find that in most cases, during the first 10 or 20 iterations z 1 decreases the energy with respect to the previous iteration, and this value depends on the inverse problems (for example, for denoising is smaller than for compressed sensing). But after at most 150 iterations the autoencoder approximation is no longer good enough and we need to perform gradient descent on z n in order to further decrease the energy. Based on these findings we set n 1 = 25 and n 2 = 150 in Algorithm 2.3 for all experiments. Note that we could also choose n 1 = n 2 = n max , since in all our experiments we observed that the algorithm auto-regulates itself, evolving from i * = 1 in the first few dozen iterations to i * = 3 when it is close to convergence. Choosing a finite value for n 1 and on a batch of 100 test images). For PSNR, higher is better and for LPIPS, lower is better. For comparison we provide the results of the baselines introduced in Section 3.1 (namely, Algorithm 1.1, CSGM [START_REF] Bora | Compressed sensing using generative models[END_REF], mCSGM (CSGM with restarts), PGD-GAN [START_REF] Shah | Solving linear inverse problems using gan priors: An algorithm with provable guarantees[END_REF] and PULSE [START_REF] Menon | PULSE: Self-Supervised Photo Upsampling via Latent Space Exploration of Generative Models[END_REF].) test images). For PSNR, higher is better and for LPIPS, lower is better. For comparison we provide the results of the baselines introduced in Section 3.1 (namely, Algorithm 1.1, CSGM [START_REF] Bora | Compressed sensing using generative models[END_REF], mCSGM (CSGM with restarts), PGD-GAN [START_REF] Shah | Solving linear inverse problems using gan priors: An algorithm with provable guarantees[END_REF] and PULSE [START_REF] Menon | PULSE: Self-Supervised Photo Upsampling via Latent Space Exploration of Generative Models[END_REF].) n 2 is only needed to ensure that i * = 3 when n → ∞, which is a necessary condition to prove statement 3 of Proposition 2.1.

Figure 5 shows the evolution of x k and G(z k ) from Algorithm 2.4 in an interpolation example. Here we can see how the exponential multiplier method in Equation (2.16) updates the values of

β k to ensure ∥G(z k ) -x k ∥ 2 ≤ ε.
Figures 6 and7 show the results of denoising, interpolation, compressed sensing, deblurring and super-resolution experiments on MNIST for different degradation levels using the proposed algorithm (JPMAP) and the baseline algorithms introduced in section 3.1. The metrics used are PSNR and LPIPS 4 [78] mean ± its standard error computed over 100 random experiments for each problem. Figure 9 displays the images of 10 representative interpolation and deblurring experiments from the hundreds of experiments summarized in figures 6 and 7.

These results show that JPMAP outperforms all other baseline algorithms in terms of PSNR and LPIPS when random restarts are not allowed. When 10 random restarts are allowed for CSGM, but not for JPMAP, then both algorithms (JPMAP and mCSGM) show a similar global performance: JPMAP tends to provide a slightly better result than mCSGM except for the most extremely ill-posed interpolation, super-resolution and compressed sensing experiments (when available measurements are less than 10% the number of pixels). In that case mCSGM outperforms JPMAP by an equally small margin. The latter case can be explained by the fact that the encoder (which is used by JPMAP but not by CSGM) struggles to generalize to images x which are very far away from M (the range of the generator). Indeed, in Section 3.4 we trained the VAE's encoder to generalize to M + n where n ∼ N (0, σ 2 DVAE I) and σ DVAE = 15/255. This value is optimal for moderately ill-posed problems, but more extreme problems may require larger values of σ DVAE or a coarse to fine scheme, where a coarse VAE (with large σ DVAE ) is used during the first few iterations and a finer VAE (with smaller σ DVAE ) is used later until convergence. Finally one may consider using random restarts for both JPMAP and CSGM for a more fair comparison.

Figure 8 performs a more detailed comparison between JPMAP and mCSGM, which also considers running times of both algorithms. For the stopping criteria used in our experiments, one run of JPMAP requires roughly as much time as mCSGM (with m = 10 restarts). In addition for moderately ill-posed problems (like interpolation of 75% missing pixels see subfigures (a) and (b)) where JPMAP's performance beats mCSGM, we can observe that JPMAP also converges much faster to that solution. For more extremely ill-posed problems (like interpolation of 90% missing pixels, see subfigures (c) and (d)) the opposite is true.

In the case of CelebA, we did not observe as much difference between JPMAP and CSGM as on MNIST. In Figure 10(a) the restorations on an interpolation problem (80% of missing pixels) are very similar to each other, but blurry. Also, although this problem is very illposed, both algorithms impressively find a solution z * very close to the code µ ϕ (x) of the ground truth image x, except for the third column where CSGM converges to a local minimum.

We hypothesize that, as CelebA is a substantially more complex dataset than MNIST, a simple model like vanilla VAE is over-regularizing the manifold of samples (underfitting problem). In particular, because of the spectral bias [START_REF] Nasim Rahaman | On the spectral bias of neural networks[END_REF] the learned manifold perhaps only contains low-frequency approximations of the true images as we can see in the reconstructions µ θ (µ ϕ (x)) of test samples (see Figure 10(b)). This may cause the posterior p Z|Y (z | y) to have fewer local minima. With more realistic generative models such as VDVAE [START_REF] Child | Very Deep VAEs Generalize Autoregressive Models and Can Outperform Them on Images[END_REF] or NVAE [START_REF] Vahdat | Nvae: A deep hierarchical variational autoencoder[END_REF], which better represent the true data manifold, we expect the objective function -log p Z|Y (z | y) to exhibit a much larger number of local minima, thus making it more difficult to optimize by a simple gradient descent scheme. In that situation the proposed JPMAP method would more clearly show its advantages.

4. Conclusions and Future work. In this work we presented a new framework to solve inverse problems with a convex data-fitting term and a non-convex regularizer learned in the latent space via variational autoencoders. Unlike similar approaches like CSGM [START_REF] Bora | Compressed sensing using generative models[END_REF], PULSE [START_REF] Menon | PULSE: Self-Supervised Photo Upsampling via Latent Space Exploration of Generative Models[END_REF] and PGD-GAN [START_REF] Shah | Solving linear inverse problems using gan priors: An algorithm with provable guarantees[END_REF] which learn the prior based on generative models, our approach is based on a generalization of alternate convex search to quasi-bi-convex functionals. This quasibi-convexity is the result of considering the joint posterior distribution of latent and image spaces. As a result, the proposed approach provides convergence guarantees that extend to a larger family of inverse problems. Experiments on denoising, interpolation, deconvolution, super-resolution and compressed sensing confirm this, since our approach gets stuck much less often in spurious local minima than CSGM, PGD-GAN or PULSE, which are simply based on gradient descent of a highly non-convex functional. This leads to restored images which are significantly better in terms of PSNR and LPIPS.

JPMAP vs related Plug & Play approaches. When compared to other decoupled plug & play approaches that solve inverse problems using NN-based priors, our approach is constrained in different ways:

(a) In a certain sense our approach is less constrained than existing decoupled approaches since we do not require to retrain the NN-based denoiser to enforce any particular property to ensure convergence: Ryu et al. [START_REF] Ryu | Plug-and-play methods provably converge with properly trained denoisers[END_REF] requires the denoiser's residual operator to be nonexpansive, and Gupta et al. [START_REF] Gupta | Cnn-based projected gradient descent for consistent ct image reconstruction[END_REF], Shah and Hegde [START_REF] Shah | Solving linear inverse problems using gan priors: An algorithm with provable guarantees[END_REF] require the denoiser to act as a projector. The effect of these modifications to the denoiser on the quality of the underlying image prior has never been studied in detail and chances are that such constraints degrade it. Our method only requires a variational autoencoder without any further constraints, and the quality and expressiveness of this prior can be easily checked by sampling and reconstruction experiments. Checking the quality of the prior is a much more difficult task for Ryu et al. [START_REF] Ryu | Plug-and-play methods provably converge with properly trained denoisers[END_REF], Gupta et al. [START_REF] Gupta | Cnn-based projected gradient descent for consistent ct image reconstruction[END_REF], Shah and Hegde [START_REF] Shah | Solving linear inverse problems using gan priors: An algorithm with provable guarantees[END_REF] which rely on an implicit prior, and do not provide a generative model.

(b) Unlike [START_REF] Ryu | Plug-and-play methods provably converge with properly trained denoisers[END_REF] which requires the data-fitting term F (x) to be strongly convex to ensure convergence, our method admits weakly convex and ill-posed data-fitting terms like missing pixels, compressed sensing and non-invertible blurring for instance.

(c) On the other hand our method is more constrained in the sense that it relies on a generative model of a fixed size. Even if the generator and encoder are both convolutional neural networks, training and testing the same model on images of different sizes is a priori not possible because the latent space has a fixed dimension and a fixed distribution. As a future work we plan to explore different ways to address this limitation. The most straightforward way is to use our model to learn a prior of image patches of a fixed size and stitch this model via aggregation schemes like in EPLL [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] to obtain a global prior model for images of any size. Alternatively we can use hierarchical generative models like in [START_REF] Karras | Progressive Growing of GANs for Improved Quality, Stability, and Variation[END_REF][START_REF] Vahdat | Nvae: A deep hierarchical variational autoencoder[END_REF] or resizable ones like in [START_REF] Bergmann | Learning Texture Manifolds with the Periodic Spatial GAN[END_REF][START_REF] Whang | Compressed Sensing with Invertible Generative Models and Dependent Noise[END_REF], and adapt our framework accordingly. map-x or map-z or joint map-x-z. In this work we explored and clarified the tight relationships between joint map-x-z estimation, splitting and continuation schemes and the more common map-z estimator in the context of inverse problems with a generative prior. On the other hand map-x estimators (which are otherwise standard in bayesian imaging) remained largely unexplored in the context of generative priors, due to the optimization challenges they impose, until the recent work of Helminger et al. [START_REF] Helminger | Blind Image Restoration with Flow Based Priors[END_REF], Whang et al. [START_REF] Whang | Compressed Sensing with Invertible Generative Models and Dependent Noise[END_REF] showed that a normalizing flow-based generative model allows to overcome those challenges and deems this problem tractable. Similarly Oberlin and Verm [START_REF] Oberlin | Regularization via deep generative models: an analysis point of view[END_REF] use Glow (an invertible normalizing flow) to compare synthesis-based and analysis-based reconstructions. Yet an extensive comparison of the advantages and weaknesses of these three families of estimators under the same prior model is still missing, and so is the link between these MAP estimators and the analysis/synthesis-based estimators in [START_REF] Oberlin | Regularization via deep generative models: an analysis point of view[END_REF]. This will be the subject of future work.

Extension to higher dimensional problems. The present paper provides a first proof of concept of our framework, on a very simple dataset (MNIST) with a very simple VAE. More experiments are needed to verify that the framework preserves its qualitative advantages on more high-dimensional datasets (like CelebA, FFHQ, etc.), and a larger selection of inverse problems.

Generalizing our proposed method to much higher dimensional problems implies training much more complex generative models which can match the finer details and higher complexity of such data. We can still use over-simplified generative models in those cases, but our preliminary experiments suggest that in that situation, not only do we obtain relatively poor reconstructions, but the objective function associated to the map-z problem presents less spurious local-minima: as a consequence our proposed joint map-x-z is overkill in that configuration, and does not present such a great competitive advantage.

The big challenge of generalizing our proposed method to much higher dimensional problems is then to train sufficiently detailed and complex generative models. And in this area VAEs traditionally lagged behind GANs in terms of quality of the generated samples, the former producing in general more blurred samples. Nevertheless some studies [START_REF] Mehdi | Assessing Generative Models via Precision and Recall[END_REF] show that VAEs and Normalizing Flows produce more accurate representations of the probability distribution. In the medium term our work should be able to benefit from recent advances in VAE architectures [START_REF] Child | Very Deep VAEs Generalize Autoregressive Models and Can Outperform Them on Images[END_REF][START_REF] Vahdat | Nvae: A deep hierarchical variational autoencoder[END_REF][START_REF] Dai | Diagnosing and Enhancing VAE Models[END_REF], and adversarial training for VAEs [START_REF] Pu | Adversarial symmetric variational autoencoder[END_REF][START_REF] Pu | Adversarial symmetric variational autoencoder[END_REF][START_REF] Zhang | Perceptual Generative Autoencoders[END_REF] that reach GAN-quality samples with the additional benefits of VAEs. These extensions are however non-trivial, since these VAEs have a huge number of parameters and they need to be retrained or fine-tuned using a denoising criterion (see section 3.4 and [START_REF] Daniel | Denoising criterion for variational auto-encoding framework[END_REF]) for our method to work properly. In addition, the latent space of the most competitive VAEs is much larger than the image space, which may reduce its regularization capabilities.

As an alternative, GAN-based generative models can be augmented with a denoising encoder network [START_REF] Donahue | Large Scale Adversarial Representation Learning[END_REF], and Normalizing Flows can also act as projectors or denoising VAEs if we split the latent space to separate the data manifold from its complement, as suggested in [START_REF] Brehmer | Flows for simultaneous manifold learning and density estimation[END_REF][START_REF] Liu | Disentangling Noise from Images: A Flow-Based Image Denoising Neural Network[END_REF]. In combination with relaxation techniques, such augmented GANs or specially tailored Flows may provide SOTA priors that fit our quasi-bi-convex optimization framework.

Towards stronger convergence guarantees under weaker conditions.. The proposed Algorithm 2.4 bears strong similarities with ADMM with non-linear constraints as introduced by Valkonen et al. [START_REF] Valkonen | A primal-dual hybrid gradient method for nonlinear operators with applications to MRI[END_REF][START_REF] Benning | Preconditioned ADMM with nonlinear operator constraint[END_REF] and analyzed by Latorre-Gómez et al. [START_REF] Latorre | Fast and Provable ADMM for Learning with Generative Priors[END_REF]. Latorre-Gómez result provides very strong convergence guarantees (linear convergence rates to a global optimum), but requires the data fitting term to be strongly convex or to satisfy a restricted strong convexity property. Our result, on the other hand, provides much weaker convergence guarantees (convergence to a stationary point), but does not require strong convexity. Further exploring these connections might hopefully lead to something closer to the best of both worlds.

In the last line we used Bayes rule to rewrite this posterior in more simple terms. However, this expression still involves the unknown conditional p Y |Z (y | z).

Let us express this maximization in terms of p Y |X (y | x).

To do so we recall the relation between the conditionals and the joint:

(B.2) p Y |Z (y | z) p Z (z) = p Y,Z (y, z) = p Z|Y (z | y) p Y (y)
We can also compute the joint distribution p Y,Z (y, z) by marginalization on a third random variable X:

p Y,Z (y, z) = p X,Y,Z (x, y, z) dx = p Y |X,Z (y | x, z) p X|Z (x | z) p Z (z) dx = p Y |X (y | x) δ(x -G(z))p Z (z) dx = p Y |X (y | G(z)) p Z (z) (B.3)
The third line follows from our graphical model Z → X → Y which implies that once we know X = x, then Z provides no additional information, therefore

p Y |X,Z (y | x, z) = p Y |X (y | x) .
The last line follows simply from the integration on x of a delta function.

p Z|Y (z | y) = 1 p Y (y) p Y |X (y | G(z)) p Z (z)
This proves the main result of this section: B.2. MAP-x. The map-x estimator is obtained by maximizing the posterior with respect to x. The generative model induces a prior on X via the push-forward measure p X = G♯p Z , which following [48, section 5] can be developed as

p X (x) = p Z G -1 (x) det S(G -1 (x)) δ M (x)
where S = ∂G ∂z T ∂G ∂z is the squared Jacobian and the manifold M = {x : ∃z, x = G(z)} represents the image of the generator G.

With such a prior p X , the x-optimization (1.2) required to obtain xmap becomes intractable (in general), for various reasons:

• the computation of det S,

• the inversion of G, and • the hard constraint x ∈ M. These operations are are all memory and/or computationally intensive, except when they are partially addressed by the use of a normalizing flow like in [START_REF] Helminger | Blind Image Restoration with Flow Based Priors[END_REF][START_REF] Whang | Compressed Sensing with Invertible Generative Models and Dependent Noise[END_REF].

Appendix C. Joint MAP-x-z, Continuation Scheme and convergence to MAP-z. The functional J 1,β introduced in Equation 1.6 can be seen from two different perspectives. From a machine learning perspective it corresponds to the joint log-posterior J 1 in the case where Σ θ (z) = 1 β I and µ θ = G, namely:

J 1,β (x, z) = 1 2 σ 2 ∥A x -y∥ 2 F (x, y) + β 2 ∥x -µ θ (z))∥ 2
H θ (x, z) = φ β (x, z)

+ 1 2 ∥z∥ 2 + C β
From an optimization standpoint it can be considered as an inexact penalisation procedure: We want to solve the constrained problem min To do so, we introduced the family of unconstrained problems min

x,z J 1,β (x, z)

and their corresponding minimizers ( xβ , ẑβ ) ∈ arg min

x,z J 1,β (x, z)

which for β = 1 γ 2 provide the map-x-z estimator.

We can show that the map-x-z estimator converges to the map-z estimator when β → ∞ (or equivalently γ → 0).

Proposition C.1. The unconstrained functional tends to the constrained functional plus the constraint:

(C.2) J 1,β (x, z) β→∞ ---→ J 1,∞ (x, z) = F (x, z) + ι x=µ θ (z) (x, z) + 1 2 ∥z∥ 2
and the unconstrained minimizers tend to the constrained minimizer as β → ∞: Proof. The pointwise convergence of φ β to ι x=G(z) as β goes to ∞ is straightforward. Let us first prove that for any sequence (β n ) n that goes to ∞, the quantity ∥ xβn -G( ẑβn )∥ goes to zero. Otherwise, for any ε > 0, there exists a subsequence (β n j ) j such that ∥ xβn -G( ẑβn )∥ > ε. In this case, for any z, one has by optimality J 1,0 (G(z), z) = J 1,βn j (G(z), z) ≥ J 1,βn j ( xβn j , ẑβn j ) > J 1,0 ( xβn j , ẑβn j ) +

(C.
β n j 2 ε 2
As a result, the nonnegative quantity J 1,0 ( xβn j , ẑβn j ) goes to -∞, which leads to a contradiction. Thus, one has x∞ = G( ẑ∞ ) for any limit point ( x∞ , ẑ∞ ) of ( xβn , ẑβn ). Assume that J 1,0 ( x∞ , ẑ∞ ) > J 1,0 (x * , z * ). Since In general, we can only prove that the limit points of the sequences generated by Algorithm 1.1 are stationary points of C.4. However, if the growth of β is sufficiently slow, then we obtain the optimality of the limit points. Indeed, given that, in Algorithm 1.1, each subloop is an exact BCD scheme, one has for any n and any j ∀ z, J 1,β k j (x n-1 , z n ) ≤ J 1,β k j (x n-1 , z)

By considering the subsequence (x n ℓ , z n ℓ ), which converges to x k j ∞ , z k j ∞ (we recall that x n ℓ and x n ℓ -1 have same limit), we can prove that ∀ z, J 1,β k j (x

k j ∞ , z k j ∞ ) ≤ J 1,β k j (x k j ∞ , z) that is, z k j
∞ is a minimizer of J 1,β k j (x k j ∞ , •). Hence, we have

∀ z, J 1,0 (x k j ∞ , z k j ∞ ) + β k j 2 ∥x k j ∞ -G(z k j ∞ )∥ 2 ≤ J 1,0 (x k j ∞ , z) + β k j 2 ∥x k j ∞ -G(z)∥ 2
Assume that β k j ∥x

k j ∞ -x∞ ∥ 2 -→ j→+∞ 0
. By letting j to ∞, we get that, for any z such that

G(z) = x k j ∞ , J 1,0 ( x∞ , ẑ∞ ) ≤ J 1,0 ( x∞ , z) + lim j→+∞ β k j 2 ∥x k j ∞ -x k j ∞ ∥ 2 = J 1,0 ( x∞ , z)
that is, ẑ∞ is a minimizer of J 1,0 ( x∞ , •) + χ C ( x∞ , •). By definition of f , this also means that ẑ∞ is a minimizer of f . However, one has to note that the growth control for β depends on the convergence speed of x k ∞ , which cannot be estimated. Algorithm 2.4 is a particular (truncated) case of Algorithm 1.1 with an adaptive choice of β that does not need to go to ∞. * )) According to the calculus above, this proves that (x * , z * ) is a stationary point of J 1,2λ * . Note that, if λ * = 0, then (x * , z * ) is a minimizer of J 1,0 . Otherwise, one has ∥x * -G(z * )∥ 2 = ε.

Hence, if we consider Algorithm 1.1 with the update rule for β k as in Algorithm 2.4 and a stopping rule saying that the iterations stop as soon as, for any given k, ∥x k ∞ -G(z k ∞ )∥ 2 ≤ ε there are two possible cases:

1. case λ * = 0: then (x k ∞ , z k ∞ ) is a solution of the constraint problem iff ∇J 1,0 (x k ∞ , z k ∞ ) = 0 (that is, x k ∞ = G(z k ∞ )); 2. case λ * > 0: unless ∥x k ∞ -G(z k ∞ )∥ 2 exactly equals ε, (x k ∞ , z k ∞ )
is not a solution of the constraint problem However, in general, (x k ∞ , z k ∞ ) is a solution of the following constraint problem min

∥x k ∞ -G(z k ∞ )∥ 2 ≤ε
J 1,0 (x, z)

with ε = ∥x k ∞ -G(z k ∞ )∥ 2 ≤ ε.
Hence, if we stop the iterations when ∥x k ∞ -G(z k ∞ )∥ 2 ≤ ε, we will get a solution of min

∥x k ∞ -G(z k ∞ )∥ 2 ≤ε
J 1,0 (x, z), ε ≤ ε which provides an error control as well.

p

  X|Y (x | y) = arg min x {F (x, y) + λU (x)}

Algorithm 2 . 1

 21 Joint posterior maximization -exact case Require: Measurements y, Autoencoder parameters θ, ϕ, Initial condition x 0 Ensure: x, ẑ = arg max x,z p X,Z|Y (x, z | y)1: for n := 0 to maxiter do 2:

  k + 1 11: end while 12: return x k , z k 3. Experimental results.

Figure 1 .

 1 Figure 1. Evaluating the quality of the generative model as a function of σDVAE. On (a) Denoising (Gaussian noise σ = 150), (b) Compressed Sensing (∼ 10.2% measurements, noise σ = 10) and (c) Interpolation (80% of missing pixels, noise σ = 10). Results of both algorithms are computed on a batch of 50 images and initialising on ground truth x * (for CSGM we use z0 = µ ϕ (x * )).

Figure 2 .

 2 Figure2. Evaluating the effectiveness of JPMAP vs CGSM as a function of σDVAE (same setup of Figure1). Without a denoising criterion σDVAE = 0 the JPMAP algorithm may provide wrong guesses z 1 when applying the encoder in step 2 of Algorithm 2.2. For σDVAE > 0 however, the alternating minimization algorithm can benefit from the robust initialization heuristics provided by the encoder, and it consistently converges to a better local optimum than the simple gradient descent in CSGM.

Figure 3 . 1 2 ∥z∥ 2 .

 312 Figure 3. Encoder approximation: (a) Contour plots of -log p θ (x|z) + 12 ∥z∥ 2 and -log q ϕ (z|x) for a fixed x and for a random 2D subspace in the z domain (the plot shows ±2Σ 1/2 ϕ around µ ϕ ). Observe the relatively small gap between the true posterior p θ (z|x) and its variational approximation q ϕ (z|x). This figure shows some evidence of partial z-convexity of J1 around the minimum of J2, but it does not show how far is z 1 from z 2 . (b) Decoded exact optimum x1 = µ θ arg maxz p θ (x|z)e 1 2 ∥z∥ 2 . (c) Decoded approximate optimum

  (a) Energy evolution, initializing with N (0, I).(b) Distance to the optimum at each iteration of (a).

Figure 4 .

 4 Figure 4. Effectiveness of the encoder approximation: We take x0 from the test set of MNIST and minimize J1(x0, z) with respect to z using gradient descent from random Gaussian initializations z0. The blue thick curve represents the trajectory if we initialize at the encoder approximation z 1 = arg min z J2(x0, z) = µ ϕ (x0). (a): Plots of the energy iterates J1(x0, z k ). (b): ℓ 2 distances of each trajectory with respect to the global optimum z * . Conclusion: Observe that the encoder initialization allows much faster convergence both in energy and in z, and it avoids the few random initializations that lead to a wrong stationary point different from the unique global minimizer.

Figure 5 .

 5 Figure 5. Evolution of Algorithm 2.4. In this interpolation example, JPMAP starts with the initialization in (a). During first iterations (b) -(d) where β k is small, x k and G(z k ) start loosely approaching each other at a coarse scale, and x k only fills missing pixels with the ones of G(z k ) (in particular the noise of y is still present). By increasing β k in (e) -(f ) we enforce ∥G(z k ) -x k ∥ 2 ≤ ε. Here we set ϵ = 3 255 2 d, that is, MSE

Figure 6 .

 6 Figure 6. Denoising, Compressed Sensing and Interpolation: Evaluating the effectiveness of Algorithm 2.3 (fixed β) and Algorithm 2.4 for different values of ϵ = α 255 2 n, with σDVAE = 15 (metrics were computed

Figure 7 .

 7 Figure 7. Deblurring and Super-resolution: Evaluating the effectiveness of Algorithm 2.3 (fixed β) and Algorithm 2.4 for different values of ϵ = α 255 2 n, with σDVAE = 15 (metrics were computed on a batch of 100

Figure 8 .

 8 Figure 8. Time/PSNR comparison between mCSGM and JPMAP. Left: Confidence intervals (for a batch of 100 random experiments) for PSNR vs computing time for both algorithms on the interpolation problem with p% of missing pixels with noise std σ = 10/255. Right: Detailed view of one of the 100 random experiments on the left. Blue lines represent m = 10 random restarts of CSGM and the orange line is the PSNR of the best z k at iteration k of mCSGM as measured by (1.4).

  (a) Results on interpolation (b) Results on deblurring.

Figure 9 .

 9 Figure 9. Experimental results on MNIST. Comparison of JPMAP with the baseline algorithms described in Section 3.1. (a) Some selected results from the interpolation experiment with 80% of missing pixels and Gaussian noise with σ = 10/255. From top to bottom: original image x * , corrupted image y, and the results computed by CSGM, mCSGM, PGD-GAN, PULSE, Algorithm 1.1 and JPMAP. (b) Same as (a) from the deblurring experiment with kernel size 3 × 3 and Gaussian noise with σ = 10/255. Conclusion: Our algorithm performs generally better than the baseline algorithms, although in some cases it falls behind mCSGM.

Figure 10 .

 10 Figure 10. (a) Some preliminary results on CelebA: 80% of missing pixels, noise std σ = 10/255. From top to bottom: original image x * , corrupted image x, restored by CSGM [6], restored image x by our framework. (b) Reconstructions µ θ (µ ϕ (x)) (even columns) for some test samples x (odd columns), showing the overregularization of data manifold imposed by the trained vanilla VAE. As a consequence, -log p Z|Y (z | y) does not have as many local minima and then a simple gradient descent yields almost the same result as JPMAP (except on third column of (a)).

  From equations (B.2) and (B.3) we can derive an expression of p Z|Y (z | y) in terms of p Y |X (• | •) and the generator G namely:

Proposition B. 1 (

 1 map-z estimator for deterministic generative models). Assume we have • a deterministic generative model where X = G(Z) and • an inverse problem characterised by the log conditional distribution log p Y |X (y | x) = -F (x, y). Then the map-z estimator is computed as xmap-z = G( ẑmap-z ) where ẑmap-z = arg max z p Y |X (y | G(z)) p Z (z) = arg min z {F (G(z), y) -log p Z (z)} . (B.4)

2 =J 1

 21 ,0 (x,z) with C = {(x, z) | x = µ θ (z)} whose solution provides the map-z estimator (C.1) (x * , z * ) ∈ arg min (x,z)∈C J 1,0 (x, z).

  3) limβ→∞ ( xβ , ẑβ ) ∈ arg min (x,z)∈C J 1,0 (x, z) = arg min x,z J 1,∞ (x, z).

J 1 ,2 ∥z∥ 2 Proof.

 12 βn ( xβn , ẑβn ) ≤ J 1,βn (x * , z * ) = J 1,0 (x * , z * ) < J 1,0( x∞ , ẑ∞ ) this leads to another contradiction.The previous result motivates Algorithm 1.1. Consider Algorithm 1.1 in the ideal case (maxiter=∞) where the internal loop converges.Proposition C.2 (Convergence of Algorithm 1.1). Let (x k ∞ , z k ∞ ) k be a sequence generated by Algorithm 1.1 when maxiter=∞. If (z k ∞ ) k is bounded, then any limit point of (x k ∞ , z k ∞ ) k is in C. Moreover, any limit point of (z k ∞ ) k is a stationary point of (C.4) f (z) = F (G(z), y) + 1 Note that, for any k, (x k ∞ , z k ∞ )is a limit point of the sequence generated by the k-th subloop in Algorithm 1.1 if it does not converge. Let (β k ) k a sequence that converges to ∞. Let k ∈ N. We consider the sequence (x k n , z k n ) n generated by ∀ n ∈ N, z k n+1 ∈ arg min z J 1,β k (x k n , z) and x k n+1 = arg min x J 1,β k (x, z k n+1 )

Proposition C. 3 (

 3 Convergence of Algorithm 2.4).Proof. Let us write the Lagrangian of the problem solved in Algorithm 2.4:∀ λ ≥ 0, L(x, z; λ) = J 1,0 (x, z) + λ (∥x -G(z)∥ 2 -ε)KKT conditions ensure that any solution (x * , z * ) of the constrained problem is associated to at least one Lagrange multiplier λ * ≥ 0 such that∂L ∂(x, z) (x * , z * ; λ * ) = 0 = ∂J 1,0 ∂x (x * , z * ) + 2λ * (x * -G(z * )) ∂J 1,0 ∂z (x * , z * ) + 2λ * (DG(z * )) * (x * -G(z

  1.1 map-z splitting Require: Measurements y, Initial condition x 0 , maxiter, k max , {β 0 , . . . , β kmax } Ensure: x = G arg max z p Z|Y (z | y) 1: for k := 0 to k max do

	2:

  good enough

	8:	done := TRUE
	9:	end if
	10:	end if
	11:

  is good enough

	17:	done := TRUE
	18:	end if
	19:	end if
	20:

In[START_REF] Laumont | On Maximum-a-Posteriori estimation with Plug & Play priors and stochastic gradient descent[END_REF] the PnP-ADMM and PnP-FBS algorithms introduced in[START_REF] Ryu | Plug-and-play methods provably converge with properly trained denoisers[END_REF][START_REF] Xu | Provable Convergence of Plug-and-Play Priors with MMSE denoisers[END_REF] are reported to converge in practice quite far beyond the conditions of the theorem, but require (to obtain optimal performance) the regularization parameter λ to be tuned to values that are far outside the region where convergence is guaranteed. Also, the performance is significantly degraded if λ is constrained to the range where convergence is guaranteed.

In another context a primal-dual optimization algorithm was proposed to solve a similar optimization problem[START_REF] Benning | Preconditioned ADMM with nonlinear operator constraint[END_REF], but this approach was not explored in the context where G is a generative model.

Code available at https://github.com/mago876/JPMAP.

MNIST images were zero-padded to 32 × 32 because LPIPS does not accept 28 × 28 images.

Acknowledgments. We would like to sincerely thank Mauricio Delbracio, José Lezama and Pablo Musé for their help, their insightful comments, and their continuous support throughout this project.

continuous support throughout this project. Funding: This work was funded by ECOS Sud Project U17E04, by the French-Uruguayan Institute of Mathematics and Interactions (IFUMI), by CSIC I+D 2018-256 (Uruguay) and by ANII (Uruguay) under Grant 11 FCE 1 2017 1 135458 and by the French Research Agency through the PostProdLEAP project (ANR-19-CE23-0027-01). Computer experiments for this work ran on a Titan Xp GPU donated by NVIDIA, as well as on HPC resources from GENCI-IDRIS (Grants 2020-AD011011641 and 2021-AD011011641R1).

Appendix A. Properties of J 1 . In this section, we establish that the objective function J 1 fulfills the assumptions required to prove the convergence of Algorithm 3, namely

• J 1 (•, z) is convex for any z;

• J 1 (•, z) has a unique minimizer for any z;

• J 1 is coercive;

• J 1 is continuously differentiable; We recall that

where

Thus, it is the sum of three non-negative terms.

A.1. Convexity and unicity of the minimizer of J 1 (•, z). Let z be fixed. Then there exists a constant C ∈ R such that ∀ x

Being the sum of two quadratic forms, J 1 (•, z) is obviously twice differentiable. Its gradient is given by

and its Hessian is

Since Σ θ (z) = γ 2 I the Hessian is positive definite (without the need to assume that A is full rank), and we have that Lemma A.1. J 1 (•, z) is strictly convex for any z.

An immediate consequence is the unicity of the minimizer of the partial function J 1 (•, z).

A.2. Coercivity of J 1 .

Lemma A.2. J 1 is coercive.

Proof. First, let us note that J 1 is the sum of three non-negative terms. If it was not coercive, then we could find a sequence (x k , z k ) → ∞ such that J 1 (x k , z k ) is bounded. As a consequence all three terms are bounded. In particular the last term ∥z k ∥ is bounded, which means that x k → ∞. From Property 1, {µ θ (z k )} and {Σ θ (z k )} are bounded for bounded {z k }. Now, from the definition of the second term of J 1 , we get that, {µ θ (z k )} and {Σ θ (z k )} being bounded and x k going to ∞ yield that H θ (x k , z k ) goes to infinity, while being bounded. This leads to a contradiction and thus proves that J 1 is coercive.

A.3. Regularity of J 1 . In the sequel we adopt the common assumption that all neural networks used in this work are composed of a finite number d of layers, each layer being composed of: (a) a linear operator (e.g. convolutional or fully connected layer), followed by (b) a non-linear L-Lipschitz component-wise activation function with 0 < L < ∞.

Therefore we have the following property:

Property 1. For any neural network f θ with parameters θ having the structure described above: There exists a constant C θ such that ∀u,

Concerning activation functions we use two kinds:

• continously differentiable activations like ELU, or • continuous but non-differentiable activations like ReLU Hence, by composition, we have that Lemma A.3. For continuously differentiable activation functions, J 1 is continuously differentiable.

Appendix B. MAP-x and MAP-z for deterministic generative models. Assume that the stochastic γ-generative model is

We now analyze the map-z and map-x estimators for the limit case when γ = 0. This is what we call a deterministic generative model, and it includes GANs for instance.

B.1. MAP-z.

By definition the map-z estimator is obtained by maximising the posterior with respect to z:

Since J 1,β k corresponds to a particular instance of J 1 , and since Algorithm 1.1 can be seen asymptotically as a particular instance of Algorithm 2.3, one can use all the results established in Proposition 2.1. In particular, the sequence (x k n , z k n ) n admits a limit point (x k ∞ , z k ∞ ) and we have

k is also bounded. Then it admits a limit point denoted ( x∞ , ẑ∞ ). Let (x

) j be a convergent subsequence of limit ( x∞ , ẑ∞ ). Let us assume that x∞ ̸ = G( ẑ∞ ). Then, there exists a > 0 and j 0 ∈ N such that

Hence, one has

which leads to a contradiction. This proves that x∞ = G( ẑ∞ ). Otherwise said, (x k ∞ -G(z k ∞ )) k goes to zero.

Since we have for any k

Note that f (z) = J 1,0 (G(z), z). One can check that f is differentiable and that ∇f (z) = (DG(z)) * ∂J 1,0 ∂x (G(z), z) + ∂J 1,0 ∂z (G(z), z)

Hence, we have proved that ∇f ( ẑ∞ ) = 0 Conclusion: If (z k ∞ ) k is bounded, any limit point of (z k ∞ ) k is a stationary point of (C.4).