A CONVEX APPROXIMATION FRAMEWORK FOR NEURAL LIKELIHOOD-BASED BAYESIAN INVERSE PROBLEMS
arXiv:2607.06252v1 [stat.ML] 7 Jul 2026
FABIAN SCHNEIDER1,2 , TAPIO HELIN2 , AND LEILA TAGHIZADEH1 .
Abstract. Many problems in science and engineering are difficult to model accurately, either due to unknown physical mechanisms, poorly quantified measurement uncertainty, or prohibitive computational costs of high-fidelity simulations. These challenges limit the applicability of classical probabilistic inference methods such as Markov chain Monte Carlo, especially in highdimensional Bayesian inverse problems. As data from scientific experiments become increasingly available, machine learning methods offer a flexible alternative to explicit parametric modelling. We study neural likelihood approximation, where the goal is to learn the likelihood function directly from data without explicit knowledge of the underlying data-generating process. A common approach trains likelihood surrogates by minimizing the Kullback–Leibler divergence between the true posterior and an approximate posterior, which is equivalent to minimizing the expected negative log-likelihood. This work improves the theoretical foundations of neural likelihood approximation by alleviating limitations of restrictive model classes: we show that, by working with un-normalized potentials and folding normalization into the training objective, the resulting learning problem is strictly convex. We show that empirical minimizers of the resulting data-driven objective converge to the true likelihood as the sample size grows. Numerical experiments for the neural likelihood approximation are conducted for a deblurring and a non-linear PDE based imaging problem.
1. Introduction In many areas of science and engineering, an accurate mathematical model of the problem of interest is difficult to obtain, or there is no canonical way to quantify the uncertainty in the measurements. Even when a suitable mathematical model exists, evaluating its discretization can be computationally expensive (for example when it requires solving one or more partial differential equations (PDEs)) which may render problems of probabilistic inference such as Markov-Chain Monte Carlo (MCMC) infeasible, especially for very high dimensional problems such as Bayesian inverse problems on function spaces. Examples of problems that fit this description include inverse problems in medical imaging [2, 15] or geophysical imaging [49], modelling in evolutionary biology [37, 50, 27], ground water models [53], transmission dynamics of bacterial infections [30, 9], financial models [3, 23] and neurophysiological problems [22]. With increasing availability of measurement data from scientific experiments, machine learning methods are becoming more and more popular over hand-tuned explicit models involving a finite amount of parameters in an approximative class. We focus on neural likelihood approximation, where the goal is to infer the likelihood function from empirical data without a mathematical model of the precise physics of the data generating process. Tuning the likelihood approximation through the KL-divergence between the posterior Date: July 8, 2026. 2020 Mathematics Subject Classification. 62G05, 65C05, 65N21, 65N75, 60B12. Key words and phrases. Neural likelihood approximation, infinite dimensional Bayesian inverse problem, datadriven methods, semiconductor devices. L. Taghizadeh acknowledges support from the Austrian Science Fund (FWF), Elise-Richter grant DOI 10.55776/V1000. Funding was also received from the European Research Council (ERC) under the Horizon 2020 research and innovation program of the European Union (Grant agreement No. 101125225). This work was supported by the Research Council of Finland (353094, 348504 and 359183). Support by the Austrian Science Fund (FWF) grant I6667-N is acknowledged. 1 Institute of Analysis and Scientific Computing, TU Wien, Vienna, Austria [email protected], [email protected] 2 School of Engineering Sciences, LUT University, Lappeenranta, Finland [email protected]. 1
2
SCHNEIDER, HELIN, TAGHIZADEH
distribution µy and an approximation µyf corresponding to the approximate likelihood has seen a lot of promising experiments in the literature [22, 34, 17, 38]. As is established (see e.g. [34, 35]), this corresponds to minimising the expected negative log-likelihood (Lemma 2.7 for details). In spite of notable achievements, the approach continues to lack a comprehensive theoretical explanation. In practice, normalization of the likelihood as a probability measure is often achieved through certain classes of parametric distributions [22, 17] or normalizing flows [34, 38] or invertible networks [35]. This is restrictive: in each case the resulting function class may be too small to contain a near-optimal likelihood. 1.1. State of the art. In the context of inverse problems, the restrictive computational cost has been addressed widely through different classes of surrogate models. One of the variants corresponds to adjusting the observational noise to compensate for modelling error [1, 16, 11] given an approximation of the forward map. Polynomial Chaos expansion [25, 26, 51] is based on the idea of finding a computationally feasible approximation of the forward map. More recently, deep learning based surrogates for the forward map have been studied extensively, see e.g. [13, 52, 20, 7] among many others. Random likelihood approximation through Gaussian processes (GP) has been studied widely in the literature including statistical properties such as consistency [44, 12]. It is applicable to either an approximation of the forward map or the likelihood directly. Notice however that often a finite amount of evaluations of the likelihood is required for the estimation. Methods in the field of simulation based inference are applicable without a given mathematical model of the forward process or likelihood. The most classical method is given by approximate Bayesian computations (ABC) [39, 24, 41], where samples from a prior distribution are passed through a simulator until a quantity similar enough to the outcome of the physical experiment is observed. To overcome limitations especially in terms of applicability when no mathematical forward model is available, approximation of probability distribution through its likelihood or a density from data is studied under the terms of neural likelihood approximation [34, 48, 22, 14] or conditional density estimation [32, 31, 38], which require only empirical data from the joint distribution of measurements and corresponding quantities of interest. Similar variants exist for neural posterior estimation [6, 14]. The statistical accuracy of neural likelihood approximation under suitable assumptions are studied in [8]. From a theoretical point of view, stability of the posterior with respect to the likelihood is studied extensively in the literature and Lipschitz continuity of the posterior with respect to a likelihood perturbation is established in KL-divergence, Hellinger and other distances [43, 42, 5]. For random surrogates, the Hellinger distance between posterior and its approximation is bounded by a moment of likelihood mismatch in [21]. 1.2. The problem setting. Let the prior measure µ on a Hilbert space H be given. Denote by L : H × Rm → R the negative log-likelihood (NLL) or potential, describing the distribution of the dependent variable (measurement) y ∈ Rm conditional on the unknown independent variable (quantity of interest) x ∈ H. By Bayes’ theorem, the posterior distribution µy of x conditional on y has a density with respect to µ given by 1 dµy (·) := exp (−L(·; y)) , (1) dµ Z(y) µ-almost everywhere with Z Z(y) =
exp (−L(x; y)) µ(dx). H
We denote by λ(dx, dy) the joint distribution of (x, y) and by π(y) the marginal distribution of y. The Kullback-Leibler (KL) divergence for measures µ ≪ ν on H is given by Z dµ DKL (µ||ν) := log (x)µ(dx). dν H
CONVEX NEURAL LIKELIHOOD APPROXIMATION
3
Notice that the KL-divergence is not symmetric, it may take different values depending on the arrangement of its arguments. The usual setup of (sequential) neural likelihood estimation can be summarized as follows (see [34] for a detailed discussion): In a first step, training data (xi , yi )N i=1 from the joint distribution λ is generated by sampling from the prior and using a simulator to compute the observation yi conditional on xi for i = 1, · · · , N . We emphasize that, since in general the noise level of the simulator is unknown, the simulator cannot serve in a likelihood function. In the next step a neural likelihood estimate qθ is fitted by estimating θ∗
(2)
:=
argmax E(x,y) log qθ (y|x) θ∈Θ
where the generated data is used to estimate the expectation. By an argument similar to Lemma 2.7 in our text, argmax E(x,y) log qθ (y|x) = argmin Ey DKL (q(y|x)||qθ (y|x)) , θ∈Θ
θ∈Θ
where q denotes the true likelihood. The class Θ can be of parameters Θ ⊂ Rd or of (normalized) functions such as autoregressive flows [33]. The approximative likelihood qθ can now be used to approximate the posterior for a given measurement y. In the sequential variant, data can be regenerated by sampling the xi from the approximative posterior conditioned on a given measurement y instead. Multiple iterations of likelihood estimation and data generation can be performed. This work is motivated by Bayesian inverse problems. Important challenges include high computational cost of the evaluation of the forward map in large scale inverse problems on function spaces, where an evaluation can correspond to solving a (system of) PDE(s). A computationally cheaper approximation of the NLL L is often necessary to determine unknown functions and allow uncertainty quantification by generating many posterior samples. Another feasible application is where no accurate mathematical formula for the forward map is available but instead limited data can be generated from a physical system. 1.3. Our contribution. As noted above, a key difficulty in neural likelihood approximation is choosing a class of potentials that is both expressive enough and normalized as a probability distribution. We address this as follows. Given a general function, we normalize it formally by dividing by its expectation, so that the normalization is absorbed into the expected KL-based loss. We show that this loss is convex in the un-normalized NLL, and we establish consistency of the approach: minimizers of the empirical loss converge to the true potential as the amount of data grows. More precisely, we consider the optimization problem argmin Ey DKL µy ||µyf , (3) f ∈F
where µyf is the posterior corresponding to the approximate NLL: (4)
dµyf dµ
(x) := R
exp (−f (x; y)) in L1 (µ) exp(−f (x; y))µ(dx) H
and f ∈ F in a suitable class of functions, see Definition 2.3 for details. Our main findings can be summarized as follows: • We remove the need to parameterize normalized likelihoods by optimizing (3) directly over a general class of un-normalized potentials. The key observation is the identity (Lemma 2.7), which reads Eπ DKL µy ||µyf = C + Eλ [f (x; y)] + Eπ [log Zf (y)], with C independent of f . Normalization is absorbed into the log-normalizer Eπ [log Zf (y)], so the objective can be evaluated from samples of the joint distribution λ without access to the true potential L. Minimizing the expected KL-divergence is thus equivalent to minimizing the expected negative log-likelihood together with this log-normalizer.
4
SCHNEIDER, HELIN, TAGHIZADEH
• In Theorem 2.8, we show that the KL-divergence Ey DKL µy ||µyf is convex as a function of f . Exploiting this convexity of course requires the class F to be convex as well. We stress that this fails for the usual normalizing-flow parametrizations [34, 38]: the set of functions f for which e−f integrates to one is in general not convex. • In Section 3, we derive an approximation of the objective (3) that relies solely on data from the joint distribution λ. This allows accurate and computationally efficient representations of the NLL to be learned directly from such data. • In Theorem 3.7 we show that, under suitable assumptions, the minimizer of the data-based objective converges to the minimizer of the population objective as the number of data points increases. • In Section 4, we present numerical experiments on two Bayesian inverse problems: a toy de-blurring problem, in which the unknown is a one-dimensional Gaussian process, and a nonlinear, PDE-based problem of estimating the doping profile of a semiconductor from voltage-current measurements. Our work is motivated by [11], where the authors show that the KL-divergence in (3) is bi-convex in the mean and covariance parameters of a Gaussian likelihood. In summary, we propose a framework that enables offline learning purely based on empirical data, and can further be used to speed up posterior sampling once measurement data is acquired. The framework is backed by theoretical convergence guarantees in the data size. It has the advantages that the optimization problem is convex and without approximations in the objective function. In many scientific problems, the unknown quantity is represented by a function or another infinite-dimensional parameter to a PDE. In particular in the inverse problems community, significant effort was employed to ensuring methods are independent of the choice of discretization [19, 18]. To this end, we formulate the optimization problem on (possibly infinite-dimensional) Hilbert spaces and show the fundamental properties in Lemma 2.7 and Theorem 2.8 in the un-discretized setting. We state the consistency property in a finite-dimensional setting since the underlying theory on bracketing is developed mostly in the setting of function classes on euclidean spaces. The rest of the paper is organized as follows: In Section 2, we state the main assumptions and use them to show how the objective function can be implemented. Further we show convexity of the objective function. In Section 3, we propose an approximate objective function based on empirical data. We show consistency of the empirical approximation. In Section 4, we implement our approach and illustrate our numerical results for two Bayesian inverse problems, namely a deblurring problem and a doping inverse problem in semiconductor devices. Finally, conclusions are drawn in Section 5.
2. The Negative log-likelihood approximation and main properties In this section, first we show that inverse problems with a Gaussian likelihood distribution satisfy certain assumptions. We define a KL-divergence based objective function for an approximative likelihood that is based on an integral in the joint distribution λ. Then, we show that this objective function is in fact convex. Finally, we define equivalence class of (likelihood) functions with corresponding norm, where two functions are equivalent if they equal up to normalization in y. It is in this norm that we can expect optimal likelihood functions to be unique. Consider a Bayesian inverse problem with a Gaussian likelihood function. More precisely, assume to have a forward map A : H → Rm and the observational noise covariance Γ ∈ Rm×m given. The prior distribution of x ∈ H before any measurement is made, is given by µ. A measurement y ∈ Rm is obtained through (5)
y = A(x) + ε,
ε ∼ N (0, Γ).
CONVEX NEURAL LIKELIHOOD APPROXIMATION
5
The NLL corresponding to (5) is given by 1 2 ∥A(x) − y∥Γ . 2 When the true NLL is not known, we are interested in estimating one of the quantities • the NLL in a class of functions F of (x, y). We denote a likelihood emulator in this class a free-form approximation. • the true forward map A from a class of forward maps A. We call an approximation in this class residual approximation. In this class, we are approximating only the forward map and assume the observational noise level to be known. • the true forward map A and the observational noise Γ, from a class of forward maps A and observational noises G, respectively. This approximation we call calibrated residual approximation. In the former case, the assumptions on the class F for this work are expressed in Assumption 2. In the latter cases, we will elaborate the necessary assumptions on the classes A and G in Examples 2.1 and 2.2 to satisfy Assumption 2. L(x, y) =
2.1. Assumptions. Let the following assumptions, which are an extension of the assumptions in [10] to a general potential, hold throughout this section. Assumption 1. The prior distribution µ on H satisfies 2
Eµ exp(C1 ∥x∥ ) < ∞.
(6) for some C1 > 0.
Assumption 2. Let F be a convex, closed set of non-negative functions and assume that the true potential L is contained in F, that is L ∈ F. Assume also that there exits C2− , C2+ > 0 and for every ε > 0 there exists δ > 0 such that for every f ∈ F and all x ∈ H, y ∈ Rm (7)
2
2
2
2
−C2− − ε ∥x∥ + δ ∥y∥ ≤ f (x, y) ≤ C2+ (1 + ∥x∥ + ∥y∥ ).
The assumption (7) is stronger than non-negativity of f . It will ensure that the marginal distribution of y under f is well defined (see Equation (9) and Lemma A.1). The following examples give convenient sets of sufficient conditions under which Assumption 2 holds for Gaussian likelihood potentials. Example 2.1. Let A be a given set of potential forward maps, and let the noise covariance Γ ∈ Rm×m be such that the following hold uniformly for all A ∈ A: (B1) (Lipschitz) There exists C3 > 0 such that for all x, x′ ∈ H ∥A(x) − A(x′ )∥Γ ≤ C3 ∥x − x′ ∥ . (B2) (A is proper) There exist R, C4 > 0 such that µ (B(0, R)) > 0 and
sup x∈B(0,R)
∥A(x)∥Γ < C4 .
Then by [10, Lemma 3.10], the class 1 2 Fresidual := f : (x, y) 7→ ∥A(x) − y∥Γ A ∈ A 2 satisfies Assumption 2. Example 2.2. Let A be a given set of potential forward maps that satisfies (B1) and (B2) in Lemma 2.1 and G a set of potential covariances such that (B3) (Bounded eigenvalues) The eigenvalues of Γ ∈ G are all uniformly bounded from below and above, that is there exist 0 < C5− , C5+ such that for all Γ ∈ G C5− ≤ γmin (Γ) ≤ γmax (Γ) ≤ C5+ , where γmin (Γ) and γmax (Γ) denote the minimum and maximum eigenvalues of Γ respectively.
6
SCHNEIDER, HELIN, TAGHIZADEH
Then the class
1 2 ∥A(x) − y∥Γ A ∈ A, Γ ∈ G 2 is contained in a set that satisfies Assumption 2. Fcalibrated :=
f : (x, y) 7→
2.2. An equivalence class for likelihood functions. Definition 2.3. To any f ∈ F , we assign a posterior distribution µyf according to dµyf dµ
(x) =
1 exp (−f (·; y)) in L1 (µ) Zf (y)
where Z (8)
Zf (y) :=
exp (−f (x; y)) µ(dx) H
with µ-almost everywhere in x. We define the marginal distribution πf (·) of y through (9)
πf (y) := R
Zf (y) . Zf (y)dy
Rm
and the joint distribution λf (dx, dy) by (10)
λf (dx, dy) = πf (dy)
e−f (x,y) µ(dx). Zf (y)
R According to Lemma A.1, ∞ > Zf (y)dy > 0, so Equation (9) is well-defined. The likelihoods exp(−L) or exp(−f ) do in general not integrate to one. Instead normalization is enforced by formally dividing Z(y), so any functions that differ by an additive constant will give rise to the same posterior distribution. We define now an equivalence class of functions FΦ , where f ∼ g for f, g ∈ F if for all x, y ∈ H × Rm f (x, y) + log Zf (y) = g(x, y) + log Zg (y). We equip FΦ with the metric ∥f − g∥Φ := ∥f (·, ·) − g(·, ·) + log Zf (·) − log Zg (·)∥L1 (H×Rm ,λ) . (11) Notice that Φ and µyf and λf are indeed independent of the choice of the representative in FΦ . We shall view the map Φ as a mapping of an equivalence class FΦ → R from now on. Lemma 2.4. For f, g ∈ FΦ , it holds that f ∼ g iff f − g depends only on y on the support of λ. Proof of Lemma 2.4. Assume f ∼ g. Then for any (x, y) ∈ supp(λ) f (x, y) − g(x, y) = log Zg (y) − log Zf (y). Since the right hand-side depends only on y so does the left hand-side. Suppose now that f = g+h, where h depends only on y on supp(λ). Then Z f (x, y) + log Zf (y) = g(x, y) + h(y) + log exp(−g(x, y) − h(y))µ(dx) H Z = g(x, y) + log exp(−g(x, y))µ(dx), H
that is f ∼ g.
□
The following lemmas ensure that ∥f − g∥Φ is finite for f, g ∈ F and will be used later. Lemma 2.5. There exists universal constants C, κ > 0 such that for every f ∈ F 2
2
Eλf eκ(∥x∥ +∥y∥ ) < C Lemma 2.6. There exists p0 > 1 such that for all p ≥ p0 and fixed h ∈ F, there exists a universal constant C > 0 and for any f, g ∈ F it holds Eπh |log Zf (y) − log Zg (y)| ≤ C ∥f − g∥Lp (λh (·,·))
CONVEX NEURAL LIKELIHOOD APPROXIMATION
7
By Lemma 2.6 the metric ∥f − g∥Φ is dominated by ∥f − g∥Lp (λ) . Since λ is sub-Gaussian by Lemma, 2.5, the norm ∥f − g∥Lp (λ) is finite. Bounding the expected difference of the normalization terms log Zf (y) will be crucial in the following. 2.3. Main properties. The first observation that we make, is that if we have access to the joint distribution λ we can compute the KL-divergence in (3) without knowledge of the mathematical expression for L, up to additive constants (independent of f ). This observation is used in practice for neural likelihood approximation in parametric classes or neural flows (see [34, 35]). Lemma 2.7. Let Assumptions 1 and 2 hold. Then Eπ DKL µy ||µyf = Eπ DKL (µy ||µ) + Eλ [f (x; y)] + Eπ [log Zf (y)] < ∞. Proof of Lemma 2.7. The decomposition is immediate from !# " !# " # " dµyf dµy dµy λ π y λ π µy (x) = E log (x) = E DKL (µ ||µ) − E log . log E E dµyf dµyf dµ We show that all quantities are finite. Since by Lemma 2.5, λ = λL is sub-Gaussian and for f ∈ F, it holds Eλ f (x, y) < ∞. By using the lower bound in (7) together with (6), 2 2 Eπ log Eµ exp(−f (x, y)) ≤ Eπ log Eµ exp C2− + C1 ∥x∥ − δ(C1 ) ∥y∥ h i 2 2 = C2− − δ(C1 ) Eπ ∥y∥ + log Eµ exp(C1 ∥x∥ ) < ∞. Finally, by using L ≥ 0, the upper bound in (7) and Jensen’s inequality to the concave function x 7→ log x Eπ DKL (µy ||µ)
≤
−Eλ [L(x, y)] − Eπ log [Eµ exp (−L(x, y))] h i 2 2 −Eπ log Eµ exp −C2+ − C2+ ∥x∥ − C2+ ∥y∥
≤
C2+ (1 + Eπ ∥y∥ ) + C2+ Eµ ∥x∥ .
=
2
2
□ Motivated by Lemma 2.7, we define an objective function Φ as Φ : FΦ → R,
(12)
f 7→ Eλ [f (x; y)] + Eπ [log Zf (y)],
where Zf (y) is given in Equation (8). The following result is motivated by an analogous result in [11] for the covariance in a Gaussian likelihood with fixed mean. We now show that Φ is convex as a function of the NLL. Theorem 2.8. Let Assumptions 1 and 2 hold. Then the function Φ given by (12) is strictly convex. Proof of Theorem 2.8. It holds for any α + β = 1, α, β > 0 and f, g ∈ F that h i α β Φ(αf + βg) = α Eλ f (x, y) + β Eλ g(x, y) + Eπ log Eµ exp (−f (x, y)) exp (−g(x, y)) ≤
α Eλ f (x, y) + β Eλ g(x, y) +α Eπ [log Eµ exp (−f (x, y))] + β Eπ [log Eµ exp (−g(x, y))]
=
α Φ(f ) + β Φ(g),
where Hölder-inequality was used. We notice further that Hölder-inequality is an equality iff there exists c(y) > 0 such that µ-almost everywhere exp (−f (x, y)) = exp (−g(x, y)) c(y), i.e. f = g in ΦF .
□
8
SCHNEIDER, HELIN, TAGHIZADEH
Remark 2.9. The convexity of Φ is a property of the objective as a functional of the negative log-likelihood f , rather than of a particular parametrization. Consequently, if f is represented by a neural network, the optimization problem in the network parameters remains non-convex. Nevertheless, strict convexity implies that the population objective admits a unique minimizer in FΦ , providing an identifiability result that underlies the consistency analysis in Section 3. Different parameterizations may therefore have many local minima in parameter space while still representing the same optimal likelihood function. 3. Consistency and data-driven methods In this section, we consider an approximation of the integral in the objective function with (finite) data. Our goal will be to show the consistency of the NLL approximation. To this end, we recall the definition of consistency and a classical theorem from [47]. For a given domain D, let ϕN : D → R be random functions and ϕ : D → R a fixed function, denote by f0 := argminf ∈D ϕ(f ) and let d be a metric on D. We say a sequence of estimators fˆN is (asymptotically) consistent, if fˆN converges to f0 in probability. Theorem 3.1 (Consistency; [47]). Assume there exists f0 ∈ D such that for every ε > 0 (13)
inf d(f,f0 )>ε
ϕ(f ) >ϕ(f0 ).
and sup ϕN (f )−ϕ(f ) → 0
(14)
f ∈D
in probability. Then any fˆN satisfying ϕN (fˆN ) ≤ inf f ∈D ϕN (f ) + oP (1) converges in probability to f0 , where oP (1) converges to zero in probability. The addition of the term oP (1) in the preceding theorem allows for incomplete minimization of the objective. For example one could choose the term 1/N , where N is the data size and stop minimization once in a 1/N neighbourhood of the minimum. Our goal is to use the preceding theorem for a suitable data-driven estimator ΦN of Φ. To this end, we first make some necessary assumptions. 3.1. Assumptions. Definition 3.2. For s > 0 integer and β ∈ R, let F be a bounded subset of n o 2 C s (Rn+m , β) := f : f (1 + ∥z∥ )β/2 ∈ C s (Rn+m ) , where C s (Rn+m ) is the space for which ∥f ∥s,∞ :=
X
∥Dα f ∥∞ < ∞
0≤|α|≤s
is finite. Assumption 3. Fix the discretization H = Rn for some ∞ > n ∈ N and let F be a bounded subset of functions in C 1 (Rn × Rm , −2). Example 3.3. Assume that the assumptions of Example 2.2 hold for A and G, let H = Rn and (B4) (Differentiability) Any A ∈ A is differentiable and there exists C6 > 0 such that ∥∇x A(x)∥ ≤ C6 ∥x∥ . Then the set Fcalibrated =
f : (x, y) 7→
1 2 ∥A(x) − y∥Γ A ∈ A, Γ ∈ G 2
is contained in a set of function that satisfies Assumptions 2 and 3.
CONVEX NEURAL LIKELIHOOD APPROXIMATION
9
3.2. Data-driven methods. Now we are ready to present a data–driven estimator ΦN of Φ. To this end, we assume that the acquisition of joint data is restricted by computational cost or physical constraints and we can only acquire a limited ensemble (xi , yi )N i=1 ∼ λ of N independent, identically distributed (iid) (supervised) training data. We further assume that data from the prior is obtained easily and in abundance so for each yi , we can get an iid (unsupervised) training M (N ) data (xi,j )j=1 ∼ µ of size M (N ) for a coercive function M : N → N. We introduce a double loop Monte-Carlo estimator for the term Λ(f ) := Eπ log Eµ e−f (x,y) , which is approximating the expectation with respect to the marginal π of y and the prior µ, by N
ΛN (f ) :=
1 X M (N ) log Zf (yi ), N i=1
M
ZfM (yi ) :=
1 X −f (xi,j ,yi ) e . M j=1
An empirical estimator ΦN of Φ is given by N
ΦN : FΦ → R,
f 7→
1 X f (xi , yi ) + ΛN (f ), N i=1
Almost sure convergence of the first summand is straight forward from the law of large numbers. Convergence properties of nested Monte Carlo are a lot more involved (see the discussions in [36, 4]). In particular, we leave the choice of an optimal function M open. The following theorem establishes convergence in probability of the nested Monte-Carlo estimator. Theorem 3.4. Let f ∈ F and Assumptions 1 and 2 hold. Further suppose (xi , yi )N i=1 ∼ λ and M (N ) N ((xi,j )j=1 )i=1 be iid. Then the convergence lim ΛN (f ) = Λ(f )
N →∞
holds in probability. 3.3. Consistency. Our goal in this section is to establish consistency of a data-driven approx∗ towards the true potential L. To do so we verify the two hypotheses of imate minimizer fN Theorem 3.1 for ϕ = Φ, ϕN = ΦN and d = ∥·∥Φ . In Theorem 3.5, we verify the condition of a well-separated minimum (13) and in Theorem 3.6 a uniform law of large numbers (14) for the empirical objective ΦN . Combining the two then yields the consistency result, Theorem 3.7, which is the main result of this section. Since existence of an exact minimizer of ΦN over FΦ cannot be guaranteed in general, we work ∗ satisfying throughout with approximate minimizers fN (15)
∗ ΦN (fN ) ≤ inf ΦN (f ) + oP (1), f ∈FΦ
where oP (1) → 0 in probability as N → ∞ (recall from the discussion after Theorem 3.1 that this allows, for instance, minimization only up to tolerance 1/N ). Theorem 3.5. Let Assumptions 1, 2 and 3 hold. Then for any f ∈ FΦ inf
∥f −L∥Φ >ε
Φ(f ) > Φ(L).
Theorem 3.6. Let Assumptions 1, 2 and 3 hold. Then sup |ΦN (f ) − Φ(f )| → 0 f ∈F
in probability. Theorems 3.5 and 3.6 verify conditions (13) and (14) of Theorem 3.1 with f0 = L (which is indeed the unique minimizer of Φ over FΦ , by Lemma 2.7 together with Lemma 2.4). Theorem 3.1 then applies directly and yields the following consistency guarantee for the approximate minimizer ∗ fN defined in (15).
10
SCHNEIDER, HELIN, TAGHIZADEH
∗ Theorem 3.7. Let fN be given by (15) and let Assumptions 1, 2 and 3 hold. Then for every ε>0 ∗ P(∥fN − L∥Φ > ε) → 0 as N → ∞,
where the symbol P(A) denotes the probability of a measurable event A under the probability measure λ. To prove Theorem 3.6, we need the following results on bracketing numbers. Definition 3.8. An ε-bracket [l, u] in LpR(λ) with p ≥ 1 with functions l ≤ u is defined as the set of all functions f such that l ≤ f ≤ u and (u − l)p dλ ≤ εp . The bracketing number N[] (ε, F, Lp (λ)) is the minimum number of ε-brackets needed to cover F . We recall the following powerful result in [29, Corollary 3] on bracketing numbers for bounded subsets of Hölder spaces. Theorem 3.9. Let s > 0 be integer, p ≥ 1 and β ∈ R. Let F be a bounded subset of C s (Rn+m , β) and 2
(1 + ∥x∥ )(γ−β)/2
(16)
<∞ p,λ
for some γ > 0; γ ̸= s. Then log N[] (ε, F, Lp (λ)) < Cε−(n+m)/ min(s,γ) . In particular, if Assumptions 3 are satisfied, the bracketing number N[] (ε, F, Lp (λ)) is finite for every p ≥ 1. (Choose any γ > 0. Since λ is sub-expontinental, the moments in (16) are finite.) Remark 3.10. It is difficult to apply classical rate-of-convergence results for M-estimators [46] to our setting. The difficulty is twofold. First, the second term in ΦN contains the nested Monte Carlo estimator ΛN (f ), which is not a linear empirical average. Although Theorem 3.4 proves convergence of ΛN (f ), it does not provide quantitative or uniform rates over FΦ . Second, the classical theory is based on localized empirical processes controlled by L2 (λ) bracketing entropy, whereas our analysis relies on Lp (λ) estimates (Lemma 2.6) for possibly large p. Extending the existing theory to this setting is left for future work. 4. Numerical demonstration In this section we consider the deblurring problem and the problem of recovering the doping profile in a semiconductor. Our goal is to test the free-form, residual and calibrated residual approximations of the negative-log-likelihood as given in Section 2. In the first example, we further compare with a GP regression of the likelihood. The GP assumes a squared exponential covariance with learnable length-scale l ∈ R and variance level σ 2 ∈ R. To emphasize the learnable parameters, we denote by f θ : Rn → R or g θ : Rn → Rm a neural-network with learnable parameters θ and by Σθ ∈ Rm×m a learnable covariance matrix. The underlying output functions and assumptions made are summarized in Table 1. Table 1. Parametrization and required data for the NLL approximations and the GP regression. Training the free-form and calibrated residual-approximation requires only joint data, residual approximation needs additionally the noise level and GP regression needs access to evaluations of the true NLL.
output required data learnable parameters
Free-form approximation f θ (x, y) (x, y)- data
calibrated residualapproximation 2 1 θ 2 g (x) − y Σθ (x, y)- data
fθ
g θ , Σθ
residualapproximation 2 1 θ 2δ 2 g (x) − y (x, y)- data and δ 2 gθ
GP regression GP(x, y) (x, y)- data and L(x, y) σ2 , l
CONVEX NEURAL LIKELIHOOD APPROXIMATION
11
4.1. Deblurring. We consider a de-blurring problem as the first example for the numerical experiments. Define the corresponding Bayesian inverse problem as following: The unknown x ∈ H is supported on the Hilbert space H = L2 (D) with D = [0, 1]. The prior is given by a GP N (0, C) with squared exponential kernel C: 1 1 2 C(x, y) = 2 exp − 2 ∥x − y∥ , x, y ∈ D, σ 2l where the parameter σ = 0.1 models the variance and l = 0.1 the length-scale of the process. The forward map A : L2 (D) → Rm blurs any given unknown x ∈ L2 (H) by computing the convolution Z (x ∗ k) (t) = x(s)k(t − s) ds [0,1]
for a given kernel function k. This convolution is then evaluated at a finite grid (t1 , · · · , td ) ⊂ Rm with m = 50: (x ∗ k)(t1 ) .. d Ax = ∈R . . (x ∗ k)(td ) We choose k to be the PDF of a Gaussian centred random variable with variance two. The measurement y ∈ Rm is obtained by corrupting Ax with additive noise: y = Ax + ε, 2
where observational noise ε ∼ N (0, δ Id ) is Gaussian with white noise statistics and noise level δ = 0.1. The free-form NLL approximation uses 4 fully connected layers, the residual approximation 5 fully connected layers and the calibrated residual approximation uses 5 fully connected layers followed by a diagonal layer that approximates the noise level. We deploy a Markov Chain Monte Carlo (MCMC) algorithm, namely the preconditioned CrankNicolson (pCN) [28, 45] algorithm to sample from the posterior distributions. In the pCN algorithm, a proposal vn given a previous sample un are generated according to p p vn = 1 − 2βun + 2βζ ζ ∼ µ, where µ is a Gaussian prior and β ∈ (0, 1) a hyperparameter. We choose β = 0.2 in this Section. The proposal vn =: un+1 is then accepted with probability pn = min (exp(L(vn , y) − L(un , y), 1) , where L is the NLL. In Figure 1, we first analyse the different NLL approximations. We demonstrate that the posterior mean of the free-form, residual and calibrated residual approximations are very close to the ground truth and contained in a given confidence interval. The free-form approximation is slightly closer to the ground truth with slightly smaller confidence band compared to the residual approximation and calibrated residual approximations. The approximations have slightly lower accuracy than using the true posterior and wider confidence bands. The calibrated residual approximation is furthest from the ground truth. Exact posterior
Free-form approximation
Residual approximation
Calibrated residual approximation
Figure 1. Posterior MCMC samples for the deblurring problem. From left to right: true posterior, free-form approximation, residual approximation and calibrated residual approximation. All methods use only 1, 000 training points.
12
SCHNEIDER, HELIN, TAGHIZADEH
In Figure 2, we illustrate that GP regression, even though it utilizes the additional information of evaluations of the true likelihood, struggles to accurately describe the likelihood especially with 1,000 training points. Even with 10,000 training points, it is less accurate than all of the NLL approximations. This experiment demonstrates that neural likelihood approximation is a highly Exact posterior
GP regression 1,000 training points
GP regression 10,000 training points
Figure 2. Posterior MCMC samples for the deblurring problem. From left to right: true posterior, GP regression with 1,000 training points and GP regression for the likelihood with 10,000 training points. data-efficient way to accurately approximate the likelihood in inverse problems. In settings where the observational noise is unknown, the free-form approximation seems preferable over a calibrated residual approximation based on this experiment. 4.2. Semiconductor devices. A non-invasive method for estimation of the doping profile inside a semiconductor is the identification by voltage-current measurements. For a thorough introduction we refer to [45] and the references therein. In short, the domain Ω = [−1, 1]2 models the crosssection of a pn-diode with Ohmic contacts ΓN = [−1, 1] × {+1} and ΓP = [−1, 1] × {−1}. Ω is separated into the sets ΩN and ΩP by a known junction line Γ. We assume that the unknown function C is given by two independent Gaussian random fields CN , CP separated by the junction line Γ: ( CN (x) x ∈ ΩN (17) C(x) = Cp (x) x ∈ Ωp . Forward model. We consider the following model based on the (close to) thermal equilibrium for the electric potential Ve given a concentration C(x) λ2 △Ve = δ 2 eVe − e−Ve − C in Ω (18)
Ve = Vbi ∇Ve · ν = 0
on ∂ΩD on ΩN .
Given Ve , we observe the quantity δ 2 eVbi ∇û · ν|ΓP , where û solves the continuity equation div (γ(x)∇û(x)) = 0 in Ω (19)
û = U on ∂ΩD ∇û · ν = 0 on ∂ΩN
and γ = eVe . We fix the parameters Vbi = 0.6, λ = δ = 1 and U = 2. Prior model. We model the spatial functions CN and Cp as independent Gaussian fields with Matern–Whittle covariance operator given ν ∥x − y∥2 21−ν ∥x − y∥2 K(x, y) = σ 2 Kν , Γ(ν) l l where Γ denotes the Gamma function, Kν the modified Bessel function and the parameters σ, ν, l denote the variance, smoothness and length scale of the process. We choose the parameters σ = 2, l = 20 and ν = 0.9.
CONVEX NEURAL LIKELIHOOD APPROXIMATION
13
Training. Training data is given by N = 10, 000 observations from the joint distribution of λ with a resolution of 100 × 100 for C and measurement resolved at 100 pixels. The PDE is solved using a finite element solver in MATLAB and with a mesh-size of 0.1. Generating the training data takes ≈ 5 hours. The boundary measurement is corrupted with additive noise of level σ = 0.01. We assume that the true forward model is unknown and test the free-form, residual and calibrated residual approximations, all of which can be employed without knowing the forward map. We use similar number of parameters for the underlying neural networks. The parametrizations utilize four convolution layers followed by four fully connected layers. The calibrated residual approximation finally has an additional diagonal layer to approximate the noise level. Training is done for 30 epochs and takes roughly 120 seconds per epoch. The computational times using MATLAB’s PDE solver and the NLL approximation are reported in Table 2. The CPU-based runtimes correspond to the same computational hardware, whereas the GPU-based NLL approximation is performed on a GPU cluster. For the approximations, we display computational runtimes for parallel evaluations, which could be used to generate many chains in parallel, which is however outside the scope of our numerical experiment. On the same hardware, the NLL approximations are ≈ 200 times faster compared with the PDE model. Table 2. Computation time in seconds for the NLL and its surrogates. Whenever possible, the 100 and 1000 evaluations are made in parallel. For the free-form, residual and calibrated residual parametrizations, computational times are essentially identical. The NLL approximations are at least ≈ 200 times faster compared to the finite elements PDE model on the same device and up to 20, 000 time faster under parallelization on a GPU device. PDE solve (CPU) single evaluation 0.7008 100 evaluations 70.0083 1000 evaluations 700.8333
NLL approx. (CPU) 0.0042 0.3972 3.6880
NLL approx. (GPU) 0.0016 0.0036 0.0331
MCMC sampling. We fix a previously unseen pair (x, y) ∼ λ and use a preconditioned CrankNicolson (pCN) sampler to sample from the posterior conditional on the measurement y. In the pCN algorithm, we use a truncated KL-expansion to sample from the Gaussian random field. For all approaches, we generate 10, 000 samples sequentially, discard the first 5, 000 as burn-in and finally thin the chain by factor five. The computational effort for posterior sampling using the NLL approximation varies due to the underlying posterior variance and acceptance rate. The computational times for the NLL approximation is < 50 minutes on CPU and < 10 minutes on a GPU cluster. In contrast, utilizing the discretized PDE model on a CPU device, posterior sampling takes 265 hours. For the NLL approximation, we require the offline computations of data generation and training, which together take less than 6 hours. In this time it would be possible to generate approximately 440 samples with the PDE model. Posterior means, biases, and variances are summarized in Figure 3. We see that the biases and variances of the posterior with the true PDE model and the free-form and residual approximations are very similar. The spatial L1 norm of the biases are given by 0.1807 (exact), 0.1839 (free-form) and 0.1954 (residual), confirming that the bias of the true PDE model is extremely close to the free-form and residual approximations. In contrast, the calibrated residual approximation fails to generate accurate posterior variance and means. This experiment demonstrates that neural likelihood approximation can lead to enormous computational acceleration in PDE-based, non-linear inverse problems. The posterior summary statistics, such as bias and variance, are highly accurate for the free-form and residual approximations. In settings where the observational noise level is unknown, the free-form approximation can recover the posterior with high fidelity. The approximation methods require only joint data for training and no evaluations of the true forward map or likelihood and are thus highly versatile.
14
SCHNEIDER, HELIN, TAGHIZADEH
True C
Exact, PDE-based NLL
Free-form approximation
Residual approximation
Calibrated residual approximation
Figure 3. Summary of MCMC samples. Top row: True quantity of interest followed by posterior mean (left), bias (middle) and posterior variance (right) of the posterior samples with the true PDE-based NLL and its approximations. The L1 norm of the biases are given by 0.1807 (exact), 0.1839 (free-form), 0.1954 (residual), 1.072 (calibrated residual). 5. Conclusions This work developed a theoretical foundation for neural likelihood approximation by overcoming limitations imposed by restrictive parametric model classes. We introduced a normalization scheme in which an unconstrained function is mapped to a valid probability density by division by its expectation, thereby embedding normalization directly into a KL-based learning objective. This eliminates the need to enforce normalization through architectural constraints while preserving
CONVEX NEURAL LIKELIHOOD APPROXIMATION
15
full flexibility of the approximative class. We showed that the resulting objective is convex in the un-normalized NLL. Moreover, we established consistency of the estimator under empirical approximation of the normalization integral, ensuring convergence to the true target in the largesample regime. We demonstrated the practical effectiveness of the proposed approach on a deblurring task and a non-linear PDE-based imaging problem, where only joint (x, y) training data is available but the forward model is assumed unknown. We studied different parametrizations of the NLL and found that a free-form neural-network based approximation can accurately approximate the NLL, even without knowledge of the observational noise level, with significant computational acceleration. A residual-based parametrization struggles when the observational noise level is not known but performs well with known observational noise. Quantification of the convergence through convergence rates remains an open problem. A further current limitation is that the consistency result is established under finite-dimensional spaces supporting the unknown and measurement due to the scope of existing empirical process theory used in the analysis. Extending these results to fully infinite-dimensional function classes remains an important direction for future work. Appendix A. Auxiliary results To prove Lemma 2.5, we show the following assertion, which is a counterpart to [10, Lemma 3.11] Lemma A.1. There exists universal constants C − , C + such that for all f ∈ F +
2
2
C − e−C2 ∥y∥ ≤ Zf (y) ≤ C + e−δ(C1 )∥y∥ . Proof. Using (7) with ε = C1 h i h i − 2 2 Zf (y) = Eµ e−f (x;y) ≤ eC2 e−δ(C1 )∥y∥ Eµ eC1 ∥x∥ . Using the lower bound of (7): i h h i + + 2 2 Zf (y) = Eµ e−f (x;y) ≥ e−C2 (1+∥y∥ ) Eµ e−C2 ∥x∥ . □ Proof of Lemma 2.5. We choose κ := min (C1 , δ(C1 )) /3. It holds h i h i1/2 h i1/2 2 2 2 2 Eλf eκ(∥x∥ +∥y∥ ) ≤ Eλf e2κ∥x∥ Eλf e2κ∥y∥ . We will now bound the terms on the right-hand side in the previous equation. Using Lemma A.1: h i h i R e2κ∥y∥2 Z (y)dy f m λf 2κ∥y∥2 πf 2κ∥x∥2 E e = E e = R R Z (y)dy Rm f R + 2κ∥y∥2 −δ(C1 )∥y∥2 C e e dy Rm R ≤ . Z (y)dy Rm f By Lemma A.1 Z
Z
Zf (y)dy ≥ C −
(20) Rm
+
2
e−C2 ∥y∥ dy > 0,
Rm
we have λf
E
h
2κ∥y∥2
e
i
R ≤
Rm
2
2
C + e2κ∥y∥ e−δ(C1 )∥y∥ dy . R + 2 C − e−C2 ∥y∥ dy Rm
By construction, 2κ < δ(C1 ) so the integral is finite. Regarding the second term, it holds by (10), (7) and (20): − h i h i −f (x,y) 2 eC2 λf 2κ∥x∥2 πf ⊗µ 2κ∥x∥2 R e e ≤ E e =E Eµ e2.5κ∥x∥ R + 2 Z (y)dy C − Rm e−C2 ∥y∥ Rm f
16
SCHNEIDER, HELIN, TAGHIZADEH
This term is indeed finite for κ ≤ C1 /3 by (6).
□
To prove Lemma 2.6, we need the following results. Lemma A.2. Let r < 1 + Cκ+ . Then there exists a C(r) > 0 such that for all f, g ∈ F 2
πg (y) < C(r). πf (y) Lr (πf ) Proof. We use the definition of πh for h ∈ {f, g} in (9) to get Z −r Z r−1 Z Z πg (y)r πf (y)1−r dy ≤ Zg (y)dy Zf (y)dy Rm
Rm
Rm
Zg (y)r Zf (y)1−r dy
Rm
We deal with the integral first. By Lemma A.1, Z Z + 2 r 1−r + r − 1−r e(−rδ(C1 )+(r−1)C2 )∥y∥ Zg (y) Zf (y) ≤ (C ) (C ) m m R ZR + 2 + r − 1−r ≤ (C ) (C ) e(r−1)C2 ∥y∥ . Rm
Since (r − 1)C2+ < κ, this term is indeed uniformly bounded by Lemma 2.5. The second factor is now dealt with by utilizing Lemma A.1: Z −r Z r−1 Z − Zg (y)dy Zf (y)dy ≤ C Rm
Rm
e
−C2+ ∥y∥2
−r Z + dy C
Rm
e
−δ(C1 )∥y∥2
r−1 dy
< ∞.
Rm
□ Lemma A.3. For q0 := 1 + 2Cκ+ r′ and any 1 < q < q0 such that for given g ∈ F and all f ∈ F 2
we have the there exists C(q) such that Eλg
dµyf
!q < C(q)
dµyg
2C +
Proof. In a first step, we introduce the Hölder exponents r := 1 + 2Cκ+ , r′ = 1 + κ2 and write 2 " ! # " ! # ! ! y q y q y q−1 y q−1 dµ dµ dµ dµ y y y f f f f = Eπf dπg Eµf Eλg = Eπg Eµg = Eπg Eµf dµyg dµyg dµyg dπf dµyg ≤
Eλf
dµyf
!r′ (q−1) 1/r′
dµyg
Eπf
πg (y) πf (y)
r 1/r .
The second term involving the integral over πf is uniformly bounded by Lemma A.2. To deal with the other term, we realize that f (x, y) ≥ 0 and hence Zg (y) ≤ 1 for all (x, y) ∈ H × Rm such that dµyf dµyg
r ′ (q−1)
′
′
=
er (q−1)(g(x,y)−f (x,y)) er (q−1)g(x,y) ′ (q−1) ′ (1−q) ≤ r r Zf (y) Zg (y) Zf (y)r′ (q−1)
≤
1 r ′ (q−1)C2+ (∥x∥2 +∥y∥2 ) r ′ (q−1)C2+ ∥y∥2 e . C− e
Equation (7) together with Lemma A.1 was used. Observe that 2r′ (q − 1)C2+ ≤ κ. Therefore ′ y r (q−1) h i dµ 2 2 f ≤ 1− Eλf eκ(∥x ∥+∥y∥ ) , Eλf y C dµg which is indeed uniformly bounded by Lemma 2.5.
□
CONVEX NEURAL LIKELIHOOD APPROXIMATION
17
Proof of Lemma 2.6. Fix in the first step y ∈ Rm . Define Λy (f ) := log Eµ e−f (x,y) . For t ∈ [0, 1] we define kt = g + t(f − g) ∈ F , since F is convex. By the fundamental theorem of calculus, Z 1 Z 1 d µ −kt (x,y) e d y y y dt E . (21) Λ (f ) − Λ (g) = Λ (kt )dt = µ −k e t (x,y) dt 0 dt 0 E By Leibniz integral rule h i d µ h −kt (x,y) i = Eµ (f (x, y) − g(x, y)) e−kt (x,y) . E e dt
(22)
The domination function G ∈ L1 (µ) in Leibniz integral rule is given by 2
2
|f (x, y) − g(x, y)| e−kt (x,y) ≤ 2C2+ (∥x∥ + ∥y∥ ) := G(x). Putting together (21) and (22), we get Λy (f ) − Λy (g)
Z 1 =
y
Eµkt [f (x, y) − g(x, y)] dt.
0
Therefore by Fubinis’ Theorem and for any 1/p + 1/q = 1 y Z 1 dµkt (g(x, y) − f (x, y) dt Eπh |Λy (f ) − Λy (g)| = Eλh dµyh 0 ! Z y q 1/q λ 1/p 1 dµ p k t ≤ E h |g(x, y) − f (x, y)| dt. Eλh dµyh 0 Notice that kt ∈ F and we can apply Lemma A.3 and fix now q < q0 and the corresponding p that satisfies 1/p + 1/q = 1. □ Appendix B. Proofs for data-driven estimates Proof of Theorem 3.4. Proof. We introduce the decomposition (23)
N 1 X M (N ) ΛN (f ) − Λ(f ) = log Zf (yi ) − log Zf (yi ) + N i=1
N
1 X log Zf (yi ) − Λ(f ) N i=1
!
and show convergence to zero in probability for the first term (I) and the second term (II) separately. (I) We define for fixed i = 1, . . . , N M (N )
ξi M (N )
We show now that ξi M (N ) 2 sup Eµ⊗π (ξ1 ) N ∈N
≤
M (N )
:= log Zf
(yi ) − log Zf (yi ).
is uniformly integrable.
2 2 log Eµ⊗π e−f (x,yi ) µ(dx) + 2 sup Eµ⊗π log
N ∈N
2 M (N ) X 1 e−f (xi,j ,yi ) . M (N ) j=1
For the first term we use (6) and (7) to find 2 2 2 log Eµ⊗π e−f (x,yi ) µ(dx) ≤ C2− + log Eµ 2eC1 ∥x∥ µ(dx) < ∞. And for the second term, we find upper and lower bounds for the logarithm-term by using (6) and Jensen’s inequality: ( PM (N ) C1 /2∥xi ,j∥2 M (N ) 1 X ≤ C2− + log M (N 1 j=1 e ) −f (xi,j ,yi ) log e C + PM (N ) 2 2 + + M (N ) j=1 ≥ −C2 − C2 ∥yi ∥ − 2 j=1 ∥xi,j ∥ . M (N )
18
SCHNEIDER, HELIN, TAGHIZADEH
By using log z ≤ z for z ≥ 1 we find 2 M (N ) X 2 1 sup Eµ C2− + log eC1 /2∥xi,j ∥ M (N ) j=1 N ∈N
≤
2 2 2(C2− )2 + 2Eµ eC1 /2∥x∥
≤
2(C2− )2 + 2Eµ eC1 ∥x∥ < ∞.
2
Through Jensen’s inequality we find 2 (N ) + M X C 2 2 4 4 2 sup Eµ⊗π (1 + C2+ ) ∥yi ∥ + ∥xi,j ∥ ≤ 2(1 + C2+ )2 Eπ ∥y∥ + 2(C2+ )2 sup Eµ ∥x∥ . M (N ) j=1 N ∈N N ∈N M (N )
M (N )
In summary, we showed that supN ∈N Eµ (ξ1 )2 < ∞ and hence ξi is uniformly integrable. M (N ) Since ξ1 → 0 converges almost surely by the law of large numbers (notice that xi,j are all iid) and the continuity of the logarithm, and is uniformly integrable, it holds (24)
Eµ⊗π
N 1 X M (N ) M (N ) ≤ Eµ⊗π |ξ1 | → 0, ξ N i=1 i
N → ∞,
PN M (N ) which implies N −1 i=1 ξi → 0 in probability, as N → ∞. (II) By the weak law of large numbers it holds that ! N 1 X log Zf (yi ) − Λ(f ) → 0 N i=1 in probability.
□
Proof of Theorem 3.5. Proof. Step 1: Optimality of L. By Lemma 2.7 applied at f = L, where µyL = µy , we obtain for any f ∈ FΦ Φ(f ) − Φ(L) = Eπ DKL µy ||µyf ≥ 0, with equality if and only if µyf = µy for π-almost every y, that is, if and only if f ∼ L in FΦ by Lemma 2.4. Hence L is the strict minimizer of Φ over FΦ . Step 2: Compactness argument. By Theorem 3.9, for every η > 0 the set F can be covered by finitely many η-brackets [lk , uk ], k = 1, . . . , K, with ∥uk − lk ∥Lp0 (λ) ≤ η. Setting mk := 12 (lk + uk ), every f in the k-th bracket satisfies ∥f − mk ∥Lp0 (λ) ≤ 12 ∥uk − lk ∥Lp0 (λ) ≤ η2 , so {m1 , . . . , mK } is a finite η2 -net for F. Since η > 0 was arbitrary, F, and hence its closure F, is totally bounded in Lp0 (λ). As Lp0 (λ) is complete and F is closed, F is complete; being complete and totally bounded, F is compact in Lp0 (λ). By Lemma 2.6, ∥·∥Φ is dominated by C ∥·∥Lp0 (λ) , so the identity map Lp0 (λ), ∥·∥Lp0 (λ) → FΦ , ∥·∥Φ is continuous, and F Φ is compact with respect to ∥·∥Φ as well. Fix ε > 0 and define Kε := {f ∈ F Φ : ∥f − L∥Φ ≥ ε}. Kε is a closed subset of the compact set F Φ , hence compact. Step 3: Attainment. If Kε = ∅, the claim is clear. Hence we assume now Kε ̸= ∅. Since Φ is Lipschitz, hence continuous, with respect to ∥·∥Φ (by |Φ(f ) − Φ(g)| ≤ ∥f − g∥Φ ), it attains its infimum on the compact set Kε : there exists f0 ∈ Kε with Φ(f0 ) = inf Φ(f ). f ∈Kε
Since ∥f0 − L∥Φ ≥ ε > 0, we have f0 ̸= L in FΦ , so by Step 1, Φ(f0 ) > Φ(L) strictly. Therefore inf
∥f −L∥Φ ≥ε
Φ(f ) = Φ(f0 ) > Φ(L).
CONVEX NEURAL LIKELIHOOD APPROXIMATION
19
□ Proof of Theorem 3.6. Proof. Notice that N[] (ε, F, Lp0 (λ)) is finite for every ε > 0 by Theorem 3.9 and for p0 given in Lemma 2.6. We adapt the proof classical proof (e.g. in [40, Theorem 3.2]). Fix ε > 0 and choose p0 finitely many ε-brackets ([lk , uk ])K < εp0 for all k ≤ K ∈ N. k=1 that cover F and with λ|uk − lk | Without loss of generality, every lk , uk satisfies (7). For any f ∈ F there exists a bracket such that f ∈ [lk , uk ]. We fix f ∈ F and the corresponding bracket. We write 1 X f (xi , yi ) Ψ(f ) := Eλ f (x, y), ΨN (f ) := N M (N ) N i h X X 1 1 ΛN (f ) := Λ(f ) := Eπ log Eµ e−f (x,y) , log exp (−f (xi,j ; yi )) N i=1 M (N ) j=i such that |Φ(f ) − ΦN (f )|
≤
|Ψ(f ) − ΨN (f )| + |Λ(f ) − ΛN (f )|.
Notice since lk (x, y) ≤ uk (x, y) λ-almost everywhere, and Ψ, Λ are increasing or decreasing, it holds for Θ ∈ {Ψ, Λ} |Θ(f ) − Θ(lk )| + |ΘN (f ) − ΘN (lk )| ≤ |ΘN (uk ) − ΘN (lk )| + |Θ(uk ) − Θ(lk )|. Thus (25)
|ΘN (f ) − Θ(f )| ≤|ΘN (f ) − ΘN (lk )| + |ΘN (lk ) − Θ(lk )| + |Θ(lk ) − Θ(f )| ≤|ΘN (uk ) − ΘN (lk )| + |ΘN (lk ) − Θ(lk )| + |Θ(lk ) − Θ(uk )|.
The right hand side depends only on the finite indices k ≤ K but not on f ∈ F. Step 1: Since ∥·∥L1 (λ) ≤ ∥·∥Lp0 (λ) it holds (26)
|Ψ(uk ) − Ψ(lk )| = Eλ |uk (x, y) − lk (x, y)| ≤ ∥uk − lk ∥Lp0 ≤ ε.
Using Lemma 2.6, the definition of the brackets, and for C independent of k: (27)
|Λ(uk ) − Λ(lk )| ≤ C ∥uk − lk ∥Lp0 = Cε
Step 2: For Θ ∈ {Ψ, Λ} and each fixed k = 1, . . . , K we have by the law of large numbers and Theorem 3.4 |ΘN (lk ) − Θ(lk )| + |ΘN (uk ) − ΘN (lk )| → |Θ(uk ) − Θ(lk )| ≤ (C + 1)ε, with convergence in probability. Combining the estimates. Since the right-hand side of (25) does not depend on f ∈ F, we obtain i X h sup |ΦN (f ) − Φ(f )| ≤ sup |ΘN (uk ) − ΘN (lk )| + |ΘN (lk ) − Θ(lk )| + |Θ(lk ) − Θ(uk )| . f ∈F
k≤K
Θ∈{Ψ,Λ}
P By Step 1, the last term satisfies supk≤K PΘ |Θ(lk ) − Θ(uk )| ≤ (C + 1)ε. By Step 2, the first two terms converge in probability to supk≤K Θ |Θ(uk ) − Θ(lk )| ≤ (C + 1)ε, since the supremum is over a finite set. Hence for any δ > 0, ! P sup |ΦN (f ) − Φ(f )| > 2(C + 1)ε + δ
−→ 0.
f ∈F
P
Since ε > 0 was arbitrary, this gives supf ∈F |ΦN (f ) − Φ(f )| − → 0.
□
20
SCHNEIDER, HELIN, TAGHIZADEH
References [1]
[2] [3]
[4] [5]
[6]
[7]
[8] [9]
[10]
[11] [12]
[13]
[14] [15] [16] [17]
[18]
[19]
S R Arridge et al. “Approximation errors and model reduction with an application in optical diffusion tomography”. In: Inverse Problems 22.1 (Jan. 2006), p. 175. doi: 10.1088/02665611/22/1/010. Simon R Arridge. “Optical tomography in medical imaging”. In: Inverse Problems 15.2 (1999), R41. Ravi Bansal and Amir Yaron. “Risks for the Long Run: A Potential Resolution of Asset Pricing Puzzles”. In: The Journal of Finance 59.4 (2004), pp. 1481–1509. doi: https : //doi.org/10.1111/j.1540-6261.2004.00670.x. Arved Bartuska et al. Double-loop randomized quasi-Monte Carlo estimator for nested integration. 2025. arXiv: 2302.14119 [math.NA]. url: https://arxiv.org/abs/2302.14119. Nada Cvetković and Han Cheng Lie. Upper and lower bounds for local Lipschitz stability of Bayesian posteriors. 2025. arXiv: 2505.23541 [math.ST]. url: https://arxiv.org/abs/ 2505.23541. Michael Deistler, Pedro J Goncalves, and Jakob H Macke. “Truncated proposals for scalable and hassle-free simulation-based inference”. In: Advances in Neural Information Processing Systems. Ed. by S. Koyejo et al. Vol. 35. Curran Associates, Inc., 2022, pp. 23135–23149. Teo Deveney, Eike H. Mueller, and Tony Shardlow. “Deep Surrogate Accelerated DelayedAcceptance Hamiltonian Monte Carlo for Bayesian Inference of Spatio-Temporal Heat Fluxes in Rotating Disc Systems”. In: SIAM/ASA Journal on Uncertainty Quantification 11.3 (2023), pp. 970–995. doi: 10.1137/22M1513113. David T. Frazier et al. The Statistical Accuracy of Neural Posterior and Likelihood Estimation. 2024. arXiv: 2411.12068 [stat.ML]. url: https://arxiv.org/abs/2411.12068. Michael U. Gutmann and Jukka Corander. “Bayesian Optimization for Likelihood-Free Inference of Simulator-Based Statistical Models”. In: Journal of Machine Learning Research 17.125 (2016), pp. 1–47. url: http://jmlr.org/papers/v17/15-017.html. Tapio Helin, Youssef Marzouk, and Jose Rodrigo Rojo-Garcia. Bayesian optimal experimental design with Wasserstein information criteria. 2025. arXiv: 2504.10092 [stat.ME]. url: https://arxiv.org/abs/2504.10092. Tapio Helin and Fabian Schneider. Optimal approximation error approach in Bayesian inverse problems. In preparation. 2026. Tapio Helin et al. “Introduction to Gaussian Process Regression in Bayesian Inverse Problems, with New Results on Experimental Design for Weighted Error Measures”. In: July 2024, pp. 49–79. isbn: 978-3-031-59761-9. doi: 10.1007/978-3-031-59762-6_3. Lukas Herrmann, Christoph Schwab, and Jakob Zech. “Deep neural network expression of posterior expectations in Bayesian PDE inversion”. In: Inverse Problems 36.12 (Dec. 2020), p. 125011. doi: 10.1088/1361-6420/abaf64. Yuga Hikida et al. Multilevel neural simulation-based inference. 2025. arXiv: 2506.06087 [stat.ML]. url: https://arxiv.org/abs/2506.06087. D.S. Holder. “Electrical Impedance Tomography: Methods, History and Applications”. In: vol. 32. Jan. 2005. isbn: 0750309520. J. Kaipio and E. Somersalo. Statistical and Computational Inverse Problems. Applied Mathematical Sciences. Springer New York, 2006. isbn: 9780387271323. Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell. “Simple and Scalable Predictive Uncertainty Estimation using Deep Ensembles”. In: Advances in Neural Information Processing Systems. Ed. by I. Guyon et al. Vol. 30. Curran Associates, Inc., 2017. Matti Lassas, Eero Saksman, and Samuli Siltanen. “Discretization-invariant Bayesian inversion and Besov space priors”. In: Inverse Problems and Imaging 3.1 (2009), pp. 87–122. issn: 1930-8337. doi: 10.3934/ipi.2009.3.87. Markku Lehtinen, Lassi Päivärinta, and Erkki Somersalo. “Linear inverse problems for generalised random variables”. English. In: Inverse Problems 5.4 (1989), pp. 599–612. issn: 0266-5611. doi: 10.1088/0266-5611/5/4/011.
REFERENCES
[20]
[21]
[22]
[23]
[24]
[25]
[26]
[27]
[28] [29]
[30]
[31] [32]
[33]
[34]
[35] [36]
[37]
21
Yongchao Li, Yanyan Wang, and Liang Yan. “Surrogate modeling for Bayesian inverse problems based on physics-informed neural networks”. In: Journal of Computational Physics 475 (2023), p. 111841. issn: 0021-9991. doi: https://doi.org/10.1016/j.jcp.2022.111841. H. C. Lie, T. J. Sullivan, and A. L. Teckentrup. “Random Forward Models and LogLikelihoods in Bayesian Inverse Problems”. In: SIAM/ASA Journal on Uncertainty Quantification 6.4 (2018), pp. 1600–1629. doi: 10.1137/18M1166523. Jan-Matthis Lueckmann et al. “Likelihood-free inference with emulator networks”. In: Proceedings of The 1st Symposium on Advances in Approximate Bayesian Inference. Ed. by Francisco Ruiz et al. Vol. 96. Proceedings of Machine Learning Research. PMLR, Dec. 2019, pp. 32–53. Ramona Maraia et al. “Bayesian synthetic likelihood for stochastic models with applications in mathematical finance”. In: Frontiers in Applied Mathematics and Statistics Volume 9 2023 (2023). issn: 2297-4687. doi: 10.3389/fams.2023.1187878. Jean-Michel Marin et al. “Approximate Bayesian computational methods”. In: Statistics and Computing 22.6 (Nov. 1, 2012), pp. 1167–1180. issn: 1573-1375. doi: 10.1007/s11222-0119288-2. Youssef Marzouk and Dongbin Xiu. “A Stochastic Collocation Approach to Bayesian Inference in Inverse Problems”. In: PRISM: NNSA Center for Prediction of Reliability, Integrity and Survivability of Microsystems 6 (Oct. 2009). doi: 10.4208/cicp.2009.v6.p826. Youssef M. Marzouk and Habib N. Najm. “Dimensionality reduction and polynomial chaos acceleration of Bayesian inference in inverse problems”. In: Journal of Computational Physics 228.6 (2009), pp. 1862–1902. issn: 0021-9991. Edward Meeds and Max Welling. “GPS-ABC: Gaussian process surrogate approximate Bayesian computation”. In: Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence. UAI’14. Quebec City, Quebec, Canada: AUAI Press, 2014, pp. 593–602. Richard Nickl. Bayesian non-linear statistical inverse problems. eng. Zurich lectures in advanced mathematics. Berlin: EMS Press, 2023. isbn: 3985470537. Richard Nickl and Benedikt Pötscher. “Bracketing Metric Entropy Rates and Empirical Central Limit Theorems for Function Classes of Besov- and Sobolev-Type”. In: Journal of Theoretical Probability 20 (June 2007), pp. 177–199. doi: 10.1007/s10959-007-0058-1. Elina Numminen et al. “Estimating the Transmission Dynamics of Streptococcus pneumoniae from Strain Prevalence Data”. In: Biometrics 69.3 (2013), pp. 748–757. (Visited on 11/19/2025). George Papamakarios. Neural Density Estimation and Likelihood-free Inference. 2019. arXiv: 1910.13233 [stat.ML]. url: https://arxiv.org/abs/1910.13233. George Papamakarios and Iain Murray. “Fast epsilon-free Inference of Simulation Models with Bayesian Conditional Density Estimation”. In: Advances in Neural Information Processing Systems. Ed. by D. Lee et al. Vol. 29. Curran Associates, Inc., 2016. George Papamakarios, Theo Pavlakou, and Iain Murray. “Masked Autoregressive Flow for Density Estimation”. In: Advances in Neural Information Processing Systems. Ed. by I. Guyon et al. Vol. 30. Curran Associates, Inc., 2017. George Papamakarios, David Sterratt, and Iain Murray. “Sequential Neural Likelihood: Fast Likelihood-free Inference with Autoregressive Flows”. In: Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics. Ed. by Kamalika Chaudhuri and Masashi Sugiyama. Vol. 89. Proceedings of Machine Learning Research. PMLR, Apr. 2019, pp. 837–848. Stefan T. Radev et al. BayesFlow: Learning complex stochastic models with invertible neural networks. 2020. arXiv: 2003.06281 [stat.ML]. url: https://arxiv.org/abs/2003.06281. Tom Rainforth et al. “On Nesting Monte Carlo Estimators”. In: Proceedings of the 35th International Conference on Machine Learning. Ed. by Jennifer Dy and Andreas Krause. Vol. 80. Proceedings of Machine Learning Research. PMLR, June 2018, pp. 4267–4276. Oliver Ratmann et al. “Using Likelihood-Free Inference to Compare Evolutionary Dynamics of the Protein Networks of H. pylori and P. falciparum”. In: PLOS Computational Biology 3.11 (Nov. 2007), pp. 1–13. doi: 10.1371/journal.pcbi.0030230.
22
[38] [39]
[40] [41]
[42] [43] [44]
[45]
[46] [47] [48]
[49]
[50] [51] [52]
[53]
REFERENCES
Jonas Rothfuss et al. “Conditional Density Estimation with Neural Networks: Best Practices and Benchmarks”. In: arXiv:1903.00954 (2019). Donald B. Rubin. “Bayesianly Justifiable and Relevant Frequency Calculations for the Applied Statistician”. In: The Annals of Statistics 12.4 (1984), pp. 1151–1172. doi: 10.1214/ aos/1176346785. Bodhisattva Sen. A Gentle Introduction to Empirical Process Theory and Applications. Lecture notes. July 2022. S. A. Sisson, Y. Fan, and Mark M. Tanaka. “Sequential Monte Carlo without likelihoods”. In: Proceedings of the National Academy of Sciences 104.6 (2007), pp. 1760–1765. doi: 10. 1073/pnas.0607208104. Björn Sprungk. “On the local Lipschitz stability of Bayesian inverse problems”. In: Inverse Problems 36.5 (Apr. 2020), p. 055015. doi: 10.1088/1361-6420/ab6f43. A. M. Stuart. “Inverse problems: A Bayesian perspective”. In: Acta Numerica 19 (2010), pp. 451–559. doi: 10.1017/S0962492910000061. Andrew M. Stuart and Aretha L. Teckentrup. “Posterior consistency for Gaussian process approximations of Bayesian posterior distributions”. In: Mathematics of Computation 87.310 (Mar. 2018). doi: 10.1090/mcom/3244. Leila Taghizadeh and Ansgar Jüngel. “Bayesian Inversion for the Identification of the Doping Profile in Unipolar Semiconductor Devices”. In: SIAM Journal on Scientific Computing 47.3 (2025), B690–B709. doi: 10.1137/24M1687042. A. van der Vaart and J.A. Wellner. Weak Convergence and Empirical Processes: With Applications to Statistics. Springer Series in Statistics. Springer, 1996. A. W. van der Vaart. Asymptotic Statistics. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 1998. Julia Walchessen, Amanda Lenzi, and Mikael Kuusela. “Neural likelihood surfaces for spatial processes with computationally intensive or intractable likelihoods”. In: Spatial Statistics 62 (2024), p. 100848. issn: 2211-6753. doi: https://doi.org/10.1016/j.spasta.2024. 100848. Hongji Wang et al. Sequential Bayesian Design for Efficient Surrogate Construction in the Inversion of Darcy Flows. 2025. arXiv: 2507.17713 [stat.ML]. url: https://arxiv.org/ abs/2507.17713. Simon N. Wood. “Statistical inference for noisy nonlinear ecological dynamic systems”. In: Nature 466.7310 (Aug. 2010), pp. 1102–1104. issn: 1476-4687. doi: 10.1038/nature09319. Dongbin Xiu. Numerical Methods for Stochastic Computations: A Spectral Method Approach. Princeton University Press, 2010. Liang Yan and Tao Zhou. “An adaptive surrogate modeling based on deep neural networks for large-scale Bayesian inverse problems”. In: Communications in Computational Physics 28.5 (Nov. 2020), pp. 2180–2205. Jice Zeng et al. “Solving high-dimensional inverse problems using amortized likelihood-free inference with noisy and incomplete data”. In: Computer Methods in Applied Mechanics and Engineering 443 (2025), p. 118064. issn: 0045-7825. doi: https://doi.org/10.1016/j. cma.2025.118064.