Bayesian Learning via Q-Exponential Process
Authors: Shuyi Li
Michael O'Connor
Shiwei Lan
Venue: NeurIPS 2023 poster (2023)
Keywords: Functional Regularization
Besov Process
$Q$-Exponential Distribution
Elliptic Contour Distribution
TL;DR: We propose a novel stochastic process named $Q$-exponential process that generalizes Gaussian process with shaper regularization and is more suitable for modeling objects with abrupt changes and sharp contrast such as images with edges.

Abstract: Regularization is one of the most fundamental topics in optimization, statistics and machine learning. To get sparsity in estimating a parameter $u\in\mathbb{R}^d$, an $\ell_q$ penalty term, $\Vert u\Vert_q$, is usually added to the objective function. What is the probabilistic distribution corresponding to such $\ell_q$ penalty? What is the \emph{correct} stochastic process corresponding to $\Vert u\Vert_q$ when we model functions $u\in L^q$? This is important for statistically modeling high-dimensional objects such as images, with penalty to preserve certainty properties, e.g. edges in the image.
In this work, we generalize the $q$-exponential distribution (with density proportional to) $\exp{(- \frac{1}{2}|u|^q)}$ to a stochastic process named \emph{$Q$-exponential (Q-EP) process} that corresponds to the $L_q$ regularization of functions. The key step is to specify consistent multivariate $q$-exponential distributions by choosing from a large family of elliptic contour distributions. The work is closely related to Besov process which is usually defined in terms of series. Q-EP can be regarded as a definition of Besov process with explicit probabilistic formulation, direct control on the correlation strength, and tractable prediction formula. From the Bayesian perspective, Q-EP provides a flexible prior on functions with sharper penalty ($q<2$) than the commonly used Gaussian process (GP, $q=2$).
We compare GP, Besov and Q-EP in modeling functional data, reconstructing images and solving inverse problems and demonstrate the advantage of our proposed methodology.

Bayesian Learning via Q-Exponential Process
Shuyi Michael O'Connor
Shiwei Lan
School of Mathematical & Statistical Sciences Arizona State University, Tempe, AZ 85287
Abstract
Regularization is one of the most fundamental topics in optimization, statisties and machine learning. To get sparsity in estimating parameter и E R", an bg penalty term, ||al|lg is usually added to the objective function. What the probabilistic distribution corresponding to such by penalty? What is the correct stochastic process corresponding to ||u[ly when we model functions u 19? This is important for statistically modeling high-dimensional objects such as images, with penalty to preserve certain properties, e.g. edges in the image. this work, we generalize the g-exponential distribution (with density proportional to) exp(-} to stochastic process named q g-exponential (Q-EP) process that corresponds the by regularization of functions. The key step to specify consistent multivariate q-exponential distributions by choosing from a large family of elliptic contour distributions. The work closely related Besov process which is usually defined in terms of series. Q-EP can be regarded definition of Besov process with explicit probabilistic formulation, direct control on the correlation strength, and tractable prediction formula. From the Bayesian perspective, Q- EP provides a flexible prior on functions with sharper penalty (q < 2) than the commonly used Gaussian process (GP, q = 2). We compare GP, Besov and Q-EP in modeling functional data, reconstructing images and solving inverse problems and demonstrate the advantage of our proposed methodology.
INTRODUCTION
Regularization on function spaces is one of the fundamental questions in statistics and machine learning. High-dimensional objects such as images can be viewed as discretized functions defined on 2d or 3d domains. Statistical models for these objects on function spaces demand regularization to induce sparsity, prevent over-fitting, produce meaningful reconstruction, etc. Gaussian process [GP 38 [2]] has been widely used as an L2 penalty (negative log-density as quadratic form) or prior on the function space. Despite the flexibility, sometimes random candidate functions drawn from GP are over-smooth for modeling certain objects such as images with sharp edges. To address this issue, researchers have proposed class of 1, penalty based priors including Laplace random field [37 34
and Besov process [31 25 [6] They have been extensively applied in spatial modeling 377, signal processing (28). imaging analysis
and inverse problems [i [5]. Figure I demonstratts an application of nonparametric regression models on functions endowed with GP, Besov and our proposed g-exponential process (Q-EP) priors respectively reconstruct blurry image a satellite, Q-EP model generates the best reconstruction, indicating its advantage over GP in modeling objects with abrupt changes sharp contrast such "edges" in image.
44 34
For these high- high-dimensional (refer to its discretization) inhomogeneous objects on d' domains D C Igd" particularly 2d images with sharp edges (d* 2), one can model them random function from
"lann @ asu
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Truth
Observadion
Gaussien
Besov
q Exporential
0.
0.2
03
60
Figure I: Image of satellite: true image, blurred observation, and reconstructions by GP, Besov and Q-EP models with relative errors 75. 19%, 21. 21.94% and 20.35% respectively.
Besov process represented by the following series for a given orthonormal basis (0r) in (D)
31 15
u: - R, u(x) E de)))()),
sd ~ (.) exp -xp(-211119
where q 2 and Y/ = x - with (inverse) variance K >> and smoothness > 0. When q q= 2 and { is chosen to be Fourier basis, this reduces to GP [16] but Besov is often used with = and wavelet basis [31] to provide "edge-preserving function candidates suitable for image analysis. Historically, (32) discovered that the total variation prior degenerates to GP prior the discretization properties in high dimensional applications. Therefore, [31] proposed the Besov prior defined as in ① and proved its discretization-inaarii property. Though straightforward, such series definition lacks a direct way to specify the correlation structure as GP does through the covariance function. What is more, once the basis {Ør} is chosen, there is no natural way to make prediction with Besov process.
mesh becomes denser and thus loses the
We propose a novel stochastic process named g-exponential process (Q-EP) to address these issues. We start with the q-exponential distribution (,)+) and generalize it to a multivariate distribution (from family of clliptic contour distributions) that is consistent to marginalization. Such consistency requires the joint distribution and the marginalized one (by any subset of components) to have the same format of density (See Section 3). We then generalize such multivariate q-exponential distribution to the process Q-EP and establish its connection and contrast to the Besov process. Note, if we view the negative log-density of the proposed distribution and process, Q-EP would impose an Ly regularization on the function space, similarly as L2 regularization given by GP whose negative log-density is quadratic form of the input variable (See Remark 2
Connection to existing works The proposed Q-EP process is related to the student-r process (TP) 41 as alternatives to GP. TP generalizes multivariate t-distribution (MVT) and is derived as scale mixture of GP. Both TP and Q-EP can be viewed as a special case of the elliptical process [ which gives the condition general elliptic distributions that can be generalized to valid stochastic process. Both papers focus on extending GP to robust models for heavier tail data, while our proposed work innovates a new Bayesian learning method on function spaces through the regularization parameter 4 (See Figure for its effect on regularization when it varies), as is usually done in the optimization. Both our proposed Q-EP and [5] are inspired by Kano's consistency result (27]. however the later focuses on a completely different process named squeezebox process. Our work on Q-EP makes multi-fold contributions the learning of functional data in statistics and machine learning:
1. We propose novel stochastic process Q-EP corresponding to the Ly regularization on function spaces.
2. For the first time we define/derive Besov process probabilistically as Q-EP with direct ways to configure correlation and to make prediction.
3. We provide flexible Bayesian inference methods based on the Markov Chain Monte Carlo (MCMC) algorithms using a white-noie noise representation for Q-EP prior models.
The rest of the paper is organized as follows. Section 2 introduces the g-exponential distribution and its multivariate generalizations. We propose the Q-EP with details in Section 3 and introduce as nonparametric prior for modeling functional data. In Section 4 we demonstrate the advantage of Q-EP over GP and Besov in time series modeling, image reconstruction, and Bayesian inverse problems (Appendix C.4) Finally we discuss some future directions in Section 5

---

2 THE Q-EXPONENTIAL DISTRIBUTION AND ITS MULTIVARIATE GENERALIZATIONS
Let start with the g-exponential distribution for scalar random variable R. It is named in [[5] and defined with the following density not in an exact form (as a probability density normalized to 1):
Ty (u) exp (
lu19).
This g-exponential distribution 2 is actually special case of the following exponential power (EP)
distribution EP(H,, 0,9) with μ
Mulu,,,9) (11/2) exp
0
where I denotes the gamma function. Note the parameter q > 0 in 33 controls the tail behavior of the distribution: the smaller the heavier tail and vice versa. This distribution also includes many commonly used ones such as the normal distribution N (u, σ²) for = 2 and the Laplace distribution L( with o 2-1/4b when 1.
How generalize it to multivariate distribution and further to stochastic process? Gomez 233 provided one possibility of a multivariate EP distribution, denoted as EP C, with the following density:
(T 2 2- x x IcI exp { [-1 (1 p(u C.9) 2)())
[ (u- 4) c ))}}}
(4)
When q reduces to the familiar multivariate normal (MVN) distribution N2 C).
Unfortunately, unlike MVN being the foundation of GP, the Gomez's 's EP distribution C,.
fails generalize to valid stochastic process because it does not satisfy the marginalization consistency as MVN does (See Section for more details). It turns out we need to seek candidates in even larger family of elliptic (contour) distributions [H, C.g):
Definition 2.1 (Elliptic distribution). A multivariate elliptic distribution ECa (5, C,g) has the follow- ing density [26]
p(u) - k. C -g-r))),, r(u) - (u (u- c-
(5)
where Ka > 0 is the normalizing constant and g(.), a one-dimensional real-valued function indepen- dent of d and ka, is named density generating function 79/.
Every elliptic (contour) distributed random vector (M, C,g) has a stochastic representation mainly due to Schoenberg [40 26),, as stated the following theorem.
Theorem 2.1. u ECA(M, C,g) if and only if
+RS
(6)
where ~~ (ddl(( uniformly distributed the unit-sphere gd+1, 1, L is the Cholesky factor of C x2
such that C = LLT, R S and R² d r(u) f(r)
- r(4)
kar - g(r).
The Gomez's EP distribution (M C,q) special elliptic distribution EC, C,g) with g(r) cxp(-tr and R ~ I(a q 3 {) [23] Not all elliptical distributions can be used to create a valid process [5) In the following, we will carefully choose the density generator in (C, C.g) define a consistent multivariate eexxoneeniaal distribution generalizable process appropriately.
3 THE Q-EXPONENTIAL PROCESS
To generalize (H, C.g) to valid stochastic process, we need to choose proper g such that the resulting distribution satisfies two conditions of Kolmogorov extension theorem [35]:

---

Exponential Power
Exxpenentall
-10
-10
-10
Figure 2: Inconsistent (Gomez's s) EP distribution EP,(M, C.q) (left) vs. consistent Q-exponential distribution q (H, C) (right). Both can be sampled using 66) with R ~ ~ d B = and R9 T (ax 2 ) respectively. Note there is significant discrepancy between the marginaliza- tion EP3 (n. C, and (M, C.q). However, the marginalization of q q-ED3 C) coincides with ( C). Empirical densities are estimated based on 10000 samples (shown as dots) for = 1.
Theorem 3.1 (Kolmogorov's Extension). For all f1....,4. ET, KEN let Vnr- be probability measures on Tidk satisfying
(K1) Mouye""oo Fots) (Fi x ... Fi)) = Vn (Fo x ... x Fg prmuatttoo
permutations 6
S(k)
(K2) : d. (Fi X Fk) V1 Jedet. ferm (FI x... x Fx X forallm
(7)
Then there exists a probability space (2, and a stochastic process {x, on Q X, : + 12" such that
Vn (Fi x ... x Fi) = PiX, E F1, X Ex
(8)
for all 4 T, k E N and all Borel sets Fr $. (K1) and (K2) are referred to as exchangeability and consistency conditions respectively.
As pointed out by Kano 27 the elliptic distribution EC, C,8) in the format of Gomez's 's EP dis- tribution with g(r) ( does not satisfy the consistency condition [also c.f. Proposition 5.1 of 23 Figure 2 (left panel) also illustrates such inconsistency numerically. However, Kano's consistency theorem [27] suggests different viable choice of make valid generalization of EC2 C. to a stochastic process [5]:
Theorem 3.2 (Kano's Consistency). An elliptic distribution consistent and only its density generator function, s(.), has the following form
g(r) 8(r)= 100=[['0
((x)) exp
} p(s)ds
15)
(9)
where p(s) is a strictly positive mixing distribution independent of d and p(s = 0) = 0.
3.1 Consistent Multivariate Q-exponential Distribution
In the above theorem 3.2 we choose p(s) = - then we have g(r) - ((-1)) exp {
9
2
which leads the following consistent multivariate g-exponential distribution q- (H C). Definition 3.1. A multivariate g-exponential distribution, denoted q DDDz(( C), has the following density
/ r(u) = (u (n- (u-uu)-c-00 (u- (10)
p(u SuluuC.q C.q) -2 ((x)) c - 1) exp
2

---

Remark 1. When q= q-EDa (M. C) reduces to MVN (a C). d= 1, if = have the density for as p(u) o ju| e differing from the original un-normalized density ag in (2) by term ju]-11 This is needed for the consistency of process generalization. Numerically, has the similar "edge-preserving" property as the Besov prior.
1, then
Remark 2. If taken negative logarithm. the density of q q-EDa in (10) yields a quantity dominated by some weighted Ln norm of u - H. i.e. 1/²² = 2 ||u μ llc. From the optimization perspective, q q-EDa. when used as a prior, imposes Ly regularization in obtaining the maximum a posterior (MAP).
Regardless of the normalizing constant, our proposed multivariate q-exponential distribution q (D, C) differs from the Gomez's s EP distribution EP₂ (H. (H. C.q) by a boxed term -)-1) As stated the following theorem, q q-ED₂ satisfies the two conditions of Kolmogorov extension theorem thus is ready to generalize to a stochastic process (See the right panel of Figure 2 for the consistency).
Theorem 3.3. The multivariate g-xxpnnential distribution both exchangeable and consistent.
Proof. See Appendix A.I
Like student-r distribution 41 and other elliptic distributions [5).
show (See Appendix A.5
that q DDD is represented as a scale mixture of Gaussian distributions for 0 < q < (27
Numerically, thanks to choice of density generator g(r) - exp { 1 2
one show
that R x (as in Appendix A.4 thus R in Theorem 2.1 can be sampled q-root of 22 random variable, which completes the recipe for generating random vector u ~ q-ED,(0,, C) based on the stochastic representation (6). This is important for the Bayesian inference detailed in Section 3.3.1 Note the matrix C in the definition (IO) characterizes the covariance between the components, as shown in the following proposition.
Proposition 3.1. If u ~ q ( [ C). then we have
2 + +
E|u] M. Cov(u)
Cidi
asd o0
(11)
dr())
Proof. See Appendix A.4
3.2 @-exponential Process as Probabilistic Definition of Besov Process
To generalize u q (O, C) to stochastic process, we need to scale it to u² =d= u so that its covariance is asymptotically finite. If u ~ q Q-ED,(0, C). then we denote u² ~ - C) following a scaled g-exponential distribution. Let 6 : L - L be kernel operator in the trace class, i.e. having eigen-pairs (z
de (x) such that ( = de (x) || = for all E N and ir(6)= Σx=1 Àg < co, Now we are ready to define the geexponential process (Q-EP) with the scaled g-exponential distribution.
Definition 3.2 (Q-EP). A (centered) q-exponential process u(x) with kernel & in the trace class, q 9-8 (0, 6), is a collection of random variables such that any finite set, u= (u(x)),. (uci),).--xx))))))) u(xa)), follows scaled multivariate q-exponential distribution, i.e. u ~~0 Lo---EDja (O,
Note, the process is defined on the d' -ddimensional space depending on the applications (x E Rd ) while d refers to the dimension of discretized process (u Rd)). While Q-EP reduces to GP when 2, q is often adopted for imaging analysis as an "edge-preserving" prior. Illustrated in Figure 3 for selected q E (0, 2), smaller y leads sharper image reconstruction with varying 4 interpolating between different regularization effects. Both Besov and Q-EP are valid stochastic processes stemming from the g-exponential distribution ny. They are both designed generalize GP to have sharper regularization (through q) but Q-EP has advantages in I) the capability of specifying correlation structure directly and 2) the tractable prediction formula.
It follows from | immediately that the covariance of the Besov process u(-) at two points x,x' IR"::
Cov(u(z),u((')) Σ (x) (x(x))
(12)
5

---

Q-EP (c=0.5)
QPEP (q=1)
O-EP (q=1.5)
Gaussian (q=2)
0.4
09
0.2
0.6
1.0
64
0.6
30 0.0
Figure 3: Image of satellite: MAP estimates by Q-EP with varying 4 parameters.
Although the smoothness and correlation strength of Besov process can be configured by proper orthonormal basis { in (12). it is less straightforward than the kernel function working for GP. On the other hand, Q-EP has more freedom on the correlation structure through (TT) with flexible choices from large class of kernels including powered exponential, Matérn and periodic where can directly specify the correlation length.
While Q-EP can be viewed as probabilistic definition of Besov, the following theorem further establishes their connection in sharing equivalent series representations.
Theorem 3.4 (Karhunen-LLovve) If u(x)) ~ 8 (0, 8) with a trace class operator & having eigen-pairs {2e, [-)] such that (d)(x) (x) [ for all & N and X 2g < 00, then we have the following series representation for uxx):
u(x)=
Σ (x)e(x);
uy :=:= J ((x))=(xx) ind q EDD (0,
(13)
where Elus] = O and Cov(ug,ur) = 2ey with Dirac function 8n = if = and 0 otherwise.
Proof. See Appendix A.5
Remark If we factor / 'he out of ug, we have the following expansion for Q-EP more comparable
to for Besov:
(xx)
(x),
udd ~
ED(O, (.)
(14)
3.3 Bayesian Modeling with Q-exponential Process
Now let us consider the generic Bayesian regression model:
y= = u(x)+++ +E, E E L~10,=
(15)
~ (du)
where L(=;0, Σ) denotes some likelihood model with zero mean and covariance E, and the mean function u can be given a prior either Besov or Q-EP. Because the definition ( in terms of expanded series, there is no explicit formula for the posterior prediction using Besov prior. By contrast, a tractable formula exists for the posterior predictive distribution for (15) with Q-EP prior Ho q- 6 (0, when the likelihood happens to be L(=;0C) = E----DDO,,, C), stated in the following theorem.
Theorem 3.5 (Posterior Prediction). Given covariates = fxit and observations y = (yil following q-ED ED in the model (15) with q- 9 prior for the same q > 0, we have the following posterior predictive distribution for u(x.) at (a) new point(s) (o x,:
((xx))y,x.x q-D( C"), === C (C+E)
C= -C.
(16)
where C (x,x), C. xx,,x.), and C.. = 8 (x..x.).
Proof. See Appendix A.6
Remark 4. From Theorem 3.5 we know that Q-EP has same predictive mean as GP. But their
predictive covariances differ by a constant (3) r +
(asymptotically 1) based on Proposition 3.1 ri3)
When the likelihood is not Q-EP, egg. multinomial, such conjugacy absent. Then we refer the following sampling method for the posterior inference.

---

3.3.1 Inference by White-Noise MCMC
We [3] to consider the pushforward (T") of Gaussian white noise Vo for non-Gaussial measures Mo =0=7 70= Vo- More specifically, we construct measurable transformation 7 - that maps standard Gaussian random variables to g-exponential random variables. The transformation based on the stochastic representation ( is more straightforward than that for Besov based on series expansions proposed by 13
%a and S ~ R9 d |zz|2 for ~ (2 La). Therefore, we have the pushforward mapping (T) and its inverse (7'11
Recall the stochastic representation (6) of u ~ q S-----EDDOA ED, C): u= RLS with R Lnd We can rewrite _
I [z)) = Lz||z|
z == - L- u
(17)
Figure B.1 illustrates that sampling with the white-noise representation (77) is indistinguishable from
sampling by the stochastic representation (6). Then we can apply such white-noise representation to dimension-innepeedden MCMC algorithms including preconditioned Crank-Nicolson (pCN) [14].
infniite-dimensional Metropolis adjusted Langevin algorithm (so-MALA) [T), infiniteddimensionall Hamiltonian Monte Carlo (co-HMC) (7], and -ifiite-dimensionall manifold MALA (ce-mMALA) g and HMC (co-mHMC)
9
See Algorithm for an implementation on pCN, hence named white-noise
pCN (wn-pCN).
3.3.2 Hyper-parameter Tuning
in the covariance function of Q-EP, e.g. variance magnitude (o²) and correlation length (p), that require careful adjustment and fine tuning. If the data process y(x) and mean u(x) are both Q-EP with the same 4, then we can have the marginal likelihood (38] as another Q-EP (c.f. Theorem 3.5). In general, when there is such tractability, hyper-parameter tuning by optimizing the marginal likelihood wingiiiollkeehoood is unavailable. However, we could impose conjugate hyper-priors on some parameters even marginalize them to facilitate the inference of them (See Appendix A.7 7 for a proposition on such conditional conjugacy for the variance magnitude o²).
As in GP, there are
To tune the correlation length (p), we could impose a hyper-prior for and sample from p(p)u). We then alternate updating and hyper- prree-pprameters in Gibbs scheme. In general, one could also use Bayesian optimization methods [33 18 for the hyper-paramedt parameter tuning.
4 NUMERICAL EXPERIMENTS
In this section, we compare GP, Besov and Q-EP by modeling time series (temporal), reconstruct- ing images (spatial) from computed tomography and solving a (spatiotemporal) inverse problem (Appendix C.4). These numerical experiments demonstrate that our proposed Q-EP enables faster convergence in obtaining a better maximum posterior (MAP) estimate. What is more, white- noise MCMC based inference provides appropriate uncertainty quantification (UQ) (by the pos- terior standard deviation). More numerical results can be found in the supplementary materials which also contain some demonstration codes. All the computer codes are publicly available at https: thub. com/lanzzi th inkinn EXP
4.1 Time Series Modeling
We first consider two simulated time series, one with step jumps and the other with sharp turnings, whose true trajectories are as follows:
4)(r) (f) = 1, E [0, 1]: 0.5,
TE (1,11.)))
2. rE (1,5,2);
0, otherwise
(t) = 1.5t, I [0, 1]: 3.5 21, (1,11.))): 3r 3,-4, re (1.5,)]
0, otherwise
We generate the time series {x} by adding Gaussian noises to the true trajectories evaluated at N 200 evenly spaced points (1:) in [0,2], that is, y. = u. (4) + & ind N(0, of (1:),, i= I.....N. .N, J,T. Let ///||=/]| = 0.015 for E [0,2] and ØT/ ||T]| = 0.01 E (0, 0.07 0007if4; E (1,2). In addition, we also consider two real data sets of Tesla and Google stock prices in 2022. See Figures (and Figures C.2 for the true trajectories (blue lines) and realizations (orange dots) respectively.
7

---

BesOv
Sphaga vernial
2.4
20
3.3
80
65
100
13
10
20
20
15
(a) Time series with sharp turnings (model fitting).
(b) Time series with turnings (prediction).
Geurs
Besov
Experentel
Gezmann
odtpprental
250
130
292 2022-00 es
2022 mna anc-a
meziti 2aazet a9z o2
mar
200 aoz
2012 anz- mis-1.c mnzuanes
aoo ynu at
220 2002 22.22200 mibiiaa
(c) Tesla stock prices in 2022 (model fitting).
(d) Tesla stock prices in 2022 (prediction).
Figure (a)()) MAP estimates by GP (left), Besov (middle) and Q-EP (right) models. (b)(d) Predictions by GP (left) and Q-EP (right) models. Orange dots actual realizations (data points). Blue solid lines are true trajectories. Black ticks indicate the training data points. Red dashed lines are MAP estimates. Red dot-dashed lines are predictions with shaded region being credible bands.
We use the above likelihood and test three priors: GP, Besov and Q-EP. For Besov, we choose the Fourier basis (t) = v2, @e())= cos (T f E N (results with other wavelet bases including Haar, Shannon, Meyer and Mexican Hat are worse hence omitted). For both GP and Q-EP, we
adopt the Matérn kernel with v 2. σ² p 0.5 and 1: C(t,t') α² 21-v w Ky)),, w
Fivj
/2v([kr-rlll/p) In both Besov and Q-EP, set q =y= 1. Figures 4a and 4c (and Figures C.2a and C.2cb compare the MAP estimates (red dashed lines). We can see that Q-EP yields the best estimates closest to the true trajectories in the simulation and the best fit to the Tesla/Google stock prices. We also investigate the negative posterior densities and relative errors, || - // ||u |l. as functions of iterations in Figure C.I Though incomparable in the absolute values. the negative posterior densities indicate faster convergence in both GP and Q-EP models than in Besov model. The error reducing plots on the right panels of subplots in Figure J.T.Iiniccat that Q-EP prior model can achieve the smallest errors. Table I compares them in terms of root mean of squared error (RMSE) and log-likee liielikkohdd (LL.).
Table 1: Time series modeling: root mean of squared errors (RMSE) and log- log-likelihood (LL) values MAP estimates by GP, Besov and Q-EP prior models.
root mean squared errors (RMSE)
log- log-likeldood (LL)
Data Sets
GP
Besov
Q-EP
GP
Besov
Q-EP
simulation (jumps)
1,2702
2.1603
1.1083
-31.4582
-89.8549
-74.0590
simulation (turnings)
1.4270
2.4556
0.9987
-39,8234
-56.7874
-87.3124
Tesla stocks
180.3769
136.8769
51.2236
-488.4558
-281.3796
-39.4070
Google stocks
44.4236
39.4809
36.8686
-386.1546
-305.0058
-265.9790
Then we consider the prediction problem. In the simulations, the last 1/8 portion and every other of the last but 3/8 part of the data points are selected for testing. The models with GP and Q-EP priors are trained on the rest of the data, as indicated by short "ticks" in Figures 4b and 4d (and Figures C.26 and C.2dp. For the Tesla/google stocks, we select every other day the first half year, every days in the 3rd quarter and every 8 days in the last quarter for training and test on the rest. They pose challenges on both interpolation (among observations) and extrapolation (at no-observation region) tasks. As we can see in those figures, uncertainty grows as the data become scarce. Nevertheless, the Q-EP yields smaller errors than GP. Note, such prediction is not immediately available for models with Besov prior.
8

---

Truth
Observation
Geussian
Besov
@Epoorential
05
00
1--e2 -010
06
Figure 5: Shepp-Logar phantom: true image, observation (sinogram), and MAP estimates by GP, Besov and Q-EP models with relative errors 68. 10%, 70.27% and 40.87% respectively.
4.2 Computed Tomography Imaging
Computed tomography (CT) is medical imaging technique used to obtain detailed internal images of human body. CT scanners use a rotating X-ray tube and a row of detectors to measure X-ray attenuations by different tissues inside the body from different angles. Denote the true imaging as function u(x) on the square unit D = [0, 1]2 taking values as the pixels. The observed data, y, (a.k.a. sinogram) are results of Radon transformation (A) of the discretized ny sensors, contaminated by noise & [
field with ng angles and
y E 8 8~~ (0,00 y YER', A P
E 1R"
In general fg < d= so the linear inverse problem is under-determined Baysian approach could fill useful prior information (e.g. edges) in the sparse data.
We first consider the Shepp-Logan phantom, a standard test image created by Shepp and Logan in 42 to model a human head and to test image reconstruction algorithms. In this simulation, we create the true image u² for a resolution of n² 128 x 128 and project at ng = 90 angles with n, 100 equally spaced sensors. The generated sinogram is then added by noise with signal noise ratio SNR Au |||||||| = 100. The first two panels Figure show the truth and the observation.
Table 2: Posterior estimates of Shepp-Logan phantom by GP, Besov and Q-EP prior models: relative error, RLE = || [ ==1l//l " of MAP (û u²") and posterior mean (ú - u) respectively, log-likelihood (LL), PSNR, SSIM and HarPPSI. Numbers in the bracket are standard deviations obtained repeating the experiments for 10 times with different random seeds.
MAP
Posterior Mean
GP
Besov
Q-EP
GP
Besov
Q-EP
RLE
0.6810
0.7027
0.4087
0.4917(6.16c-7)
0.4894(3.53c-5)
0.4890(4.796
LL
-1.55c+6
-1.54c+6
-1.57c+5
-5.21c45(8.47)
-.80c+5(116.34)
-4.6c+5(307799
07.77
PSNR
15.5531
15.2806
19.9887
18.3826(1.09) (09c-5)
(8.4226(6.27c-4)
(8.44))(
SSIM
0.4028
0.3703
0.5967
0.5561(3.92e-)
0.5535(2.38e-4) (.2403)5.26e-4)
HaarPSI
0.0961
0.0870
0.3105
0.3126(1.52e-8)
0.3126(3.36e-4)
0.3122(3.06e-))
Note, the computation involving full sized (d x dd) kernel matrix C for GP and Q-EP is prohibitive. Therefore, we consider its Mercer's expansion (12) with Fourier basis for a truncation the first L = 2000 items. [ that while GP and Besov models reconstruct very blurry phantom images, the Q-EP prior model produces MAP estimate of the highest quality. For each of the three models, we also apply wn-pCN to generate 10000 posterior samples (after discarding 5000) and use them to reconstruct u (posterior mean or median) and quantify uncertainty (posterior standard deviation).
Sabler/summrrizz the errors relative MAP (u*) and posterior mean (77 respectively, -|@- u² I1 [M""]] (with ù being u² or w), log-lik ligelikeldoood (LL), and several quality metrics in imaging analysis including the peak signal-to-noise ratio (PSNR) 20). the structured similarity index (SSIM) [46), and the Haar wavelet-based perceptual similarity index (HaarPSI) [39]. Q-EP attains the lowest error and highest quality scores most cases. In Figure C.3 we compare the uncertainty by these models. seems that GP has uncertainty filed with more recognizable shape than the other two. However, the posterior standard deviation by GP is much smaller (about 1% of that with Q-EP) compared with the other two. Therefore, this raises a red flag that GP could be over-confident about less accurate estimate.

---

Truth
Observation
Geussian
Besov
Exporential
0.
00
04
80
C4
06
Truth
Observv jon
Gauasion
Deruov
q Expomential
5.3
0.2
03 00
50
15.40 669
04
54
02
24
66
Figure 6: CT of human head (upper) and torso (lower): true image, observation (sinogram), and MAP estimates by GP,, Besov and Q-EP models with relative errors 29.99%, 22.41 and 22.24% (for head) and 26.1 %,, 21.77% and 21.53% (for torso) respectively.
Finally, we apply these methods to CT scans of a human cadaver and torso from the Visible Human Project [ These images contain n² 512 x 512 pixels and the sinograms obtained with m 200 angles and n, = 512 sensors. The first two panels of each row Figure show a highly calibrated CT reconstruction (treated "truth") and the observed sinogram. The rest three panels illustrate that both Besov and Q-EP models outperform GP in reconstructions, as verified in the quantitative summaries in Table C.2
Figure C.4 indicates that GP tends to underestimate the uncertainty.
In these CT reconstruction examples, we observe larger discrepancy of performance between Besov and Q-EP in the low-dimensional data-sparse application Sehepp-Logaae phantom at resolution n = 128 x 128 with na = 90 angles and п, = 100 sensors) compared with the high-dimensional data-intensive applications (two human body CTs at resolution n² 512 512 with næ 200 angles and n, 512 sensors). This may be due to the truncation Mercer's 's kernel representation (12) and different rates of posterior contraction
2
We will explore them in another journal paper.
5 CONCLUSION
In this paper, we propose the g-exponential process (Q-EP) as a prior on L functions with a flexible parameter q > control the degree of regularization. Usually, q = adopted capture abrupt changes or sharp contrast in data such as edges in the image the Besov prior has recently gained popularity for. Compared with GP, Q-EP can impose sharper regularization through q. Compared with Besov, Q-EP enjoys the explicit formula with more control on the correlation structure as GP. The numerical experiments in time series modeling, image reconstruction and Bayesian inverse problems demonstrate our proposed Q-EP is superior in Bayesian functional data modeling.
In the numerical experiments of current work, manually grid-search for the optimal hyper- parameters. The reported results are not sensitive to some of these hyper-parameter such as the variance magnitude (σ²) and the correlation length (p) but may change drastically to others like the regularity parameter (V) and the smoothness parameter (s). In future, we will incorporate hyper-priors for some of those parameters and adopt hierarchical scheme overcome such shortcoming. We plan to study the properties such as regularity of function draws of Q-EP and the posterior contraction, and compare the contraction rates among GP, Besov and Q-EP priors (2] 22 2]. Future work will
also consider operator based kernels such as graph Laplacian [6
29
Acknowledgments and Disclosure of Funding
SL is supported by NSF grant DMS-211 2134256.
10

---

References
[1] The visible human project.
[2] Sergios Agapiou, Masoumch Dashti, and Tapio Helin, Rates of contraction of posterior distributions based p-exponential priors. Bernoulli, 2763):1616 1642, 2021.
[3] D. F. Andrews and C. Mallows. Scale mixtures of normal distributions. Journal the Royal Statistical Society. Series (Methodological), 36(1):99-102, 1974.
Noor Awad, Neeratyoy Mallik, and Frank Hutter. DEHB: Evolutionary hyberband for scalable, robust and efficient hyperparameter optimization. In Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence International Joint Conferences on Artificial Intelligence Organization, aug 2021.
[5] Maria Bänkestad, Jens Sjölund, Jalil Taghia, and Thomas Schön. The elliptical processes: a family of fat-tailed stochastic processes. 03 2020.
[6] Johnathan M. Bardsley. Applications of a nonnegatively constrained iterative method with statistically based stopping rules to ct. pet, and spect imaging. Electron. Trans. Numer. Anal., 38:34 43, 2011.
[7] A. Beskos, F. J. Pinski, J. M. Sanz-Serna, and A. M. Stuart. Hybrid Monte-Carlo on Hilbert spaces. Stochastic Processes and their Applications, 121:2201-2230, 2011.
[8] Alexandros Beskos. A stable manifold MCMC method for high dimensions. Statistics & Probability Letters, 90:46-52, 2014.
[9] Alexandros Beskos, Mark Girolami, Shiwei Lan, Patrick E. Farrell, and Andrew M. Stuart. Geometric MCMC for infiniteddimensional inverse problems. Journal of Computational Physics, 335, 2017.
[10] Alexandros Beskos, Gareth Roberts, Andrew Stuart, and Jochen Voss. MCMC methods for diffusion bridges. Stochastics and Dynamics, 1003)339-3500 2008.
[11] Alessandro Buccini, Mirjeta Pasha, and Lothar Reichel. Linearized krylov subspace bregman iteration with nonnegativity constraint. Numerical Algorithms, 87(3):1 sep 2020.
[12] Stamatis Cambanis, Steel Huang, and Gordon Simons. On the theory of elliptically contoured distributions. Journal of Multivariate Analysis, 11(3):368-383 1981.
[13] Victor Chen, Matthew M. Dunlop, Omiros Papaspiliopoulos, and Andrew M. Stuart. Dimension-robusst meme in bayesian inverse problems. 03 2018.
[14] Simon L Cotter, Gareth O Roberts, AM Stuart, and David White. MCMC methods for functions: modifying old algorithms to make them faster. Statistical Science, 28(3):22-4466 2013.
[15] Masoumeh Dashti, Stephen Harris, and Andrew Stuart. Besov priors for bayesian inverse problems. Inverse Problems and Imaging, 6(2):183 122:183-200, may 2012.
[16] Masoumeh Dashti and Andrew M. Stuart. The Bayesian Approach to Inverse Problems, pages 311-4311 1-428. Springer International Publishing, Cham, 2017.
[17] Matthew M. Dunlop, Dejan Slepčev, Andrew M. Stuart, and Matthew Thorpe. Large data and zero noise limits of graph-baser eenissuperiieed learning algorithms. Applied and Computational Harmonic Analysis, 49(2)65556677 2020.
[18] Stefan Falkner, Aaron Klein, and Frank Hutter. BOHB: Robust and efficient hyperparameter optimization at scale. In Jennifer Dy and Andreas Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 1437-1446. PMLR, 10-15 15 Jul 2018.
[19] K. Fang and Y.T. Zhang, Generalized Multivariate Analysis. Science Press, 1990.
[20] Osama S. Faragallah, Heba El- Hi-Hoseny, Walid El-Shafai, Wael Abd El- El-Rahman, Hala S. El-Sayed, EI- Sayed M. El- El-Raaaie, Fathi E. Abd El El-Samie, and Gamal G. N. Geweid. A comprehensive survey analysis for present solutions of medical image fusion and future directions. IEEE Access, 9: 11371, 2021.
[21] Subhashis Ghosal, Jayanta K. Ghosh, and Aad W. van der Vaart, Convergence rates of posterior distributions. The Annals of Statistics, 28(2). apr 2000.
[22] Subhashis Ghosal and Aad van der Vaart. Fundamentals of nonparametric bayesian inference. 2017.

---

[23] E. Gómez. M.A. Gomez-Vilegge and J.M. Marin. A multivariate generalization of the power exponential family of distributions. Communications in Statisticy Theory and Methods, 27(3):589-600, jan 1998.
[24] A. P. Dawid J. M. Bernardo, J. O. Berger and A. F. M. Smith. Regression and classification using gaussian process priors. Bayesian Statistics, 6:475-501, 1998.
[25] Junxiong Jia, Jigen Peng, and Jinghuai Gao. Bayesian approach to inverse problems for functions with variable-index besov prior. Inverse Problems, 3288000 2016.
[26] Mark E. Johnson. Multivariate Statistical Simulation, chapter Elliptically Contoured Distributions, pages 106-124, Probability and Statistics. John Wiley & Sons, Ltd, 1987.
[27] Y. Kano. Consistency property of elliptic probability density functions. Journal of Multivariate Analysis, 51(1):139 147, 1994.
[28] Tomasz J. Kozubowski, Krzysztof Podgórski, and Igor Rychlik. Multivariate generalized laplace distri- bution and related random fields. Journal of Multivariate Analysis, 113:59-72, 2013. Special Issue on Multivariate Distribution Theory in Memory of Samuel Kotz.
[29] Shiwei Lan. Learning temporal evolution spatial dependence with generalized spatiotemporal gaussian process models. Journal of Machine Learning Research, (3(25))1-553. 2022.
[30] Shiwei Lan, Shuyi Li, and Babak Shahbaba. Scaling up bayesian uncertainty quantification for inverse problems using deep neural networks. SIAM Journal of Uncertainty Quantification, 2022. to appear.
[31] Matti Lassas, Eero Saksman, and Samuli Siltanen. Discretizatiooninvarian hayesian inversion and besov space priors. Inverse Problems and Imoging, 311):87-122 2009.
[32] Matti Lassas and Samuli Siltanen. Can one use total variation prior for aige-preserving Bayesian inversion? Inverse Problems, 20(5):1537, 2004.
[33] Lisha Li, Kevin Jamieson, Giulia DeSalvo, Afshin Rostamizadeh, and Ameet Talwalkar. Hyperband: A novel bandit-based approach to hyperparameter optimization. J. Mach. Learn. Res., 18(1):6765-688 68166, jan 2017.
[34] Felix Lucka. Fast markov chain monte carlo sampling for sparse bayesian inference in high- nigh-dimensional inverse problems using 11-type priors. Inverse Problems, 28(12):12501,, nov 2012.
[35] Bernt Øksendal. Stochastic Differential Equations. Springer Berlin Heidelberg, 2003.
[36] N. Petra and G. Stadler. Model variational inverse problems governed by partial differential equations. Technical report, The Institute for Computational Engineering and Sciences, The University of Texas at Austin., 2011.
[37] Krzysztof Podgórski and Jörg Wegener. Estimation for stochastic models driven by laplace motion. Commanications in Statistics Theory and Methods, 40(18):3281 3302, sep 2011.
[38] Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning. The MIT Press, 2005.
[39] Rafael Reisenhofer, Sebastian Bosse, Gitta Kutyniok, and Thomas Wiegand. A haar wavelet-basee perceptual similarity index for image quality assessment. Signal Processing: Image Communication, 61:33-43, 2018.
[40] J. Schoenberg, Metric spaces and completely monotone functions. Annals of Mathematics, 39:811 1938.
[41] Amar Shah, Andrew Wilson, and Zoubin Ghahramani. Student-t Processes as Alternatives to Gaussian Processes. In Samuel Kaski and Jukka Corander, editors, Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics, volume 33 of Proceedings of Machine Learning Research, pages 877-885, Reykjavik, Iceland, 22-25 Apr 2014. PMLR.
[42] I. A. Shepp and B. F. Logan. The fourier reconstruction of a head section. IEEE Transactions on Nuclear Science, 21(3):11-33 1974.
[43] Andrew M Stuart. Inverse problems: a Bayesian perspective. Acta Numerica, 19:451-559, 2010.
[44] Simopekka Vänskä, Matti Lassas, Samuli Siltanen, and Rolf Insitute. Statistical x-ray tomography using empirical besov priors. International Journal of Tomography and Statistics, 11, 06 2009.
12

---

[45] Umberto Villa, Noemi Petra, and Omar Ghattas. hIPPYlib: An extensible software framework for large-scale inverse problems govemndd by PDEs; part i: Deterministic Inversion and Linearized Bayesian Inference. ACM Transactions on Mathematical Software, 47(2):1-34, jun 2021.
[46] Zhou Wang. Alan C Bovik, Hamid R Sheikh, and Eero Simoncelli. Image quality assessment: from error visibility to structural similarity. IEEE transactions on image processing, 13(4);600- 612, 2004.
[47] MIKE WEST. On scale mixtures of normal distributions. Biometrika, 74(3):646-648, 09 1987.
13

---

Supplement Document for "Bayesian Learning via Q-Exponential Process"
A PROOFS
A.1 Proof of Theorem
Proof. First we prove the exchangeability of (H, C) with general (non-identity) covariance matrix C= (Iidda for some kernel function E. It actually holds for all elliptic distributions including MVN. Their densities contain the essential quadratic form r(u) =-C-'u which is invariant under any permutation of coordinates.
Denote - -.- :,,,""' " , " '.. Myy] 1" Without loss of generality, only need to show r(u) invariant by switching two coordinates, say, 4 ++ Denote u' = "' "My,'" - " u4gl' . Switching & and I leads to a different covariance matrix C' obtained by switching both i-th and j-th and columns simultaneously in C. If we denote the elementary matrix Eij as derived from switching i-th and j-th rows of the identity matrix I. Then we have
E,ju, E,,
Note Ei is idempotent, i.c. E Therefore
E,,
Next, the consistency directly follows from Kano's consistency Theorem 3.2 with our choice of gr). The proof is hence completed.
A.2 Theorem of Q-EP as a mixture of Gaussians
Theorem A.1. Suppose q (D,(EDO for 0 q 2, then there exist random variable V 0 and a standard normal random vector Z. ~ (2, independent of each other such that 7./V.
Proof. Based on
3 it suffices to show (- ((r) 2 0 for all & E N. Observe that
g'(r) = (g 1) 1)4-11 4 +x exp
y $ 0 when S 2. Denote (-d) (-g)) =
(x 9-11/11
1112 11 1=) - where the coefficients of polynomial pt are all non-negative. Then
have
(-ar) [d)
k+1
8(r) = d
-aa0----a(((-111100 13/2 , 9
r - -1 (1 11/2 , expl
expl-
where Pk+1 1 (r(3-11/2 1 being the term in the square bracket has all positive coefficients because the powers (3 1)/2 and -1 appear coefficients in 4. ( 11/2 11) and are both negative. The proof is completed by induction.
A.3 Proposition distribution of r(u)
The following proposition determines the distribution of / r(u) as g-root gamma (also chi-squared) distribution thus gives a complete recipe for generating random vector (0, based on the stochastic representation (6).
Proposition A.1. Ifu ~ q- (-EDa(0, C),
then we have
r + [(
R T B 2 Zd: and E[R(] 2
~ d as d + VkEN (18)
Proof. With out chosen g(r), the density of r becomes
14
f(r) x -1, 1) exp
exp
2
2
14

---

A change of variable r + r2 yields the density of R = rå that can be recognized as the density of %i
On the other hand, since V [( = B 2), we have: 1
E| / vd fvr)dr
v ++ 2- exp
1 2")
dv
r(2)
[ ()
d
- I(5)
~2
de
where we use T(x+ a) ~ [(x)x*
as + 0o with x and a = A when d 00
A.4 Proof of Proposition 3.1
Proof. By Theorem 2.6.4 in 99 for q (-EDa(H,C) = EC,(H, C,g) with our chosen & we know
Eluj = f and Cov(u)- = (E[R*) (C)n(C)) C. It follows by letting = 2 in Proposition A.1 and using
the similar asymptotic analysis.
A.5 Proof of Theorem
Proof. Note we can approximate (x(x) (D) with simple functions (x(x)
E KD (x,(,,,))))) (x) where
D s are measurable subsets of D and XD, (x) = if E D, and 0 otherwise. By the linear combination property of elliptic distributions [c.f. Theorem 2.6.3 in 19], ñ; SD (x))x 4-ED(O,e e) with
28 2 +)
c= o E[u?]] to be determined. Note Og =
- comes from Proposition 3.1 and the
ar())
r(g
scaling u* d2 uu in Definition 3.2 We have ay (())
as d oo, Taking the limit oc, we have ue In (xx 02 (x) ED(O,c). In general, by the similar argument we have Cov(us, ue) = E[uuur]
J
(x)xdx'
= . 8 (xx,))))()))) (xxxx = ( = 2.3
Thus completes the proof.
A.6 Proof of Theorem 3.5
Before proving Theorem 3.5 we first prove the following lemma based on the conditional elliptic
distribution 122 19).
Ci1
C12
Lemma A.1. !fu = (u,, u2) ~ (-D C) with μ =
and
u Rdd u, e IR4
[2
C21 C2
for 1,2 and di + d2 d. then have the following conditional distribution
un un q- DD, (M 12 C11-2),
M μ4 C12 C (n2 M2), C112-11 C12C C21
Proof. This directly follows from [Corollary 5 of Theorem 5 in 12 or [Corollary 3 of Theorem 2.6.6 in 19 for q-ED,(H, ECC, (H, C,g) with our chosen g.
Now we prove the Theorem 3.5
Proof. By the linear combination property of the elliptic distributions 26 ED(O, C+ Then based on the consistency, we have the joint distribution
19
have y ~
C+E
Σ
C, C..
ED
0.
[(xx)
CT
Therefore, the conclusion follows from Lemma A.I
15

---

A.7 Proposition of Conditional Conjugacy for Variance Magnitude (o*)
Proposition A.2. If we assume proper inverse taeee-ggammaa prior for the variance magnitude such that u aD oo), and σ ~ r- 1-rrraaa (a. then we have
o u ( (', d'= d'=atz d
(u (1-1) Co (u-H
B' B
(19)
Proof. Denote r0- = (u- m1) cT' (u 4.. We can compute the joint density of u and σ²
p(u, =
- 9 (2x) $|oo - - o
(2-1)
of
ro
y
Ba
exp
-
o 4
2
[(a)
9
< x< (a+ (a+
exp
o
B B
2
By identifying the parameters for σ we recognize that 6 u another inverse -aaema with parame- ters o' and B' as given.
B ALGORITHM
Algorithm 1 White-noise Preconditioed Crank-Nicolson (wn-pCN) for Q-EP Prior Models
1: Fix B E (0, 1]. Choose initial state (0 E TRd.
2: for k k=0, 0,... K do
3: Propose 2(8) = (1 B2) + BB', o' ~ ((.I).
4: Set Set ((++)) (k)) with acceptance probability
( L(T min { { L(T(
5: else set (L+1) (4).
6: end for
C ADDITIONAL EXPERIMENTAL RESULTS
In this section, we present some additional numerical experimental results that cannot be included in the main text due to the page limit.
First, numerically verify the equivalence between the stochastic representation (6 and the white- noise representation (17) of q-ED, random variable in Figure B.I More specifically, we generate 10000 samples using each of these representations and illustrate in Figure B.I that the two samples yield empirical marginal distributions (Id and 2d) close enough to each other.
C.1 Time Series Modeling
For modeling the simulated time series and stock prices, include the optimization trace of negative (log)- (oot)rrostrrio densities and relative errors for the two simulations and two stocks prices in Figure As commented in the main text, these plots show that Q-EP model can converge faster to lower errors compared with GP and Besov models.
Next, we compare MAP estimates by GP, Besov and Q-EP models in Figure C.2a for simulated time
series with step jumps and in Figure C.2c for the Google stock prices 2022. We also investigate
the prediction results by GP and Q-EP in these two examples in Figures C.2b C.2d Table C.1 summarizes the RMSE of estimated stock prices by the three models and standard deviation repeating the experiments 10 times independently.
16

---

mooch-rap
X
Figure B.1: Comparison in sampling q ED using the stochastic representation 66) (organge) and the white-noiss representation (17) (blue). Numerical results show their sampling distributions are indistinguishable. Empirical densities are estimated based 10000 samples (shown dots).
boce
391
390
Eopoaree
43
155
30
(a) Time series step jumps.
(b) Time series with sharp turnings.
10
Dandoon
Exrizion
(c) Tesla stock prices in 2022.
(d) Google stock prices in 2022.
Figure C.I: Negative posterior densities (left) and errors (right) as functions of iterations in the BFGS algorithm used to obtain MAP estimates. Early termination is implemented if the error falls below some threshold or the maximal iteration (1000) is reached. Relative errors are compared against truth in the simulation and the actual data in the Tesla stock.
C.2 Computed Tomography Imaging
In the problem of reconstructing human head and torso CT images, Table C.2 compares GP, Besov and Q-EP models in terms relative error (RLE), log-likelihood (LL), and imaging quality metrics including PSNR, SSIM and HamPSI. In most cases, Q-EP outperforms, achieves comparable scores with the other two methods.
Lastly, Figures and C.4 that the posterior standard deviations estimated by wn-pCN using GP model could misleading because seemingly more recognizable shape deludes the fact that they are about two orders of magnitude smaller in value compared with the other two models. This implies that GP might underestimate the uncertainty present in the observed sinograms in the CT imaging analysis.
17

---

2.2
70
05
0.
00
03
25
0.0
35
(a) Time series with step jumps (model fitting).
(b) Time scries with jumps (prediction).
140
satt
sad
2022 es-n
2022 08-)
ant anan mnte tiat
aide eiz anda osst
ait mnte maruii
(c) Google stock prices in 2022 (model fitting).
(d) Google stock prices in 2022 (prediction).
Figure C.2: (c)(c) MAP estimates by GP (left), Besov (middle) and Q-EP (right) models. (b)() Predictions by GP (left) and Q-EP (right) models. Orange dots are actual realizations (data points). Blue solid lines are true trajectories. Black ticks indicate the training data points. Red dashed lines are MAP estimates. Red dot-dashed lines are predictions with shaded region being credible bands.
Gaussian
Besov
q-Exponential
06
0.0200
0.8
omm
0.029
94
0.1
62
0.0050
0.8
Figure C.3: Shepp-Logan phantom: uncertainty field (posterior standard deviation) given by GP, Besov and Q-EP models. GP tends to underestimate the uncertainty values (about I 5% of that with Q-EP).
Gaussian
Besov
E-Cxooeeniall
20a
26.9
688
0.29
35
37.5
011
0.15
78
0.2
0.05
500
5.
23
Gaussian
Besoav
q Expoorential
3.9
125
20.8
04
0.75
75
5.100
58
0.0
0.6
0.0
0,000
0.0
0.2
0.6
0.0
02
0.6
0.0
00
Figure C.4: CT of human head (upper) and torso (lower): uncertainty field (posterior standard deviation) given by GP. Besov and Q-EP models. Note GP tends to underestimate the uncertainty values (about 1% of that with Q-EP).
18

---

Table C.I: Posterior estimates Tesla and Google stock prices by GP, Besov and Q-EP prior models: 22. Results are repeated 10 times with different random seeds.
RMSE:
Tesla
Google
GP
Besov
Q-EP
GP
Besov
Q-EP
RMSE std(RMSE)
171.8515
90.3086
83.8130
20,4095
25.2012
18.3597
1.8018
1.1478
2.6949
0.7115
1.9698
0.9617
Table C.2: MAP estimates for CT of human head and torso by GP, Besov and Q-EP prior models: relative error, RLE == û u |//||* of MAP (ú u²), log-likelihood (LL), PSNR, SSIM and HamPSI.
Head
Torso
GP
Besov
Q-EP
GP
Besov
Q-EP
RLE LL
0.2999
0.2241
0.2224
0.2611
0.2177
0.2153
-4.05e+5
-1.12e+4
-1.17e+4
-3.30e+5
-3.86e+3
-4.37e+3
PSNR
24.2321
26.7633
26.8281
23.6450
25.2231
25.3190
SSIM
0.7010
0.7914
0.8096
0.5852
0.6983
0.6982
HaarPSI
0.0525
0.0593
0.0587
0.0666
0.0732
0.07190
C.3 Noisy/Blurry Image Reconstruction
Next we consider reconstructing a (128 x 128 pixels) image of satellite shown on the leftmost of Figure 1 from blured observation next it. The image itself be viewed function u(x) on the square unit D [0, 1J2 taking values as the pixels. When evaluating u(x) the discretized domain, u(x) becomes matrix of size 128 128, which further be vectorized
E Fd with d = 1282. The true image, denoted as u², is blurred by applying motion blur point spread function [IS and adding 5% Gaussian noise. The actual observation, y(x), can be written as in the following linear model:
y(x) = Au(x) +E, E~ N (0,
De/I|Au|| where A E /R/dd is the blur motion PSF with / d and de |||| = 5%. Note, the blurring effect in the observed image (the second from left of Figure is mainly due to the PSF operator A, the small Gaussian noise.
Gaussian
Gaussian
negative posterior
Besov
Besov
100
q-Exponential
g-Expnnential
relative error
105
0.7
50
100
150
2000
250
300
530
400
1
50
100
150
200
2.5
300
330
400
iteration
iteration
Figure C.5: Image of satellite: negative posterior densities (left) and errors (right) as functions of iterations in the BFGS algorithm used obtain MAP estimates. Early termination is implemented the error falls below some threshold or the maximal iteration (1000) is reached.
We compare the reconstructions by MAP estimate in Figure The output by GP is blurry and close the observed image, which means that GP does "de-noise" much. The result by Besov much better than GP due to the L1 regularization but it is still not sharp enough. We can see that the Q-EP prior model produces the reconstruction of the highest quality. [ngureejdemonstratte
the effect of q >> 0: the smaller 4, the more regularization and hence sharper reconstruction. We also compare
their negative posterior densities and relative errors, || l u²1 /||="]], ||. Figure C.5 The Q-EP prior model yields the smallest error among all the three models.
19

---

Timer
met
2.0
031
(a) True initial condition (top left), and the solutions a(x,7) different time points
(b) Spatiotemporal observations at 80 selected locations
(color dols) across different time points.
Truth
Gaussian
Besov
Q-Exponential
0.6
00
0.
34
66
0.
1.0 60
62
0.
06
0.5
1.0
0.
02
0.4
06
0.8
10
00
02
04
06
10
Figure C.7: Advection-diffusion inverse problem: true initial condition un and posterior mean estimates by GP, Besov and Q-EP models.
C.4 Advection-Diffusion Inverse Problem
Finally, we consider a Bayesian inverse problem governed by a time-dependent advection-diffusic equation [36] (3) that can be applied to heat transfer, pollution tracing, etc. The inverse problem involves inferring an unknown initial condition 40 L (22) from spatiotemporal point measurements {y(x,.t)) as
y(x,1) = 9 ( (x.1), (x.1) ~ N (0,
The forward mapping & &u maps the initial condition uy pointwise spatiotemporal obser- vations of the concentration field u(x,1) through the solution of the following advection-diiftusion equation [36 45]:
1 Re
KAU+ Vu = 0 in Q x (0,T)
Av p Vy 0 in Ω
u(-00) 4g in Ω
V v 0 in Ω
KVH π 0, on aΩ (0,7)
v g, on aΩ
where [0, 1 is bounded domain shown in Figure C.6a K 10 is the diffusson coefficient, and T > 0 is the final time. The velocity field is computed by solving the following steady- Navier-Stokes equation with the side walls driving the flow [36]. Here, p is the pressure, and Re is the Reynolds number, which is set to 100 in this example. The Dirichlet boundary data g R² is given by g = e2 = (0, 1) on the left wall, g = e2 on the right wall, and g everywhere else.
To generate data, we set the true value parameter rip in C.4 as n 0.5 A cxp(-1000 (x0 0.35) (y- 0.7) illustrated in the top left panel of Figure C.6a which also shows few snapshots of the solutions u(x,t) other time points
regular grid mesh of size 61 61. Spatiotemporal
observations (y(X;,;) (/) are collected of two inner boxes (See Figure C.6a and also Figure C.6b
80 selected locations {x:) around the boundary
across / = 16 time points {tz) evenly
20

---

forward prediction (x-0.325, y-0.01)
forward prediction (x-0.249, y-0.225)
forward prediction (X-0.219. y.0...00
0.15
0.15
0.20
0.10
0.10
0.15
0.05
0.05
0.10
0.00
0.05
0.05
-0.05
-0.10
Truth
Truth
0.00
Truth GausSian
-0.15
-0.15
-0.05
-0.20
@EEpponeniaal
@Exponential
q-Exponential
-0.20
0
5
Figure C.8: Advection-difftsion inverse problem: comparing forward predictions, 9 (x,1,) based on the GP (blue dashed lines) and Q-EP (orange dot-dashed lines) prior models at three selective locations (0.325,
(0.249, 0.2.299,.225) and (0.249, 0.350). Blues dots are observations.
distributed between 1 and 4 seconds (thus denoted as Øu) with noise variance I = / where On = 0.5 max du, i.e. y(x,,,)) = 9 (u]) ni = u(x,(f)) + Mij-
To solve the inverse problem of finding the initial condition in the Bayesian framework 43 16 we impose Un with GP, Besov and Q-EP priors respectively and seek the posterior p(uply). For GP and Q- EP, we adopt covariance kernel, 6 (8. yA) - defined through the Laplace operator A, where governs the variance of the prior and y/ controls the correlation length [16] 30). We set y 1 and 0=8 8 in this example. For Besov, we adopt 2d Fourier basis of the format Φ.j = cos(xix)) cos π jx2) and truncate the series (1) for the first L. = 1000 terms.
We apply wn-pCN to this challenging nonlinear inverse problem with high dimensionality (3413) of spatially discretized и at each time f. Figure C.7 compares the posterior mean estimates of uo given by these three models. Because the truth (leftmost) has clear edge its cutoff by 0.5, Q-EP is more appropriate than GP and it indeed generates better estimate closer to the truth. Figure C.8 plots the prediction of forward mapping at few selective locations on the left side of lower inner box by 3 (xx4)) D. I (1²7)) with u(s) p(uoly). Compared with GP, Q-EP predicts the solution path closer to the truth 9 (4)) (x,t.) where the observations see more dynamical changes. More importantly. Q-EP provides proper UQ with credible bands wide enough to include the true trajectories. On the other hand, the posterior estimates by GP come with much narrower error bands that miss the truth. Again, we observe GP prior model being overconfident about less accurate estimates.
21


Reviews:
[{"review_id": "WbSTCIWe0l", "reviewer": "Reviewer_V1qe", "summary": "The authors propose the 'q-exponential process', a stochastic process interpretation of Lq function regularization, which can induce sparsity in the solutions to optimization problems and can be used as a functional prior for Bayesian applications in time series regression, image reconstruction, and other applications.  ", "questions": "It may be interesting to further unpack Remark 3 regarding the similar mean & covariance for Q-EP and GP, for example showing how this affects predictions for Q-EP and GP in a tractable toy example.", "limitations": "The authors note the need for grid search over hyperparameters such as regularity parameter and smoothness parameter.    The method is tested only for q=1 and it would be interesting to explore other values.", "rating": 6, "confidence": 1, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The proposed q-exponential process has several important benefits:\n- enforcing sparsity or more sharpness in function space\n- control of correlation structure, similarly to GP\n- flexibility in choice of kernel\n- conjugacy for posterior prediction with appropriate likelihoods, in contrast to the Besov case", "weaknesses": "The paper could provide more background on the problem setting and Besov processes, as the current introduction is difficult to follow.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the suggestion on more background about Besov processes. In addition to the existing introduction which includes its mathematical definition, we will elaborate more on its implication and applications on imaging analysis.\n\nWe also appreciate the reviewer's advice on expanding Remark 3. In the revision, we will restore more mathematical details that were in the initial draft but hidden from the submission for brevity. However, we have two tractable toy examples in section 4.1 for which we compare the predictions for Q-EP and GP in Figure 3(b) and in Figure C.2(b) as well.\n\nRegarding the limitation, we emphasized in the manuscript that $q=1$ is the case often opted in (Conclusion). We did explore the effect of regularization parameter $q$ for a spectrum of different values in Figure C.6 and mentioned that briefly in the paragraph \"Connection to existing work\" in the introduction. We will make it more explicit with proper emphasis."}, {"review_id": "SdGPeOLZxh", "reviewer": "Reviewer_Yr66", "summary": "This paper proposed a new random process prior that corresponds to estimating parameters with $\\ell_q$ penalty. The process, named Q-EP, can be used to provide a shaper penalty than the standard Gaussian process. Empirical experiments show the practical use case for Q-EP.", "questions": "What is your insight when comparing the Bayesian models with various denoising works in deep learning? I believe they are very different in many aspects.", "limitations": "I do not see any potential negative societal impact of this paper.", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The derivation of Q-EP looks solid to me. I am not sure about the novelty of the work. \n\nThe experimental results look convincing, though I do not know if there is a standard benchmark.", "weaknesses": "It is hard for people unfamiliar with the field to understand the paper. In particular, the authors do not give a fair amount of text to explain the background. It is hard to relate the abstract math in the introduction with the examples provided (e.g. Figure 1). There is no clear statement claiming the connection between Q-EP with the $L_p$ regularization. The Bayesian model is not introduced until page 6. In sum, I suggest a major rewriting of the paper, potentially with a different organization.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer's critics. As mentioned in the introduction, the novelty lies in the first probabilistic definition of Besov process (which is widely used in imaging analysis and Bayesian inverse problems) with explicit specification of correlations and tractable prediction formula.\n\nFor nonparametric regression with novel priors, mean squared error (MSE) (or root MSE, RMSE) and log-likelihood (LL) (See Table 1) are standard measures to compare, as included in the paper by Professor Zoubin Ghahramani's group [37], the paper by Bankestad et~al [4], and the seminar book by Rasmussen and Williams [Gaussian Processes for Machine Learning, 2006]. We also added multiple standard quality metrics in imaging analysis (PSNR,SSIM, HaarPSI as in Table 2). They all support our claimed numerical advantages for q-EP. While we are open to any other metrics the reviewer might suggest, we are also willing to defend against \"limited evaluation\", given 4 time series examples, 3 CT image reconstructions and one Bayesian inverse problem regarding an advection-diffusion equation reported in this paper.\n\nRegarding the connection between q-EP and $L_q$ regularization, we kind of assumed for granted and that might have caused some confusion. It is similar to the relationship between GP and $L_2$ regularization: negative density of Gaussian distribution yields the $L_2$ regularization term usually added to the objective function, that is, $\\frac{1}{2}\\Vert u-\\mu\\Vert^2$. In q-EP scenario, that $L_q$ regularization term is $\\frac{1}{2} r(u)^{\\frac{q}{2}}$ with $r(u)=\\langle u-\\mu, \\mathcal{C}^{-1}(u-\\mu)\\rangle$ plus some other comparatively smaller term $\\log r$. We will add explicit explanation and more background relating the math to Figure 1 to the main text.\n\nRegarding the structure, we admit that we spent some space on explaining the marginalization consistency -- we did that to emphasize its importance and  the univariate q-exponential distribution may fail to generalize to a valid stochastic process without care, as happened in other literatures -- we think that is part of our novelty. The Bayesian model is discussed after the new prior q-EP is fully introduced. This is in the same spirit of Rasmussen and Williams' GP book, which does not jump to the model at beginning.\n\nWe thank the reviewer for this good question. First, the paper introduces a nonparametric modeling tool with the novel q-EP prior imposing more regularization than GP. It can be applied to imaging analysis but not limited to image denoising. Second, one of the advantages of the proposed Bayesian models over the majority of optimization based deep learning techniques is the uncertainty quantification (UQ) (refer to Figures C.3 and C.4). UQ is of scientific interest and is the natural byproduct of Bayesian approaches however not the main focus of many vanilla versions of denoising works in deep learning. Last but not the least, these two approaches can interact and evolve to new methods. For example, the q-EP prior can provide a new loss function other than MSE for training the denoising neural network (NN). On the other hand, there might exist certain format of NN whose limiting behavior mimics q-EP prior, similar to the relationship between DNN and GP [Neal 1994a, Lee et al 2018 ICLR]. This is also an interesting direction the authors would like to pursue."}, {"review_id": "NteNsTuZnW", "reviewer": "Reviewer_av1p", "summary": "As a generalization of GP, it is important to construct a stochastic process (prior) that can express various degrees of smoothness. As such process, the Besov processes have been proposed but they are defined in the form of a series expansion, and the corresponding probability distribution is not given in an easy-to-handle form.\nIn this work, the q-exponential process is proposed, which addresses the drawback of the conventional Besov process. A probability distribution called the q-exponential distribution (q-ED) that facilitates sampling from the Besov process is proposed. It is extended to multivariate in a consistent manner in the sense of marginalization. q-ED is expressed as an elliptic distribution, which also satisfies the exchange rate, thus satisfying the conditions of Kolmogorov's extension theorem, and the existence of the corresponding stochastic process is guaranteed. The KL expansion shows that the q-exponential process (q-EP) corresponding to the q-ED can be defined and has a series representation almost equivalent to the Besov process.\nThe q-EP has the strong advantage that the posterior predictive distribution can be constructed by MCMC, and Bayesian regression can be performed.", "questions": "Q1: In my understanding, the proposed q-EP is essentially equivalent to the Besov process, and the major difference is its explicit control of the covariance structure. If this is the case, what is the main reason for the performance gap between Besov and q-EP in experiments, particularly for image reconstruction? If not, please reveal the difference between Besov proc. and q-EP. In particular, does q-EP includes Besov proc as a special case?\n\nQ2: Related to the problem mentioned in 6.Weakness, I'm not sure the experimental setting in subsection 4.1 is fair. For the Besov process, the Fourier basis is chosen and fixed. It is plausible but the inferior performance of the Besov process could be contributed to the bad choice of the basis. Please validate this experimental setting in more detail.\n\nQ3: I don't understand the meaning of the second sentence in \"Introduction\". How the \"High-dimensional objects\" can be viewed as \"evaluation of proper functions\", and what is the \"proper function\"?\n\nQ4: What is the difference between \"d^{\\star}\" in line 34 and \"d\" in other places?", "limitations": "I'm concerned that the scale mixture expression is not possible for all q. \nSee, e.g., \nM. West, ``On scale mixtures of normal distributions'', Biometrika (1987).\nI found that the range 0 < q < 2 is specified in the Appendix, but should be explicitly stated in the body of the paper to make the applicability of the methodology clear.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "It is an important contribution from the viewpoint of statistical modeling to define an explicit probability distribution as q-ED and to provide a sampling method for Besov processes described in series representation. q-ED is firstly defined as a one-dimensional probability distribution and then is extended to a multivariate distribution paying attention to consistency. The method of extending to a stochastic distribution and the method of performing MCMC inference (posterior distribution calculation) based on stochastic representation is standard and reasonable in the context of SDE simulation (e.g., a scaled mixture of normals). \nThe proposed q-EP will have wide applicability and some successful application results are presented in the paper.", "weaknesses": "Minor issue:\nIt is claimed that \"we have less control on the correlation strength once the orthonormal basis is chosen\". It is true but when we are free to choose the orthonormal basis, I feel the conventional Besov process has equivalent or even better freedom for specifying the correlation structure to the proposed q-EP. I acknowledge that specifying the correlation structure via the covariance function is handy, but in terms of the degree of freedom, it is fair to consider the case we can freely choose the basis.\n\nIt is claimed in line 218 that \"codes will be publicly available at to-be-released\", hence currently I cannot evaluate the reproducibility of the experimental results.\n\nVery minor presentation issues:\n- Eq.(5) and elsewhere, p(u) should be treated as the function of u, not r. I mean, for example in Eq.(5), \np(u) = k_d |C|^{-1/2} g(r(u)).\n\n- F_i in Theorem 3.1. should be defined.\n- \\mathbb{R}^k in Eq.(7) should be \\mathbb{R}^{n}\n\n- References should be properly described. For example, ref [1] does not have a link to the project, and ref [12] lacks bibliographic information.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for supporting the contribution and the potential impact of q-EP. We agree with the reviewer that the correlation strengthen of Besov process can be configured through the choice of basis functions $\\{\\phi_\\ell(x)\\}$, as spelled out in Equation (12). But it is just less straightforward than specifying in the kernel function $C(\\cdot,\\cdot)$ as usually done in Gaussian process. We will reword the relevant sentences to reflect this.\n\nThe Github repository hosting all the codes is currently private but will be made public when the work gets published. However, along with the submission, we provided a zip file in the supplementary materials containing Python codes for demo and an example of reconstructing blurred image of satellite. All the codes in the supplement are anonymous and accompanied with a `readme` help document so the reviewer should be able to reproduce some results by following the instruction. We hope the reviewer can understand and find it helpful.\n\nWe thank the reviewer for careful reading and we will fix all the minor issues raised in \"Weaknesses\". We also appreciate all the questions and now answer them one by one:\n\nQ1) Yes, the statement about the relationship between q-EP and Besov process is right. As highlighted in the introduction, q-EP can be viewed a probabilistic definition of Besov process with explicit specification of correlation strength and tractable prediction formula. Their representation equivalence is further elaborated in Theorem 3.4 and Remark 2. Due to their difference in the mathematical format (probabilistic distribution vs series representation), they behave differently in numerics. We found q-EP is superior than Besov in relatively lower-dimensional cases (See Figures 3,4 and Tables 1,2) and they become more similar when the dimensions go much higher and dimension reduction (partial eigendecomposition of C) is implemented (See Figure5 and Table C.2). We agree that it is an interesting phenomenon and will investigate it in future work.\n\nQ2) We tried multiple wavelet bases including Hard, Shannon, Meyer and Mexican Hat etc. but none of them generated better result. We could include these reconstruction results in supplementary materials in the revision.\n\nQ3) Sorry for the confusion. What we tried to say is that each image can be viewed as a function defined on a bounded domain whose values are the pixels. We will revise this sentence to avoid confusion.\n\nQ4) $d^\\star$ in line 34 is the dimension of space where subjects of interest are defined. For example, in our numerical examples, $d^\\star=1$ for the time series and stocks and $d^\\star=2$ for various CT images. While $d$ in other places refers to discrete dimension of the processes. For example, $d=N=200$ is the number of time points for the time series and $d=n^2=128\\times 128$ is the image size for CT images. We will clarify them in the revision.\n\nAlso thanks for pointing out the condition $0<q<2$ for the scale mixture result only mentioned in the appendix. We will make it explicit in the main text when revising the paper."}, {"review_id": "1nypuGP5Td", "reviewer": "Reviewer_Krqg", "summary": "Motivated by the correspondence between Gaussian Process priors and ridge regularization for non-parametric regression problems, the authors in this paper develop a stochastic process prior, namely the $\\textbf{Q-exponential (Q-EP) process}$, which can correspond to $\\ell_q$-regularization. Specifically, by starting from multivariate  $q$-exponential distributions, the authors verify a Kolmogorov Consistency criterion to eventually develop the aforementioned process. Subsequently, the authors further provide justifications how the (Q-EP) process process allows a more flexible modeling of covariance kernels compared to parallel Besov processes also designed for similar tasks and derive a posterior predictive formula. Lastly, the benefits of the prior are demonstrated through numerical examples pertaining to problems in functional data analysis, image reconstruction, and solving inverse problems.", "questions": "There is substantive literature on \"minimax optimal\" posterior contraction over standard function spaces (e.g. Sobolev, Besov etc.) results for suitable designed Gaussian process priors  -- can the authors discuss whether there are some set ups of non-parametric regression problems where some version of the priors developed here offers optimal posterior contraction over some function spaces?", "limitations": "None noted.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Offers a detailed probabilistic description and construction of a stochastic process prior for non-parametric Bayesian regression which (i) corresponds to $\\ell_q$-regularization; (ii) is developed from easy to understand first principles; and (iii) allows easy computation of posterior predictive distributions. ", "weaknesses": "Theoretically, it was not discussed in which type of non-parametric problems and function spaces one might get better rates of posterior contraction using the developed process type priors compared to Gaussian process priors. Most of the theory is basic probabilistic calculations and quite stratighforward.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for raising good points in the \"Weaknesses\" and \"Questions\" sections. Q-EP is proposed as a nonparametric prior for flexible Bayesian models including regression, classification, density estimation, inverse problems, etc. The motivation is to impose more regularization than Gaussian process (GP) on function spaces while enjoying the similar tractability (on correlation and prediction) as GP. The function space is assumed to be $L^q$ and we require the kernel function $\\mathcal{C}$ is a trace-class (having summable eigenvalues) so the quadratic $r(u)=\\langle u, \\mathcal{C}^{-1}u\\rangle$ and hence the process is well-defined.\n\nWe are aware of Professor Ghosal, Professor Van Der Vaart and their collaborators' seminal works on posterior contraction theory on non-parametric Bayesian models. We actually have been working on similar theories for q-EP and a spatiotemporal Besov extension that relies on q-EP. However, given the limited space of 9 pages, we decided to focus on the introduction of q-EP, its application to nonparametric modeling and demonstration of its advantages over GP when modeling subjects with abrupt changes or sharp contrast (such as edges in image). We leave the rigorous treatment of contraction theory to a journal submission.\n\nWe thank the reviewer for raising such a good question on comparing the contraction rates between q-QP and GP and we will investigate it in our ongoing work."}, {"review_id": "VYtjx66iMs", "reviewer": "Reviewer_n1kr", "summary": "This paper proposes a generalization of Besov processes to higgher dimensions that satisfies the stochastic process constraints that previous works could not satisfy. This is achieved by finding the right radius density function so that the corresponding elliptic distribution satisfies Kolmogorov's extension theorem. They then demonstrate how one may do Bayesian modeling and inference with this Q-EP process and provide experiments demonstrating their results.", "questions": "* How should one choose q? Should q=1 always be the default starting point?", "limitations": "The authors have adequately addressed the limitations.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "* Very interesting stochastic process family that I can see having a lot of use cases\n* Compelling experimental results showing the efficacy of their method", "weaknesses": "* The experimental results use pretty small images. How well do these methods scale up to larger images and high-dimensional domains more generally?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the nice comments and the generous support.\n\nWe have three CT examples in Section 4.2: the Shepp–Logan phantom is of size $128\\times 128$ and the other two human body parts CT images are of size $512\\times 512$. They are pretty standard sizes for images in machine learning. Note, the discrete dimension $d=128^2$ or $512^2$. We need to specify covariance matrices of size $d\\times d$ for GP and q-EP which are extremely large ($512^2\\times512^2=68,719,476,736$). Therefore, we did dimension reduction by partial Eigen-decomposition (taking the first $L=2000$ largest eigenvalues) using randomized algorithms (Joel A. Tropp, ACM 204 CalTech). If the reviewer refers \"larger images\" to super-resolution satellite images, the methodology should still work in theory, but care needs to be taken when implemented in powerful GPUs or TPUs.\n\nAs commented in the conclusion, $q=1$ is often adopted for q-EP (and Besov process as well) to impose more regularization as opposed to $q=2$ which corresponds to GP. We studied the effect of the parameter $q$ in Figure C.6 in the supplementary materials by varying $q\\in (0,2]$. We found that smaller $q$ led to sharper reconstruction of the blurred image of satellite. In general, we could do cross validation to choose $q$ or even impose a hyper-prior on $q$ for a fully Bayesian treatment, which is expected to be more complicated and will be a good future direction."}]

Author rebuttal:
We thank all the anonymous reviewers for their careful reading, constructive advices and critics. All the reviews acknowledge the novelty of the proposed q-EP as a nonparametric prior and its potential impact in statistics and machine learning applications. Some demand clarification and presentation improvement, and others suggest theoretic exploration. All their comments and opinions are highly appreciated. Below we have made one-to-one responses and hope the reviewers and program/area chairs can find them useful when evaluating the paper. We will be happy to discuss more with reviewers in the discussion stage. Thanks!

Meta-review: Strengths:
- The paper presents a novel Q-exponential process that brings several advantages to nonparametric statistical modeling.
- The methodological derivation is solid, and the application of the method appears to be successful in particular contexts.

Weaknesses:
- Multiple reviewers mention that the paper is not easily accessible due to insufficient background explanation, particularly about the problem setting and Besov processes. 
- Concerns are raised about the theoretical depth of the paper, especially in comparison to Gaussian process priors.
- Reviewers ask for further clarity on the relationship between Q-EP and Besov processes, implications when comparing with Bayesian and deep learning models, and the settings where the proposed method might be advantageous. These areas need to be addressed.

Overall, the paper is expected to be a strong and impactful contribution to the field, and recommended for acceptance.  Howerver, the authors should address the raised concerns in the final version.
Decision: Accept (poster)