Conditional score-based diffusion models for Bayesian inference in infinite dimensions
Authors: Lorenzo Baldassari
Ali Siahkoohi
Josselin Garnier
Knut Solna
Maarten V. de Hoop
Venue: NeurIPS 2023 spotlight (2023)
Keywords: score-based generative models
diffusion models
inverse problems
bayesian inference
infinite dimensions
TL;DR: We propose a method to learn the posterior distribution in infinite-dimensional Bayesian linear inverse problems using amortized conditional score-based generative models.

Abstract: Since their initial introduction, score-based diffusion models (SDMs) have been successfully applied to solve a variety of linear inverse problems in finite-dimensional vector spaces due to their ability to efficiently approximate the posterior distribution. However, using SDMs for inverse problems in infinite-dimensional function spaces has only been addressed recently, primarily through methods that learn the unconditional score. While this approach is advantageous for some inverse problems, it is mostly heuristic and involves numerous computationally costly forward operator evaluations during posterior sampling. To address these limitations, we propose a theoretically grounded method for sampling from the posterior of infinite-dimensional Bayesian linear inverse problems based on amortized conditional SDMs. In particular, we prove that one of the most successful approaches for estimating the conditional score in finite dimensions—the conditional denoising estimator—can also be applied in infinite dimensions. A significant part of our analysis is dedicated to demonstrating that extending infinite-dimensional SDMs to the conditional setting requires careful consideration, as the conditional score typically blows up for small times, contrarily to the unconditional score. We conclude by presenting stylized and large-scale numerical examples that validate our approach, offer additional insights, and demonstrate that our method enables large-scale, discretization-invariant Bayesian inference.

Conditional score-based diffusion models for Bayesian inference in infinite dimensions
Lorenzo Baldassari University of Basel
Ali Siahkoohi Rice University aliskOrice.o edu
Josselin Garnier
Ecole Polytechnique, IP Paris josselin garmierropolyyeechiiuuuu....
lorenzo. baldassari Ounibas.ch
Knut Sølna University of California Irvine ksolna@math. werttthh..ii...o uci. edu
Maarten V. de Hoop Rice University
mvd2erice. edu
Abstract
Since their initial introduction, score-basee diffusion models (SDMs) have been suc- cessfully applied solve variety of linear inverse problems in finite-dimensional vector spaces due to their ability to efficiently approximate the posterior distribu- tion. However, using SDMs for inverse problems in infiiieedimmnnsionll function spaces has only been addressed recently, primarily through methods that learn the unconditional score. While this approach is advantageous for some inverse problems, it is mostly heuristic and involves numerous computationally costly for- ward operator evaluations during posterior sampling. To address these limitations, we propose theoretically grounded method for sampling from the posterior of Bayesian linear inverse problems based on amortized condi- tional SDMs. In particular, prove that one of the most successful approaches for estimating the conditional score in finite dimensions the conditional denoising estimator can also be applied in infinite dimensions. A significant part of our analysis dedicated to demonstrating that extending to the conditional setting requires careful consideration, the conditional score typically blows up for small times, contrarily to the unconditional score. We conclude by presenting stylized and large-scale numerical examples that validate our approach, offer additional insights, and demonstrate that our method enables large-scale, discretization invariant Bayesian inference.
mnfinie-dimennioe
SDMs
Introduction
Inverse problems seek to estimate unknown parameters using noisy observations or measurements. One of the main challenges is that they are often ill-posed. problem ill-posed if there are no solutions, or there are many (two more) solutions, or the solution is unstable in relation to small errors in the observations [1]. A common approach transform the original ill -ooep problem into well-posed one is to formulate it as least-squares optimization problem that minimizes the difference between observed and predicted data. However, minimization of the data misfit alone negatively impacts the quality of the obtained solution due to the presence of noise in the data and the inherent nullspace of the forward operator [2, 3). Casting the inverse problem into a Bayesian probabilistic framework allows, instead, for a full characterization of all the possible solutions [4-6]. The Bayesian approach consists of putting a prior probability distribution describing uncertainty the parameters of interest, and finding the posterior distribution over these parameters [7]. The prior must be chosen appropriately in order to mitigate the ill- ill-posedness of the problem and facilitate computation of the posterior. By adopting the Bayesian formulation, rather than finding single solution to the inverse
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

problem (e.g., the maximum a posteriori estimator [8]), distribution of solutions- the posterior- is finally obtained, whose samples are consistent with the observed data. The posterior distribution can then be sampled to extract statistical information that allows for uncertainty quantification [9].
Over the past few years, deep learning-based methods have been successfully applied to analyze linear inverse problems in Bayesian fashion. In particular, recently introduced score-based diffusion models (SDMs) [10] have become increasingly popular, due to their ability of producing approxi- mating samples from the posterior distribution [11, 12]. An SDM consists of a diffusion process, which gradually perturbs the data distribution toward tractable distribution according to a prescribed stochastic differential equation (SDE) by progressively injecting Gaussian noise, and a generative model, which entails a denoising process defined by approximating the time-reversal of the diffusion. Crucially, the denoising stage is also diffusion process [13] whose drift depends on the logarithmic gradients of the noised data densities- the scores scores-which are estimated by Song et al. [10] using neural network. Among the advantages of SDMs over other deep generative models is that they produce high-quality samples, matching the performance of generative adversarial networks [14], without suffering from training instabilities and mode-collapse [10, 15J. Additionally, SDMs are not restricted to invertible architectures like normalizing flows [16]. which often limits the complexity of the distributions that can be learned. Finally, and most importantly to the scope of this work, SDMs have demonstrated superior performance in a variety of inverse problems, such as image inpainting [10, 17], image colorization [10]. compressing sensing, and medical imaging [12, 18].
In the aforementioned cases, SDMs have been applied by assuming that the data distribution of interest is supported on finite- finie-dimennional vector space. However, in many inverse problems, especially those governed by partial differential equations (PDEs), the unknown parameters to be estimated are functions (e.g., coefficient functions, boundary and initial conditions, or source functions) that exist in a suitable function space, typically an infinit-diimessioaal Hilbert space. The inverse heat equation or the elliptic inverse source problem presented in [19] are typical examples of ill-posed inverse problems that are naturally formulated in infinie-diimensinall Hilbert spaces. In addition to these PDE-based examples, other interesting cases that are not PDE-based include geometric inverse problems (e.g., determining the Riemann metric from geodesic information or the background velocity map from travel time information in geophysics involving singular integral operators [21]. A potential solution for all of these problems could be to discretize the input and output functions into fintte-dimensional vectors and apply SDMs to sample from the posterior. However, theoretical studies of current diffusion models suggest that performance guarantees do not generalize well on increasing dimensions [22-24]. This is precisely why Andrew Stuart's guiding principle to study a Bayesian inverse problem for functions-"avoit discretization until the last possible moment" i5]-is critical to the use of SDMs.
[20]) and inverse problems
Motivated by Stuart's 's principle, this work we define a conditional score in the infinite-dimensiona setting, critical step for studying Bayesian inverse problems directly in function spaces through SDMs. In particular, we show that using this newly defined score as a reverse drift of the diffusion process yields generative stage that samples, under specified conditions, from the correct target conditional distribution. We carry out the analysis by focusing on two cases: the case of a Gaussian prior measure and the case of a general class of priors given as density with respect to a Gaussian measure. Studying the model for a Gaussian prior measure provides illuminating insight, not only because it yields an analytic formula of the score, but also because it gives full characterization of SDMs in the infiniteddimensional setting, showing under which conditions we are sampling from the correct target conditional distribution and how fast the reverse SDE converges to it. It also serves as a guide for the analysis in the case of a general class of prior measures. Finally, we conclude this work by presenting, in Section 6, stylized and large-scale numerical examples that demonstrate the applicability of our SDM. Specifically, we show that our SDM model (i) is able to approximate non-Gaussian multi-modal distributions, a challenging task that poses difficulties for many generative models [25]; (ii) is discretization iiniionnnivaaiatt,, property that is consequence of our theoretical and computational framework being built on the -iieene-immesionall formulation proposed by Stuart [5]; and (iii) is applicable to solve large-scale Bayesian inverse problems, which we demonstrate by applying to a large-scale problem in geophysics, i.e., the linearized wave-equation- based imaging via the Born approximation that involves estimating 256x256-dimension unknown parameter.
Related works Our work primarily motivated by Andrew Stuart's comprehensive mathematical theory for studying PDE-governed inverse problems in Bayesian fashion [5]. In particular, we are
2

---

interested in the infiniteddimnnsoone dimensional analysis [7, 26], which emphasizes the importance of analyzing PDE-governed inverse problems directly in function space before discretization.
Our paper builds upon a rich and ever expanding body of theoretical and applied works dedicated to SDMs. Song et al. [10] defined SDMs integrating both score-based (Hyvärinen [27]; Song and Ermon [17]) and diffusion (Sohl-Dickstein et al. [28]; Ho et al. [29]) models into a single continuous-time framework based on stochastic differential equations. The generative stage in SDMs is based on a result from Anderson [13] proving that the denoising process is also a diffusion process whose drift depends on the scores. This result holds only in vector spaces, which explains the difficulties to extend SDMs to more general function spaces. Initially, there have been attempts to project the input functions into a dimens-ammnnional feature space and then apply SDMs (Dupont et al. [30]; Phillips et al. [31]). However, these approaches are not discretization invariant. It is only very recently that SDMs have been directly studied in function spaces, specifically fininte-dimensionll Hilbert spaces. Kerrigan et al. [32] generalized diffusion models to operate directly in function spaces, but they did not consider the time-continuous limit based on SDEs (Song et al. [10]). Dutordoir et al. [33] proposed denoising diffusion generative model for performing Bayesian inference of functions. Lim et al. [34] generalized score matching for trace-class noise corruptions that live in the Hilbert space of the data. However, as Kerrigan et al. [32] and Dutordoir et al. [33], they did not investigate the connection to the forward and backward SDEs as Song et al. [10] did in finite dimensions. Three recent works, Pidstrigach et al. [24]. Franzese et al. [35] and Lim et al. [36], finally established such connection for the unconditional setring. In particular, Franzese et al. [35] used results from
SDEs theory (Föllmer and Wakolbinger [37]; Millet et al. [38])
close to Anderson [13].
Among the mentioned works, Pidstrigach et al. [24] is the closest to ours. We adopt their formalism to establish theoretical guarantees for sampling from the conditional distribution. Another crucial contribution comes from Batzolis et al. [39], we build upon their proof to show that the score can be estimated by using a denoising score matching objective conditioned on the observed data [17, 40J. A key element in Pidstrigach et al. [24], emphasized also in our analysis, is obtaining the expected square norm of the score that needs
estimate on to be uniform in time. We explicitly compute the expected square norm of the conditional score in the case of Gaussian prior measure, which shows that a uniform in time estimate is not always possible in the conditional setting. This is not surprising, given that the singularity in the conditional score noise vanishes is well-known phenomenon in finite dimensions and has been investigated in many works, both from a theoretical and a practical standpoint [41, 42]. In our paper, we provide set of concrete conditions to be satisfied to ensure uniform estimate in time for a general class of prior measures in infinite dimensions.
Pidstrigach et al. [24] have also proposed a method for performing conditional sampling, building upon the approach introduced by Song et al. [12] in a finite-dimensionl setting. Like our approach, their method can be viewed as contribution to the literature
likelihood-frie ffree, simulation- inference [43, 44). Specifically, the algorithm proposed by Pidstrigach et al. [24] relies on a projection- type approach that incorporates the observed data into the unconditional sampling process via proximal optimization step to generate intermediate samples consistent with the measuring acquisition process. This allows Pidstrigach et al. [24] to avoid defining the conditional score'. While their method has been shown to work well with specific inverse problems, such as medical imaging [12], it is primarily heuristic, and its computational efficiency varies depending on the specific inverse problem at hand. Notably, their algorithm may require numerous computationally costly forward operator evaluations during posterior sampling. Furthermore, their implementation does not fully exploit the discretizaatonniivarian since they employ UNet to parametrize their score, limiting the evaluation of their score function to the training interval. The novelty of our work is then twofold. First, we provide theoretically grounded guarantees for an approach that not heuristic and can be implemented such that it not constrained to the grid on which we trained our network. As result, we show that we effectively take advantage of the iiscreizationniivaraance invariance property achieved by adopting the dnfiiit-dimmesiionll
property achieved by studying the problem in infinite dimensions
We note that, in newer version of their paper submitted to arXiv on October 3, 2023 (four months after submission aiv and NeurIPS 2023). Pidstrigach al. [24] abandoned the projection-type approach. Instead, they invoke the conditional score function to perform posterior sampling and solve an inverse problem in a similar fashion to ours (we refer to Section 8.3 of their paper for details). However, they still do not address the well- posedposedeess of the forward-revent reverse conditional SDE and the singularity of the conditional score, and their implementation is based on UNets and, thus, is not cattreeizaaionnnvvriiiaa
3

---

formulation proposed by Stuart [5]. Second, we perform discretization-invanni Bayesian inference by learning an amortized version of the conditional score. This done by making the score function depending the observations. As result, provided that we have access to high-quality training data, during sampling we can input any new observation that we wish to condition on directly during simulation of the reverse SDE. In this sense, our method is data-driven, as the information about the forward model is implicitly encoded in the data pairs used to learn the conditional score. This addresses critical gap in the existing literature, as the other approach using siinittedimeesionall SDM resorts to projections onto the measurement subspace for sampling from the posterior-a method that not only lacks theoretical interpretation but may also yield unsatisfactory performance due to costly forward operator computations. There are well-documented instances in the literature where amortized methods can be a preferred option in Bayesian inverse problems [45-50) as they reduce inference computational costs by incurring an offline initial training cost for a deep neural network that is capable of approximating the posterior for unseen observed data, provided that one has access to a set of data pairs that adequately represent the underlying joint distribution.
Main contributions The main contribution of this work is the analysis of conditional SDMs in infiniteddimensioonal Hilbert spaces. More specifically,
We introduce the conditional score in an ciminiedi setting (Section 3).
We provide comprehensive analysis of the forward-reverse conditional SDE framework in the case of a Gaussian prior measure. We explicitly compute the expected square norm of the conditional score, which shows that uniform in time estimate is not always possible for the conditional score. We prove that as long as we start from the invariant distribution of the diffusion process, the reverse SDE converges to the target distribution exponentially fast (Section 4).
We provide a set of conditions to be satisfied to ensure a uniform in time estimate for a general class of prior measures that are given as density with respect to a Gaussian measure. Under these conditions, the conditional score-use usedd as a reverse drift of the diffusion process in SDMs-yieldd generative stage that samples from the target conditional distribution (Section 5).
We prove that the conditional score can be estimated via a conditional denoising score matching objective in infinite dimensions (Section 5).
We present examples that validate our approach, offer additional insights, and demonstrate that our method enables large-scale, inppeiizaationnnvviienn Bayesian inference (Section 6).
2 Background
Here, we review the definition of unconditional score-based diffusion models (SDMs) in infinite- dimensional Hilbert spaces proposed by Pidstrigach al. [24], as we will adopt the same formalism to define SDMs for conditional settings. We refer to Appendix A for a brief introduction to key tools of probability theory in function spaces.
Let Adaata be the target measure, supported on a separable Hilbert space (H, (.. Consider a forward anfiniteddimensione diffusion process (X:)tejo.r) for continuous time variable t e [O,T). where Xo is the starting variable and X, its perturbation at time t. The diffusion process is defined by the following SDE:
1 X.dt + VaWW. dX, X
where C: H + H is fixed trace class, positive-definite symmetric covariance operator and W. is a Wiener process on H. Here and throughout the paper, the initial conditions and the driving Wiener processes in (1) are assumed independent.
The forward SDE evolves Xo ~ Po towards the Gaussian measure (0, C') as t - oc. The goal of score-based diffusion models is to convert the SDE in (1) to a generative model by first sampling X (0, C),, and then running the correspondent reverse-time SDE. In the dimen--imessioall case, Song et al. [10] show that the reverse-time SDE requires the knowledge of the score function V loog pr (Xt), where (X,) )) the density of the marginal distribution of X, (from now on denoted P.) with respect to the Lebesgue measure. In inmeniiediimensionall Hilbert spaces, there is no natural analogue of the Lebesgue measure (for additional details, see [5I]) and the density is thus no longer well defined. However, Pidstrigach et al. [24, Lemma 1] notice that, in the mnite-iimensional setting

---

where H = =R, the score can be expressed as follows:
CV2 lo) = <""
(x-e""
=
(2)
fort >> O. Since the eighthhand side of the expression above is also well-defined in infinite dimensions,
Pidstrigach et al. [24] formally define the score as follows:
Definition 1. In the infinite ciiitt-ddmnnsionnl setting, the score or reverse drift defined by
S(,,z) == -(1 -1---(1)-- (1= (x-
-12 = =1).
Assuming that the expected square norm the score is uniformly bounded in time, Pidstrigach et al. [24, Theorem 1] shows that the following SDE
1 Z, S(T Zi) + VCdW. dz,
Zo Pr.
(4)
the time-reversal of (1) and the distribution of Zr is thus equal to Pao proving that the forward- reverse SDE framework of Song et al. [10] generalizes to the infinit-ddimensional setting. The reverse SDE requires the knowledge of this newly defined score, and one approach for estimating it is, similarly to [10], by using the denoising score matching loss [40]
E [|s(t..x)) + (1 (1- - -x) x0]]
(5)
where 5(1, X,) is typically approximated by training a neural network.
The conditional score in infinite dimensions
Analogous to the score function relative to the unconditional SDM in infinite dimensions, we now define the score corresponding to the reverse drift of an SDE when conditioned on observations. We consider setting where Xo is an H-valued random variable and H is an space. Denote by
Hilbert
Y = AXo B,
noisy observation given by n linear measurements, where the measurement acquisition process is represented by a linear operator A H - R", and B (O, CB) represents the noise, with Cg n x n nonnegative matrix. Within a Bayesian probabilistic framework, solving (6) amounts to putting an appropriately chosen prior probability distribution Ho on Xg. and sampling from the conditional distribution of Xo given Y
To the best of our knowledge, the only existing algorithm which performs conditional sampling using infinite-dimensional diffusion models on Hilbert spaces is based on the work of Song et al. [12]. The idea, adapted to infinite dimensions by Pidstrigach et al. [24], is to incorporate the observations into the unconditional sampling process of the SDM via a proximal optimization step to generate intermediate samples that consistent with the measuring acquisition process. Our method relies instead on utilizing the score of introduce in this work. We begin by defining the conditional score, by first noticing that, in finite dimensions, we have the following lemma:
SDMs conditioned on observed data, which we
Lemma 1. In the function for t > as
setting where = RR., we can express the conditional score
CV log -1- (x 1/2 X = y, X, = =,)
Since the right-hand side of (7) is well-defined in infinite dimensions, by following the same line of thought of Pidstrigach et al. [24] we formally define the score as follows:
Definition 2. In the infinite- dimensi-imen setting, the conditional score is defined by
S(1, = -(1 -11---->-
(XolY = y, v, v.X.= =i))
(8)
Remark 1. It is possible to define the conditional score in infiniteddimnsiione Hilbert spaces by resorting to the results of [37, 38], see Appendix C.I.
5

---

For Definition 2 to make sense, we need to show that if we use (8) as the drift of the time-reversal of the SDE in (1) conditioned on y, then will sample the correct conditional distribution of Xo given Y = y in infinite dimensions. In the next sections, we will carry out the analysis by focusing on two cases: the case of a Gaussian prior measure (0, Cp), and the case where the prior of Xo is given as a density with respect to Gaussian measure, i.e.,
dpo du
Xo po.
(zo)
1 = (O,
(9)
E. [X,]]
where Cj is positive and trace class and @ is bounded with E, < to0.
4 Forward-reverse conditional SDE framework for Gaussian prior measure
We begin our analysis of the forward-reverse conditional SDE framework by examining the case where the prior of Xo is Gaussian measure. This case provides illuminating insight, not only because it is possible to get an analytic formula of the score, but also because offers a full characterization of SDMs in the adfnite-dimensionall setting, showing under which conditions we are sampling from the correct target conditional distribution and how fast the reverse SDE converges to it. We also show that the conditional score can have a singular behavior small times when the observations are noiseless, in contrast with the unconditional score under similar hypotheses.
We assume that + 0 in (9). All distributions in play Gaussian:
Xo ~ NO,,
(10)
X,
(11)
XolY ~ N(MXY,
(12)
X ~
(13)
where C, + CB) and C. = Cp. theorem [52], there exist (n₂) in [0, +00) and an orthonormal basis (1z) in H such that CRuy = fyy Vj. We consider the thinite-dimensione case with μj > 0 Vj. We assume that C₁ is trace class so that Σ, M4 < +00. We assume that the functions Uy are eigenfunctions of C and we denote by Xÿ the corresponding eigenvalues.
+ By Mercer
We assume an observational model corresponding to observing finite-dimensional subspace of H spanned by (y(1)....., (y((n) corresponding to gh Vy((k), I, ,....,,,, n. where gj E H is such that (Af), (gj.f); f). We denote I(n) {n(1)...... .t.).....ntt m())). We assume moreover CB o In. Let Z
be the solution of reverse-ting -imese SDE:
dz. = Z,dt + S(T I, Ze, + VCaW,, Zo ~ XT y,-
(14)
We want show that the reverse SDE have just formulated in (14) indeed constitutes reversal of the stochastic dynamics from the forward SDE in (I) conditioned on y. To this aim, we will need the following lemma:
Lemma 2. We define z(i) = (oy, Z). pb5) >g//s for all j. We also define g(5) My for
j E I(") and yاا) = 0 otherwise, and q(3) = p3//0z la/ozz/oo for E I(m) and qld) = 0 otherwise. Then we can write for all j
az -
dt 1 -1n dt + Vydwo,
(15)
with WG) independent and identically distributed standard Brownian motions,
e'pl5) + g())) - I+ (e (1 +())² n (1)))()) = 2
e 2/12)0g0)) (1055) (1) 1 (e' - 1))00) + gb))
(16)
Proof. The proof a Gaussian calculation. relies on computing (S, vz), which yields an analytic formula. See Appendix B.
O
Lemma 2 enables us discuss when we are sampling from the correct target conditional distribution XutY ~ NM,Y, C.). We can make a few remarks:

---

In the limit T + oo, we get (( (T -) + 1122 and y( (T - t) +
. If j & I(n) then we have the same mode dynamics as in the unconditional case. Thus we
sample from the correct target distribution if I' is large if we start from 2/6) NN(0 x()) for
x( = + 21 - ---), which is the distribution of x = =Xr, U2) given Y = y.
j E I(m) and we start from zo²1 ~ N(26)) >
), then we find z/-1 ~ N xX)) with
2/' -100)
ex
My
(E (1-1+ 1)p (1 + g(o)
(1 (1+ - 1)p10 -1--
g(
g()))s
+ 1+
(17)
=10) T//2
00)9(0)
(1- 1+ (e 1)p()(1 (1 + -110--11)0111 g(3)) (1-TTO0-11000
=0)
T 1+ (e - (1 q011-
20
+ 1 g(0)
(18)
The distribution of (5) (Xo, Uz) given Y y (1)(19001 (00//( ((1 + 90)), 2 /1 + g()). As
3/ + y0) /(1 + (1)) and Enj + 11/(11 + 900) as T + +o0, this shows that we sample from the exact target distribution (the one of X, given Y y) for T large.
If we start the reverse-time SDE from the correct model
--122)10)(1)
IS uye I
1+ 1(96) + 21
(19)
g()
then indeed 2I"² N the same distribution as Xo given Y = y, which is the exact target distribution. We can show similarly that -r- has the same distribution as X2 given Y y for any E [0.T].
(11401)))))))000)) ) g)))). This shows that, for any T, Zr has
In the case that dB that we observe the mode values perfectly for j I(m), then
e²/2
(1))
(1) 2 - e'
(lw)) (1)
(20)
e
and indeed lim,TT - gوو) a.s. Indeed the t-1 singularity at the origin drives the process to the origin like in the Brownian bridge.
Our analysis shows that, as long as we start from the invariant distribution of the diffusion process, we are able to sample from the correct target conditional distribution and that happens exponentially fast. This proves that the score of Definition 2 is the reverse drift of the SDE in (14). Additionally, the analysis shows that the score is uniformly bounded, except when there is no noise in the observations, blowing up near ( = 0. Remark 2. Note that, for g0 = -, we obtain the unconditional model:
dzz")= = - (T dt + ViydW(1), with (t) 1
E e'pli)
(21)
1+ (e' 1))(()
=C the square expectation of the norm and the Lipschitz constant of the score are uniformly bounded in time: sup (n.0 = 1/2.
Proposition 1. The score is S,,, = (t, (2, (W)e, oe (t,,)
=
(1(T))(T - t) 1/2) xli) + plos)(T (T - t)y(0) and it satisfies
e' g()))
A²
El| S(I, X l = y) - ET+ (e' - () + g()) Pj
(22)
Proof. The proof is a Gaussian calculation given in Appendix B.
L
x3 which is equal to E, Ag In the unconditional setting, we have Et||S(t, X))31 = 1+(e--1)pth My when C C Cp.- It indeed uniformly bounded in time. In the conditional and noiseless setting (aB = 0), we have Et][S(t, Xx. w) d
y] = EjEzo IpL My Ejezo T which blows up as 1/t as t - 0. This result shows that the extension of the based diffusion models to the conditional setting is not trivial.
7

---

5 Well-posedness for the reverse SDE for a general class of prior measures
We are now ready to consider the case of general class of prior measures given density with respect to a Gaussian measure. The analysis of this case resembles the one of Pidstrigach et al. [24] for the unconditional setting. The main challenge is the singularity of the score for small times, an event that in the Gaussian case was observed in the noiseless setting. In this section we will provide a set of conditions be satisfied by & in (9), so that the conditional score is bounded uniformly in time. The existence of this bound is needed to make sense of the forward-reverse conditional SDE, and to prove the accuracy and stability of the conditional sampling.
We start the analysis by recalling that, in the infinit-ddimensional case, the conditional score is (8). It easy to get first estimate:
St Xt. y) 1 Y = yl s (1 - (1-c-))-T(()),
(23)
The proof follows from Jensen inequality and the law of total expectation, see Appendix C. Note that (23) is indeed an upper bound of (22) since Tr(C)- Σ3 A3-
Note that the bound (23) is also valid for the unconditional score S(t, = -(1 - -
Xx = =i). We can observe that the upper bound (23) blows up in the limit of small times. (X e-//2
We can make a few comments:
The bound (23) is convenient for positive times, but the use of Jensen's 's inequality results in very crude bound for small times. As shown in the previous section, we know that there exists a bound (21) for the unconditional score in the Gaussian case that is uniform in time.
The singular behavior as 1/t at small time is, however, not artificial. Such a behavior is needed in order to drive the state to the deterministic initial condition when there are exact observations. This behavior has been exhibited by (20) and (22) in the Gaussian case when og = 0. This indicates that the following assumption (24) is not trivial in the conditional setting.
For Definition make sense in the more general case where the prior of Xo is given as density with respect to a Gaussian measure, we will need to make the following assumption.
Assumption 1. For any y E R", have
sup E[]iS(t, Xt, V Y = v] < oc. tE(O.T]
(24)
We are now ready to state the analogous result to Pidstrigach et al. [24, Theorem 1]. Proposition 2. Under Assumption 1, the solution of the reverse-time SDE
dZ 1 Zedt + S(T 1, Zr, dt + VCdW,, Zo XT|Y = y:
(25)
satisfies 7t Y y-
Proof. Given Assumption 1, the proof follows the same steps as the one given in [24] for the unconditional score. See Appendix C for the full proof.
Assumption I is satisfied under some appropriate conditions. In the following proposition, we provide set of conditions that ensure the satisfaction of this assumption. It shows that it is possible to get an upper bound in (23) that is uniform in time provided some additional conditions are fulfilled.
Proposition 3. We assume that Cy in (9) and C in (1) have the same basis of eigenfunctions (vy) and we define x(") = (Xe, uz) and St) x. y) = (S(t, x. z,y), 0z) so that in (1) S(t, x. y) = S0 y)u,- We assume an observational model as described in Section 4 and that the (1 q))) are uniformly bounded with respect j and that C is of trace class. We make modified version of assumption in (9) as follows. We assume that 1) the conditional distribution of Xo given Y- y is absolutely continuous with respect to the Gaussian measure μ with Radon- Nikodym derivative proportional to Hap(-WWroo v)); 2) we have ( v) = E, (x) .y).
( = =(00, vz); for (260), ,v) (ipp-v())(s00 (o(bi,,)) v)) we have
1
{1000000)) n)1 K, (x(0), w)) (00 .v) S Lio0i0 200),, (26)
K
where K and L do not depend on j. Then Assumption I holds true.
8

---

Proof. The proof is given in Appendix C.
To use the new score function of Definition 2 for sampling from the posterior, we need to define a way to estimate it. In other words, we need to define a loss function over which the difference between the true score and neural network se st 2t, y) is minimized in 0. A natural choice for the loss function is
(t-0110.))...3--CC11 (t, v2, so 24, 2
(27)
however it cannot be minimized directly since we do not have access to the ground truth conditional score {t,, Ft, y). Therefore, in practice, different objective has to be used. Batzolis al. [39] proved that, in finite dimensions, denoising score matching loss can be used:
E (1
[ [I[V In p(pp
.., v)|l"
(28)
This expression involves only V.. log p(ra) (2)) which can be computed analytically from the transi- tion kernel of the forward diffusion process, also in infinite dimensions. In the following proposition, we build on the arguments of Batzolis et al. [39] and provide proof that the conditional denoising estimator a consistent estimator of the conditional score in infinite dimensions.
Proposition 4. Under Assumption I, the minimizer in 0 of
E (, (xxx= [I] -1--- e e1/2x0) so Tz, v)
(29)
the same as the minimizer of
Ex, [|SS(t, (IISSt ot, v)
(30)
The same result holds if we add U(,,T) T) in the expectations.
Proof. The proof combines some of the arguments of Batzolis et al. [39] and steps of the proof of Lemma 2 in [24], see Appendix C.
Remark 3. A statement of robustness can be written as in [24, Theorem 2].
6 Numerical experiments
To put the presented theoretical results into practice, we provide two examples. The first stylized example aims at showcasing (i) the ability of our method in capturing nontrivial conditional distribu- tions; and (ii) the discretization-invariann example, we sample from the posterior distribution of linearized seismic imaging problem in order to demonstrate the applicability of our method to large-scale problems. both examples, order to enable learning in function spaces, parameterize the conditional score using Fourier neural operators [53]. Details regarding our experiment and implementation? are presented at Appendix D.
property of the learned conditional SDM. In the second
Stylized example Inspired by Phillips et al. [31], we define the target density via the relation ao = ay² + E with E ~ [(1, 2), a ~ 1], and y E [-3,3]. Figure la illustrates the samples Ty evaluated on a fine y grid. After training (details in Appendix D). sample the conditional distribution on uniformly sampled grids between [-3, 3,, each having 20 to 40 grid points. Figures 1b and Ic show the predicted samples for grid sizes of 25 and 35, respectively. The marginal conditionals associated with in the bottom row of Figure 1 indicates the ground truth density, and colored estimated densities correspond to different discretizations of the horizontal axis. The visual inspection of samples and estimated densities indicates that our approach is indeed discretizationninvaiaa
- 1.0, 0.0, 0.5 are shown in Figures Id-1f, respectively. The gray shaded density
Linearized seismic imaging example In this experiment, we address the problem of estimating the short- wavelegth component of the Earth's 's subsurface squared-slownes: model (i.e., seismic image; cf. Figure 2a) given surface measurements and a long- waveeengg,,,, smooth squared-slowness model (cf. Figure 2b). Following Orozco al. [54], in order to reduce the high dimensionality of surface measurements, we apply the adjoint of the forward operator, the Born scattering operator, to the measurements and use the outcome (cf. Figure 2c) instead of measured data condition the SDM. After training, given previously unseen observed data, we use the SDM to sample 103 posterior
Code reproduce results can be found at hutp://gihhb..co//iiiiakkooiei/

---

True fuactions
Predicted functions
Prodacted functions
grdd eton 25
grit sitw
(a)
(b)
(c)
1.00
y 0.50
0.6
denely
Roscction
True deratty
demilt
04
0.0 -1
(d)
(e)
(1)
Figure 1: A depiction of the method's discretization invariance. Top row displays ground truth (a) and predicted samples (functions) on uniformly sampled grid with (b) 25 and (c) 35 grid points. Bottom row shows conditional distribution marginals for (d) y =- -1.0, (e) y 0.0, and (f) y 0.5.
0.0
Ground truth tunknown/ image
L55
0.0
Smooth bockground model
k1pm
Reverse migrated
0.5
0.5
13
3.0
0.0
1.0
2.0
2.0
2.5
1.0
2.5
0.5
3.0
3.0
3.0
0.
Hostennia dalunce Suza)
Hostenntall disuace Oaz)
Horizootal dolance Dant)
(a)
(b)
0.0
Conditional ostimate
0.0
Pointwiisn standard devintion
x18
9.0
0.0
Prodiction arror
0.5
1.0
0.5
1.0
0.5
-005
2.0
2.0
0.6
1.0 2.5
-15
(d)
Figure 2: Seismic imaging and uncertainty quantification. (a) Ground-truth seismic image. (b) Background squared-slownen (c) Data after applying the adjoint Born operator. (d) Conditional (posterior) mean. (e) Pointwise standard deviation. (f) Absolute error between Figures 2a and 2d.
samples estimate the conditional mean (cf. Figure 2d), which corresponds to the mininum--variane estimate [55], and the pointwise standard deviation (cf. Figure 2c), which we use to quantify the uncertainty. As expected, the pointwise standard deviation highlights areas of high uncertainty. particularly in regions with complex geological structures- such as near intricate reflectors and areas with limited illumination (deep and close to boundaries). We also observe strong correlation between the pointwise standard deviation and the error in the conditional mean estimate (Figure 2f), confirming the accuracy of our Bayesian inference method.
Conclusions
We introduced a poooretically-goounde method that is able to perform conditional sampling in infinite- dimensional Hilbert (function) spaces using score-based diffusion models. This is a foundational step in using diffusion models to perform Bayesian inference. To achieve this, we learned the dimente-dimeniionall score function, as defined by Pidstrigach et al. [24], conditioned on the observed data. Under mild assumptions on the prior, this newly defined score used as the reverse drift of the diffusion process-yiele generative model that samples from the posterior of linear inverse problem. In particular, the well-known singularity in the conditional score for small times can be avoided. Building these results, we presented stylized and large-scale examples that showcase the validity of our method and its pnccretizationinnaaianne theoretical and computational framework being built on difnitt-dimensionall spaces.
a property that is consequence of our
10

---

Acknowledgments
JG was supported by Agence de l'Innovation de Défense AID via Centre Interdisciplinaire d' E'tuddes pour la Défense et la Sécurité CIEDS (project 2021 PRODIPO). LB, AS, and MVdH acknowledge support from the Simons Foundation under the MATH + X program, the Department of Energy under grant DE-SC0020345, and the corporate members of the Geo-Mathematical Imaging Group at Rice University. was supported by Air Force Office of Scientific Research under grant FA9550-22-1-005 22-1-0176 and the National Science Foundation under grant DMS-2308333
References
[1] Hadamard. Lectures the Cauchy's Problem in Linear Partial Differential Equations. Yale University Press, 1923.
[2] R. C. Aster, B. Borchers, and C. H. Thurber. Parameter estimation and inverse problems. Elsevier, 2018.
[3] K. Ito and B. Jin. Inverse Problems. World Scientific, 2014.
[4] M. Lehtinen, L. Paivarinta, and E. Somersalo. Linear inverse problems for generalised random variables. Inverse Problems, 5(4):599, 1989.
[5] A. M. Stuart. Inverse problems: Bayesian perspective. Acta Numerica, 19:451-35 115,, 2010.
[6] A. Tarantola. Inverse problem theory and methods for model parameter estimation. SIAM,
2005.
[7] B. T. Knapik, A. W. van der Vaart, and J. H. van Zanten. Bayesian inverse problems with Gaussian priors. The Annals of Statistics, 39(5):2626 2657, 2011.
[8] A. Bennett. Inverse Modeling of the Ocean and Atmosphere Cambridge University Press, 2002.
[9] A. Stuart. Uncertainty quantification in Bayesian inversion. 2014 SIAM Conference on Uncertainty Quantification, 2014.
[10] Y. Song, Sohl-Dickstein, Kingma, Kumar, S. Ermon, and B. Poole. Score-based generative modeling through stochastic differential equations. In International Conference on Learning Representationss 2021.
[11] B. Kawar, G. Vaksman, and M. Elad. SNIPS:
Solving noisy inverse problems stochastically. In Advances Neural Information Processing Systems, volume 34, pages 21757-21769, 2021.
[12] Y. Song, L. Shen, L. Xing, and S. Ermon. Solving inverse problems in medical imaging with score-based generative models. In International Conference on Learning Representations, 2022.
[13] B. D. Anderson. Reverse-time diffusion equation models. Stochastic Processes and their Applications, 12(3)311-326, 1982.
[14] Dhariwal and A. Nichol. Diffusion models beat GANs on image synthesis. In Advances in Neural Information Processing Systems, volume 34, pages 8780-8794, 2021.
[15] D. Gnaneshwar, B. Ramsundar, D. Gandhi, R. Kurchin, and V. Viswanathan. Score-based generative models for molecule generation, 2022. URL https://arxiv.org/abb.. 04698. Preprint.
[16] H. Lee, J. Lu, and Y. Tan. Convergence of score-based generative modeling for general data distributions. In Proceedings of The 34th International Conference on Algorithmic Learning Theory, volume 201, pages 946-985, 2023.
[17] Y. Song and S. Ermon. Generative modeling by estimating gradients the data distribution. In Advances in Neural Information Processing Systems, volume 32, 2019.

---

[18] A. Jalal, M. Arvinte, G. Daras, E. Price, A. G. Dimakis, and J. Tamir. Robust compressed sensing MRI with deep generative priors, Advances Neural Information Processing Systems, 34:14938-14999 14954, 2021.
[19] M. Dashti and A. M. Stuart. The Bayesian approach to inverse problems, pages 311 311-428. Springer International Publishing, 2017.
[20] G. Uhlmann and A. Vasy. The inverse problem for the local geodesic ray transform. Inventiones Mathematicae, 205(():83-120, 2016.
[21] A. Dynin. Inversion problem for singular integral operators: C* -"-approach. Proceedings the National Academy of Sciences, 75(10):4668 46700, 1978.
[22] S. Chen, S. Chewi, Li, Y. Yi, Salim, and Zhang. Sampling is as easy as learning the score: Theory for diffusion models with minimal data assumptions. In The Eleventh International Conference on Learning Representations, 2023.
[23] V. D. Bortoli. Convergence of denoising diffusion models under the manifold hypothesis. Transactions on Machine Learning Research, 2022.
[24] Pidstrigach, Y. Marzouk, S. Reich, and S. Wang. Infinite-dimensionall diffusion models for function spaces, 2023. URL https: :tap:::///araae //arxiv...vv/////
http.ttts./aai...or///aa/aa www//or//////////oo/o 2 10130v1. Preprint.
[25] S. Arora, A. Risteski, and Y. Zhang. Do GANs learn the distribution? Some theory and empirics. In International Conference on Learning Representations, 2018.
[26] S. Lasanen. Measurements and infinite-dimeensionall statistical inverse theory. Proceedings in Applied Mathematics and Mechanics, 7: 7:1080101-1080102 2007.
[27] A. Hyvärinen. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(24))695-709, 2005.
[28] J. Sohl-Dickstein, E. Weiss, N. Maheswaranathan, and S. Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In Proceedings of the 32nd International Conference on Machine Learning, volume 37, pages 2256-2265, 2015.
[29] Ho, A. Jain, and P. Abbeel. Denoising diffusion probabilistic models. In Advances in Neural Information Processing Systems, volume 33, pages 2020.
[30] E. Dupont, H. Kim, S. M. A. Eslami, D. J. Rezende, and D. Rosenbaum. From data to functa: Your data point is function and you can treat it like one. In Proceedings of the 39th International Conference on Machine Learning, volume 162, pages 3694-5725 2022.
[31] A. Phillips, Seror. M. Hutchinson, V. Bortoli, A. Doucet, and E. Mathieu. Spectral diffusion processes. In NeurIPS 2022 Workshop on Score-Based Methods, 2022.
[32] G. Kerrigan, J. Ley, and P. Smyth. Diffusion generative models in infinite dimensions. In Proceedings of The 26th International Conference Artificial Intelligence and Statistics, volume 206, pages 3538-9563 2023.
[33] V. Dutordoir, A. Saul, Z. Ghahramani, and F. Simpson. Neural diffusion processes. In Proceedings of the 40th International Conference on Machine Learning, volume 202, pages 0090-9012, 2023.
[34] J. H. Lim, N. B. Kovachki, R. Baptista, C. Beckham, K. Azizzadenesheli, J. Kossaifi, V. Voleti, J. Song, K. Kreis, J. Kautz, C. Pal, A. Vahdat, and A. Anandkumar. Score-based diffusion models in function space, 2023. URLhttps::/ httpt:t/t///aaee www.aawwwww...//////0/0000 20000/000000000..... Preprint.
[35] G. Franzese, G. Corallo, S. Rossi, M. Heinonen, M. Filippone, and P. Michiardi. Continuous- time functional diffusion processes. In Thirty-seventh Conference on Neural Information Processing Systems, 2023. URL
12

---

[36] S. Lim, E. Yoon, T. Byun, T. Kang, S. Kim, K. Lee, and S. Choi. Score-based generative modeling through stochastic evolution equations hilbert spaces, In Thirty-seventee Conference on Neural Information Processing Systems, 2023. URL. htw.:: httaap:///ppereeeee.eeeete id=GrElRvXnEj.
[37] H. Föllmer and A. Wakolbinger. Time reversal of infiniteddimensiinall diffusions. Stuchastic Processes and their Applications, 22(1):59-77 1986.
[38] A. Millet, D. Nualart, and M. Sanz. Time reversal for infiniteddimension diffusions. Probability Theory and Related Fields, 82(3):315 337,, 1989.
[39] G. Batzolis, J. Stanczuk, C.-B. Schönlieb, and C. Etmann. Conditional image generation with score-based diffusion models, 2021. URL
http://aaxxivv.r/ab..r.///// 116.//111.1.166. Preprint.
[40] Vincent. A connection between score matching and denoising autoencoders. Neural Compu- tation, 23(7):1661-1674 2011.
[41] D. Kim, S. Shin, K. Song, W. Kang, and I.-C. Moon. Soft truncation: A universal training technique of score-based diffusion model for high precision score estimation. In Proceedings of the 39th International Conference on Machine Learning. volume 162, pages 1 1201-11228, 2022.
[42] Dockhorn, A. Vahdat, and K. Kreis. Score-based generative modeling with critically-aamie Langevin diffusion. In International Conference on Learning Representations, 2022.
[43] K. Cranmer, J. Brehmer, and G. Louppe. The frontier of simulation- inference. Proceedings of the National Academy of Sciences, 1eees,,117444 117744):30055-30062 2020.
[44] A. Lavin, D. Krakauer, H. Zenil, Gottschlich, T. Mattson, J. Brehmer, A. Anandkumar, Choudry, K. Rocki, A. G. Baydin, C. Prunkl, B. Paige, O. Isayev, E. Peterson, P. L. McMahon, J. Macke, K. Cranmer, J. Zhang, H. Wainwright, A. Hanuka, M. Veloso, S. Assefa, S. Zheng, and A. Pfeffer. Simulation intelligence: Towards new generation of scientific methods, 2022. Preprint.
URL https:: /ttp://aaivvveeea org/ abaaaass//a/ 2212/222222.2.....
[45] R. Baptista, B. Hosseini, N. B. Kovachki, and Y. Marzouk. Conditional sampling with monotone GANs: From generative models likelihood inference, 2023. URL https /arxiv.. org/abs/ rrg/abs/2006 06755. Preprint.
[46] Y. Kim, S. Wiseman, A. Miller, D. Sontag, and A. Rush. Semi-amortized variational autoen- coders. In Proceedings of the 35th International Conference on Machine Learning, volume 80, pages (678-2687,, 2018.
[47] J. Kruse, G. Detommaso, U. Köthe, and R. Scheichl. HINT: Hierarchical invertible neural transport for density estimation and Bayesian inference. Proceedings of the AAAI Conference on Artificial Intelligence, (319))81118111 2021.
[48] S. T. Radev, U. K. Mertens, A. Voss, L. Ardizzone, and U. Köthe. BayesFlow: Learning complex stochastic models with invertible neural networks. IEEE Transactions on Neural Networks and Learning Systems, 33(4):1452-146) 2022.
[49] A. Siahkoohi and F. J. Hermann. Learning by example: Fast reliability-awar seismic imaging with normalizing flows. In First International Meeting for Applied Geoscience & Energy Expanded Abstracts, pages 1580-1585, 2021.
[50] A. Siahkoohi, G. Rizzuti, Orozco, and F. Herrmann. Reliable amortized variational inference with physics-based latent distribution correction. Geophysics, 8873):R297-R322 2023.
[51] G. Prato. An Introduction Infinite-Dimmnsiional Analysis. Springer, 2006.
[52] J. Mercer and A. R. Forsyth. XVL. Functions of positive and negative type, and their connection the theory of integral equations. Philosophical Transactions of the Royal Society of London. Series A, Containing Pupers of a Mathematical or Physical Character, 209(441-458): 446,, 1909.
13

---

[53] Z. N. B. Kovachki, K. Azizzadenesheli, B. liu, K. Bhattacharya, A. Stuart, and A. Anand- kumar. Fourier neural operator for parametric partial differential equations. In International Conference Learning Representations, 2021.
[54] R. Orozco, A. Siahkoohi, G. Rizzuti, T. van Leeuwen, and F. J. Herrmann. Adjoint operators enable fast and amortized machine learning based Bayesian uncertainty quantification. In Medical Imaging 2023: Image Processing, volume 12464, page 12464IL, 2023.
[55] B. D. Anderson and J. B. Moore. Optimal Filtering. Prentice-Hall, Englewood Cliffs, NJ, 1979.
[56] D. Kingma and J. Ba. Adam: method for stochastic optimization, 2017. URL https: //ar //axxiv.. brrg/abs/10 1412. 6980. Preprint.
[57] G. Lambaré, J. Virieux, R. Madariaga. and S. Jin. Iterative asymptotic inversion in the acoustic approximation. Geophysics, 57(9): 1338-1554 1992.
[58] Nemeth, C. Wu, and Schuster. Least-squares migration of incomplete reflection data. Geophysics, 64(1):208-22) 221, 1999.
[59] Veritas. Parihaka 3D marine seismic survey Acquisition and processing report. Technical Report 3460, New Zealand Petroleum & Minerals, Wellington, 2005.
[60] WesternGeco, Parihaka 3D PSTM final processing report. Technical Report 4582, New Zealand Petroleum & Minerals, Wellington, 2012.
[61] M. Louboutin, M. Lange, F. Luporini, N. Kukreja, P. A. Witte, F. J. Herrmann, P. Velesko, and G. J. Gorman. Devito (V3.1.0): An embedded domain-specific language for finite differences and geophysical exploration. Geoscientific Model Development, 12(3):1665-188 2019.
[62] Luporini, M. Louboutin, M. Lange, Kukreja, P. Witte, J. Hückelheim, C. Yount, P. Kelly, F. J. Herrmann, and G. J. Gorman. Architecture and performance of Devito, system for automated stencil computation. ACM Transactions on Mathematical Software, 46(1), 2020.
14

---

A Probability measures on infnite-dimensional Hilbert spaces
In this section, we briefly present some fundamental notions related to probability measures on infinit-dimensionall spaces, specifically separable Hilbert spaces (H, (--)). There is abundant literature on the subject. For more details we refer to Stuart [5], Pidstrigach et al. [24], Kerrigan et al. [32], Prato [51] and references therein.
A.1
Gaussian measures on Hilbert spaces
Definition 3. Let (S, IP) be probability space. A measurable function X Ω H is called a Gaussian random element (GRE) if for any h E H, the random variable (h, X) has a scalar Gaussian distribution.
Every GRE has a mean element m E H defined by
X (())d(ww)
Jo
and linear covariance operator C:: H: - H defined by
Ch J (h, (())
(m, h)m, vh E H.
We denote X ~ (mm, C²) for a GRE in H with mean element m and covariance operator C. It can be shown that the covariance operator of a GRE is trace class, positive-defintt and symmetric. Conversely, for any trace class, dositive-definite and symmetric linear operator C H - H and every m E H, there exists GRE with (m,, C). This leads us to the following definition:
Definition 4. If X is GRE, the pushforward of through X, denoted by Ixx, is called Gaussian probability measure on H. We will write Px (m,, C)).
Let (m, C). We can make few remarks:
1) For any E H, we have (h, X) ~ ((hh, m), (Ch, h))-
2) C' is compact. By Mercer theorem [52] there exists (X₂) and an orthonormal basis of eigenfunctions
(v)) such that A 2 and Cvy- = Vj. We consider the infnite-dimensional case in which Aラ >
Vj.
3) Suppose m = 0 (we call the Gaussian measure of X centered). The expected square norm X is
given by
x X -
Σ uy x
(Cez UC) = ~ Tr(C), Σ 1=1
j=1
which is finite since C is trace class.
A.2 Absolutely continuous measures and the Feldman- -ajek theorem
Here we introduce the notion of absolute continuity for measures.
Definition Let μ and be two probability measures on H equipped with its Borel o-algebra B(H)). Measure μ is absolutely continuous with respect (we write μ. < v) if μ(Σ) 0 for all Y B(H) such that v(Σ)
Definition 6. If & and < μ then μ and v are said to be equivalent and we write μ v. If μ and v are concentrated on disjoint sets then they are called singular; in this case we write μ 1
Another notion that will be used throughout the paper the Radon Radon-Nikodye Nikodym derivative.
Theorem 1. Let a and be two measures ((H) and be o-finite. If μ. < then there exists a v-measurable function f on H such that
u(A') J fdu, VA' E B(H).
Furthermore, f is unique v-a.e. and is called the Radon-Nikodym derivative of i with respect to v. is denoted by du/dv.
15

---

Remark 4. In the paper, we will sometimes refer to f as the density of with respect to v.
We finally able state the Feldman-Hajek theorem in its general form. Theorem 2. The following statements hold.
1. Gaussian measures μ (m,, C1), (m,, C2) are either singular or equivalent.
2. They are equivalent if and only if the following conditions hold:
(i) u and y have the same Cameron-Martin space Ho - -
(ii) mi -m2 E Ho-
(iti) The operator (CT C2 C2/2 102
Hilbert-Schhiic operator on the closure
Ho.
3. If n and v are equivalent and Ci = C2 = C, then v-a.s. the Radon-Nikodym derivative du/do given by
du (h) (h)) dy
where (h) = (C-1/2(m - m2). C-1/2(h - (m-m)) _/11c-1/2 - mm E H.
A.3 Bayes' theorem inverse problems
Let H and K be separable Hilbert spaces, equipped with the Borel o-algebra, and A : H K measurable mapping. We want to solve the inverse problem finding X from Y, where
Y A(X)= + B
and B E K denotes the noise. We adopt a Bayesian approach to this problem. le (X,Y) E HxK be a random variable and compute X X|Y. We first specify (X, Y) as follows:
1) Prior: X ~ fo measure on H.
2) Noise: ~ ?/o measure on K, with B independent from X.
The random variable Y|X is then distributed according to the measure Tae the translate of noo by A(X). We assume that The < 1/0. Thus for some potential v H K + R.
dn (y) V(a.))
dro
The potential V. y) satisfying the above formula is often termed the negative log likelihood of the problem. Now define vo be the product measure V/ P4 Vi- We can finally state the following infimiteddimensional analogue of the Bayes' theorem.
Theorem 3. Assume that H K R is vo measurable and define
7(y)= /
((.ww) dpo-
dpo-
Then
dpol- ( = (x) Z(y)
1 (w..))
dji
where Pol- = y) is the conditional distribution of X given Y y.
B Proofs of Section 4
B.1 Proofs Lemma and Proposition
We assume that Cp in (9) and C in (1) have the same basis of cigenfunctions (oy) and that Caul μ Crj = Xyn Vj. We define x(i) (X,, uz). = (v. vz) and S61 x, y) = (S(1, x. y), 0y) so that in (1) Stt, = E, (t) x, v)u.. We assume j I(n), so that we consider a mode corresponding to an observation. We then have
dx(0)
1 xx dt Vx
2
16

---

with WC) standard Brownian motions which are independent for the different modes j. Note also that with CN and C having the same basis of eigenfunctions the system modes is diagonalized so that the x processes are independent with respect to mode j. both for the observed and un-observed modes. Thus we have
x ( = v yu) = x(6) toBm x = e (>)(1-)-))h")
for !/ (1) independent standard Gaussian random variables. We then seek
2 (.v) -[ | =, y y =[x x = ,(),,(()) y(t) = ylli),
which in this Gaussian setting is the L projection of x onto x(") and y(i). Thus we can write To (3-z) =xx + by(1) with (a, b) solving
Bl + bY x ( B((ax,) + +b xx²)xxx 0,
which gives
2/2
1 (E -1 1 + (el 1)pt0( + g(ā))*
a 1 (e - 1) (1 + g()))
for pbt) >3/11500 g(5) = /
We then get in view of (8)
-TT (et - 1)p( + q()) xb1)+ + T+ - (T+ 1) + QOT)) 10) (31)
e1 9(5)
e122(610())
(t) x, - (-
T+
Note that with some abuse of notation we then have
(t,,,), St),t 1(1))
which is important since then also the time reversed system diagonalizes. We remark that for an unobserved mode we get by a similar, but easier, calculation
xt, = (1 e'p(i) 1 (e' - 1)p()) x901.
which simply corresponds to setting dB = co in (31).
Consider next E[S (t, x, v) Note first that
[[x(()) = y) (g))--122) Var[x(() Y = v] vVar(X(() (i) = y] + >;(1 -e
(1+ y(u)
(e' +g)))
1 gio)
We then easily check that the score is conditionally centered B(s (t, Xx, Y) Y y] 0 and we then get
c'pl)) (())
Y y)
Var(x. Y y)
I+ (et- 1)p00> q(o))
c11,(0(0))31++ ( 9())
1+ (ef-1)p00) + q(i))'
which gives Proposition upon summing over the mode index j. where we define g() = 0 for the unobserved modes.
17

---

C Proofs of Section 5
C.1 Discussion about an alternative approach
The following lemma is a complementary result related to Remark 1. It shows that we can actually derive the expression of the score from the results contained in Millet et al. [38]. The result is powerful. but requires the verification of technical conditions.
Lemma 3. Under the conditions stated Proposition the score is defined by
U(,,U,t) = -- (x- -er-e-1/3
X Xx =x, Y
(32)
and the time reversed diffusion takes the form in (14).
Proof.
Define
x( =U1; X) for Cuty AzUz-
Then
dx = 1 x x dt + xxdWW() for W = (oj, W) 2
(33)
and where we assume that C' is of trace class. This is then an infinite dimensional system of the type considered in Millet et al. [38]. We proceed to verify some conditions stated in Millet et al. [38]: (i) the coefficients of the system (33) satisfy standard growth and Lipschitz continuity conditions (assumption (H1, H4) satisfied); (ii) the coefficients depend on finitely many coordinates (assumption (H2) satisfied); the system time independent and diagonal (assumption (H5) satisfied). Moreover
define z05) (II,.... ...,...;... f-- 1, ...+,,...,, then the law of x given x-()) has for 1 >> 0 density (r(2""1X(()) = 800) Y = y) with respect to Lebesgue measure and so that for to > 0 and each j: J B|00,00 log(p, (og(p,(>))))
y) co. Then it follows from Theorems 3.1 and 4.3 in Millet et al. [38] that the time reversed problem is associated with the well-posed martingale problem defined by the coefficients in (14) for the score being:
a A ai
(p((x()))( = xGi),, Y = v)
(uy, 8(x,,,)))
Pe(z(0) X( $6),, Y y)
with the convention that the right hand side is null on the set (2)()))) = ز6ز,, Y = y) = -].
It then follows for t > 0
(Ug, S(),
X a Az (pr(z(0) x( = x00, x(x)) = 2 Y =
f duo(aa)) (1) 1x(1) $0), Y = y)
pe(x00))(() = +00, Y = 4)
We then get
y. f)) = duo(lo = _(0),x Y =y)
(1)
= x00) x( = x =
Y = =n)
(
x
1- c-
(x( 1x(0 +00,, y)
1 - e (5)
dpo(=0)) ( [r+ = y.yYv =)
x0
6-
18

---

C.2 A preliminary lemma
The following lemma is the equivalent of [24, Lemma 3]. is used in the forthcoming proof of Proposition 2.
Lemma 4. In the finite dinenniame setting E RD, we have any 0 S s ssts
V log Pt,y (t. = e(t-,)) V s = =r),
where By is the expectation with respect the distribution of Xo and W given Y y and pt,y is the pdf of Xt under this distribution.
Proof. We can write
(x, = /
J Pe,y (2. (xela)d
where P (- (atx)) is the pdf of X, given Y and X, = xs. is, in fact, equal to the pdf of X given X, = Zg. which is the pdf of the multivariate Gaussian distribution with mean exp(-(tt - s)/2)zz and covariance (1 eexp(-(t s)))C. Therefore
(x)) = S (esu(E))tts(eell)
We can then deduce that
1 (pt.y(x1)= D/2) e- (t- det(C))1/2 / dxspey((s)
JRO
x Vz, exp
2(1 (1-s)) (ot exp( - 8)/2)ar.) - -xpp
- s)/2)xr.)
c(c-5)/2
(2T) - e (t- ((--))) det(C))1/2 20 drspar((1s)
1
X Vz, exp
2(1 - -e-(t-3)) (+t exp((((-- 8)/2).))))
C (xx expl-(( ())-((())2)))))) s)/2)ars)
(2T) (t- s1/1/2 det(()/1 / dx. (Pa,((t))
c(1-0)/2
RD
X exp
2(1 e-((-s)) (I) - expl- s)/2),)) (², - expl-(t- -8) =))2)(z.))
[e(1-0)12 / J dx. v., (xx) (o
which gives
(pt = (1-5)) (rjo( (x.,) V log logpss((n dd...
( get Pay(z.)
Using again that Pelawv (.12r=) Pt)o(110r, and Pala (1)
(pt,((2)) (1-0)//2
RD (r Ds. Vog Po (xx )x.
= (1-s\\2
((..)).0 I log
JDD
(pn) and (xal(t)))))) = P(annolentt))) Pewl
Since V log
we get that
v logg (og = (1-0)
Pale.y x V log Pa Da (I )xc.
RD
= (110)) V log (X.)) = xz).
19

---

C.3 Proof of Proposition 2
The proof adapts the one of [24] to the conditional setting. The only difference is that the expectation is Ey, which affects the distribution of Xo but not the one of W. Moreover, Lemma 4 shows that the key to the proof (the reverse- time martingale property of the fiiteeddimensional score) is still valid. Here By is the expectation with respect to the distribution of Xo and W given Y = y.
To prove Proposition 2, are left to show that the solution of the reverse-time SDE
dZp
2 Lt + S(T t, Zz, y) + VCdW,, Zo ~ XT =
(34)
satisfies ZT Xo y. We recall that X, is the solution to the SDE
dX, 2 X,dl + VdW,, Xo ~ P4-
We first notice that X, is given by the following stochastic convolution:
('" (-(t=s)/
=
For pD the orthogonal projection on the subspace H spanned by U....., (the eigenfunctions of C), XP- = (X,) are solutions to
1 x, dt + (CC)dw,
dxp =
2
where
CD = PP'" W. W',
We define X. = XM xp. Then
[ e -00001/2 (CD:M JaWD
x D x..M = e 1/2
where the superscript D M indicates the projection onto span (UD+1..... UMM. holds that
By sup
|| x PXPMII [i] S 2e EEylEy 1 2(1 c
Σ
4
imD41
D -> 00, where used Doob's 1.2 inequality bound the stochastic integral. Therefore (XA" is a Cauchy sequence and converges to X, in (.) Consequently, the distribution of X.N given y y converges to the distribution of X, given Y y N +oc.
Recall that
S(t, -(1 -1- E, e-1/2 Xo Xe).
and recall that
C Vogo (1 -1t----))-ppEEp e 1/2 |xP]. In particular, due to the tower property of the conditional expectations, CDV log pl (xP) EylSStt, v) x.).
Since, by Assumption 1,
[IS v), < 00; the quantities [X, -1/2 XPP bounded in (P(PP,) and will converge to the limit, [St, Xt, y) Xi] = Xe, y), by the Martingale convergence theorem. We get rid of the
projection pD by
(1 c B, [I| log D (xP) - S(t, Xt, 1
EE][| Eyy [X, - e 1/2 Xa XPI Ey,(X, - -t/2 xo X.]]2A]
1/2
< EyllEElXx c- x/ xxI E,IX, - -1/220 x
+ [(I PD PD)TE, [X, - 1/2 t/2 x.P111+]
< [E -
- e- 1/2 Xo XP] - E,,,,,
[Xe - - -112 Xo X.]1
+ [[][[[PP[[ -e -12 xx)
20

---

The first term vanishes due to our previous discussion. The second tern vanishes since
/ (t- VCdW,,1
Ey - (X, e Xp) = PD)
= [I - pP) / e (1-s))/2 VCdW CdWWalit)
(1 (1-c-)
Σ 23 +0
D+1
as - Cc.
D We now make use of the fact that V log Pa, y rection by Lemma 4. We therefore get sequence of continuous Li-bounded Martingales converging to stochastic process. Since the space of continuous Lf-bounded martingale is closed and pointwise convergence translates to uniform convergence, get that S is Le-bbunded martingale, with the convergence of of C log ply to S being uniform in time.
is a square-integrable Martingale in the reverse-time di-
We have that
2, - 40
Zods Zods- 2 0
1
-
S'' cD (2)) = V(CD)W.
=
Since all the terms the left-hand side converge in L?, uniformly in t, so does the right-hand side. Using again the closedness of the spaces Martingales and Levy's characterization of Wiener process, we find that V(CP) converges to VCW. Therefore
Z, Zo + 2 Z + S[ S(t, Ze,9) VCW.
Therefore, Ze is indeed solution to (34) and Zz ~ Xo Y = y. Using uniqueness of the solution we then conclude that this holds for any solution Z,.
C.4 Proof of (23)
(t, x,, = y| -1- -E x = X.1li =y
= (1 - E
--1-)//2 = X. Xe]1 = yl
(1 "1 EE|
(1-9)/2 VCaWall2 = x, = vl
(1 2EIl
(1- 1/ CCW..Il/ y]
(1 ") TE|
(1 )) TT(C).
C.5 Proof of Proposition
Note that with the assumptions in Proposition with C₁ and C having the same basis of eigen- functions and the separability assumption the Radon-Nikodym derivative for the modes, the system for the modes again diagonalizes. However, in this case the (conditional) distribution for x is non-Gaussian in general and the change of measure with respect to the Gaussian measure characterized by w) (z))), y). We let the superscript g denote the Gaussian case with e = 1, then
21

---

we have:
5ti (t, x.
--(1-e-1) -000)00000---0)0000 - --11/2(Xx0 X, =x, =x
(09))
x= vl) =
1 == ( x: (15))(()))) (x00 y)
((:))
x(6) -02/2 / xo
=- -1-16- --
dro
1/ (xbi))
(1)
1 (1.9)
(I0
x21) 2x01) (0))) + x01) cl/2,1))
= (1 -
e-22
6
doo
R
m
(1)
(x(1))
= ا (t, v. (t, 26),,,0),)) v) + R (t, y),
for ((10 the mode score in the Gaussian case given in (31) and with
11 (19)) 1
(z.(6) ((60))(G))()11)) (z:10) .y)
(t,
(J) 1
(2li1)
11
,,)
-1- 2/2
(16)00) (0o + xl5) x(0) (2)10 (zo + xl0e - - (x10) .1) dro- 10.0//1y
Tu
Mat by (x(ا))
We have for Φx the centered Gaussian density second moment
JR (xp0 we ve- t/=)) (ps (v xy ( )d
Tt, =
My Su (p0 we 1/2))
(v x3
(1) (ob (v. (v.y)/ /)/
yidu
for 1 - = pg/(1 (()) and () = -(و) / (1 + 9io) for gاو = (v, 0z).
Here ,1))) Py are respectively the mean and variance of x given y and where we used the parameterization set forth in Section B.1. We then have
(TO) v)1 K2,
lim T6) (t, 1. t40
We moreover have
JR 2 (o - xoe 1/2)0 (1)) (2o - (1) )dro
v)1 SE 5
- -
JR ( -E- 112)0 (/) (u zÿ 216) )oh(1))))))) (v, y)du
Sx 260 (000)p1(000) -age =00e )0² 11y
(ro ydoo
S e 112 L
Se (x21((100) we vee-t/2j0 (6 (v - (1) dd
We then find
limsup 18 / (1) 410
(t, -y)l S ( LK,
and, moreover
RR00 (t, 4)15
(1)) - x222
(e'-1) (1+ (e'- 1)pli)(1 (e'--i))) q()))}²
2/
1 1) q()) + Aj
(1()))2
Consider in the Gaussian case in (22) mode so that (1 + (0)) t oo and X, fixed, then the contribution of this mode to the score norm blows up in the small time limit. The situation with
22

---

pl) + q(o)) t oo would happen for instance in limit of perfect mode observation so that JB + 0 and thus g t oc. Indeed in the limit of small (conditional) target mode variabilty relative to the diffusion noise parameter the score drift becomes large for small time drive the mode the conditional target distribution. We here thus assume p + g is uniformly bounded with respect to mode (j index), moreover, that C' is of trace class. We then find that Assumption is satisfied with the following bound
L Xe, Ve. =
sup xE(0,7]
rye"" (ktp5)(1 + g(0) + 22ge (LK)² (1+3 (10) + g0)
< 2
We remark that in the case that we do not have uniform bound the p + y it follows from (23) that the rate of divergence of the expected square norm of the score at most t- as t 0 with C of trace class.
C.6 Proof of Proposition 4
We start from (30):
E (1, xr, - s0 vr, v) Ex, [1|5(1, m) (1
+ ((x.,)
w) v.x.. (S((,, oe, v), 80 (1, vt. v))
From Definition 1 we have
E (a (S(t, Fr, v). solt, x1,,,,,,))) 21,
=-(1 -1-- (1--et) t) Ea,, ((x, ( 22 e 2/2 1
/ro)) so so 2e, 2)) Xx x =x. ( De e 1/2 2/2200 so xe, y)
=-(1 --- Ea B. ((x
--(1 - t 'E (o ( (I.. -/ so so(1,2e111)
We obtain that
Ex [1S((t, - so(t, v)
=B+ Bl [II (1 - 1001-e-t/2.p x122 s0, 2))))
with
B Em.y []] Ir, wi]l - B(xoz)) - --- -
- t/2 zo) 1
that does not depend on O. Since L(X,/Xo To,Y = y) L(XilXo = xo) we finally get that
(11(), [1,, sa(1, vt, 1)
B + Re
[I] -11-6-1--2 -1/22 so Ft. v) 11
D Numerical experiments: details and additional results
In this section, provide additional details regarding numerical experiments. In both experi- ments, we parameterize the conditional score sp Tr, y) using discretizationninvaaian Fourier neural operators [FNOs; 53]. This parameterization enables mapping input triplets (t, zt, y) to the score conditioned on y at time t. Once trained-by minimizing the objective function in equation (28) with respect d-we the FNO approximation to the conditional score sample realizations of the conditional distribution by simulating the reverse-timg SDE in equation (14).
D.1 Stylized example
In this example, the conditional distribution that we approximate is defined using the relation
23

---

Predicted functions
Predicted functions
Prodacted functions
grsl Jon 20
grdd eton 30
grid stzw 40
-2
(a)
(b)
(c)
Figure The proposed method's discretization invariance. Predicted samples (functions) on uniformly sampled grid with (a) 20, (b) 30, and (c) 40 grid points.
Tu ay² E.,
(35)
where S ~ 2) and a a u[-1,1) Here, I refers to the Gamma distribution, and 1 denotes the uniform distribution over the set {-1, Having explicit expression characterizing the conditional distribution allows us to easily evaluate the obtained result through our method-as opposed to needing use baseline method, e.g., Markov chain Monte Carlo.
Training data The discretization invariance of our model enables us to use training data that live on varying discretization grids. We exploit this property and simulate training joint samples (zo, y) by evaluating the expression in equation (35) over values of y that are selected as nonuniform grids over the domain |-3 3] that contain 15-50 grid points.
Architecture In this example, the FNO comprises (i) a fully connected lifting layer that maps the direc-dimensional inclur-inellding the timestep the y value, and the corresponding æ value-for each grid point to 128-dimensional lifted space; (ii) five Fourier neural layers, as introduced in [53], which contain pointwise linear filters applied the five lower Fourier modes; and (iii) two fulle-connected layers, separated by ReLU activation function, that map the 128-dimensional lifted space back the conditional score scalar) for each grid point.
Optimization details To train the FNO, we minimized the objective function in equation (28) using 2 x 10 training steps. At each step, we simulated a batch of 512 training pairs selected from grid with varying numbers of discretization points. We utilized the Adam stochastic optimization method [56] with an initial learning rate of 10-3, which decayed to 5 10-4 during optimization, following a power-law rate of - -/33. Regarding the diffusion process, we followed the approach outlined in Ho et al. [29] and employed standard Gaussian noise with linearly increasing variance for the forward dynamics, which was discretized into 500 timesteps. The initial Gaussian noise variance was set to 10-4 4 and linearly increased over the timesteps until it reached x 10-2. The training hyperparameters were chosen by monitoring the validation loss over 1024 samples. The training process took approximately minutes on a Tesla V100 GPU device. For further details, please refer to our open-source implementation on GitHub.
illustrates the predicted samples associated with the remaining testing grid sizes, whose densities were shown in Figures Id-1f.
Additional results Figure 3
D.2 Linearized seismic imaging example
The inverse problem we are addressing in this example involves the task of estimating the short- wavelength component of the Earth's 's unknown subsurface squared-slowness model using measure- ments taken the surface. This particular problem, often referred as seismic imaging, can be recast a linear inverse problem when we linearize the nonlinear relationship between surface measurements and the qquared-slowness model, as governed by the wave equation. its simplest acoustic form, the linearization with respect to the slowness model-around a background smooth squared slowness model results in a linear inverse problem for the estimation of the true seismic image 8mn* using the following forward model,
d; Jmmo q: jm)mm= + en ~ p(e),
(36)
where d= = [d represents collection of n, linearized surface measurements, i.e., data wherein the zero-order term of Taylor's series has been subtracted, and J J(mo, q)) denotes the lincarized
24

---

Born scattering operator, which is defined in terms of the source signature qi and the background squared -lownessss model mg. Due to noise and linearization errors, the above expression contains the term 6,. In addition to this noise term, the forward operator has non-trivial nullspace due to the presence of shadow zones and finite-aperture data [57,58]. To tackle this challenges, we set seismic imaging into Bayesian framework and learn the associated posterior distribution via our proposed method.
Training data We generated synthetic data by applying the Born scattering operator to 4750 2D seismic images, each with dimensions of 3075 in x 5120 m. These images were extracted from the Kirchhoff migrated Parihaka-3D dataset, which contains seismic images obtained by imaging the data collected in New Zealand [59, 60J. We parameterize the Born scattering operator using a background squared-slwwness slowness model (recall Figure 2b). The data acquisition geometry involves 102 sources with a spacing of 50 m, each recorded for two seconds via 204 receivers spaced at 25 m located on top the image. The source wavelet used was a Ricker wavelet with central frequency of 30 Hz. To replicate a more realistic imaging scenario, we add hand-limited noise to the the data, obtained by filtering white noise with the source function. To create training pairs, we first simulated noisy seismic data for all the 2D seismic images based on the aforementioned acquisition design. Subsequently, we reduce the dimensionality of seismic data by applying the adjoint of the Born scattering operator to the data. We use Devito [61, 62] for the wave-equation based simulations.
Architecture In this example, the FNO is composed of (i) fully connected lifting layer that maps the five-dimensional vector-including the timestep t, the two spatial coordinates, the data value (y. obtained after dimensionality reduction by applying the adjoint of the forward operator), and the corresponding x value-for each grid point to a 32-dimensional lifted space; (ii) four Fourier neural layers, as introduced in [53]. which contain pointwise linear filters applied to the 24 lower Fourier modes; and (iii) two cully-connected layers, separated by a ReLU activation function, that map the 128-dimensional lifted space back to the conditional score for each grid point.
Optimization details We train the FNO according to the objective function in equation (28) with the Adam [56] stochastic optimization method with batch size 128 for 300 epochs. We use an initial stepsize of 2 x 10m3 decaying to 5 x 10 during optimization with a power-law rate of -1/3. We use similar diffusion process as the previous example. The training hyperparameters are chosen by monitoring the validation loss over 530 samples. The training process takes approximately two hours and 15 minutes on a Tesla V100 GPU device. For further details refer to our open-source implementation on GitHub.
Additional results Figures 4-7 illustrate more results regarding using the SDM to sample from the posterior distribution of multiple seismic imaging problem instances where the ground truth images (see Figures 4a-7a) are obtained from the test dataset. In each plot, we use the dimensionality-reduee data as the conditioning input to the FNO (see Figures 4b-7b). Through the SDM, we obtain 103 posterior distribution samples and use them to estimate the conditional mean (see Figures 4c-7c) and the pointwise standard deviation among samples (see Figures 4d-7d), with the former serving as a measure of uncertainty. In all cases, the regions of significant uncertainty correspond well with challenging-to-ineee image sections of the model, which qualitatively confirms the accuracy of our Bayesian inference method. This observation becomes more apparent in Figures 4g-7g, displaying two vertical profiles with 99% confidence intervals (depicted orange-colored shading) for each experiment, which demonstrate the expected trend of increased uncertainty with depth. Furthermore, we notice that the ground truth (indicated by dashed black lines) mostly falls within the confidence intervals for most areas. We also consistently observe a strong correlation between the pointwise standard deviation and the error in the conditional mean estimate (see Figures 4e-7e), which further asserts the accuracy of our method in characterizing the posterior distribution in this large-scale Bayesian inference problem. To prevent bias from strong amplitudes in the estimated image, we present the normalized pointwise standard deviation divided by the envelope of the conditional mean in Figures 4f-7f. These visualizations provides an amplitude-ing indeedddeddeeeeeett assessment of uncertainty. highlighting regions of high uncertainty at the onset and offset of reflectors (both shallow and deeper sections). Additionally, the normalized pointwise standard deviation underscores uncertainty in areas of the image where there are discontinuities in the reflectors (indicated by black arrows), potentially indicating the presence of faults.
25

---

Image
Reverse time migrated Image
Profle 0.2 km
Proftle at 4.9 km
0.5
Conditionel
Ground-truch Image
05
%:9 confideace interval.
wwwm  10 /w www. MA
0.0
-15
Hnviorrtil destnnee
Elortanbbl ducwers
(a)
(b)
tean
0.0
Polstwise standard devistion
xte RD
1.0
100
0.5
65
RD
Degith
1.0
-15
Hocionatal datamoo
0.0
Prediction ernor
vipe
0.5
1.
1.0
1.00
1.5
06
0.6
0.30
1000
-1000
-1000
Perurbation
Parturbation
Figure 4: Seismic imaging and uncertainty quantification. (a) Ground-truth seismic image. (b) Data after applying the adjoint Born operator (known as the reverse-time migrated image). (c) Conditional (posterior) mean. (d) Pointwise standard deviation. (e) Absolute error between Figures 4a and 4c. (f) Normalized pointwise standard deviation by the envelope of the conditional mean. (g) Vertical profiles of the ground-trute image, conditional mean estimate, and the 99% confidence interval at two lateral positions in the image.
26

---

Image
Reverse time migrated Image
Profle 0.2 km
Proftle at 4.9 km
0.5
Conditionel meen
Ground-truch
05
%39 confideace Interval
0.0
-15
Elortanbal ducwers
(a)
(b)
tean
0.0
Polstwise standard deviation
xte
9D
1.0
1.0
0.5
65
0.
Degit
1.0
-15
2.0
0.0
Prediction error
kipet
0.5
1.
1.00
1.5
06
0.6
030
1000
-1000
-1000
Perurbation
Parturhation
Figure 5: Seismic imaging and uncertainty quantification. (a) Ground-truth seismic image. (b) Data after applying the adjoint Born operator (known as the reverse-time migrated image). (c) Conditional (posterior) mean. (d) Pointwise standard deviation. (e) Absolute error between Figures 5a and 5c. (f) Normalized pointwise standard deviation by the envelope of the conditional mean. (g) Vertical profiles of the ground-trute image, conditional mean estimate, and the 99% confidence interval at two lateral positions in the image.
27

---

Image
Reverse time migrated Image
Profle at 0.2 km
Proftle at 4.9 km
0.5
Conditional
Ground-truch Image %39 confideace Interval
0.
0.0
Hnviorrtal destnnee
Elortanbbl ducwers
(a)
(b)
tean
0.0
Polstwise standard deviation
xte
9D
1.0
1.0
0.5
65
0.
Deptts
0.5
1.0
5.2
-1
2.0
datanee
Hocionatal
0.0
Prediction ernor
vipet
0.5
1.0
1.00
1.5
06
0.6
0.30
1000
-1000
-1000
Perurbation
Porturbation
Figure 6: Seismic imaging and uncertainty quantification. (a) Ground-truth seismic image. (b) Data after applying the adjoint Born operator (known as the reverse-time migrated image). (c) Conditional (posterior) mean. (d) Pointwise standard deviation. (e) Absolute error between Figures 6a and 6c. (f) Normalized pointwise standard deviation by the envelope of the conditional mean. (g) Vertical profiles of the ground-trute image, conditional mean estimate, and the 99% confidence interval at two lateral positions in the image.
28

---

Image
Reverse time migrated Image
Profie 0.2 km
Proftle at 4.9 km
0.5
Conditionel
Ground-truch Image
05
%:9 confideace Interval
0.0
1o00 M Mmm Wi
-15
Invioattal destnnee
Flortanbal ducwere
(a)
(b)
tean
0.0
Poistwise si standard devistion
xte
9D
1.0
1.0
0.5
85
0.D
Degtt
55
1.0
-1
2.0
Hocionatal datano
www
0.0
Prediction ernor
1107
0.5
1.0
1.00
06
0.6
0.30
1000
-1000
1000
-1000
Perurbation
Parturhation
Figure 7: Seismic imaging and uncertainty quantification. (a) Ground-truth seismic image. (b) Data after applying the adjoint Born operator (known as the reverse-time migrated image). (c) Conditional (posterior) mean. (d) Pointwise standard deviation. (e) Absolute error between Figures 7a and 7c. (f) Normalized pointwise standard deviation by the envelope of the conditional mean. (g) Vertical profiles of the ground-trute image, conditional mean estimate, and the 99% confidence interval at two lateral positions the image.
29


Reviews:
[{"review_id": "njfJ0YF6XH", "reviewer": "Reviewer_pjkQ", "summary": "The author extends score based diffusion from finite dimensional processes to separable Hilbert space processes.\nThey demonstrate that on a non-linear toy data set that the method can work.", "questions": "-When you in  section 4. let $(Af)_i=(v_k,f)$ where $v_k$ is a eigenvector, is it not very obvious what you get? I mean the problem is a infinite series of independent processes and the data don't couple the processes. Does it not mean that you are back to the finite dimensional case (since if $v_j \\notin A$ then for that $j$ the processes is equivalent to the prior?\n\n-In the application setting I don't get the prior on $X_0$ and when I look in the code it looks one has used $x_0$ as just a line?\n- What is the operator $C$ in this example?\n\nvery very minor:\nDo you really have to write infinite dimensional in each sentence?", "limitations": ".", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper reads very well and is easy to follow. It is important to study what happens in general separable Hilbert spaces as many algorithms and models breaks down in the infinite dimensional setting. And even though you always do real application in the finite dimensional setting it is still important because if the method works in a separable Hilbert space setting it will not break down when you increase the precision your finite dimensional discretization.\n\n", "weaknesses": "The numerical experiment is not very convincing. The results for this toy problem are not very impressive. \nI suspect that if one runs multiple parallel version of a  Crank–Nicolson algorithm one would get better result.\nIt also would be interesting to see how the numerical method scale in practice when you increase the resolution of the grid.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the Reviewer for their critical feedback.\n\nWe are happy to clarify our manuscript in response to the Reviewer's questions. \n\n1. **Numerical experiment not very convincing.** We acknowledge your remark, along with those of the other reviewers, regarding the need for improvement in our toy example. In our global response, we have included a PDF that contains additional experiments that will be incorporated in the paper. We have addressed your specific concern by incorporating results indicating the discretization invariance of our method by sampling the toy conditional distribution over grids with varying discretization. These preliminary results anticipate that we are able to sample the bi-modal non-Gaussian conditional distribution (Figure 1 in the attached PDF). Furthermore, we have included a large-dimensional example that involves learning an amortized approximation to the posterior distribution of a wave-equation-based inverse problem over a $256 \\times 256$ dimensional unknown, illustrating the applicability of our method to large-scale problems (Figure 2 in the PDF). \n\n   We sincerely hope that these updates could lead to an improvement in your assessment of the paper.\n\n2. **When in Section 4 you let $(Af)_i = (v_k,f)$ where $v_k$ is a eigenvector, is it not very obvious what you get?** The Gaussian setting of Section 4 makes it possible to carry out explicit and detailed calculations, because it is indeed possible to work $j$ by $j$. The main purpose of Section 4 is to show that the extension of the score-based diffusion models to the conditional setting is not trivial, but possible, in the infinite-dimensional setting. As a byproduct, it also shows that this extension in the finite-dimensional setting also requires some conditions because the blow-up of the score that we exhibit also holds in the finite-dimensional setting with noiseless data. So even the finite-dimensional setting is not obvious~! The singularity in the conditional score as noise vanishes is indeed a well-known phenomenon in the finite-dimensional setting. It has important implications. For a discussion about the difficulties this singularity causes in practice and numerical methods for approximating the score at small times in score-based models in the finite dimensional case, we refer to [1, 2]. As for our paper, we anticipate that Assumption 1 is reasonably easy to satisfy. However, we must be mindful of the possibility of encountering a blow-up of the score under the conditions described by Reviewer hg8t. To address this concern, we will add a remark in our paper.\n\n3. **$C$ and prior on $x_0$ in the toy example.** Incorporating the suggestions we have received, we have made the decision to extend and improve Section 6. Consequently, we will be updating our implementation. The revised Section 6 will showcase an extended numerical experiment, demonstrating not only the success of our method in approximating a bi-modal non-Gaussian conditional distribution but also its discretization-invariance in practice as the grid resolution is varied (Figure 1 in the PDF included in the global response). Moreover, we will include a new experiment in the appendix, demonstrating the applicability of our method to a large-scale inverse problem in geophysics, specifically linearized seismic imaging via the Born approximation (Figure 2 in the PDF), that involves estimating a $256 \\times 256$-dimensional unknown parameter.\n\n   To address your other questions, in the toy example currently provided in the paper, we do not employ any explicit prior information, and in our amortized variational inference approach the prior is implicitly learned during training using the provided dataset of joint samples. We acknowledge that the notation used in Section 6 might have been misleading, as our focus in the current example in the paper is on learning the conditional distribution of $y$ given $x_0$, while in the rest of the paper it was the opposite (we fixed the notation in the PDF included in the global response). As for the operator $C$, we are utilizing the finite-dimensional approximation outlined in [3, Appendix G]. Specifically, we employ [4, Equation 11] as the projected equation.\n\n3. **Too many \"infinite-dimensional\" in the paper.** Thank you for the feedback. We will work on reducing the redundancy in the writing.\n\n**References**\n\n[1] D. Kim, S. Shin, K. Song, W. Kang, and I.-C. Moon, *Soft Truncation: A Universal Training Technique of Score-based Diffusion Model for High Precision Score Estimation*, 2022.\n\n[2] T. Dockhorn, A. Vahdat, and K. Kreis, *Score-Based Generative Modeling with Critically-Damped Langevin Diffusion*, 2022.\n\n[3] J. Pidstrigach, Y. Marzouk, S. Reich, and S. Wang, *Infinite-dimensional diffusion models for function spaces*, arXiv:2302.10130.\n\n[4] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole, *Score-based generative modeling through stochastic differential equations*, ICLR 2021."}, {"review_id": "PluZ681Bdk", "reviewer": "Reviewer_hg8t", "summary": "Score-based diffusion models are successful in solving inverse problems in a finite-dimensional setting, but infinite-dimensional diffusion models needs to be constructed with care, as the definitions for Lebesgue measures and densities become less clear. The authors extends the work of Pidstrigach et al. [33] of unconditional score matching in infinite dimensions to a conditional setting. ", "questions": "I have a few questions regarding the general aspects of Bayesian inference in infinite dimensions, as I am not familiar with the exact formulation. \n- Observational model in the paper occupies a finite-dimensional subspace. Is there a scenario where one cannot find an orthonormal basis such that the observation $y$ only spans a finite subspace?\n- While I understand Assumption 1 _can_ be satisfied under certain conditions given by the prior measure's Radon-Nikodym derivatives, but does this assumption suffer from finite training data? If we think about extreme settings with only a few training data, the score matching essentially tries to memorize these noiseless data points, causing an automatic violation of the assumption. ", "limitations": "- Line 15 typo: it should read \"extension of ... to the *conditional* setting\". \n- Minor citation error on Line 72: The seminal score matching paper [17] has Hyvärinen as the sole author. \n- Line 293 typo: “… proposition _on_ such set of conditions”", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "I recommend a strong accept for the paper because of its theoretical soundness and its approachable presentation in its explanation. \n\nThe paper presents a theoretically elegant and principled approach to solving inverse problems in infinite dimensions, as it is guided by Stuart’s principle and does not involve projection to finite vector spaces and discretization when unnecessary. The paper also analyzes a general scenario with prior distributions absolutely continuous w.r.t. the Gaussian measure, and propose analogous results to Pidstrigach et al. [33]. ", "weaknesses": "I cannot identify a specific point of weakness that has to be addressed. One can argue against its simple experiment but I see a proof-of-concept experiment sufficient for this paper. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the Reviewer for their positive feedback. \n\nWe are happy to clarify our manuscript in response to the Reviewer's questions. \n\n1. **Finite-dimensional observational model.** Indeed, if the number of observations is finite, we may think that the observations \nonly span a finite subspace. There are, however, instances where considering infinite-dimensional measurements can prove advantageous. This is particularly relevant when aiming to demonstrate the robustness of theoretical results, such as the asymptotic behavior of sample errors, regardless of the method used to discretize measurements. Such situations often arise when dealing with data being functions observed on a dense array. Andrew Stuart provides an example in [1, Section 3.5]: the inverse problem is determining the initial condition for the heat equation, given noisy observation of the solution at a positive time.\n\n2. **About Assumption 1.** This is a good remark. Indeed we may deal with a blow up of the score under such conditions and the results  in Section 4 make this explicit in the case with a Gaussian prior.  However, this situation should only occur with noiseless data. As soon as data are noisy, we believe Assumption 1 should be reasonably easy to satisfy.\n\n3. **Typos.** Thank you for noticing the typos. We will fix them.\n\n**References**\n\n[1] A. Stuart, *Inverse problems: a Bayesian perspective*, Acta Numerica 19 (2010), 451-559."}, {"review_id": "Ta858UjdMk", "reviewer": "Reviewer_WxYS", "summary": "This paper proposed a method to deal with inverse problems in infinite dimensions using conditional-score-based models. Specifically, they propose to directly learn the posterior distribution in infinite-dimensional Bayesian linear inverse problems using amortized conditional SDMs. Moreover, this paper also discussed the robustness of the learned distribution against perturbations of the observations.  A numerical experiment is conducted to validate the efficiency.", "questions": "Some additional questions:\n\n1. What is the difference between (2) and (3)? \n\n2. The experiment part only considers a simple low-dimensional problem, is it possible to add some results of large-scale problems?", "limitations": "See above.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. This paper proposed an interesting method to deal with infinite-dimensional Bayesian linear inverse problems.\n\n2. It provides a detailed analysis of the forward-reverse conditional SDE framework in the case of a Gaussian prior measure.\n\n3. It provides a set of conditions to ensure a uniform in-time estimate for a general class of prior measures.", "weaknesses": "1. Regarding the introduced definition of the conditional score and the result that the conditional score can be estimated via a conditional denoising score matching, it seems that they are straightforward extensions of the unconditional case. For score-based models, the key is to learn a general distribution using score matching, regardless of whether it is a conditional distribution or a non-conditional distribution. I mean, there is actually no fundamental difference between a conditional distribution and non-conditional distribution, i.e., score-based matching applies to them or other general distributions equally.\n\n2. From my understanding, and also as suggested in the abstract, this paper only focuses on infinite-dimensional Bayesian linear inverse problems, rather than arbitrary inverse problems. The title of this paper is kind of misleading.  On the other hand, for linear inverse problems, the score of the likelihood term is relatively easy to obtain, compared to directly training a conditional score network. Please correct me if I am wrong.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the Reviewer for their feedback.\n\nWe are happy to clarify our manuscript in response to the Reviewer's questions. We hope that this could lead to an increase in their score.\n\n1. **Conditional SDMs.** Various approaches have been proposed in the literature for dealing with conditioning, both in finite and infinite-dimensional cases. In [1], an effective approach involves subspace projection during time-reversed diffusion. The subsequent noise-level-dependent update yields favorable application outcomes. In [2] three methods are discussed from theoretical and practical viewpoints in finite dimensions. It is indeed not obvious what is the best way to enforce the constraint associated with observations and what the role of the noise level is, and how to incorporate the constraints from a computational perspective. These challenges are particularly acute in the infinite-dimensional case (see next answers). Here we remark that we agree that a Bayesian and a conditional score perspective is a natural approach. It is important however to understand how the conditioning affects the score and the training and the role of noise. Here we have shed light on this important challenge by considering the Gaussian context discussed in Section 4 where we get explicit insight about the conditional score. This reveals in particular the singular behavior that one may observe for the time reversed diffusion in case of low noise observations. Moreover, we generalized the results in [2] regarding a conditional denoising estimator in infinite dimensions. This result is important for efficient training and provides a theoretical underpinning for efficient conditional generation and further research. Indeed in the context of inverse problems efficient implementation of the conditioning is the central challenge. \n\n2. **Title.** Regarding the potentially misleading title, we don't know if it can be modified at later stages. While we agree that we could have been more vocal in the title about the specific problems we are addressing, we remark that the vast majority of papers concerning SDMs and inverse problems (IPs) focus on linear IPs. Addressing nonlinear IPs is a challenging task, and the limited subset of papers that delve into more general cases typically highlight this in their titles from the outset. However, we will make a concerted effort to emphasize in the abstract and throughout the text that we are exclusively considering linear IPs.\n\n3. **Score of the likelihood.** Contrarily to methods that incorporate the gradient of the log-likelihood in order to sample from the posterior, ours does not assume the knowledge of the forward model, as only data pairs are used to learn the score. We remark that, while using the log-likelihood in SDMs linear IPs looks straightforward, it becomes in fact analytically intractable in terms of DMs, due to their dependence on time. That's why existing approaches resort to projections onto the measurement subspace [3, Section 5.4]. However, we remark that Pidstrigach's implementation [3, Section 6] utilizes UNet to parameterize the score, limiting the evaluation of the score function to the training interval only. In contrast, our method is not limited to the grid on which we initially train our network. In other words, we truly incorporate the advantages of the infinite-dimensional approach. Additionally, we remark that the projection-type methods are primarily heuristic and suffer from instability when dealing with ill-posed IPs. Recently, some workarounds have been proposed to tackle these issues in finite dimensions [4, 5], but we remark that our method, together with not assuming any knowledge of the forward model and operating directly in infinite dimensions, is also designed to offer an alternative to data-specific inference. In fact, projection-type methods can involve costly forward operator computations. Our method, instead, learns an amortized version of the conditional score and, by doing so, addresses a critical gap in the literature involving infinite-dimensional SDMs. There exist cases in which amortized methods can be a preferred option for Bayesian IPs (see line 119 and accompanying references in our paper).\n\n3. **Difference between (2) and (3).** In finite dimensions, the drift of the reverse stochastic differential equation (SDE) involves the score function $\\nabla \\log p_t$, where $p_t$ represents a density with respect to a Lebesgue measure. However, in infinite-dimensional Hilbert spaces, this density is no longer well-defined (the Heine-Borel theorem does not hold in such spaces). As a result, the left-hand side of equation (2) cannot be interpreted literally in infinite dimensions. In equation (3), $S(t,x)$ is thus defined formally, leveraging the fact that the right-hand side of equation (2) is well-defined in infinite dimensions. This implies that we need to demonstrate that $S(t,x)$ truly represents the score function, i.e., the drift of the reverse SDE. Analyzing the forward-reverse conditional SDE framework is an important contribution of our paper.\n\n4. **Is it possible to add results on large-scale problems?** Yes, see the PDF included in our global response. Figure 2 illustrates a new experiment, which we will include in the appendix, demonstrating the applicability of our method to a large-scale inverse problem in geophysics, i.e. linearized seismic imaging via the Born approximation with a $256 \\times 256$-dimensional unknown parameter.\n\n**References**\n\n[1] Y. Song et al, *Solving inverse problems in medical imaging with score-based generative models*, 2022.\n\n[2] G. Batzolis et al, *Conditional image generation with score-based diffusion models*, 2021.\n\n[3] J. Pidstrigach et al, *Infinite-dimensional diffusion models for function spaces*, 2023.\n\n[4] H. Chung et al, *Diffusion posterior sampling for general noisy inverse problems*, 2023.\n\n[5] Y. Wang et al, *Zero-shot image restoration using denoising diffusion null-space model*, 2022. "}, {"review_id": "hBRUoA9gm0", "reviewer": "Reviewer_Q2NB", "summary": "This paper mathematically examines linear inverse problems in infinite dimensional vector spaces. Particularly, it is proved that the conditional denoising estimator is a consistent estimator of the conditional score in infinite dimension. ", "questions": "Gaussian process (GP) is a popular method for inverse problem in Hilbert space. I wonder if the current problem would have a certain connection to GP, in particular, in the case of Gaussian prior discussed in section 4. ", "limitations": "Descriptions on the motivation to examine the inverse problem in Hilbert space are lacking.  In what practical situations do inverse problems in Hilbert space come out? For instance, GP is widely used for Bayes optimization, which is an optimization scheme for black-box functions. Additional writing about possible applications of the inverse problem in Hilbert space would make the paper more attractive. ", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The consistency of the conditional denoising estimator in infinite dimensional vector spaces is mathematically shown. For a specific case of Gaussian prior, the forward-reverse SDEs are solved exactly, which shows an exponentially fast convergence in the reverse SDE. A sufficient condition for the success of the score-based diffusion model framework is presented for the infinite dimensional version.  ", "weaknesses": "The considered inverse problems are not of the plug-and-play type, which would limits the practical utility. Numerical example is limited to a one-dimensional toy model. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the Reviewer for their positive feedback.\n\nWe are happy to clarify our manuscript in response to the Reviewer's remarks and questions.\n\n1. **Practical utility against plug-and-play-type approach.** While plug-and-play methods are indeed very popular, we would like to emphasize that we deliberately chose not to focus on them, first and foremost because we wanted to provide theoretically-grounded guarantees for an approach that is not heuristic. Furthermore, it is worth noticing that current implementation of plug-and-play approaches using diffusion models for conditional sampling in the infinite-dimensional setting [1, Section 6] does not fully exploit the ``discretization-invariance'' property achieved by studying the problem in infinite dimensions. Pidstrigach employ a UNet to parametrize their score, which restricts the evaluation of their score function to the training interval. Consequently, when dealing with a specific inverse problem on a new grid, they would need to gather new training data on that grid, retrain the neural network, and consequently limit the use of the score to that specific grid. In contrast, our implementation allows us to handle new discretizations without requiring additional training data, as demonstrated in the new experiment that we intend to include in the final version of the paper (see Figure 1 in the attached PDF for preliminary results). Therefore, we believe that our approach is not only theoretically grounded but also offers a broader applicability for our method.\n\n2. **Numerical example is limited to a one-dimensional toy model.** We acknowledge your remark, along with those of the other reviewers, regarding the need for\nimprovement in our toy example. In our global response, we have included a PDF that\ncontains additional experiments, illustrating the applicability of our method to a large-scale inverse problem in geophysics, specifically linearized sesmic imaging via the Born\napproximation (Figure 2 in the attached PDF) with a $256 \\times 256$-dimensional unknown parameter. We will add these experiments in the paper.\n\n3. **Connection to Gaussian Process.** Gaussian process regression is indeed related to the calculations presented in section 4 which is in the Gaussian framework. GP has been used to solve inverse problems because of - at least - two fundamental reasons. First it naturally arises when data is contaminated by noise, and additive Gaussian  noise is the most simple noise model. Second it allows the use of simple but powerful theorems (such as Gaussian conditioning theorem). Moreover, in the infinite-dimensional setting, the GP approach is possible via the use of Gaussian measures (and appropriate theorems such as Feldman-Hajek theorem). In our paper we leverage these aspects of Gaussian analysis to gain novel and fundamental insights about the conditional score. Moreover, we emphasize the robustness of these insights  by introducing a prior that is absolutely continuous with respect to Gaussian measure. \n\n4. **Possible applications of inverse problems in Hilbert spaces.** Bayesian approach to inverse problems makes it possible to deal with under-determined and/or noisy inverse problems by an appropriate prior modelling. In the infinite-dimensional PDE context, this prior should sample functions in the suitable space functions, which are very often Hilbert spaces. For instance, the inverse heat equation (how to determine the initial condition of a heat or convection-diffusion equation from noisy measurements) or the elliptic inverse problem (how to determine the source of an elliptic equation from noisy measurements) presented in [2] are concrete examples of noisy and under-determined inverse problems and they are naturally formulated in Hilbert spaces. We will mention these examples in the paper.\n\n**References**\n\n[1] J. Pidstrigach, Y. Marzouk, S. Reich, and S. Wang, *Infinite-dimensional diffusion models for function spaces*, arXiv:2302.10130.\n\n[2] M. Dashti and A. M. Stuart, *The Bayesian approach to inverse problems*, in Handbook of Uncertainty Quantification, Springer, 2017, pp. 311-428."}, {"review_id": "dj01uysHci", "reviewer": "Reviewer_XSPZ", "summary": "- The study proposes a method to learn the posterior distribution in infinite-dimensional Bayesian linear inverse problems using amortized conditional Score-based Diffusion Models (SDMs). This extends conditional SMDs into the infinite-dimensional function space setting, as existing conditional SDMs have previously only dealt with finite-dimensional vector spaces (noting also that _unconditional_ SDMs have recently been extended to infinite dimensional vector spaces by Pidstrigach et. al).  This leads the way for applications in, for example PDEs, where the unknown parameters to be estimated take the form of functions. \n- The key technique underlying their approach is to define the _conditional_ score in an infinite dimensional setting, extending the method of Pidstrigach et. al, who defined the _unconditional_ score in the infinite dimensional setting. \n- Using their definition of the conditional score in infinite dimensions, this allows them to avoid having to solve a potentially expensive proximal optimization step, as was done by Pidstrigach et. al. \n- The authors then provide a comprehensive theoretical analysis of the use of their conditional score in SDMs and show: \n  - How this newly defined score is used as a reverse drift of the diffusion process, which leads to a generative model that samples from the correct target conditional distribution under certain conditions.\n  - That as long as you start from the invariant distribution of the diffusion process, the reverse SDE converges to the target distribution exponentially fast \n  - By explicitly computing the expected square norm of the conditional score, they shows that a uniform in time estimate is not always true for the conditional score. This leads them to provide a set of conditions to be satisfied to ensure a uniform in time estimate for a general class of prior measures.\n  - That the conditional score can be estimated via a conditional denoising score matching objective in infinite dimensions.\n  - That the conditional denoising estimator is a consistent estimator of the conditional score in infinite dimensions.\n  - That unlike the unconditional score, for noiseless observations the conditional score blows up as T->0\n- They present a small toy experiment that validates their approach by demonstrating the applicability of their method in approximating non-Gaussian multi-modal distributions.\n\n\n\n", "questions": "- What are some other applications of this work besides PDEs? Can you include a few of them in the introduction so the reader doesn't think applications of your work is limited only to PDEs? \n\n- In the conclusion the authors say their method \"is able to perform conditional sampling directly on infinite dimensional Hilbert spaces.\" Don't Gaussian processes also allow you to perform conditional sampling directly on infinite dimensional Hilbert spaces? What are the differences between your approach and Gaussian processes? When would I want to use one vs. the other? Could a GP be used as a baseline in your Figure 1? If so it would be very interesting to see how it compares. \n\n- The authors also say \"We also show that the conditional score can have a singular behavior at small times when the observations are noiseless, in contrast with the unconditional score under similar hypotheses.\". Is this just a limiting phenomenon of theoretical interest, or do you expect this to cause difficulties in practice? Are there any applications where we should expect observations to be noiseless or does that assumption never hold in practice? \n", "limitations": "- There does not appear to be a \"Limitations and Broader Impacts\" statement in this work. ", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "I will address strengths and weaknesses across the four dimensions (Originality, Quality, Clarity, Significance) below. ", "weaknesses": "**Originality:** \n- Are the tasks or methods new? \n  - Yes, the authors make a novel contribution to the quickly-growing diffusion model literature, by showing how to extend conditional SDMs into the infinite-dimensional vector space, thus opening the door to a wider range of applications (PDEs, etc)\n- Is the work a novel combination of well-known techniques? Is it clear how this work differs from previous contributions? \n  - The authors are very clear about: \n     1. How they are starting with the recently-proposed framework of Pidstrigach et al\n     2. The specific point where they deviate from, and then extend, Pidstrigach's work (Definition 2, eq. 8)\n  - Further, they provide a comprehensive analysis of the use of their defn 2 in SDMs in sections 4 and 5. \n  - However I should state that I am not adequately familiar with the mathematical techniques deployed in this paper to check their technical claims for correctness or novelty, so please differ to another reviewer with more expertise in this area. \n- Is related work adequately cited?\n  - Yes. \n\n**Quality:** \n- Is the submission technically sound?\n  - It appears to be, although as previously stated I am not sufficiently versed in their mathematical techniques to be 100% sure. \n- Are claims well supported (e.g., by theoretical analysis or experimental results)? \n  - Theoretically yes, experimentally no. I realize this is a theory paper and don't necessarily expect full-scale experiments on PDEs, but in Section 6 I expected to see at least 2 obvious baselines which were foreshadowed in the body of the text but not experimentally validated on the toy experiments. These are:\n    1. How does their method compare to the crude, discretization-based approach discussed in lines 56 - 61:\n        > A straightforward solution may be to discretize the infinite-dimensional input and output function spaces into finite-dimensional vectors, and apply SDMs to learn the posterior. Yet theoretical studies of current DMs suggest that performance guarantees do not generalize well on increasing dimensions [7, 9, 33]. This is precisely why Stuart’s guiding principle to study a Bayesian inverse problem for functions— “avoid discretization until the last possible moment” [41] — is more than ever critical to the use of SDMs. \n    A primary motivation for their method is \"Stuart's Principle\", which says to avoid discretization until the last possible moment. I would have liked to see experimentally why this principle is so important on the simple toy examples they have provided.\n    2. How does their method compare to the approach proposed by Pidstrigach for conditional sampling, which uses a proximal optimization step? This is a second primary motivation for their approach (and appears in the abstract) - their method is more performant because it avoids solving an optimization problem at each timestep. As far as I can tell it has not been experimentally validated that their approach is faster or more performant than the baseline method of Pidstrigach. Evidence showing their approach either gets better samples, or get samples of the same quality but more efficiently, is needed, given that a primary motivation for their approach is that the baseline method of Pidstrigach may be too computationally costly because of their use of proximal optimization. \n- Are the methods used appropriate? \n  - Yes\n- Is this a complete piece of work or work in progress? \n  - Yes, modulo the missing baselines discussed above. I believe such baselines are needed to consider this a complete piece of work, given how \"Stuart's Principle\" and \"avoid proximal optimization\" play a key role in the storyline and motivation for their technique. \n- Are the authors careful and honest about evaluating both the strengths and weaknesses of their work?\n  - Yes. In particular I found their observation, that unlike the unconditional score, for noiseless observations the conditional score blows up as t->0, particularly interesting. However in the conclusion I would have liked to read more about the authors' reflections on the strengths/weaknessness/future directions of their approach, both from a technical standpoint and from the standpoint of potential downstream applications of this work. \n\n**Clarity:** \n- Is the submission clearly written? \n  - I found all writing up until section 4 relatively easy to follow. I began to get lost around sections 4 and 5 and couldn't follow the math, but I attribute this largely to not being comfortable with stochastic differential equations. I could still follow the high-level plot in these sections, but will defer to other reviews to evaluate the technical claims. \n- Is it well organized? \n  - Yes \n- Does it adequately inform the reader? \n  - Yes \n\n**Significance:**\n- Are the results important? \n  - Yes, although I would have liked the authors to motivate the applications of their approach more thoroughly. I believe they only listed PDEs as an example for why you would want to use this approach but surely there are more applications than just PDEs, no? \n- Are others (researchers or practitioners) likely to use the ideas or build on them? \n  - Yes\n- Does the submission address a difficult task in a better way than previous work? \n  - Theoretically it appears so, but baselines are needed to experimentally validate these claims. \n- Does it advance the state of the art in a demonstrable way?\n  - Theoretically it appears so, but baselines are needed to experimentally validate these claims. \n-  Does it provide unique data, unique conclusions about existing data, or a unique theoretical or experimental approach?\n  - Yes, they build on the work of Pidstrigach in a novel way. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the Reviewer for their feedback.\n\nWe are happy to clarify our manuscript in response to the Reviewer's questions. We hope that this could lead to an improvement in their assessment of the paper.\n\n1. **Discretization-based approaches and Pidstrigach's procedure.** While it would have been interesting to extend the section dedicated to numerical experiments with explicit comparison with other methods, we want to point out that the limitations of crude discretization-based approaches are well-documented in the literature on Bayesian inverse problems. Consequently, we made a deliberate decision to avoid focusing on them experimentally. Moreover, in the literature on score-based diffusion models (SDMs), there are already instances where results show degradation when generalizing to infinity [1, Figure 5]. \n\n    Regarding the comparison between our method and the one by Pidstrigach [2], we want to stress that our method is not necessarily more performant$-$that was never our claim. It really depends on the task to be solved. Our aim was to frame the comparison between the two methods through the lens of the more general comparison between case-specific inference and amortized inference. Amortized methods *can* be a preferred option in Bayesian inverse problems. In our paper, we cite some works where these examples are provided (line 119).\n\n   Given that there is abundant literature on both topics---the importance of discretizing as late as possible and amortized methods---we believe that it is more interesting, since this is a theory paper, to focus on the fact that while Pidstrigach's method for conditional sampling is heuristic, ours *provably* samples from the posterior and is discretization-invariant. Our main contribution is to offer theoretical guarantees for sampling from the posterior, with the additional appeal of providing a method for practitioners seeking discretization-invariant amortized DMs. In this sense, it is important to note that our method stands in comparison to Pidstrigach's one. In fact, their implementation is essentially finite. Pidstrigach use a UNet to parametrize their score [2, Section 6], which restricts its evaluation to the training interval. If dealing with a specific inverse problem on a new grid, they would need to gather new training data, retrain the NN, and limit the use of the score to that specific grid. In contrast, our approach allows us to move away from the initial grid, effectively taking advantage of the discretization-invariance property. This flexibility ensures a broader applicability of our method. Your comment helped us recognizing that we haven't stressed this difference enough, so we will add a remark and a new numerical experiment (see Figure 1 in the PDF included in our global response).\n\n   Finally, we agree with your remark that our toy example needs some improvement to better fit the storyline. In the PDF we included preliminary results, namely Figure 1, illustrating that the method can handle different discretizations. This serves to underscore the importance of Stuart’s principle. Additionally, we incorporated an experiment on a large-scale inverse problem in geophysics, specifically linearized seismic imaging via the Born approximation with a $256 \\times 256$-dimensional unknown parameter, to enhance the comprehensiveness of our results (see Figure 2).  We will add these experiments in the paper.\n\n2. **Concluding remarks.** Given the theoretical nature of the paper, summarizing a series of details from the proofs and presenting them in a self-contained section is challenging. Nonetheless, we recognize the value in discussing the strengths, weaknesses, and future directions of our method. Therefore, we intend to address this aspect by adding a few remarks throughout the text incorporating, among others, Reviewer hg8t's second question.\n\n3. **Applications beyond PDEs.** While applications of our work to PDE-based inverse problems are the most natural ones, they are by no means the only ones. In the non-PDE class we may think at geometric inverse problems (e.g. how to determine the Riemann metric from geodesic information, or the background velocity map from travel time information in geophysics) [3] or inverse problems for singular integral operators [4]. We will cite these examples in our paper.\n\n4. **Can a GP be used as a baseline in Figure 1?** In the Gaussian framework of Section 4 this would make sense but in the general framework of Sections 5-6 this does not seem possible. It is clear that the conditional distribution in the example of Section 6 is strongly bimodal so GP regression does not seem appropriate to address this example.  \n    \n5. **Singularity of the conditional score.** The singularity in the conditional score as noise vanishes has been investigated in finite dimensions. It is not merely a theoretical phenomenon but has practical implications. For a discussion about the difficulties this singularity causes and numerical methods for approximating the score at small times in finite dimensions, we refer to [5, 6]. As for our paper, we anticipate that Assumption 1 is reasonably easy to satisfy. However, we must be mindful of the possibility of encountering a blow-up of the score under the conditions described by Reviewer hg8t. To address this concern, we will add a remark.\n\n**References**\n\n[1] A. Phillips et al, *Spectral diffusion processes*, 2022.\n\n[2] J. Pidstrigach et al, *Infinite-dimensional diffusion models for function spaces*, 2023.\n\n[3] G. Uhlmann and A. Vasy, *The inverse problem for the local geodesic ray transform*, 2016.\n    \n[4] A. Dynin, *Inversion problem for singular integral operators:* $C^*$*-approach*, 1978. \n\n[5] D. Kim et al, *Soft Truncation: A Universal Training Technique of Score-based Diffusion Model for High Precision Score Estimation*, 2022.\n\n[6] T. Dockhorn et al, *Score-Based Generative Modeling with Critically-Damped Langevin Diffusion*, 2022."}]

Author rebuttal:
We thank the Reviewers for their valuable and constructive feedback. 

Based on your comments, we have taken significant steps to enhance our Section 6. We are now including additional experiments that demonstrate the applicability of our method to large-scale problems and showcase its discretization invariance. **Preliminary results can be found in the attached PDF, which includes relevant figures for your reference.** 

In accordance with your suggestions, we have outlined the following improvements: 

1. **First Experiment.** Taking advantage of the additional content page for the camera-ready version, if our paper is accepted we will extend Section 6 by incorporating an experiment demonstrating the discretization invariance of our proposed method by sampling a bi-modal non-Gaussian conditional distribution over grids with varying discretizations. In particular, we present results related to sampling the posterior over various discretizations of the $[-3, 3]$ domain. Figure 1 in the attached PDF showcases the outcomes of this experiment by displaying the samples and the marginal conditional distributions at $y=-1, 0, 0.5$ (we used the relation $X_0= a y^2 + \epsilon$). We enhanced robustness against different discretizations by training the score-based model on data residing on nonuniform grids containing $15$ to $50$ grid points. No changes were made to other hyperparameters compared to the original experiment. 

2. **Second Experiment.** We will add a new experiment in the appendix to demonstrate the applicability of our method to a large-scale inverse problem in geophysics, specifically linearized seismic imaging via the Born approximation (Figure 2 in the PDF). More specifically, the problem we address involves estimating the short-wavelength component of the Earth's unknown subsurface squared-slowness model (refer to Figure 2a in the uploaded PDF), using measurements collected at the surface. This inverse problem, known as *seismic imaging*, can be formulated as a linear inverse problem by linearizing the nonlinear relationship between recorded data and the squared-slowness model, as governed by the wave equation. In its simplest acoustic form, the linearization around a background squared slowness model (illustrated in Figure 2b of the uploaded PDF) leads to an inverse problem for estimating the true seismic image (depicted in Figure 2a).

   Given the high dimensionality of the observed data, we summarize it by projecting the data back into seismic image space using the adjoint Born operator, as shown in Figure 2c. This process leads to the reverse-time migrated image, which we utilize instead of data to train our conditional score-based model. We carried out the training for $300$ epochs using $4750$ pairs of true seismic images and associated reverse-time migrated images from the 3D Parihaka real dataset, with each image being of size $256 \times 256$. We used a batch size of $128$ and an initial learning rate of $2 \times 10^{-3}$, which decays to $5 \times 10^{-4}$ over the epochs following a power-law rate of $-1/3$.

   Post-training, for a new seismic image (refer to Figure 2a in the uploaded PDF), we simulate seismic data and compute the reverse-time migrated image using the adjoint Born operator. The trained conditional score-based model is then employed to sample from the posterior distribution of the seismic image, given the reverse-time migrated image. We draw $1000$ samples to compute the conditional (posterior) mean estimator, visualized in Figure 2c of the PDF. These samples are also utilized to calculate pointwise standard deviations (Figure 2d) as a measure of uncertainty.

   As anticipated, the pointwise standard deviation highlights areas of high uncertainty, particularly in regions with complex geological structures$-$such as near intricate reflectors and areas with limited illumination (deep and close to boundaries). The regions of significant uncertainty correspond well with challenging-to-image sections of the model. This observation becomes more apparent in Figures 2h and 2i, displaying two vertical profiles with $99$\% confidence intervals (depicted as orange-colored shading), which demonstrate the expected trend of increased uncertainty with depth. Furthermore, we notice that the ground truth (indicated by dashed black lines) largely falls within the confidence intervals for most areas.

   We also observe a strong correlation between the pointwise standard deviation and the error in the conditional mean estimate (Figure 2e), confirming the accuracy of our Bayesian inference method. To prevent bias from strong amplitudes in the estimated image, we present the normalized pointwise standard deviation divided by the envelope of the conditional mean in Figure 2f. This visualization provides an amplitude-independent assessment of uncertainty, highlighting regions of high uncertainty at the onset and offset of reflectors (both shallow and deeper sections). Additionally, the normalized pointwise standard deviation underscores uncertainty in areas of the image where there are discontinuities in the reflectors (indicated by black arrows), potentially indicating the presence of faults.

With these updates and the additional improvements we are implementing based on our responses in the rebuttal, we sincerely hope to have addressed the concerns raised by the Reviewers. If our responses satisfied your concerns, we kindly ask that you consider revisiting your score accordingly. Please let us know if there is anything else that we can do or clarify to enhance the quality of this paper. Once again, thank you for your feedback.

Meta-review: The reviewers agree that is paper is an interesting and useful contribution to the literature, especially after the discussions with the authors. Please reread the reviews and your responses carefully, and make the necessary changes in the final submission.
Decision: Accept (spotlight)