Variational Gaussian processes for linear inverse problems
Authors: Thibault Christophe RANDRIANARISOA
Botond Szabo
Venue: NeurIPS 2023 poster (2023)
Keywords: Linear inverse problems
Gaussian processes
Variational inference
Inducing variables
Asymptotics
Contraction rates
TL;DR: We derive minimax contraction rates results for variational Gaussian process posteriors for linear inverse problems.

Abstract: By now Bayesian methods are routinely used in practice for solving inverse problems. In inverse problems the parameter or signal of interest is observed only indirectly, as an image of a given map, and the observations are typically further corrupted with noise. Bayes offers a natural way to regularize these problems via the prior distribution and provides a probabilistic solution, quantifying the remaining uncertainty in the problem. However, the computational costs of standard, sampling based Bayesian approaches can be overly large in such complex models. Therefore, in practice variational Bayes is becoming increasingly popular. Nevertheless, the theoretical understanding of these methods is still relatively limited, especially in context of inverse problems.

In our analysis we investigate variational Bayesian methods for Gaussian process priors to solve linear inverse problems. We consider both mildly and severely ill-posed inverse problems and work with the popular inducing variable variational Bayes approach proposed by Titsias [Titsias, 2009]. We derive posterior contraction rates for the variational posterior in general settings and show that the minimax estimation rate can be attained by correctly tunned procedures. As specific examples we consider a collection of inverse problems including the heat equation, Volterra operator and Radon transform and inducing variable methods based on population and empirical spectral features.

Variational Gaussian Processes For Linear Inverse Problems
Thibault Randrianarisoa Department of Decision Sciences Bocconi University
Botond Szabo
Department of Decision Sciences
Bocconi University
via Roentgen 1, 20136, Milano, MI, Italy
via Roentgen 1, 36, Milano, MI, Italy
botond. szaboouniboccon........
Abstract
By now Bayesian methods are routinely used in practice for solving inverse prob- lems. In inverse problems the parameter or signal of interest is observed only indirectly, an image of given map, and the observations are typically further corrupted with noise. Bayes offers a natural way to regularize these problems via the prior distribution and provides a probabilistic solution, quantifying the remain- ing uncertainty in the problem. However, the computational costs of standard, sampling based Bayesian approaches Therefore, practice variational Bayes is becoming increasingly popular. Never- theless, the theoretical understanding of these methods is still relatively limited, especially in context of inverse problems. In our analysis we investigate variational Bayesian methods for Gaussian process priors solve linear inverse problems. We consider both mildly and severely ill-posed inverse problems and work with the popular inducing variables variational Bayes approach proposed by Titsias (57]. We derive posterior contraction rates for the variational posterior in general settings and show that the minimax estimation rate can be attained by correctly tunned procedures. As specific examples we consider a collection of inverse problems including the heat equation, Volterra operator and Radon transform and inducing variable methods based on population and empirical spectral features.
be overly large in such complex models.
Introduction
In inverse problems we only observe the object of interest (i.e. function or signal) indirectly, through a transformation with respect some given operator. Furthermore, the data is typically corrupted with measurement error or noise. In practice the inverse problems are often ill-posed, i.e. the inverse of continuous. Based on the level of ill-posedness ill-posed cases. The ill-posedness of the problem prevents us from simply inverting the operator it would blow up the measurement errors in the model. Therefore, overcome this problem, regularization techniques are applied by introducing a penalty term in the maximum likelihood approximation. Standard examples include generalized Tikhonov, total variation and Moore-Penross estimators, see for instance 4 9
the operator
distinguish mildly and severely
18 56] or a recent survey 3 on data-driven methods for
solving inverse problems
An increasingly popular approach to introduce regularity to the model is via the Bayesian paradigm, for instance 3 12
29 54 and references therein. Beside regularization, Bayesian methods provide probabilistic solution to the problem, which can be directly used to quantify the remaining uncertainty of the approach. This is visualised by plotting credible sets, which are sets accumulating prescribed percentage of the posterior mass. For computing the posterior typically MCMC algorithms used, however, these can scale poorly with increasing sample size due the complex structure of the likelihood. Therefore, in practice often alternative, approximation methods are used. Variational
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Bayes (VB) casts the approximation of the posterior into an optimization problem. The VB approach became mereasingly popular scale up Bayesian inverse problems. see for instance the recent papers 31 42] and references therein. However, until recently these procedures were considered black box methods basically without any theoretical underpinning. Theoretical results are just starting to emerge [2] 46 63 64, 66], but we still have limited understanding of these procedures in complex models, like inverse problems.
[22
In our analysis we consider Gaussian process (GP) priors for solving linear inverse problems. For Gaussian likelihoods, due to conjugacy, the corresponding posterior has an analytic form. Nevertheless, they are applied more widely, in non-conjugate settings as well. However, training and prediction even in the standard GP regression model, scales as O (n*) ) (or O (n²) for exact inference in the recent paper [16] leveraging advances in computing hardware) and (n²),. respectively, which practically limits GPs to a sample size n of order 10t. Therefore, in practice often not the full posterior, but an approximation is computed. Various such approximation methods were proposed based on some sparse or low rank structure, see for instance [13 14 33 49; 50 51 52 57. Our focus here lies on the increasingly popular inducing variable variational Bayes method introduced in 57 58
In our work we extend the inducing variable method for linear inverse problems and derive theoretical guarantees for the corresponding variational approximations. More concretely we adopt a frequentist Bayes point-of-view in our analysis by assuming that there exists a true data generating functional parameter of interest and investigate how well the variational posterior can recover this object. We derive contraction rates for the VB posterior around the true function both in the mildly and severely ill-posed inverse problems. We then focus on two specific inducing variable methods based on the spectral features of the prior covariance kernel. We show that for both methods if the number of inducing variables are chosen large enough for appropriately tunned priors the corresponding variational posterior concentrates around the true function with the optimal minimax estimation rate. One, perhaps surprising aspect of the derived results is that the number of inducing variables required to attain the optimal, minimax contraction rate is sufficiently less in the inverse setting than in the direct problem. Therefore, inverse problems can be scaled up at higher degree than standard regression models.
Related literature. The theory of Bayesian approaches to linear inverse problems is now well established. The study of their asymptotic properties started with the study of conjugate priors
26 27 28) before addressing the non-conjugate case [25] 45] and rate-adaptive priors By now we have a good understanding of both the accuracy of the procedure for recovering the true function and the reliability of the corresponding uncertainty statements. The theory of Bayesian non-linear inverse problems is less developed, but recent years have seen increasing interest in the topic, see the monograph (36] and references therein. Some algorithmic developments for variational Gaussian approximations in non-linear inverse problems, and applications MCMC sampling, can be found in 440
[
20
26
55
The inducing variable approach for GPs proposed by [57] 58 has been widely used in practice. Recently, their theoretical behaviour was studied in the direct, nonparametric regression setting. In [ it was shown that the expected Kullback-Leible divergence between the variational class and posterior tends to zero when sufficient amount of inducing variables were used. Furthermore, optimal contraction rates (37] and frequentist coverage guarantees [88 59 60) were derived for several inducing variable methods. Our paper focuses on extending these results to the linear inverse setting.
Organization. The paper is organized follows. In Section 2 we first introduce the inverse regression model where we carry out our analysis. Then we discuss the Bayesian approach using GPs and its variational approximations in Sections 2.1 and 2.2, respectively. As our main result derive contraction rates for general inducing variable methods, both in the mildly and severely ill- poo-p cases. Then in Section 2 we focus two specific inducing variable methods based on spectral features and provide more explicit results for them. We apply these results for a collection of examples, including the Volterra operator, the heat equation and the Radon transform in Section 3 Finally, we demonstrate the applicability of the procedure in the numerical analysis of Section and conclude the paper with discussion in Section 5 The proof of the main theorem together with technical lemmas and additional simulation study are deferred to the supplementary material.
Notation. Let C, c be absolute constants, independent of the parameters of the problem whose values may change from line to line. For two sequences (a,)) and (b,,) of numbers, an S b, means that there
2

---

exists a universal constant such that an I cb,, and we write a, x bn if both an X bn and b, S an hold simultaneously. We denote by an < b, if lan/bal tends zero. The maximum and minimum of two real numbers a and b are denoted by V and a A b, respectively. We use the standard notation 8.j j=j- For m 2 1, we note 444 the set of positive-deffnite matrices of size m x m.
2
Main results
In our analysis we focus on the non-parametric random design regression model where the functional parameter is observed through a linear operator. More formally, we assume to observe i.id. pairs of random variables (I,, Yi)==...... satisfying
Y (Afo) (1ra) Zi,
Z, N(O, 3z cid G i= i=1,...,,
(1)
where fo E (2),; for some domain T TRd and measure l on T, the underlying functional parameter of interest and A: [1) + (22;; (e(A')), G), for the measure G == Au on X, is known, injective, continuous linear operator. In the rest of the paper we use the notation P/o and E for the joint distribution and the corresponding expectation, respectively, of the data (X,Y) (xi, I.......... Furthermore, we denote by Ex, Eyjx, Pylx the expectation/distriuuiti G and the law of (Yi)s given the design respectively. Finally, for simplicity we take o² in our computations.
under
In the following, denoting A" the adjoint of A, we assume that the self-adjoint operator A: L2(T:;) + (2(T: possesses countably many positive eigenvalues (k])) with respect to the eigenbasis (ez) (which is verified if A is compact operator for instance). We remark that (g5), defined by Ae,- Kygj is an orthonormal basis of (2; G). We work the ill ill-posed problem where Kj + O, the rate of decay characterizing the difficulty of the inverse problem.
Definition 1. We say the problem is mildly ill-posed problem of degree p 0 Of Ki = j. has polynomial decay. In the severely ill-posed problem, the rate we consider is exponential, Ky x 55 for e > 0,p 2 1, and p is the degree of ill-posedness once again.
In nonparametrics it is typically assumed that fo belongs to some regularity class. Here we consider the generalized Sobolev space
H = {
(2T:N)):
1 < x 11 = -EMRI..
(2)
for some 3 O. We note that the difficulty in estimating fo from the data is twofold: needs to deal with the observational noise, which is a statistical problem, as well as to invert the operator A. which comes from inverse problem theory. As result of the ill-posedness of the problem, recovering fo from the observations may suffer from problems of unidentifiability and instability. The solution these issues is incorporate some form of regularization the statistical procedure. The Bayesian approach provides natural way to incorporate regularization into the model via the prior distribution on the functional parameter. In fact penalized likelihood estimators can be viewed as the maximum a posteriori estimators with the penalty term induced by a prior. For example Tikhonov type regularizations can be related to the RKHS-norm of Gaussian Process prior, see 35 for more detailed discussion.
2.1 Gaussian Process priors for linear inverse problems
We focus on the Bayesian solution of the inverse problem and exploit the Gaussian likelihood structure by considering conjugate Gaussian Process (GP) priors on f. A GP GP (n(-), ()) is a set of random variables (f(() E T), such that any finite subset follows Gaussian distribution. The GP is described by the mean function and a covariance kemel k(t,t'). We consider centered GPs priors (i.e. take 1) = 0). Then the bilinear, symmetric nonnegative-deffite function k: T x T + determines the properties of the process (e.g., its regularity). In view of the linearity of the operator A the corresponding posterior distribution is also a Gaussian process. The mean and covariance function of the posterior is given by
+ KIAS (KASAS + +z2 y, (t, 8) + k(t,s) KIAS (KATAS+ 1z2 KA
(3)
3

---

where = (y1......u.) Af (=.())))==1.... KASAs - -EAAAf E R with En denoting the expectation with respect to the GP prior II, KTAS = (EILA(E)))())) i1 E R", see the supplement for the detailed derivation.
Due to the closed-form expressions for the posterior and the marginal likelihood, well as the simplicity with which uncertainty quantification may be produced, GP regression has gained pop- ularity [44]. Furthermore, the asymptotic frequentist properties of posteriors corresponding to GP priors in the direct problem, with A taken to be the identity operator, is well-established by now. Optimal contraction rates and confidence guarantees for Bayesian uncertainty quantification were derived in the regression setting and beyond, see for instance [10 39 47
53
48
61
62 65 and
references therein. In the following, we say that En is an po-poster contraction rate for the posterior II
[X YJ [ if for any Mn +
E f: II foll 2 M.C. x,Y] + 0.
In our analysis we consider covariance kernels with eigenfunctions coinciding with the eigenfunctions of the operator A",, i.e we take
k(t, = E, (55)(0ee(())
(4)
where (Aj)y denote the corresponding eigenvalues. The asymptotic behaviour of the corresponding posterior has been well investigated in the literature both in the mildly and severely ill-posed inverse
problems. Rate optimal contraction rates and frequentist coverage guarantees for the resulting credible
sets were derived both for known and unknown regularity parameters [15] 26 27 28 [55. These
results were further extended for other covariance kernels where the eigenfunctions do not exactly match the eigenfunctions of the operator A, but in principle they have to be closely related, see 20 25 45).
However, despite the explicit, analytic form of the posterior given in (3 and the theoretical underpin- ning, the practical applicability of this approach is limited for large sample size n. The computation of the posterior involves inverting the n-by-n matrix K KASAS + L²Ln, which has computational complexity Onn3). Therefore, in practice often not the true posterior, but scalable, computationally attractive approximation is applied. Our focus here is the increasingly popular inducing variable variational Bayes method introduced in (57
2.2 Variational GP for linear inverse problems
In variational Bayes the approximation of the posterior is casted as optimization problem. First tractable class of distributions Q is considered, called the variational class. Then the approximation T is computed by minimizing the Kullback-Leible divergence between the variational class and the true posterior, i.e.
V arg infeeo KL (QIII IX,Y
There is a natural trade-off between the computational complexity and the statistical accuracy of the resulting approximation. Smaller variational class results in faster methods and easier interpre- tation, while more enriched classes preserve more information about the posterior ensuring better approximations.
In context of the Gaussian process regression model (with the operator A taken to be the identity), (58) proposed a low-rank approximation approach based inducing variables. The idea is to compress the information encoded in the observations of size n into m so called inducing variables. We extend this idea for linear inverse problems. Let consider real valued random variables (UI,.... Um E Le (II), expressed as measurable linear functionals of f and whose prior distribution is HI, In view of the linearity of u, the joint distribution of (f, u) is Gaussian process, hence the conditional distribution flu denoted by II(- u)), is also a Gaussian process with mean function and covariance kernel given by
t + Ktn K-lu and (t, + k((,s) Kru K Kus:
respectively, where K., = ((())) E T and Kun En(uu)) E RR.x Then the posterior is approximated via probability measure V, "T",B("")) by V [[[:|u]d (u), which is

---

absolutely continuous against II and satisfies dV (f) (
dv,,
(u(f)) We note that the variables u dI.
were first considered to be point evaluations of the GP prior process before these ideas were extended interdomain inducing variables, c.g. integral forms of the process 30 [58.
Taking u = =W(H (H Σ) as a multivariate Gaussian, the corresponding v x V u is a Gaussian process, with mean and covariance functions
t - Kiu SteKKutte Kudthe and (t,s) + k(t, - Ktn (Kuu -E) K)K
(5)
Letting n. and E be the free variational parameters, the variational family is taken as
@ = (v V = (u En), My E R", En E 5
consisting of "m-sspare" Gaussian processes.
By similar computations as those from 337. it can be shown that IT |X, YJ is equivalent to any element of Q (they are mutually dominated) so that the KT. divergence always finite and there exists V corresponding to the minimizer & of KL (VI]II[-|X,, Y]). Furthermore, we have
dvv
dtI (5)
dvr dll
(u)
exp(
(
/E (Y; Af(x)))) (TT())))
20 -202
Σm (Y, Karlayaa K kul))²
x erp
20
where KAs(a.j= (f,,)() One can observe that the parameters of the variational approxi- mations are
14 o 2un (Kuu (Kuu+o + KuA -1 KUATY,
(6)
Mn =Kuu (Kuu + o KuAF Kuus
for KuAf = RT = Enu(AP)T the m x matrix whose jth column is K Af(ajj and (Af)T
Then the explicit form for the variational posterior V/" can be attained by
plugging in the parameters (6) into the variational mean and covariance function (5) We also define
QAJAS KT K Kun KuAf-
We investigate the statistical inference properties of the above variational posterior distribution Vp" More concretely we focus how well the variational approximation can recover the underlying true functional parameter fo of interest in the indirect, linear inverse problem (]. We derive contraction rate for * both in the mildly and severely ill-posed inverse problem case. Furthermore, we consider both the standard exponential and polynomial spectral structures for the prior, i.e. we assume that the cigenvalues of the prior covariance kernel satisfics either XÚ X xj xj Ej" or X3 X j--1-2 2n for some 0, >> O. Finally, in view of (37), introduce additional assumptions on the covariance kemel of the conditional distribution of f u ensuring that the variational posterior is not too far from the true posterior in Kullback-Leibh divergence.
Theorem 1. Let's assume that fo E H and l1,5ll 5 for 3 0,7 2.
In the mildly-ill posed problem where Kj x j p > 0, if x j j-1-2a for a > and
(a 3) + p 3/2 + 27, the posterior contracts at the rate gim
n2 1420
2. In the severely ill-posed problem where K2 X e - 0.p 8/v
2 1, 1f = -3"" for
a > the posterior contracts at the rate elm log
Furthermore, if there exists a constant C. independent of n. such that
Ex |KAAAS QAsAsl SC, and ExTr (KTTAKAJA QASAS) Cne2
where En = n 1-20420 1, and En = n. -c/(5+22)) log 8/ptca/(Et22) in 2 V contracts around fo at the rate ehth for the mildly and severely ill-posed problems i.e.
aAS4p
E V [r: |1
- Joll 2 Mn inv
- 0, M₁ + oc.
5

---

Proof. We provide the sketch of the proof here, the detailed derivation of the theorem is deferred to the supplementary material. In a first step, we derive posterior contraction rates around Afo in empirical Z--norm under fixed design. In particular, we obtain an exponential decay of the probability expectation in the form
By E nt E i=1 - AAo) 2 Mne | X, X,Y]LA, LA, S Ce
(8)
eootp for arbitrary Mn + DO, where En = n 142 in the mildly and en =n ETIE log टाश in the severely ill-posed problems and A, is an event on the sample space IK" with probability tending to one asymptotically. In the mildly ill-posed case, this follows from results in [2].. while additional care is needed the severely ill-posed case. As a second step, we go back to the random design setting. We show, using concentration inequalities and controlling the tail probability of GPs in the spectral decomposition, that the empirical and population L2 -ooms are equivalent on large enough event. This implies contraction rate with respect to the || || (1(2))))) -oom around Afo- similarly to (8) In the third step, using the previous result on the forward map, we derive contraction rates around fo- To achieve this apply the modulus of continuity techniques introduced in [25 Notably, we extend their ideas to infinite Gaussian series priors in the severely ill- ill-posed case as well. Since in all these steps we can preserve the exponential upper bound for the posterior contraction (on large enough event), we can apply Theorem 5 of [46), resulting in contraction rates for the VB procedure. It requires control of the expected KL divergence between these two distributions, which follows from our assumptions on the expected trace and spectral norm of the covariance matrix of Af fu, see Lemma 3 in [37 for the identity operator A.
We briefly discuss the above results. First of all, the L2(T: ()) ccontracciionn rate of the true posterior, to the best of our knowledge, wasn't derived explicitly in the literature before, hence it is of interest in its own right. Nevertheless, the main message is that the variational posterior achieves the same contraction rate the true posterior under the assumption (7). Note that in the mildly ill-posed inverse problem case for eigenvalues > > j-188 (i.e. taking a = 3), the posterior contracts with the minimax rate n (0(1-22+2p) Note that the d- d-dmensiional case directly follows from this result when one defines the regularity class (2) and ill ill-posedness (Definition] I with B/d and p/d which would imply the rate -9/(4+2422p), Similarly in the severely-ill posed case one can achieve the minimax logarithmic contraction rate. Furthermore, the choice of the eigenvalue structure in the theorem was done for computational convenience, the results can be generalised for other choices of Xy as well. Though we considered the random variables u fixed as we do not optimize them above, they could conceivably be considered as free variational parameters and selected at the same time as y and Σ..
In the next subsection we consider two specific choices of the inducing variables, i.e. the population spectral feature method and its empirical counter part. We show that under sufficient condition on the number of inducing variables condition (7) satistied implying the contraction rate results derived in the preceding theorem.
2.3 Population and empirical spectral features methods
We focus here on two inducing variables methods, based on the spectral features (i.e. eigenspectrum) of the empirical covariance matrix K KASAS = An and the corresponding population level covariance operator (z, y) + En)
We start with the former method and consider inducing variables of the form
Ug Σ" Af(a), -= m,
where v, (o)...., e3) is the eigenvector of KAfAf corresponding to the jth largest eigenvalue Pj of this matrix. Similarly to the direct problem studied in [7], this results in (KArAF); = -pjoj. (KAru); = p, QASAs = ΣØ= P P3 and KAFAF QAFAF E P3V The computational complexity of deriving the first m eigenvectors of KASAf is O(n2)). This still quadratic in which sets limitations to its practical applicability, but it can be computed for arbitrary

---

choices of the prior covariance operator and map A. We also note that this choice gives the optimal rank-m upproximation Q Af of K KAfAf and it was noted in g that gives the minimiser of the trace and norm terms in
The second inducing variables method is based the eigendecomposition of covariance kernel (x. + En Let us consider the variables
" , ((())))())d((()), j= i=1,...;m.
(10)
Again, by extending the results derived in the direct problem [8] to the inverse setting, this results in (KAfAf)is = Agngdis, (KAfu).s = XyK QAIAS = Σ5= Azh H and KAJAS QAFAS = -i=++ )g Ø0 where Φ, = tational complexity of this method is O(nm²), which is substantially faster than its empirical counter part. However, requires the exact knowledge of the eigenfunctions of the prior covariance kernel, and therefore in general has limited practical applicability.
The compu-
Corollary 1. Let's assume that So E H1 l195lloo j for 3 0,7 2 0 and in the
1. mildly-ill posed case K2 x j.; take prior eigenvalues >j x j j-1-2a for some or > 0, (a A 8) + > 332 + 2%, number of inducing variables mn 2 n 0+2p+22) and denote by S +720
2. severely ill ill-posed case Kj > ==-5" take prior eigenvalues X3 x xj E5
- for a 2
& >> O. number of inducing variables m. 2 (5 + 2c) "log and introduce the notation
En = log 3/P n.
Then both for the population (if = O in 7 and empirical spectral features variational methods the corresponding variational posterior distribution contracts around the truth with the rate e i.e.
Ejo V f |If foll (((=))) 2 Ma + 0, Mn - 0c.
Remark 1. In the mildly ill-posed inverse problem taking a = 3 results in the minimax contraction
rate for my 2 142p4t5 Note that it is substantially less compared to the direct problem with O, hence the computation even faster in the inverse problem case.
3
Examples
In this section we provide three specific linear inverse problems as examples. The Volterra (integral) operator and the Radon transformations are mildly ill-posed, while the heat-equation is severely ill-posed inverse problem. We show that all cases by optimally tunning the GP prior and including enough inducing variables, the variational approximation of the posterior provides (from minimax perspective) optimal recovery of the underlying signal fo-
3.1
Volterra operator
First, let us consider the Volterra operator, A: L22 [2 - L2 (0, 1] satisfying that x Af(x) I f(s)ds, A f(x) = f f(s)ds.
(11)
The eigenvalues of A'A and the corresponding eigenbases are given by x3 1/2)-2,r= ey(x) 2 cos ((j ( 1/2)mx), (x)x) sin ((j 1/2) respectively, see (21). Therefore the problem is mildly ill-posed with degree p = and these bases are uniformly bounded, ie. sup, ||e V lلg co. The following lemma is then direct application of Corollary
(j
Corollary 2. Consider the Volterra operator in ① and assume that E H ffor some B 1/2. Set the eigenvalues in (4) >j --- Then the variational posterior V resulting from either the empiricul or population spectral features inducing variable methods achieves the minimax contraction rate if the number of inducing variables exceeds m. 2 n nz428, i.e. for arbitrary M,, - 00
Eso V [IV-- follza(o.1) 2 M,,nnM,nBB0 (B/((+43))) +.
7

---

3.2 Heat equation
Next let us consider the problem of recovering the initial condition for the heat equation. The heat equation is often considered as the starting example in the PDE literature and, for instance, the Black-Scholes PDE can be converted to the heat equation well. We consider the Dirichlet boundary condition
u(z,t) t), i u(z,t) t) 0x2
0²
u(x =)(x), u(0,t) u(L,t) 0,
(12)
Ot
for u defined on [0. 1] x (0,7). T O. For y L2(0, 1], u(x, 1) v2 Miee ff " sin(jax)) Therefore, if A D + D. with D = {f E 12] f(0) f(1) 00) such that, for μ f. Af(x) responding singular-values and singular- functions of the operator A are Ki = e j-'"x and e²(x) (x) = v2 sin( Therefore it is severely ill-posed problem with = 2 and c= n TT We also note that sup, ||e V |19 2o < co. This problem has been well studied both in the freeuentist 6 [9] 32] and Bayesian setting (28 [54]. Then, by direct application of Corollary D we can provide optimality guarantees for the variational Bayes procedure in this model as well.
with Mj = v2 .o ((s))in(jsssds.
u(2, T). then the cor-
Corollary 3. Consider the heat equation operator A above in the linear inverse regression model ① and assume that Jn HP for some B > O. Furthermore, set the eigenvalues X3 = j a 120 0, &> > 0 in 4 Then the variational approximation viz resulting from either of the spectral 122 n achieves the minimax
features inducing variables method with my Z (E 6 TTT) - 11/2 log-/2 contraction rate, i.e. for arbitrary Mn. + DO
Es 40 [| foll 10 2 Mn log 882 mg +0. 0.
3.3 Radon transform
Finally, we consider the Radon transform [24], where for some [Lebesgue--squar--nneeaab function f: : D + I defined on the unit disc D = { x R² ||~r||2 S 1}. we observe its integrals along any line intersecting D. we parameterized the lines by the length s E [o, 1] of their perpendicular from the origin and the angle [0.2x) of the perpendicular to the x-axis, we observe
/ f(8 cos - -ssin Φ sin & + cos Φ)dl,
Aʃ(s,) x)
(13)
2 1-
where ()) E S [000, 2T). The Radon transform is then map from A: (2
L2 (S; G), where y is -- times the Lebesgue measure and dG'(s, Φ) 2II dsdo. Then A is a bijective, mildly ill-posed linear operator of order = 1/4. Furthermore, the operator's singular value decomposition can be computed via Zernike polynomials Z (degree m, order k) and Chebyshev polynomials of the second kind Um (cos 0) sin ((m + 1)0) sin 0 I + 1, see [24]. Translating to the single index setting, we get for some functions 1, m N + N satisfying m(j) Vs and [(])1 S m(j), that
(r 0) = vm(j)= + 12!h6) 20))
(s( U Um() (s)e
polar coordinates are used on D. Therefore, we have sup, (lle v ||g llx) / < xc. Then, by directly applying Corollary to this setting the minimax contraction rate.
show that the variational Bayes method achieves
Corollary 4. Consider the Radon transform operator 13 in the inverse regression model (T) and let us take fo E H., 3 > 9/4. Taking polynomially decaying eigenvalues Xj X j-1125, the empirical spectral features variational Bayes method achieves the optimal minimax confraction rate if mn Z n. (3/22+20), i.e. for any M,, - 50
E "= [I foll (2(1;)) M,n -0/3322322228) - 0.
4 Numerical analysis
We demonstrate the approximation accuracy of the variational Bayes method on synthetic data. We consider here the recovery of the initial condition of the heat condition 3.2) which is a severely
8

---

ill-posed. In the supplement we provide additional simulation study for mildly ill ill-posed inverse problems as well. We set the sample size 8000, take uniformly distributed covariates on [0, 1), and let
0.4 (5xr), j odd,
Fo(() - -22 oj a;j00aaa (1+8) sin(jnt), Cj 2.5 + 2 v2 j even,
for 3 = The independent- observations are generated as Y, ~ (N(Af((z))) (2;), 1), depending on the solution of the forward map Afo after time T 10-2
We consider the prior with X = -Jj for { 10 view of Corollary the optimal number
of inducing variables is (g 2T² 1/2 log 1/2 6. We consider the population spectral feature method described (10) and plot the variational approximation of the posterior for m = and = 3 inducing variables in Figure We represent the true posterior mean by solid red and the upper and lower pointwise 2.5% quantiles by dashed red curves. The true function is given by blue and the mean and quantiles of the variational approximation by solid and dotted purple curves, respectively.
Observe that with 6, see left part of Figure I the variational approximation results in similar 95% pointwise credible bands and posterior mean as the true posterior, providing an accurate approximation. Also note that both the true and the variational posterior contain fo most of the points, indicating frequentist confidence validity of the set. At the same time, by taking a factor of two less inducing points, i.e. frequentist coverage, but suboptimally large posterior spread, see the second plot in Figure
3, the credible
will be overly large, resulting in good
The computations were carried out with 2,6 GHz Quad-Core Intel Core i7 processor, The compu- tation of the exact posterior mean and covariance kernel on a grid of 100 points took over half an hour (in CPU time), while the variational approximation was substantially faster, taking only 50.5 ms, resulting in a 3.68 10t times speed.
posterior mean
Puintwise
Variational mean (m=6)
Puintwise (variationa)
Variational mean (m=3)
(varistional)
Figure 1: True and variational posterior means and credible regions for Gaussian series prior (sine basis) the initial condition μ fo of the heat equation 12 for (left) or (right) inducing variables from method 100
A more extensive numerical analysis is available in the appendix, considering the application of our method to the settings of Sections and as well. We conduct these experiments several times and compare the average Mean Integrated Squared Error (MISE), see appendix A, and compute time for different choices m. We observe that all our examples, while increasing m results in longer computation, the MISE does not improve after threshold close the one presented in results.

---

Therefore, it is sufficient to include as many inducing variables as we considered in Corollary to obtain better performance. More than that would would only increase the computation complexity. In the Appendix, we also provide a literature review and some justifications of how relevant these problems are in practice
5 Discussion
We have extended the inducing variables variational Bayes method for linear inverse problems and derived asymptotic contraction rate guarantees for the corresponding variational posterior. Our theoretical results provide a guide for practitioners on how to tune the prior distribution and how many inducing variables to apply (in the spectral feature variational Bayes method) obtain minimax rate optimal recovery of the true functional parameter of interest. We have demonstrated the practical relevance of this guideline numerically on synthetic data and have shown that using less variables results in highly suboptimal recovery.
In our analysis have considered priors built the singular basis of the operator A. In principle our results can be extended to other priors as well, until the eigenbasis of the covariance operator is not too different from the basis of the operator A. This, however, would complicate the computation of the Kullhack -Leibler divergence between the variational family and the posterior, resulting in extra technical challenges. In this setting the empirical spectral features method seems practically more feasible, especially, if the eigenbasis of the covariance kernel is not known explicitly. In the literature several different types of inducing variable methods were proposed, considering other, practically more relevant approaches of interest. Furthermore, extension to other type of inverse problems is also feasible. For instance the deconvolution problem, when fo is convoluted with rectangular kernel, the eigenvalues given by the SVD are the product of polynomially decaying and oscillating part and the "average degree" of sol-posedsess does not match the lower and upper bounds [23]. Extension to non-linear inverse problem is highly relevant, as these problems tend to be computationally even more complex, but very challenging. One possible approach is to linearize the problem and take its variational approximation. Finally, it is of importance to derive frequentist coverage guarantees for VB credible sets. Our approach cannot directly be extended for this task. However, in the direct case, for some special choices of inducing variables, frequentist coverage guarantees were derived using kernel ridge regression techniques [38] (60). This result, although computationally somewhat cumbersome, in principle can be extended to the inverse setting as well. One last drawback of our results is that the priors we consider are non non-adaptative in the mildly ill-posed case. Minimax contraction rates are attainable only if the covariance eigenvalues smoothness B. While this is not an issue the severely ill-posed case in our results, we keep the study of adaptation for future works as it is much more involved question.
properly tuned, given the
Funding. Co-funded by the European Union (ERC, BigBayesUU, project number: (01041064) Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.
References
[1] AGAPIOU, S., LARSSON, S., AND STUART, A. M. Posterior contraction rates for the Bayesian approach to linear ill-posed inverse problems. Stochastic Process. Appl. 123, 10 (2013), 3828-3860.
[2] ALQUIER, P., AND RIDGWAY, J. Concentration of tempered posteriors and of their variational approxima- tions. Ann. Statist, 48, 3 (2020), 1475-1497.
[3] ARRIDGE, S.. MAASS, P., ÖKTEM. O., AND SCHÖNLIEB, C.-B. Solving inverse problems using data-driven models. Acta Numerica 28 (2019), 1-174.
[4] BERTERO, M. Regularization methods for linear inverse problems. In Inverse problems (Montecatini Terme, 1986), vol. 1225 of Lecture Notes in Math. Springer, Berlin, 1986, pp. 52-112.
[5] BISSANTZ, N., HOHAGE, T., MUNK, A., AND RUYMGAART, F. Convergence rates of general regu- larization methods for statistical inverse problems and applications. SIAM Numer. Anal. 45, 6 (2007), 2610-2636.
[6] BISSANTZ. N., AND HOLZMANN, H. Statistical inference for inverse problems. Inverse Problems 24, 3 (2008), 034009.
10

---

[7] BRIOL, F.-X.. OATES. C. J.. GIROLAMI. M., OSBORNE, M. A., AND SEIDINOVIC, D. Probabilistic integration. Statistical Science 34, (2019). 1-22.
[8] BURT, D., RASMUSSEN, C. E., AND VAN DER WILK, M. Rates of convergence for sparse variational Gaussian process regression. Proceedings of 36th International Conference Machine Learning (09-15 Jun 2019), K. Chaudhuri and R. Salakhutdinov, Eds., vol. 97 Proceedings of Machine Learning Research, PMLR, pp. 862-871.
[9] CANDÈS, E. J., AND WAKIN, M. B. An introduction to compressive sampling. IEEE signal processing magazine 25, 2 (2008), 21-30.
[10] CASTILLO, I. Lower bounds for posterior rates with Gaussian process priors. Electron. J. Stat. 2 (2008). 1281-1299.
[11] CAVALIER, L. Nonparametric statistical inverse problems. Inverse Problems 24, 3 (2008), 034004, 19.
[12] COTTER, S, L., DASHTI, M., AND STUART, A. M. Approximation of Bayesian inverse problems for PDEs. SIAM Numer. Anal. 48, (2010), 322-345.
[13] CSATó, L. Gaussian processesciterative sparse approximations. If you have discovered material in AURA which unlawful e.g. breaches copyright, (either yours or that a third party) or any other law, including but not limited to those relating to patent, trademark, confidentiality, data protection, obscenity, defamation, libel, then please read our Takedown Policy and contact the service immediately. March 2002.
[14] CSATÓ, L., AND OPPER, M. Sparse On-Line Gaussian Processes, Neural Computation 14. (03 2002). 641-668.
[15] FLORENS, J.-P., AND SIMONI, A. Regularizing priors for linear inverse problems. Econometric Theory 32, 1 (2016), 71-121.
[16] GARDNER, J., PLEISS, G., WENNERGER, K. Q., BINDEL, D., AND WILSON, G. Gpytorch: Blackbox matrix- matrix gaussian process inference with gpu acceleration. Advances in Neural Information Processing Systems (2018), S. Bengio. H. Wallach, H. Larochelle, K. Grauman, N. R. Garnett, Eds., vol. 31, Curran Associates, Inc.
and
[17] GHOSAL. S., AND VAN DER VAART, A. Convergence rates of posterior distributions for non-i.i.d. observations. Anm. Statist. 35, 1 (2007), 192 223.
[18] GOLDENSHLUGER, A., AND PEREVERZEV, S. V. On adaptive inverse estimation of linear functionals in Hilbert scales. Bernoulli 9, (2003), 783-807.
[19] GOLUBEV, G. K., AND KHAS' MINSKII, R. Z. statistical approach differential equations. Problemy Peredachi Informatsii 35, 2 (1999), 51-66.
inverse problems for partial
[20] GUGUSHVILL, S., VAN DER VAART, A., AND YAN, Bayesian linear inverse problems in regularity scales. Annales de l'Instiut Henri Poincaré, Probabilités et Statistiques 56, 3 (2020), 2081 2107.
[21] HALMOS, P. R. A Hilbert space problem book, vol. 19. Springer Science & Business Media, 2012.
[22] HEGDE, P., YILDIZ, Ç., LÄHDESMÄKI, H., KASKI, S., AND HEINONEN, M. Variational multiple shooting for bayesian odes with gaussian processes. In Uncertainty in Artificial Intelligence (2022), PMLR, pp. 790-799.
[23] JOHNSTONE, 1. M., AND RAIMONDO, M. Periodic boxcar deconvolution and Diophantine approximation. Ann. Statist. 32,5 (2004),
[24] JOHNSTONE, I. M., AND SILVERMAN, B. W. Speed of estimation in positron emission tomography and related inverse problems. Ann. Statist. 18, 1 (1990), 251-280.
[25] KNAPIK, B., SALOMOND, J.-B. A general approach posterior contraction in nonparametric inverse problems. Bernoulli 24, 3 (2018), 2091 2121.
[26] KNAPIK, SZABÓ. T., VAN DER VAART, A. W., AND VAN ZANTEN, H. Bayes procedures for adaptive inference in inverse problems for the white noise model. Probab. Theory Related Fields 164, 3-4 (2016), 771-813.
[27] KNAPIK, B. T., VAN DER VAART, A. W., AND VAN ZANTEN, J. H. Bayesian inverse problems with Gaussian priors, Ann. Statist. 39, (2011). 2657.

---

[28] KNAPIK, B. T., VAN DER VAART. A. W., AND VAN ZANTEN, J. Н. Bayesian recovery of the initial condition for the heat equation. Communications in Statisticx-heeory and Methods 42, 7 (2013), 1294- 1313.
[29] LAW, K., STUART, A., AND ZYGALAKIS, K. Data assimilation. Cham, Switzerland: Springer 214 (2015), 52.
[30] LÁZARO-GREDIEE GREDILA,, M., AND FIGUEIRAS- -VVDAL,, A. Inter-domain gaussian processes for sparse inference using inducing features. In Advances in Neural Information Processing Systems (2009). Y. Bengio, D. Schuurmans, J. Lafferty, C. Williams, and A, Culotta, Eds., vol. 22, Curran Associates, Inc.
[31] MAESTRINI, L.. AYKROYD, R. G., AND WAND. M. P. A variational inference framework for inverse problems. arXiv preprint S.KIv:2103.09900 (2021).
[32] MAIR, B. A., AND RUYMGAART, F. H. Statistical inverse estimation in hilbert scales. SIAM Journal on Applied Mathematics 56, 5 (1996), 1424-1444
[33] MATTHEWS, A. G. G., HENSMAN, J., TURNER, R., AND GHAHRAMANI, 7. On sparse variational divergence between stochastic processes. Proceedings of the 19th In-
methods and the kullback-leine ternational Conference on Artificial Intelligence and Statistics (Cadiz, Spain, 09-11 1 May 2016), A. Gretton and C. C. Robert, Eds., vol. 51 of Proceedings of Machine Learning Research, PMLR, pp. 231-239.
[34] MENG, R., AND YANG, X. Sparse gaussian processes for solving nonlinear pdes. arXiv preprint arXiv:2205.03760 (2022).
[35] NASHED, M. Z., AND WAHBA, G. Generalized inverses in reproducing kemel spaces: An approach to regularization of linear operator equations. SIAM Journal on Mathematical Analysis 5, 6 (1974), 974-987.
[36] NICKL, R. Bayesian non-linear statistical inverse problems. Lecture Notes ETH Zurich (2022).
[37]
ZANTEN, H. Contraction for sparse variational approximations in gaussian process regression. Journal of Machine Learning Research 23, 205 (2022). 1-26.
NIEMAN, D., SZABO, B.,
[38] NIEMAN, D., SZABO, B., AND VAN ZANTEN, H. Uncertainty quantification for sparse spectral variational approximations in gaussian process regression. arXiv preprint arXiv:2212.1103 (2022).
[39] PATI, D., BHATTACHARYA, A., AND CHENG, G. Optimal Bayesian estimation random covariate design with a rescaled Gaussian process prior. J. Mach. Learn. Res. 16 (2015), 2837-2851.
[40] PINSKI, F. J., SIMPSON, G., STUART, A. M., AND WEBER, H. Algorithms for kullback-leibler approximation of probability measures in infinite dimensions. SIAM Journal on Scientific Computing 37, 6 (2015), A2733-A2757.
[41] PINSKL, F. J., SIMPSON, G., STUAKT, A. M., AND WEBER, Н. Kullback-leibler approximation for probability measures on infinite dimensional spaces. SIAM Journal on Mathematical Analysis 47, 6 (2015). 4091-4122.
[42] POVALA, J., KAZLAUSKAITE, I., FEBRIANTO, E., CIRAK, F., AND GIROLAMI, M. Variational bayesian approximation of inverse problems using sparse precision matrices. Computer Methods in Applied Mechanics and Engineering 393 (2022), 114712.
[43] QUIÑONERO CANDELA, J., AND RASMUSSEN, C. E. A unifying view of sparse approximate gaussian process regression. J. Mach. Learn. Res. 6 (dec 2005), 1939-1959.
[44] RASMUSSEN, C. E., AND WILLIAMS, C. K. Gaussian processes for machine learning. Adaptive Computation and Machine Learning. MIT Press, Cambridge, MA, 2006.
[45] RAY, K. Bayesian inverse problems with non-conjugate priors. Electronic Journal of Statistics none (2013), 2516 2549.
[46] RAY, K., AND SZABÓ, B. Variational bayes for high-dimensional linear regression with sparse priors. Journal of the American Statistical Association 117, 539 (2022), 1270-1281.
[47] ROUSSEAU, J., AND SZABO, B. Asymptotic behaviour of the empirical Bayes posteriors associated to maximum marginal likelihood estimator. The Annuls of Statistics 45, 2 (2017), 833 865.
[48] ROUSSEAU, J., AND SZABO, B. Asymptotic frequentist coverage properties of Bayesian credible sets for sieve priors. The Annals of Statistics 48, 4 (2020), 2155 2179.
12

---

[49] SEEGER, M. Bayesian gaussian process models: Pac-bayesian generalisation ertor bounds and sparse approximations. Tech. rep., University of Edinburgh, 2003.
[50] SEEGER, M. Pac-hayesian generalisation eтror bounds for gaussian process classification. J. Mach. Learn. Res. null (mar 2003), 33-2699.
[51] SEEGER, M. W., WILLIAMS, C. K. I., AND LAWRENCE, N. D. Fast forward selection to speed up sparse gaussian process regression. In Proceedings the Ninth International Workshop on Artificial Intelligence and Statisties (03- Learning Research, PMLR, pp. 354-261. Reissued by PMLR on 01 April 2021.
2003), C. M. Bishop and B. Frey, Eds., vol. R4 of Proceedings of Machine
[52] SNELSON, E., AND GHAHRAMANI, Z. Sparse gaussian processes using pseudo-inpuls. In Advances in Neural Information Processing Systems (2005), Y. Weiss, B. Schölkopf, and J. Platt, Eds., vol. 18, MIT Press.
[53] SNIEKERS, S., AND VAN DER VAART, A. Adaptive Bayesian credible sets in regression with Gaussian process prior. Electronic Journal of Statistics 9, 2 (2015), 2475 2527.
[54] STUART, A. M. Inverse problems: a bayesian perspective. Acta numerica 19 (2010), 451-559.
[55] SZABÓ. B.. VAN DER VAAKT, A. W., AND VAN ZANTEN, H. Frequentist coverage adaptive nonpara- metric bayesian credible sets. Ann. Statist. 43, 4 (2015), 1391-1428. (with discussion).
[56] TIKHONOV, A. Solution incorrectly formulated problems and the regularization method. Soviet Math. 4 (1963), 1035-1038.
[57] TITSIAS, M. Variational learning of inducing variables sparse gaussian processes. Proceedings of the Twelth International Conference on Artificia Intelligence and Statistics (Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA, 16-18 Apr 2009). D. van Dyk and M. Welling, Eds., vol. 5 of Proceedings of Machine Learning Research, PMLR, pp. 567-574.
[58] TITSIAS. M. K. Variational model selection for sparse gaussian process regression. Report, University of Manchester, UK (2009).
[59] TRAVIS, L., AND RAY, K. Pointwise uncertainty quantification for sparse variational gaussian process regression with a brownian motion prior. NeurlPS 2023 (2023).
[60] VAKILI, S., SCARLETT, J., SHIU, D.-S., AND BERNACCHIA, A. Improved convergence rates for sparse approximation methods in kernel-based learning. In International Conference on Machine Learning (2022), PMLR, pp. 21960-21983.
[61] VAN DER VAART, A., AND VAN ZANTEN. H. Information rates of nonparametric Gaussian process methods. Mach. Learn. Res. 12 (2011), 2095-2119.
[62] VAART, A. W., AND VAN ZANTEN, Н. Rates of contraction of posterior distributions based on Gaussian process priors. Ann. Statist. 36, (2008), 1435-1463.
[63] WANG, Y., AND BLEI. D. M. Prequentistcoossitennnyy 527 (2019), 1147-1161.
consistency of variational Bayes. Amer. Statist. Assoc. 114.
[64] YANG, Y., PATI, D., AND BHATTACHARYA, A. o-variational inference with statistical guarantees. Ann. Statist, 48, (2020), 886-905.
[65] YANG, Y.. AND TOKDAR, S. T. Minimax-optinal nonparametric regression high dimensions. The Annals of Statistics 43, (2015), 652 674.
[66] ZHANG, F., AND GAO, C. Convergence rates of variational posterior distributions. The Annals of Statistics 48, 4 (2020), 2180 2207.
13


Reviews:
[{"review_id": "t2UgdBzIlJ", "reviewer": "Reviewer_rih9", "summary": "The authors propose an inducing-point variational Bayes approximation for linear inverse problems. They study the frequentist properties of the arising variational posterior approximation. Specifically, they prove contraction rates for the variational posterior, both in the mildly and severely ill posed case. Moreover, they illustrate the practical relevance of their results with a number of synthetic data experiments.", "questions": "In 2020, three papers were published in the Annals of Statistics providing a comprehensive and general study of the frequentist properties of variational posteriors. Two of those papers are cited by the authors as references [2] and [56]. Another one I would add to this list is “$\\alpha$-VARIATIONAL INFERENCE WITH STATISTICAL GUARANTEES” by Yang, Pati, Bhattacharya. I wonder if the authors would explain whether the results of those papers could be applied to the variational GP approximations for linear inverse problems, as considered in the authors’ submitted work. If they could, then I would like the authors to give reasons why they decided to study this problem directly rather than applying the theory provided previously in those three papers. If they couldn’t, I would like the authors to briefly describe why.\n\nI also have a few other questions and suggestions. In Definition 1 the authors say that a problem is mildly ill-posed if $\\kappa_j$ has exponential decay. Do the authors mean polynomial decay? Also in Theorem 1, only at the end of the statement do the authors define what it means for a posterior (or a variational posterior) to contract at a certain rate. I would include a definition of what it means to contract at a certain rate before stating Theorem 1. Moreover, I’m not sure if the notation $S_{++}^m$ used at the bottom of page 4 and the notation $\\Pi_u$ used at the top of page 5 are introduced anywhere in the paper before they are used.", "limitations": "A limitation of this work concerns its practical usability. For example, the authors’ results hold under the assumption that the model is well-specified. It is hard to expect that such an assumption will be satisfied in real-life applications and for real data sets. I believe that the authors should acknowledge this limitation and perhaps comment on possible extensions to instances of model misspecification.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well-written and its objectives and scope are clear. It provides an important contribution. As far as I am aware, frequentist properties of variational posteriors arising as inducing-point approximations in linear inverse problems have not been studied explicitly before. Such properties and contraction rates are useful in that they illustrate how the choice of the prior and the number of inducing points in the approximation impacts the statistical properties of the variational posterior. In that sense, they can provide useful guidance for the user.", "weaknesses": "I am not exactly sure if the comparison to the existing literature provided by the authors is substantial enough – I address this among my questions to the authors. ", "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  his/her interest in our work and the suggestions.\n\n\n## Literature review\n\nWhile the three papers on variational Bayes published in AoS in 2020 are indeed all cited in our\narticle [2,54,56], they consider direct problems, while our focus is on ill-posed inverse problems where the standard approaches can not be applied directly.\nThe challenge is that the operator \\mathcal{A} does not possess a continuous inverse hence one natural way to connect the direct and inverse problems \nis through modolus of continuity [22]. The modulus of continuity measure quantitatively the continuity of the operator (on a subset) of the function space, \nsee (26) for the precise definition in the Supplement. To have a meaningful bound and hence a good control on this quantity one has to consider appropriate\n metrics both in the \"inverse\" $L_2(\\mathcal{T},\\mu)$ and \"forward\" $L_2(\\mathcal{X},G)$ function spaces. A natural and (relatively) easy to handle choice is\n the $L_2$-norm in both cases. A further difficulty of using standard testing based methods (on which the general approaches are typically built on) is that\n they provide contraction rates in the direct, forward model with respect to the Hellinger distance or empirical $L_2$-norm, none of which are particularly \nappropriate to derive meaningful bounds for the modulus of continuity. Finally, we would like to note that for Gaussian processes (GP) a more specialised,\n powerful machinary was developed in the literature based on GP concentration functions [52]. We advocate that stating the conditions in terms of the \ncorresponding GP concentration inequality provides more insights than general approaches presented in the previous papers. This considerably simplifies\n the treatment of the true posterior distribution. Then to connect the true posterior with its variational approximation one can use a simple formula based\n on the duality of Kullback-Leibler divergence. In our case the most convenient formulation was derived in [40], but very similar techniques were used in\n the earlier papers [2,54,56] as well. \n\n\n## Typos and other minor comments:\n\nThank you for pointing out this typo, it is indeed mildly ill-posed, we have corrected it. Based on your suggestion, we have now included the definition\n of a posterior contraction rate before stating our main results. We have also updated the notation section with $S_{++}^m$ denoting the set of $m\\times m$\n positive-definite matrices and $\\Pi_\\mathbf{u}$ the induced prior distribution of the inducing variables $\\mathbf{u}$.\n\n## Misspecification\n\nWe agree that studying the robustness of our methodology in case of misspecification is an interesting and important problem. In case of misspecified models\n the posterior typically contracts around the Kullback-Leibler projection of the true distribution to the considered misspecified class. This has been well\n established in parametric and semi-parametric frameworks, but in non-parametric models, let alone inverse problems (to the best of our knowledge) no formal general\n result is available. Our focus here was on deriving the first contraction rate result for variational Bayes in case of ill-posed inverse problems,\n hence we leave this extension to future research."}, {"review_id": "qyXCdlJrXG", "reviewer": "Reviewer_TgVX", "summary": "This theoretical paper addresses the problem of inferring a function $f$ when\n\n$y(x_i) = \\[Af\\](x_i) + \\epsilon_i$ with $i=1, ...,n $\n\nand $A$ is suitably defined linear operator with $\\epsilon_i$ being i.i.d Gaussian noise. In words this is a linear inverse problem.\n\nA (standard) Bayesian approach is adopted by putting a Gaussian process prior on $f$ with a carefully chosen covariance function.\n\nThen an interdomain variational inducing point approach is adopted. All this so far is standard and pre-existing.\n\nThe contribution of the paper is to derive posterior contraction rates for this setup under certain assumptions and to apply it some specific examples with different choices of the the linear operator $A$. Specifically these are a volterra operator, a diffusion equation and a Radon transform. ", "questions": "I have no questions other than the ones mentioned above.", "limitations": "There are no negative impacts emerging from this work that I can see.\n\nLimitations are adequately discussed in the submission for the most part. The \"weaknesses\" section of this review has further discussion.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This is an important and challenging topic. The standard regression case has been studied by Nieman et al. 2022 which corresponds to taking the $A$ to be the identity in this paper. The extension is certainly of interest since the Bayesian approach to inverse problems is a thriving and important area. \n\nThe required extension is treated well by the authors. The particular challenges of introducing an operator $A$ are clearly described and handled. The results derived require significant technical heavy lifting. \n\nI was not able to go through the proof in detail but the results and approach look sensible and I suspect it to be correct. \n\nThe examples equations chosen were interesting and informative. \n\nThe numerical simulations, if brief, were nonetheless adequate in a paper of this type and compelling.\n\nLimitations were communicated where they arose. ", "weaknesses": "**The paper requires strong assumptions on the kernel of the Gaussian process.** To some extent these are forced by theory since a very mispecified covariance function would not contract. But, as the authors clearly concede, there is a more general set of covariance functions that would work than the ones considered.\n\n**The two methods of choosing the inducing variables are somewhat idealized.** For a mildly ill-posed problem and the empirical spectral feature method the complexity is $O(N^2 M)$. Since $M$ the number of inducing points typically needs to be an increasing function of $N$ the number of data points then this is not much of an improvement on the exact inference which is $O(N^3)$ in practice and $O(N^{2.3})$ in theory. The authors do mention this. Alternatively the population method of choosing inducing variables uses $O(NM^2)$ but requires oracle knowledge of the eigensystem of $A^*A$. Admittedly the authors do show some interesting examples where they are able to do this. The required $M$ is more benign for the strongly ill-posed system and correspondingly the heat equation example given is quite compelling.  \n\n**The inducing inputs are not optimized** This is known to make a big difference in the more practical part of the literature particularly in higher input dimensions. It is also one of the things that distinguishes the Titisias variational approach from the Deterministic Training Conditional (DTC) approach using the nomenclature of this paper https://www.jmlr.org/papers/volume6/quinonero-candela05a/quinonero-candela05a.pdf. \n\n**The literature review is not complete.** \n\nThere are two significant papers from Pinski et al. Specifically: \n\n*Kullback--Leibler Approximation for Probability Measures on Infinite Dimensional Spaces. Pinski, Simpson, Stuart and Weber. 2015*\n\n*Algorithms for Kullback--Leibler Approximation of Probability Measures in Infinite Dimensions. Pinski, Simpson, Stuart and Weber. 2015*\n\nThese papers are early, rigorous and are in a similar area to the submission. The method is distinct since they work principally with the inverse covariance operator. Importantly they also address the challenging non-linear case. The type of guarantee derived is also different since they do not derive contraction directly. Instead in the second they use their approximating for an MCMC proposal and gain asymptotic guarantees that way. Nonetheless these papers should be discussed.\n\nA historical note. The work of *Burt et al. 2019* is excellent and highly impactful but perhaps a bit too much is attributed to it here. It provides *rates* of convergence to the posterior but there was some work on convergence before this. This is largely because convergence follows straightforwardly from the variational formulation. *Titsias 2009* (the longer version) Section 2.3 has some discussion of improvement in terms of KL. This was formulated in terms similar to those used in the present submission by *Matthews et al. AISTATS 2015 On Sparse variational methods and the Kullback-Leibler divergence between stochastic processes.* The AISTATS paper mentioned also references several earlier foundational papers and the present submission would benefit from mentioning them too. For instance the early work authored by *Opper*, *Seeger* and *Csato* amongst others. None of these existing papers removes the main contribution of the current work. \n\n**Smaller points/ typos:**\n\nFor the mildly ill-posed case described in Definition 1 shouldn't the decay be polymonial rather than exponential? \n\nOn line 346 I believe the units need to be milliseconds for the description to make sense. So I think the typo is that \"50.5 m\" should be \"50.5 ms\". ", "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 carefully reading our paper and for the helpful comments.\n\n## Assumptions on the GP kernel\n\nWe thank the referee for this excellent comment. It is indeed standard in the literature to model covariance functions in accordance to the linear\n operator $\\mathcal{A}$. More general priors were also considered e.g. [1,39] but the covariance structure in all these cases were closely\n related to the SVD bases of the operator $\\mathcal{A}$. For example, wavelet priors can provide a feasible alternative in several examples [39]. \nA general understanding, however, is still missing even for the true posterior, let alone its variational approximation. This poses an interesting,\n but very challenging and technical research topic, going beyond the scope of our paper.\n\n\n## Idealized choice of inducing variables\n\nWe agree that in the mildly ill-posed case, the improvement from $\\mathcal{O}(N^3)$ to $\\mathcal{O}(NM^2)$ in the population and to $\\mathcal{O}(N^2M)$ in\n the empirical spectral features method (where the SVD basis of $\\mathcal{A}$ is not available) is somewhat limited. However, we argue, that even this\n relatively mild improvement allows us to analyse substantially bigger dataset in practice (and speeds up the computations as demonstrated in our numerical analysis).\nWe completely agree with the referee, that in the case of severely ill-posed problems the improvement is substantial, allowing basically linear computational complexity.\n\n## Inducing points are not optimised\n\nYes, we completely agree, that inducing points methods are highly important in practice. However, their analysis is considerably more challenging \nand to the best of our knowledge no posterior contraction rates were derived even in case of the direct problem. Therefore we leave this important \ntopic to future work and thank again for the suggestion.\n\n\n## Literature review and typos\n\nWe are particularly grateful to the reviewer for the two references from Pinsky et al. which are indeed closely related to our problem.\n While their objectives are different from ours, we now cite them as suggested by the reviewer, to make the literature review more complete.\n We have also included the reference to Seeger, Csato and Opper in our revised manuscript, thanks for pointing it out.\n\nThank you for reporting these typos, we have corrected them in the revised manuscript.\n"}, {"review_id": "mLD0IC05Lq", "reviewer": "Reviewer_VbCA", "summary": "This paper analyses posterior contraction rates in Variational GP for mildly and severely ill-posed linear inverse problems. The authors prove that under certain conditions variational posterior achieves same contraction rates as the true posterior. Moreover, they show that given enough inducing points (in spectral domain) the above-mentioned condition holds.\nAs specific examples, the minimal number of inducing points needed for optimal contraction and contraction rate is derived for two mildly ill-posed problems (Volterra operator, Radon transformations) and one severely ill-posed problem (heat equation).", "questions": "Suggestions:\n1. Increase readability of figure 1. The legend way too small.", "limitations": "The limitations are adequately addressed.", "rating": 5, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "Analysis of posterior contraction rates improves our theoretical understanding of VGP for linear inverse problems.\n\nThe presentation is well structured.\n\nThe contribution is clear and well argued for.", "weaknesses": "The experimental part in the main manuscript is very limited, with most experiments taken out to the appendix.  \n\nI believe the paper is more suited for a journal format. It is very technical and requires a thorough review which is hardly possible in a conference format.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your interest in our paper. We acknowledge the brievity of the experimental part in the main paper due to the space limitations.\n Please note that our work besides proposing a new methodological approach also provides theoretical guarantees in general settings and considers\n three specific examples. This theoretical part is arguably the main contribution of our work. Therefore, we have decided to move most of the simulations\n to the supplement in order to comply with the paper format. Nevertheless, we hope that the numerical results (although mostly moved to the supplement)\n carry the intended message and illustrate our theoretical findings well.\n\nWe agree with the reviewer that our paper is well-suited for a (good) Stats/ML journal. However, we have decided to submit our work to NeurIPS, because\n this is one of the most prestigous machine learning conference receiving a lot of visibility from a broad audience with different backgrounds. We note\n that variational Bayes and inverse problems are highly relevant for various community, including Machine learning, statistics and numerical analysis.\n  Therefore, in view of the multidisciplinary aspect of our work, we hope it would foster the links between these topics, and we deemed NeurIPS the best\n venue to do so. As for the theoretical aspect of our work, we argue that theory is necessary when we seek to use ML techniques in a responsible way,\n as it provides a justification for their use in practice and improves the explainability of the approach. Please also note that it is not uncommon\n to publish theoretical and highly technical papers in NeurIPS and we made a big effort in rigorously deriving and presenting our proofs.\n\nThank you for pointing out the readability issue with Figure 1, we have addressed it in our revision.\n"}, {"review_id": "jSbqKDK3x4", "reviewer": "Reviewer_KBMK", "summary": "This paper presents a theoretical analysis of the posterior contraction rates for the variational posterior in the context of linear inverse problems. Specifically, the variational posterior is obtained using the sparse Gaussian processes approach, with inducing variables learned through KL divergence minimization. The authors investigate three types of operators: the Volterra operator, the heat equation, and the Radon transform, each exhibiting varying degrees of ill-posedness. The numerical analysis focuses on synthetic toy problems, and the authors demonstrate that in this case, even a small number of inducing points can yield highly accurate results, achieving significant speed-up.", "questions": "Figure 1 shows that using m=3 is insufficient, while m=6 yields highly accurate results for the studied case, while this is a single case not quite informative. I am wondering whether it is possible to establish a \"phase-transition curve\" that can inform practitioners about the required number of inducing points for various types of inverse problems and different numbers of measurements. \n\nIf the results are not understood correctly, it might leads to poor choice on the number of inducing points,  and over-or-under estimate of the posterior uncertainty. This could be harmful in mission-critical applications such as medical imaging. How to prevent this happen? \n\nCould you please further explain why the posterior uncertainty being over-estimated while in traditional variational approximation the uncertainty usually got underestimated? How well is the uncertainty calibrated?", "limitations": "The authors do have discussed to what extend the presented results can be extended to nonlinear cases and other types of inverse problems.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "[Originality] I am not up-to-date with theoretical results with frequentist-Bayes analysis, and how far the presented results are based on or deviates from the monograph Nickl, Richard, 2023.  \n\nNickl, Richard. \"Bayesian non-linear statistical inverse problems.\" (2023).\n\n[Quality] The theoretical results seems sound. \n\n[Clarity] The paper is well-written; however, the use of jargon limits its accessibility to a narrower audience.\n\n[Significance] The inducing-point of Gaussian processes is widely used, so the theoretical results concerned with it should be of interest to the community working on inverse problems.", "weaknesses": "[Quality] The experimental evaluation is limited to a synthetic toy example only.  \n\n[Significance] The authors motivate the variational Bayes approach for computational reasons. However, recent developments of scalable GP methods can reduce the complexity of exact GP from O(n^3) to O(n^2) and effectively leverage the modern computing hardware such as GPU (Gardner, Jacob, et al.  2018). In this regard, one can afford many more inducing points in practice. \n\nGardner, Jacob, et al. \"Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration.\" *Advances in neural information processing systems* 31 (2018).\n\nThe paper is a theory paper, but its practical significance appears to be lacking. It would greatly enhance the paper's value if the authors could at least mention some real-world applications in which the studied inverse problem commonly arises. For instance, discussing the application of the Radon transform in computed tomography (CT) could serve as a relevant and illustrative example.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the referee for his/her interest in the paper and the constructive comments and suggestions.\n\n## Theoretical literature on frequentist-Bayes analysis\n\nThe monograph from Richard Nickl is multiscale analysis applied to non-linear inverse problems. This is a completely differnt techniques basically deriving simultaneous semi-parametric Bernstein-von Mises theorems leading to a Bernstein-von Mises type of result on a weak function space. For the investigated linear-inverse problems it is substantially simpler to follow our strategy deriving contraction rates in the direct problem, then use modulus of continuity and connect the true posterior with its variational approximation.\n\n\n## Experiments\n\nWhile our numerical analsys indeed considers only synthetic datasets, we argue that the selected examples are highly relevant and often form the baselines\n in the literature. For example, the heat equation is often consdiered as the starting example in the PDE literature and for instance the Black-Scholes PDE\n can be converted to the heat equation as well. In addition, the Radon transform is a mathematical technique with various applications, particularly in the \nfield of medical imaging and image analysis. It is used to analyze and transform data from the spatial domain to the Radon domain, which provides a different\n perspective on the data that can be useful for specific tasks. Inverting the Radon transform has found a lot of applications where lower-dimensional integrals\n of the inside of an object are more readily available than the object of interest itself. For instance, in Computed Tomography (CT) Imaging, X-ray measurements\n are taken from different angles around a patient, and the Radon transform is used to reconstruct a cross-sectional image (slice) of the patient's body.\n This helps doctors visualize internal structures and diagnose various medical conditions. Medical Single Photon Emission Computed Tomography (SPECT) is a nuclear\n medicine imaging technique that uses gamma-ray detectors to generate 3D images of the distribution of radioactive tracers within a patient's body, applying the\n Radon transform in the image reconstruction process. [*Barrett*, 1984,  *Natterer and Wübbeling*, 2001, *Ramlau and Scherzer*, 2019, *Ambartsoumian and Quinto*, 2020]. \nFurther applications include Seismic Imaging, Geophysical imaging, Radar and Sonar Imaging, Particle Tracking and Material Science and Crystallography, which we \nall discuss in details (with references) in the supplement of our manuscript.\n\n## Recent developments on scalable GPs\n\nWe sincerely thank the reviewer for bringing the paper on GPytorch to our knowledge. While their method improves the scalability of both exact and sparse approximate\n computations of posterior GPs, this should not encourage us to include many more inducing variables in our approach. Indeed, one of the main results of our paper is \na sufficient lower bound on the number of variables we need.  Adding more won't improve substantially the quality of the inference as demonstrated by our theory and numerical analysis.\nNevertheless, the combination of scalable approaches will allow us to tackle larger data sets and/or less regular models more efficiently. Please also note that our approach \nprovides substantial speed up in case of severely ill-posed inverse problems (where almost linear computational complexity applies) and for highly irregular functions.\n\n## Phase-transition curve\n\nWe thank the reviewer for this nice idea. In our general response we have addressed this comment in details.\n\n\n## Uncertainty quantification\n\nWe thank the reviewer for this really interesting question. The observation that usual variational methods tend to under-estimate the posterior uncertainty\n is based on the study of mean-field variational method or multimodal posterior distributions estimated by a unimodal one, see e.g. \"Bishop: Pattern recognition and machine learning\". However,\nour case does not fall into these categories. The inducing variable GP approximation has two components. Conditioned on the inducing points its law is coming\n from the GP prior. Then the inducing variables are fitted to best approximate the posterior (in KL divergence). Therefore, one can view this GP approximation as a\n combination of the prior and the posterior. Since the posterior is tunned to have the correct variance, while the prior has variance of O(1)\n (it is not incluenced by the data), this results in a conservative (possibly overly large) approximation of the uncertainty in contrast to the common belief originating from the mean-field approximation.\n"}, {"review_id": "MCA1xoN1Bw", "reviewer": "Reviewer_Sbdt", "summary": "This paper analyzes posterior contraction rates for the inducing variables variational Bayes approximation for Gaussian process. In particular, mildly and ill-posed inverse problem settings, which are characterized by profiles of eigenvalues of measurement operators, are considered. For the both cases, the posterior contract rates are assessed theoretically. The usefulness of the results are shown for applications to three example problems. ", "questions": "What happens if the target function does not belong to the assumed functional space?", "limitations": "Yes. ", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The posterior contract rates are assessed analytically for the inducing variable variational Bayes approximation of Gaussian process. This is useful to determine the optimal number of inducing variables in the approximation. ", "weaknesses": "The analysis relies on the assumption that the functional space where the target function belongs is known. However, there are many practical situations where such an assumption does not hold. The analysis is not applicable for such situations. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the summary and the question on adaptivity, which is of high importance in practice.\n As it is common in the literature, our first article focuses on the non-adaptive case when the function space \nis known and derived guarantees in this setting. The next, natural, follow up question is to address adaptation \nto the unknown function class. Please note that for the variational Bayes procedure adaptation to unknown function classes\n has been considered only for specific examples in the literature focusing mainly on mean-field VB methods. These results\n can not be directly applied to the inducing variable Gaussian Process  framework (where the variational class is not mean-field)\n even in the direct case, let alone in the more complicated ill-posed inverse problem setting. Since this is a very challenging\n question going beyond the scope of our already quite technical and long paper, we have left it to future work. "}]

Author rebuttal:
# General rebuttal:

We would like to thank the reviewers for carefully reading our manuscript, for acknowledging the importance of the topic
 ("important and challenging topic" - TgVX) and that the examples are " interesting and informative" (TgVX).
We are happy to hear that they found the paper "clear and well argued" (VbCA), "well-written" (KBMK, rih9), and that it
 "provides an important contribution" (rih9). We would also like to thank the reviewers for their constratice comments
 and suggestions which lead to a further improved version of the manuscript.

We highlight here the main changes. Below, we address the comments of each referee point-by-point and provide a more detailed answer to their questions.
- We have corrected the typos and provided the few missing definitions/concepts. We also discuss the strength and limitations of our results in more details.
- We have provided an extended description of the pracical application of the considered examples following the suggestion of KBMK.
- Following the suggestion of KBMK we have compared numerically the empirical and theoretical "phase-transition curves". 
We have conducted a numerical experiment for the Voltera operator, please find the figure attached. We computed the (logarithm
of the) fraction of the mean integrated squared error (MISE) corresponding to the true and variational posterior means (we used 20 repetitions to
 empirically approximate these quantities). We have considered n ranging from 100 to 10000 and m from 1 to 17. We have also plotted the phase transition
 curve (white line) coming from our theoretical analysis (see the PDF file attached). One can note that the theoretical curve closely resembles the curve where the phase transition occurs
in the empirical study. We believe that this provides a further evidence that our theoretical results indeed turn out to be a useful practical guide.


Meta-review: The reviewers agrees that the paper provides new theoretical insights into the tuning of sparse variational Gaussian process models for solving linear inverse problems, which arise frequently in scientific applications.
Decision: Accept (poster)