Deconditional Downscaling with Gaussian Processes
Authors: Siu Lun Chau
Shahine Bouabid
Dino Sejdinovic
Venue: NeurIPS 2021 Poster (2021)
Keywords: statistical downscaling
spatial disaggregation
deconditioning
Gaussian processes
kernel methods
conditional mean embeddings
conditional expectation
two-stage regression

Abstract: Refining low-resolution (LR) spatial fields with high-resolution (HR) information, often known as statistical downscaling, is challenging as the diversity of spatial datasets often prevents direct matching of observations. Yet, when LR samples are modeled as aggregate conditional means of HR samples with respect to a mediating variable that is globally observed, the recovery of the underlying fine-grained field can be framed as taking an "inverse" of the conditional expectation, namely a deconditioning problem. In this work, we propose a Bayesian formulation of deconditioning which naturally recovers the initial reproducing kernel Hilbert space formulation from Hsu and Ramos (2019). We extend deconditioning to a downscaling setup and devise efficient conditional mean embedding estimator for multiresolution data. By treating conditional expectations as inter-domain features of the underlying field, a posterior for the latent field can be established as a solution to the deconditioning problem. Furthermore, we show that this solution can be viewed as a two-staged vector-valued kernel ridge regressor and show that it has a minimax optimal convergence rate under mild assumptions. Lastly, we demonstrate its proficiency in a synthetic and a real-world atmospheric field downscaling problem, showing substantial improvements over existing methods.


Deconditional Downscaling with Gaussian Processes
Siu Lun Chau University of Oxford
Shahine Bouabid University Oxford
Dino Sejdinovic" University of Oxford
Abstract
Refining low-resolution (LR) spatial fields with high-resolution (HR) information, often known as statistical downscaling, is challenging as the diversity of spatial datasets often prevents direct matching of observations. Yet, when LR samples are modeled as aggregate conditional means of HR samples with respect to a mediating variable that is globally observed, the recovery of the underlying fine-grained field can be framed as taking an "inverse" of the conditional expectation, namely deconditioning problem. In this work, we propose a Bayesian formulation of deconditioning which naturally recovers the initial reproducing kernel Hilbert space formulation from Hsu and Ramos [[]. We extend deconditioning a downscaling setup and devise efficient conditional mean embedding estimator for multiresolution data. By treating conditional expectations as inter-domain features of the underlying field, posterior for the latent field can be established as solution the deconditioning problem. Furthermore, we show that this solution can be viewed as two-staged vector-valued kernel ridge regressor and show that it has minimax optimal convergence rate under mild assumptions. Lastly, we demonstrate its proficiency in synthetic and a real-world atmospheric field downscaling problem, showing substantial improvements over existing methods.
Introduction
Spatial observations often operate at limited resolution due to practical constraints. For example, remote sensing atmosphere products [2]
provide measurement of atmospheric properties such as cloud top temperatures and optical thickness, but only at low resolution. Devising methods refine low-resolution (LR) variables for local-scale analysis thus plays crucial part our understanding of the anthropogenic impact on climate
When high high-resolution (HR) observations of different covariates are available, details can be instilled into the LRR field for refinement. This task is referred to as statistical downscaling or spatial disaggregation and models L.R observations as the aggregation of unobserved underlying HR field. For example, multispectral optical satellite imagery [6] 7 typically comes at higher resolution than atmospheric products and can be used refine the latter.
Statistical downscaling has been studied various forms, notably giving it a probabilistic treatment [ 9,10 D1 22 [3] in which Gaussian processes (GP) [4 are typically used in conjunction with [I5] to recover the underlying unobserved HR field. Our approach follows this line of research where we do not observe data from the underlying HR groundtruth field. On the other hand, deep neural network (DNN) based approaches [7] [88 study this problem from different setting, where they often assume that both HR and L.R matched observations available for training. Then, their approaches follow standard supervised learning setting in learning mapping between different resolutions.
sparse variational formulation
Indicates equal contribution
Department of Statistics, Oxford, UK, OX1 3LB, www..chwuuusssee...... wtw.s.oo.tttttee...... ox.ac.uk, shahine.bouabide ieaati.dssttt.oooo www.s.o..ac.....oooooo akku dino.sejding mnnosseddiovviee ox.ts.x..c..k>>
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

However, both lines of existing methods require access to bags of HR covariates that are paired with aggregated targets, which in practice might be infeasible. For example, the multitude of satellites in orbit not only collect snapshots of the atmosphere at different resolutions, but also from different places and at different times, such that these observations are not jointly observed. To overcome this limitation, we propose to consider more flexible mediated statistical downscaling setup that only requires indirect matching between LR and HR covariates through a mediating field, We assume that this additional field can be easily observed, and matched separately with both HR covariates and aggregate LR targets. We then use this third-party field to mediate learning and downscale our unmatched data. In our motivating application, climate simulations [99 20 21 based on physical science can serve as a mediating field since they provide a comprehensive spatiotemporal coverage of meteorological variables that can be matched to both LR and HR covariates.
Formally, let ' = {x(1), ..... x!")} C X' be a general notation for bags of HR covariates, f - R the field of interest we wish recover and the LR aggregate observations from the field f. We suppose that "x and ž are unmatched, but that there exists mediating covariates y. 7 E 3/ such that ('x, y) are jointly observed and likewise for (Ú) B) as illustrated in Figure 1 We assume the following aggregation observation model = Ex [x = 5] + with some noise E. Our goal in mediated statistical downscaling is then to estimate ʃ given (ox,, v) and (ū, 2), which corresponds to the deconditioning problem introduced in 1
Motivated by applications in likelihood-Fie free inference and task- transfer regression, Hsu and Ramos [T] first studied the decondi- tioning problem through the lens of reproducing kernel Hilbert space (RKHS) and introduced the framework of Deconditional Mean Embeddings (DME) as its solution.
In this work, we first propose a Bayesian formulation of decon- ditioning that results into a much simpler and elegant way to arrive to the DMFE-based estimator of Hsu and Ramos (1. using the conditional expectations of f. Motivated by our mediated statistical downscaling problem, we then extend deconditioning to multi-resolution setup and bridge the gap between DMES and existing probabilistic statistical downscaling methods [9. By placing a GP prior on the sought field f. we obtain poste- rior distribution of the downscaled field a principled Bayesian solution to the downscaling task on indirectly matched data. For scalability, we provide tractable variational inference ap- proximation and an alternative approximation to the conditional mean operator (CMO) (22 to speed up computation for large multi-resolution datasets.
Figure The LR response (blue) and the bag HR covariates f (green) are unmatched. The down- scaling is mediated through bag- level LR covariates and ğ (or- ange).
From theoretical stand point, we further develop the frame- work of DMEs by establishing it as a two-staged vector-valued regressor with natural reconstruction loss that mirrors Grünewälder et al. (23 work conditional mean embeddings. This perspective allows us to leverage distribution regression theory from [24 (5 and obtain convergence rate results for the deconditional mean operator (DMO) estimator. Under mild assumptions, we obtain conditions under which this rate is a minimax optimal in terms of statistiall-computtatinnaal efficiency.
Our contributions are summarized as follows:
We propose a Bayesian formulation of the deconditioning problem of Hsu and Ramos M We establish its posterior mean DME-based estimate and its posterior covariance gauge of the deconditioning quality.
We extend deconditioning to a multi-resolution setup in the context of the mediated statistical downscaling problem and bridge the gap with existing probabilistic statistical downscaling methods. Computationally efficient algorithms are devised.
We demonstrate that the DMO estimate minimises two-staged vector-valued regression and derive its convergence rate under mild assumptions, with conditions for minimax optimality.

---

We benchmark our model against existing methods for statistical downscaling tasks in climate science, on both synthetic and real-world multi-resolution atmospheric fields data, and show improved performance.
Background Materials
2.1 Notations
Let Y,Y be pair of random variables taking values in non-empty X' and respectively, Let k: X' x X - R and :: : x ) - I be positive definite kernels. The closure of the span of their canonical feature maps k2 = k(r,) -) and lg := E(y.) -) for x E X and y y respectively induces RKHS H4 S Ex and He R endowed with inner products (,, and (t-
We observe realizations 'D1 = (4r,, M/ } of bags br, = {-/ ) from conditional distribution Px|Y=E,> with bag-level covariates y/y sampled from Py. We concatenate them into vectors x [ [ and y = 8/
yN
For simplicity, our exposition will use the notation without bagging - i.. where D1 = [zyz Mj و when the generality of our contributions will be relevant to the theory of deconditioning from Hsu and Ramos [1] in Sections 2] [ We will come back a bagged dataset formalism Sections 3.3
5 which corresponds to our motivating application of mediated statistical downscaling.
With an abuse of notation, we define feature matrices by stacking feature maps along columns as Eyw] and we denote Gram matrices as Kxx
2x = [Riz,
kizg and My == In d 4x = [,[2,, ",),,r,, and Lyy = y = [e(y), iillis The notation abuse (-)T())
a shorthand for the elementwise RKHS inner products when clear from the context.
Let Z denote the valued random variable stemming from the noisy conditional expectation of some unknown latent function X + R. as Z [x e. We suppose one observes
another set of realizations D2 = {ūg د from Pyz,, which is sampled independently from Dj.
Likewise, we stack observations into vectors y := [5 define feature map My [Ci
jas] 2 == [ ΣM]] and
L.
2.2 Conditional and Deconditional Kernel Mean Embeddings
Marginal and Joint Mean Embeddings Kernel mean embeddings of distributions provide a pow- erful framework for representing and manipulating distributions without specifying their parametric form [22 26). The marginal mean embedding of measure Px defined as μ. == Ex kxl E HA and corresponds to the Riesz representer of expectation functional f + Ex [x)] It can hence be used to evaluate expectations Ex [xX) = (ʃ, fx)x.. the mapping Px + px is injective, the kernel k is said to be characteristic 277 property satisfied for the Gaussian and Matérn kernels on Rd
28 In practice, Monte Carlo estimator fax NN E liz provides an unbiased estimate of
Extending this rationale to embeddings of joint distributions, we define Cry = Eyler @ ly] He & He and Cxy == Ex [kx @ ly] E Hs & Hr, which can be identified with the cross- covariance operators between Hilbert spaces Cry : He + Hr and Cxy He Ha. They correspond to the Riesz representers of bilinear forms (g. g') + Cov(g(Y),, gYY) = (g. Cyyg')e and (f.g) Cov(f(X), = (f. 9) As above. empirical estimates
obtained as EEmeg Cy. @ ly and ČxY } P 1/ EN kg, @ ly Again, notation abuse (.)(-T is shorthand for element-wise tensor products when clear from context.
Cyy By Dy
Conditional Mean Embeddings Similarly, one can introduce RKHS embeddings for conditional distributions referred to Conditional Mean Embeddings (CME). The CME of conditional prob- ability measure IP PXiy== is defined as pX\Y= == Ex [ = y] E Hy. As introduced by Fukumizu et al. (27). is common formulate conditioning in terms of Hilbert space op- erator CxlY He + H4 called the Conditional Mean Operator (CMO). CxiY satisfies by definition Cx pX|Y= and C Ex(f(X)\Y -I, e He. where CT de-
3

---

notes the adjoint of Cxly. Plus, the CMO admits expression Cxjy = CxrCrr- provided by E |Range(Crr)) Vy E y 27). Song et al. [30] show that nonparametrie empirical form of the CMO writes
Čx ( + NAIN)
(1)
where X > 0 is some regularisation ensuring the empirical operator is globally defined and bounded.
As observed by Grünewälder et [23], since Cxly defines mapping from He to Has it can be interpreted as the solution to a vector-valued regression problem. This perspective enables derivation
of probabilistic convergence bounds on the empirical CMO estimator 23 25
Deconditional Mean Embeddings Introduced by Hsu and Ramos ① as new class of embeddings, Deconditional Mean Embeddings (DME) are natural counterparts of CMEs. While CME μ'|Y==
Hs. allows to take the conditional expectation of any / E Hr through inner product [x = v] (J, =x] the DME denoted MX=xJY E Hr solves the inverse problem and allows to recover the initial function of which the conditional expectation was taken, through inner product (Ex[S(X) = =). = fz)).
The associated Hilbert space operator, the Deconditional Mean Operator (DMO), is thus defined as the operator Dxlr : HA - He such that D Xx = J.VS E He- It admits an expression in terms of CMO and cross-covariance oper- ators Dxly = (rx) (CxIxCrr(Cxxx)
provided ly Range(Cry)) and kn E Eango(CxiyCrrCC , vy E y and Vo E X. A nonparametrie empirical es- timate of the DMO using datasets Dí and D2 as deseribed above, is given by DXIY = Vy (A A + MMMM) A I where > 0 is a regularisation term and A = (Lyy XAI "Lyg can be interpreted as a mediation operator. Applying the DMO to expected responses ž, Hsu and Ramos [ are able to recover estimate of as
f(x) = k(z, x) (A KxxA MeIM) -1
(2)
Note that since separate samples ý can be used to estimate Cry, this naturally fits mediating variables setup where x and the conditional means à are not jointly observed.
3 Deconditioning with Gaussian processes
In this section, we introduce Conditional Mean Process (CMP), stochastic process stemming from the conditional expectation of a GP. We provide characterisation of the CMP and show that the corresponding posterior mean of its integrand is a DME-based estimator. We also derive in Appendix B]a variational formulation of our model that scales to large datasets and demonstrate its performance in Section
For simplicity, we put aside observations bagging in Sections 3.1 and 3.2 our contributions being relevant the general theory of DMEs. We return to a bagged formalism in Section 3.3 and extend deconditioning to the multiresolution setup inherent to the mediated statistical downscaling application. In what follows, X' is a measurable space, y Borel space and feature maps k- and Ey are Borel-measuree functions for any X, y. All proofs and derivations of the paper are included the appendix.
3.1 Conditional Mean Process
Bayesian quadrature 14 31 32 is based on the observation that the integral of a GP with respect to some marginal measure is Gaussian random variable. We start by probing the implications of integrating with respect to conditional distribution Px and considering such integrals as functions of the conditioning variable y. This gives rise to the notion of conditional mean processes.
Definition 3.1 (Conditional Mean Process), Let f ~ GP(m, k) with integrable sample paths, i.c. fx [f| dFx < xo a.s. The CMP induced by f with respect to Px y is defined as the stochastic
'the slightly unusual notation px=a|y is taken from Hsu and Ramos D and is meant to contrast the usual conditioning y

---

process (g(1) y E y} given by
( = Exxs(x) = y| | f(x) dP x(I==((()))
(3)
By linearity of the integral, it is clear that g(y) is a Gaussian random variable for each y E y. The sample paths integrability requirement ensures g is well-defined almost everywhere. The following result characterizes CMP as a GP on y/
Proposition 3.2 (CMP characterization). Suppose Ex [x) < oo and Ex [||xx | oo and let YX' ~ Pxy. Then g is Gaussian process 9 ~ GP(v, q) a.s. specified by
v(y) (m) y) Ay. y' E y. Furthermore, q(y.3') = XXIIY=y-XXII--> (
((,,0') Exxx (X(X, X') y. Y' v/]
(4)
a.s.
Intuitively, the CMP can be understood GP on the conditional means where its covariance 4(y. y/ is induced by the similarity between the CMEs at and 1/. Resorting to the kernel defined on 3 we can reexpress the covariance using Hilbert space operators as ( yy) = (Cxly Ey. Cx &y )a- natural nonparametric estimate of the CMP covariance thus comes using the CMO estimator from ①I
as {y, 5²) yv. (yy NAI ) Kxx (Ly NAIN) [(y, y'). When
Hk. the mean
function can be written as vvyy = [X =y: m) for which we can also use empirical estimate i(y) =
Eyy.y y) yy NAIN) Φ m. Finally, one also quantify the covariance between the CMP g and its integrand f. i.e. (or( (vy) = Ex [k(x,, X) = y]. Under the same assumptions as Proposition 3.2 this covariance be expressed using mean embeddings, i.e. Cov(f(x), g(v))
(kis (xX(Y=y)> and admits empirical estimate k(x, x) (Ly + NAIN)) E(y,y).
3.2 Deconditional Posterior
Given independent observations introduced above, D, = (x, and D2 = (ý, 2), we may now con- sider an additive noise model with CMP prior on aggregate observations ž (g(()), o'I)) Let Qyy q(y.3 y) be the kernel matrix induced by q on ý and let Y := Cov(f(x), 2) o Cxpy be the cross-covariance between f(x and 2. The joint normality of f(x) and ã gives
f(x)
mm(x)
KXX r
y. yy.- ~ N
Qyy + +
(5)
Using Gaussian conditioning, GP field f given the aggregate observations ă corresponding to y. The latter can naturally be degenerated if observations are paired, i.e. y = y. This formulation can be seen example of the inter-domain GP 33 where we utilise the observed conditional means ă as inter-domain inducing features for inference of f.
then readily derive the posterior distribution of the underlying
Proposition 3.3 (Deconditional Posterior). Given aggregate observations & with homoscedastic noise o², the deconditional posterior of f is defined as the Gaussian process flž gP ka) where
ma(x) m(x) + k Cxly y Qgys + Im) -(v()))), ka(x, ,²) = k( x') - k.j Cx v QQy + o I' - Kai
(6)
(7)
Substituting terms by their empirical forms, can define nonparametric estimate of the ma as ma(x) m(x) + k(æ;, x) A I'IM)
(IM))1 - (?-D())))
(8)
which, when m = 0, reduces to the DME-based estimator in (2 by taking the noise variance o² as the inverse regularization parameter. Hsu and Ramos [1] recover a similar posterior
expression in their Bayesian interpretation of DME. However, they do not link the distributions of f and its CMP, which leads much more complicated chained inference derivations combining fully Bayesian and MAP estimates, while we naturally recover it using simple Gaussian conditioning.
Likewise, empirical estimate of the deconditional covariance is given by
ka( = k(x, A(x, x) AQ) + + L'L)) ' AT(xx'').
(9)
5

---

Interestingly, the latter can be rearranged to write as the difference between the original ker- nel and the kernel undergoing conditioning and deconditioning steps (a = k(a,=²)
(Ka, D Dx This can be interpreted measure of reconstruction quality, which degen- erates in the case of perfect deconditioning, i.e. DxIr = Idd.-
3.3 Deconditioning and multiresolution data
Downscaling application would typically correspond to multiresolution data, with bag dataset 1=-1- In this setup, the mismatch in size between vec- "Di {("z,,,;)) ,, where "x, {ar p
(a[")) ... TN = and [M1 ... un) prevents from readily applying ①1 to estimate the CMO and thus infer the deconditional posterior. There is, however, a straight- forward approach to alleviate this: simply replicate bag-level covariates yj to match bags sizes ni-
tor concatenations
Although simple, this method incurs a O({E)= n) cost due to matrix inversion in
Alternatively, since bags Ly are sampled from conditional distribution Px|Y= unbiased Monte Carlo estimators of CMEs are given by px|Y=
Ei 1gg. Let My
".x Ax denote their concatenation along columns. We can then rewrite the cross- covariance operator Cxv = [Ex [v] and hence take + My Vy an estimate for Cxy. Substituting empirical forms into Cx = CxiyCry and applying Woodbury identity, we obtain alternative CMO estimator that only requires inversion of N x N matrix. We call it Conditional Mean Shrinkage Operator and define as
$ *ÜxIV := :y My + +NINN)
(10)
This estimator can be seen as generalisation of the Kernel Mean Shrinkage Estimator [34 to the conditional case. We provide Appendix modifications of (8) and (9) including this estimator for the computation of the deconditional posterior.
C
4 Deconditioning as regression
obtain DMO-based estimate for the posterior mean of flž. When the estimate gets closer to the exact operator, the uncertainty collapses and the Bayesian view meets the frequentist. It is however unclear how the empirical operators effectively converge in finite data size. Adopting an alternative perspective, we now demonstrate that the DMO estimate can be obtained as the minimiser of two-staged vector-valued regression. This frequentist turn enables to leverage rich theory of vector-valued regression and establish under mild assumptions a convergence rate on the DMO estimator, with conditions fulfill minimax optimality in terms of statistiaal-compuuationnn efficiency. In the following, we briefly review CMO's vector-valued regression viewpoint and construct an analogous regression problem for DMO. We refer the reader [35 for a comprehensive overview of rector-valued RKHS theory. As for Sections ann general theory of DMEs and we hence put aside the bag notations.
In Section 3.2
[..2) this section contributes to the
Stage 1: Regressing the Conditional Mean Operator As first introduced by Grünewälder et al. 23 and generalised to infinite dimensional spaces by Singh et al. (25]. estimating C is equivalent solving vector-valued kernel ridge regression problem in the hypothesis space of Hilbert-Schmidt operators from Hr to Hr, denoted as HS (H. He). Specifically, we may consider the operator-valued kernel defined over H. as ['(f, f') == (f. f)x Idx, We denote Hr the He-valued RKHS spanned
by I' with norm || ||r, which can be identified to HS(H,, H. [ Singh et al. (25) frame CMO
regression as the minimisation surrogate discrepancy (C) == Ex,y [[]kx C Ey [12], to which they substitute empirical regularised version restricted to Mr given by
N
(C
Σ IR - tya - A|C
CE HE
A>0
(11)
N
This H4 Ha-valued kernel ridge regression problem admits closed-form minimiser which shares the same empirical form as the CMO, i.e. C arg minceH. (C) 23 25
'Hr Span (I, Ha, E He} Span {f @ f f E Hr, h E He) Hs & He 2 HS(Ha, Hz)

---

Stage 2: Regressing the Deconditional Mean Operator The DMO on the other hand is defined as the operator Dxly= H4- + 'He such that vf E H.. Dx D = Since deconditioning corresponds to finding a pseudo-inverse to the CMO, it natural to consider a reconstruction objective Ea(D) := Ey | DCx x fr Introducing a novel characterization of the DMO, we propose to minimise this objective in the hypothesis space of -chbert- operators HS(H,, He) which identifies to Hr. As per above, we denote ČxjY the empirical CMO learnt in Stage and we substitute the loss with an empirical regularised formulation on Hr
M
Ea(D) == M
1
Σ 114 DDXXXX l12 +cI DE Hr
E> 0
(12)
==1
Proposition 4.1 (Empirical DMO vector-valued regressor). The minimiser of the empirical reconstruction risk is the empirical DMO, Le. Dx argmin (a(
Since requires to estimate the CMO first, minimising 12 can be viewed two-staged vector value regression problem.
25 this perspective enables us to state Convergence results Following the footsteps of 24
the performance of the DMO estimate in terms of asymptotic convergence of the objective & As in Caponnetto and De Vito [36,, we must restrict the class of probability measure for Pxy and Py to ensure uniform convergence even when Hg is infinite dimensional. The family of distribution considered is general class of priors that does not assume parametric distributions and is parametrized by two variables: > controls the effective input dimension and e E]1, 2] controls functional smoothness. Mild regularity assumptions are also placed on the original spaces X, their corresponding RKHS Hgg Hyg and the vector-valuee RKHS Hr. We discuss these assumptions in details in Appendix [ Importantly, while Hgg can be infinite dimensional, we nonetheless have to assume the RKHS Hi is finite dimensional. In further research, we will relax this assumption.
Theorem 4.2 (Empirical DMO Convergence Rate). Denote Dry = arg minpeur Ea(D). Assume assumptions stated in Appendix Dare satisfied. In particular, let (b, c) and (0, c') be the parameters
of the restricted class of distribution for Py and IP Xx respectively and let EJO, 1] be the Hölder alc'=1) M Re'=1)
continuity exponent in Hr. Then, if we choose N 2t,, N the following result,
where > 0, we have
If a S
b(++)) then Ea(Dxly) EEa(De, O(M .ा with M न bot
bc+1) Ifa bet1 then &a (Dxxx) ED(D)y) = (MM with M पहfτ
This theorem underlines a trade-off between the computational and statistical efficiency with respect to the datasets cardinalities N = |D MM = |Dzl and the problem difficulty b, b, c'. c'. For S b(c+1) smaller means less samples from D2 at fixed M and thus computational savings. But it also hampers convergence, resulting in reduced statistical efficiency. At Mc-1)
bc+11
bo-1 < 2. convergence rate is minimax computaaionall-saaisiiiie efficiency optimal, i.e. convergence rate is optimal with smallest possible M. We note that at this optimal, N > -M and which means less samples are required from D2. a 2 b(+++))
be+11 does not improve the convergence rate but only increases the size of Dj and hence the computational cost it bears.
Deconditional Downscaling Experiments
We demonstrate and evaluate our CMP-based downscaling approaches on both synthetic exper- iments and a challenging atmospheric temperature field downscaling problem with unmatched multi-resolution data. We denote the exact CMP deconditional posterior from Section 3 as CMP, the CMP using with efficient shrinkage CMO estimation as S-CMP and the variational formulation as VARCMP. They are compared against VBAGG [9] which we describe below and GP regres- sion [4] baseline (GPR) modified to take bags centroids as the input. Experiments are implemented in PyTorch [37 38) all code and datasets are made available and computational details are provided Appendix

---

Variational Aggregate Learning VBAGG is introduced by Law et al. [9] as variational aggregate learning framework disaggregate exponential family models, with emphasis on the Poisson family, We consider its application to the Gaussian family, which models the relationship between aggregate targets 2, and bag covariates { J. by bag-wise averaging of GP prior on the function of interest. In fact, the Gaussian VBAGG corresponds exactly to special case of CMP on matched data, where the bag covariates are simply one hot encoded indices with kernel {(j.j') 8(j.j") where ó the Kronecker delta. However, VBAGG cannot handle unmatched data as bag indices do not instill the smoothness that is used for mediation. For fair analysis, we compare variational methods VARCMP and VBAGG together, and exact methods CMP/S-CMP to an exact version of VBAGG, which we implement and refer as BAGG-GP.
5.1 Swiss Roll
The scikit-learn [39] swiss roll manifold sampling function allows to generate a 3D manifold of points E IS mapped with their position along the manifold E R. Our objective will be to recover for each point x by only observing / at an aggregate level. In the first experiment, we compare our model to existing weakly supervised spatial disaggregation methods when all high righrresolution covariates matched with a coarse-level aggregate target. We then proceed to withdraw this requirement in a companion experiment.
5.1.1 Direct matching
Experimental Setup Inspired by the experimental setup from Law et al. 19, we regularly split space along height B times as depicted in Figure 2 and group together manifold points within each height level, hence mixing together points with very different positions on the manifold. We obtain bags of samples {( f
where the j"h bag contains n, points "x, = (z) 1 and
their corresponding targets Łú {t( (8) ) We then con- struct bag aggregate targets by taking noisy bag targets average 2, = - E (() + Ey, where Ey ~ NO(,o²). We thus obtain matched weakly supervised bag learn- ing dataset D = {(0z,, 24)] Since each bag corre- sponds to height-level, the center altitude of each height split B/ E R is a natural candidate bag-level covariate that informs on relative positions of the bags. We
Figure 2: Step 1: Split space regularly along height. Step 2: Group points into height-level bags. Step 3: Average points targets into bag-level aggregate targets.
augment the above dataset as = {( Mj =5) ==) Using these bag datasets, we now wish to downscale aggregate targets 2; to recover the unobserved manifold locations {t, and be able to query the target at any previously unseen input æ.
Models We use zero-mean prior on f and choose Gaussian kernel for h and l. Inducing points location is initialized with K-means++ procedure for VARCMP and VBAGG such that they spread evenly across the manifold. For exact methods, kernel hyperparameters and noise variance o² are learnt on by optimising the marginal likelihood. For VARCMP, they are learnt jointly with variational distribution parameters by maximising evidence lower bound objective. While CMP-based methods can leverage bag-level covariates y/j, baselines are restricted to learn from D 40 is used in all experiments.
Adam optimiser
Results We test models against unobserved groundtruth (t,) )m by evaluating the root mean square error (RMSE) to the posterior mean. Table shows that CMP, S-CMP and VARCMP outperform their corresponding counterparts i.e. BAGG-GP and VBAGG, with statistical significance confirmed by a Wilcoxon signed-rank test in Appendix E Most notably, this shows that the additional knowledge on bag-level dependence instilled by é reflected even setting where each bag is matched with an aggregate target.
8

---

Table 1: RMSE of the swissroll experiment for models trained over directly and indirectly matched
datasets: ; scores averaged over 20 seeds and s.d is reported * indicates our proposed methods.
Matching
CMP
S-CMP*
VARCMP*
BAGG-GP
VBAGG
GPR
Direct
0.33
0.25 2004
0.18
0.60
0.22
0.70+aas
Indirect
0.80
1.05
0.87+aor
1.13
1.46+034
1.04
5.1.2 Indirect matching
Experimental Setup We now impose indirect matching through mediating variable y/- We ran- domly select N [ bags which we consider to be the N first ones without loss of generality and split D into D, = [ My " and D2 = {(ūs, )= = {(UN+j, NN =N+;)]Ø=1 such that no pair of covariates bag "ry and aggregate target 2, are jointly observed.
Models CMP-based methods are naturally able to lear from this setting and are trained by inde- pendently drawing samples from D1 and D2. Baseline methods however require bags of covariates to be matched with an aggregate bag target. To remedy this, we place a separate prior ~ GP(0, E) and fit regression model 2 ((فق) ) + Ej over D2. We then the predictive posterior mean augment the first dataset as Dj = { (Lz,, Big(,, 3 This dataset can then be used to train BAGG-GP, VBAGG and GPR.
Results For comparison, we use the same evaluation as the direct matching experiment. Table ] underlines anticipated drop in RMSE for all models, but we observe that BAGG-GP and VBAGG suffer most from the mediated matching of the dataset while CMP and VARCMP report best scores by substantial margin. This highlights how using separate prior g to mediate Dj and D2 turns out to be suboptimal in contrast to using the prior naturally implied by CMP. While it is more computationally efficient than CMP, we observe relative drop in performance for S-CMP.
5.2 Mediated downscaling of atmospheric temperature
Given the large diversity of sources and formats of remote sensing and model data, expecting directly matched observations is often unrealistic [41]. For example, two distinct satellite products will often provide low and high resolution imagery that can be matched neither spatially nor temporally 2 3 5). Climate simulations [19 20 21
on the other hand provide a comprehensive coarse resolution coverage of meteorological variables that can serve as a mediating dataset.
For the purpose of demonstration, we create an experimental setup inspired by this problem using Coupled Model Intercomparison Project Phase 6 (CMIP6) [19] simulation data. This grants us access to groundtruth high-resolution covariates to facilitate model evaluation.
Experimental Setup We collect monthly mean 2D atmospheric fields simulation from CMIP6 data [4] 43 for the following variables: air temperature at cloud top (T), mean cloud top pressure (P), total cloud cover (TCC) and cloud albedo (or). First, we collocate TCC and a onto HR latitude- longitude grid of size 360x to obtain fine grained fields (latitude HR longitude altitude HR TCC ARR augmented with static HR surface altitude field. Then we collocate P and onto LR grid of size 21 x42 to obtain coarse grained fields (latitude. longitude pll. T..)). We denote by B the number of low resolution pixels.
Our objective is to disaggregate TLR to the HR fields granularity. We assimilate the jth coarse temperature pixel an aggregate target 2j := TLR corresponding to bag j. Each bag includes HR (6) )i = (Iatitude HR(() longitude HR((). altitude HR(S) HR(6) TCCHR()) .9)
HR()) ()) To emulate unmatched observations, we randomly select = _ of the bags ('x, ", and keep the opposite half of LR observations [ such that there is no single aggregate bag target that corresponds to one of the available bags. Finally, we choose the pressure field PLR as the mediating variable. We hence compose a third party low resolution field of bag-level covariates y; : (latitude],,, longitude/,R PIR) which can separately be matched with both above sets to obtain datasets D1 {(x,, 33) and D2 [(ū, 2
covariates fzz {
B-N
[(N + 2x+5)] B-N
9

---

Figure 3: Left: High-resolution atmoshperic covariates used for prediction; Center-Left: Observed now-resolution temperature field, grey pixels are unobserved; Center Unobserved high-resolution groundtruth temperature field; Center-Right: VARCMP deconditional posterior mean; Right 95% confidence region size on prediction; temperature values are in Kelvin.
280
260
240
TCC
at
Observed T
Groundtruth Deconditional Mean
95% Confidence region
Table 2: Downscaling similarity scores of posterior mean against groundtruth high resolution cloud top temperature field; averaged over seeds; we report s.d. "!": lower better ; "f": higher better.
Model
RMSE
MAE 4
Com:
SSIM
VARGPR
8.02:+0.28
5.55
0.831
0.212:#0.011
VBAGG
8.25+0015
5.82+t.11
0.801
0.182:#0.044
VARCMP
7.40
5.34
0.848
0.212
Models Setup We only consider variational methods to scale to large number of pixels. VARCMP naturally able learn from indirectly matched data. We use Matémn-1.5 kernel for rough spatial covariates (latitude, longitude) and Gaussian kernel for atmospheric covariates (P, TCC, a) and surface altitude. & and are both taken sums of Matérn and Gaussian kernels, and their hyperparameters are learnt along with noise variance during training. A high-resolution noise term also introduced, with details provided in Appendix E Inducing points locations are uniformly initialized across the HR grid. We replace GPR with an inducing point variational counterpart VARGPR [13]. Since baseline methods require matched dataset, we proceed as with the unmatched swiss roll experiment and fit GP regression model with kernel é on D2 and then use its predictive posterior mean to obtain pseudo-targets for the bags HR covariates from D
Results Performance is evaluated by comparing downscaling deconditional posterior mean against original high resolution field THR available CMIP6 data 43 which we emphasise never observed. We use random Fourier features [44] approximation of kernel & to scale kernel evaluation to the HR covariates grid during testing. As reported in Table VARCMP substantially outperforms both baselines with statistical significance provided in Appendix F the reconstructed image with VARCMP, plots for other methods are included the Appendix E The model resolves statistical patterns from HR covariates into coarse resolution temperature pixels, henceforth reconstructing faithful HR version the temperature field.
Discussion
introduced a scalable Bayesian solution the mediated statistical downscaling problem, which handles unmatched multi-resolution data. The proposed approach combines Gaussian Processes with the framework of deconditioning using RKHSs and recovers previous approaches as its special cases. We provided convergence rates for the associated deconditioning operator. Finally, demonstrated the advantages over spatial disaggregation baselines in synthetic data and a challenging atmospheric fields downscaling problem.
In future work, exploring theoretical guarantees of the computationally efficient shrinkage formulation in multi-resolution setting and relaxing finite dimensionality assumptions for the convergence rate will have fruitful practical and theoretical implications. Further directions also open up quantify uncertainty over the deconditional posterior since it is computed using empirical estimates of the CMP covariance. This may be problematic the mediating variable undergoes covariate shift between the two datasets.
10

---

Acknowledgments
SLC is supported by the EPSRC and MRC through the OxWaSP CDT programme EP/L0167 10/1. SB supported by the EU's Horizon 2020 research and innovation programme under Marie Skłodowska- Curie grant agreement No 860100. DS is supported in partly by Tencent Al lab and partly by the Alan Turing Institute (EP/N510129/1).
References
[1] Kelvin Hsu and Fabio Ramos. Bayesian Deconditional Kernel Mean Embeddings. Proceedings of Machine Learning Research. PMLR, 2019.
[2] I.. A. Remer, Y. J. Kaufman, D. Tanré, S. Mattoo, D. A. Chu, J. V. Martins, R.-R. Li, C. Ichoku, R. C. Levy, R. G. Kleidman, T. F. Eck, E. Vermote, and B. N. Holben. The MODIS Acrosol Algorithm, Products, and Validation. Journal of the Atmospheri Sciences, 2005.
[3] S. Platnick, M.D. King, S.A. Ackerman, W.P. Menzel, B.A. Baum, J.C. Riedi, and R.A. Frey. The MODIS cloud products: algorithms and examples from Terra. IEEE Transactions on Geoscience and Remote Sensing, 2003.
[4] Graeme L Stephens, Deborah G Vane, Ronald Boain, Gerald G Mace, Kenneth Sassen, Zhien Wang, Anthony J Illingworth, Ewan J O' comnor, William B Rossow, Stephen L Durden, et al. The CloudSat mission and the A-train: new dimension of space- based observations of clouds and precipitation. Bulletin of the American Meteorological Society, 2002.
[5] William L. Barnes, Thomas Pagano, and Vincent Salomonson. Prelaunch characteris- ties of the Moderate Resolution Imaging Spectroradiometer (MODIS) on EOS-AMI. IEEE Transactions on Geoscience and Remote Sensing, 1998.
[6] Zbyněk Malenovsky, Helmut Rott, Josef Cihlar, Michael E. Schaepman, Glenda Garefa-Santor Richard Fernandes, and Michael Berger. Sentinels for science: Potential of Sentinel-1, -2, and -3 missions for scientific observations of ocean, cryosphere, and land. Remote Sensing of Enviromment, 2012.
[7] Edward J. Knight and Geir Kvaran. Landsat-8 operational land imager design. Remote Sensing of Environment, 2014.
[8] Yivan Zhang, Nontawat Charoenphakdee, Zhenguo Wu, and Masashi Sugiyama. Learning from aggregate observations. In Advances in Neural Information Processing Systems, 2020.
[9] Leon Ho Chung Law, Dino Sejdinovic, Ewan Cameron, Tim C.D. Lucas, Seth Flaxman, Katherine Battle, and Kenji Fukumizu. Variational learning aggregate outputs with Gaussian processes. In Advances Neural Information Processing Systems, 2018
[10] Oliver Hamelijnck, Theodoros Damoulas, Kangrui Wang, and Mark A. Girolami. Multi- resolution multi-task Gaussian processes. In Advances in Neural Information Processing Systems, 2019.
[11] Fariba Yousefi, Michael Thomas Smith, and Mauricio A. Álvarez. Multi-task learning for aggregated data using Gaussian processes. In Advances in Neural Information Processing Systems, 2019.
[12] Yusuke Tanaka, Toshiyuki Tanaka, Tomoharu Iwata, Takeshi Kurashima, Maya Okawa, Yasunori Akagi, and Hiroyuki Toda. Spatially aggregated Gaussian processes with multivariate areal outputs. Advances in Neural Information Processing Systems, 2019.
[13] Arto Klam Ville Tanskanen, Krista Longi Non-Linearities Gaussian Processes with Integral Observations. IEEE international Workshop on Machine Learning for Signal, 2020.
[14] C Rasmussen and C Williams. Gaussian Processes for Machine Learning, 2005.
[15] Michalis K. Titsias. Variational learning of inducing variables in sparse gaussian processes. In AISTATS, 2009.
[16] Anna Vaughan, Will Tebbutt, J Scott Hosking, and Richard E Turner. Convolutional conditional neural processes for local climate downscaling. Geoscientific Model Development Discussions, pages 1-25, 2021.
11

---

[17] Brian Groenke, Luke Madaus, and Claire Monteleoni. Climalign: Unsupervised statistical downscaling of climate variables via normalizing flows. In Proceedings of the 10th International Conference Climate Informatics, pages 60-66, 2020.
[18) Thomas Vandal, Evan Kodra, Sangram Ganguly, Andrew Michaelis, Ramakrishna Nemani, and Auroop R Ganguly. Deepsd: Generating high resolution climate change projections through single image super-resolution. In Proceedings of the 23rd acm sigkdd international conference on knowledge discovery and data mining, pages 1663-1672, 2017.
[19] Veronika Eyring, Sandrine Bony, Gerald A Meehl, Catherine A Senior, Bjorn Stevens, Ronald J Stouffer, and Karl E. Taylor. Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization. Geoscientific Model Development, 2016.
[20]
Gregory M. Flato. Earth system models: an overview. WIREs Climate Change, 2011.
[21] Marko Scholze, J. Icarus Allen, William J. Collins, Sarah E. Cornell, Chris Huntingford, Manoj M. Joshi, Jason A. Lowe, Robin S. Smith, and Oliver Wild. Earth system models. Cambridge University Press, 2012.
[22] Krikamol Muandet, Kenji Fukumizu, Bharath Sriperumbudur, and Bernhard Schölkopf. Kernel mean embedding of distributions: review and beyond. arXiv preprint arXiv:1605.099
2016.
[23] Steffen Grünewälder, Guy Lever, Luca Baldassarre, Sam Patterson, Arthur Gretton, and Mas- similano Pontil. Conditional Mean Embeddings as Regressors. In Proceedings of the 29th International Coference on International Conference on Machine Learning, 2012.
[24] Zoltán Szabó, Bharath K. Sriperumbudur, Barnabás Póczos, and Arthur Gretton. Learning theory for distribution regression. Journal of Machine Learning Research, 2016.
[25] Rahul Singh, Maneesh Sahani, and Arthur Gretton. Kernel instrumental variable regression. In
Advances in Neural Information Processing Systems, 2019.
[26] Le Song, Kenji Fukumizu, and Arthur Gretton. Kernel embeddings of conditional distributions:
A unified kernel framework for nonparametric inference in graphical models. IEEE Signal Processing Magazine, 2013.
[27] Kenji Fukumizu, Francis R Bach, and Michael Jordan. Dimensionality reduction for supervised learning with reproducing kernel hilbert spaces. Journal of Machine Learning Research, 2004.
[28] Kenji Fukumizu, Arthur Gretton, Xiaohai Sun, and Bernhard Schölkopf. Kernel measures of conditional dependence. In Advances in Neural Information Processing Systems, 2008
[29] Bharath K. Sriperumbudur, Kenji Fukumizu, Arthur Gretton, Bernhard Schölkopf, and Gert R. G. Lanckriet. On the empirical estimation integral probability metries. Electronic Journal of Statistics, 2012.
[30] Le Song, Jonathan Huang, Alex Smola, and Kenji Fukumizu. Hilbert space embeddings of conditional distributions with applications to dynamical systems. In Proceedings of the 26th Annual International Conference on Machine Learning. 2009.
[31] FM Larkin. Gaussian measure in Hilbert space and applications in numerical analysis. The Rocky Mountain Journal of Mathematics, 1972.
[32] François- François-Xavier
Xavier Briol, Chris J Oates, Mark Girolami, Michael A Osborne, Dino Sejdinovic, et al. Probabilistic integration: A role in statistical computation? Statistical Science, 2019.
[33] Tim GJ Rudner, Dino Sejdinovic, and Yarin Gal. Inter-domain deep Gaussian processes. In International Conference on Machine Learning, 2020.
[34] Krikamol Muandet, Bharath Sriperumbudur, Kenji Fukumizu, Arthur Gretton, and Bernhard Schölkopf. Kernel mean shrinkage estimators. Journal of Machine Learning Research, 2016.
[35] Vern Paulsen and Mrinal Raghupathi. An introduction to the theory of reproducing kernel Hilbert spaces. Cambridge university press, 2016.
[36] Andrea Caponnetto and Ernesto De Vito. Optimal rates for the regularized least-squares algorithm. Foundations of Computational Mathematics, 2007.
[37] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary Vito,, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy,
12

---

Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. PyTorch: An Imperative Style, High-Performance Deep Learning Library. Advances in Neural Information Processing Systems 32. 2019.
[38] Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. Gpytorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration. In Advances in Neural Information Processing Systems, 2018.
[39] Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg. J. Vanderplas, A. Passos, D. Couraaeea, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 2011.
[40] Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization. In 3rd International Conference on Learning Representations, ICLR, 2015.
[41] D. Watson-Parris, N. Schutgens, N. Cook, Z. Kipling, P. Kershaw, E. Gryspeerdt, B. Lawrence, and P. Stier. Community Intercomparison Suite (CIS) 11.4.:: a tool for intercomparing models and observations. Geoscientific Model Development, 2016.
[42] Malcolm Roberts. MOHC HadGEMS-GC311-H model output prepared for CMIP6 High- ResMIP hist- 1950. Earth System Grid Federation, 2018 v20180730. doi: 10.22033/ESGF/
CMIP6.6040.
[43] Aurore Voldoire. CNRM-CERFACS CNRM-CM6-1-HE I-H model output prepared for CMIP6 HighResMIP hist-1950, Earth System Grid Federation, 2019 v20190221. doi: 10.22033/ESGF/ CMIP6.4040.
[44] Ali Rahimi, Benjamin Recht, et al. Random features for large-scale kernel machines. In Advances Neural Information Processing Systems, 2007.
[45] Olav Kallenberg. Foundations of Modern Probability. Springer, 2002.
[46] Charles A Micchelli and Massimiliano Pontil. On learning vector valued functions. Neural computation, 2005.
[47] Ho Chung Law, Peilin Zhao, Leung Sing Chan, Junzhou Huang, and Dino Sejdinovic. Hyper- parameter learning via distributional transfer. In Advances in Neural Information Processing Systems, 2019.
[48] Ingo Steinwart and Andreas Christmann. Support Vector Machines. Springer Publishing Company, Incorporated, 2008.
13


Reviews:
[{"review_id": "xM1Ew_0GK_4", "reviewer": "Reviewer_KjUw", "summary": "The authors consider the problem of refining low-res (LR) spatial observations using high-res (HR) covariates. \n\nBasically, the problem is as follows. The analyst observes LR outcome Z, HR covariate X, and LR covariate Y from the process Z=E[f(X)|Y]+e. The analyst wishes to recover the function f.\n\nIn more detail: the analyst observes two datasets: D1 and D2. D1 consists of N observations (x_j, y_j) where x_j may be a bag.  D2 consists of M observations (y_j, z_j). The function f is modelled as a function in an RKHS (with Bayesian or frequentist approach).\n\nThe authors characterize the Bayesian version of this problem and derive the GP posterior. The authors characterize the frequentist version of this problem and prove minimax optimal finite sample rates.\n", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Originality: It seems that the estimator previously appears in Hsu and Ramos. Prop 3.2 and Prop 3.3 appear to be new characterizations of the Bayesian approach. Prop 4.1 demonstrates that the estimator is the solution to a 2-stage (frequentist) regression. Theorem 4.2 relates the problem considered in this paper to the nonparametric instrumental variable problem in order to match symbols with an existing theorem.\n\nQuality: The theoretical results are complete and thorough. Why is H_\\ell finite dimensional? This seems to be unnecessary when appealing to the quoted theorem. There seems to be a small clash between the theory and empirics. In theory, k is characteristic and ell is finite. In empirics, k and ell are sums of Guassian and Matern kernels, so they are infinite dimensional and not characteristic.\n\nClarity: It took me a long time to understand the problem set up: what is observed, what is the goal, and what will be estimated. Figure 1 helps, and should be larger. Please improve the introduction to this learning problem. Theorem D.4 should be in the main text since otherwise its assumptions are hidden.\n\nSignificance: It is difficult for me to assess the significance of this paper in the spatial observation resolution literature. Within the RKHS literature, the theoretical contributions are modest: the estimator and main theorem are essentially from previous papers, but the Bayesian characterization and connection to nonparametric instrumental variable regression are new and insightful. I think these should be emphasized more. If there is a strong significance for practitioners, I will raise the score.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "What assumptions could be violated, and how would that impact conclusions and policy recommendations?", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "dDcgl4foe2w", "reviewer": "Reviewer_GawS", "summary": "The paper proposes a novel Gaussian process approach for settings in which we have two datasets, D_1={(x_i, y_i)}_i=1 ^N and D_2=((y_i,z_i))_i=1^M  that are linked by the mediating variable Y. The datasets are not matched, and we want to learn a function from X to Z. In contrast to existing work, the algorithm generalizes to settings in which the resolution of the two datasets differ.", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper comes along with a theoretical analysis of the proposed approach and a careful experimental validation. In addition, the paper is accomodated by a github repository that looks well documented and easy to run. While the paper seems to be theoretically sound, I found parts of it quiet confusing. However, this might also be attributed to my rather scarce background knowledge on conditional mean embeddings.\n\nThe title and abstract of the paper suggest that the main innovation of the paper is a new GP formulation that allows for different resolutions. However, for large parts of the paper, this seems to not be very relevant. In fact, the down-scaling only enters in Section 3.3 and in the experiments and is not even part of the main notation.\n\nThe major contribution of the paper seems to me is more the definition of the conditional mean process which allows for mediating variables in general. Interestingly,  one arrives at a GP formulation including analytical solutions for mean and covariance. This was surprising to me because chaining two Gaussian processes together and marginalizing out the intermediate variable leads in general to non-Gaussian distributions. \n\nThe posterior mean recovers a solution that has been published previously using the deconditional mean embedding [19].\nFor ease of readibility, I would suggest to move the complete deconditional mean embedding discussion in Section 4 and not scattered between Section 1, Section 2.2 and Section 4. It is not strictly needed for understanding the conditional mean process and distracts from it. \n\nI also wonder if the paper would not be better split into two separate contributions - one with the general theoretical framework and one with the specific down-scaling application in mind. For me, the two parts seem to be a bit disconnected.\n\nIn the experimental section, I would suggest to add a second metric that takes the uncertainty estimates of the GP into account. After all, this is one of the main benefits of the GP-based method. It would also be interesting to compare the uncertainty estimates between the proposed solution and the GP in [19] for the special case that the granularity levels are the same on both datasets.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The societal impact paragraph is sufficient.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "aJyU0b_gQBN", "reviewer": "Reviewer_sfJ5", "summary": "This paper presents a GP model for refining coarse-grained spatial data, which can be used for handling aggregated data that are unmatched spatially and temporally. The authors also develop the two-stage regression algorithm for downscaling. The effectiveness of the proposed model is demonstrated using synthetic and real-world datasets.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The task of recovering spatial fields from aggregated data is interesting and important. This paper is well-written and technically sound. The use of Conditional and Deconditional Kernel Mean Embeddings to the downscaling problem is novel. The experiments on atmospheric temperature data are encouraging. Accordingly, my opinion tends to accept; but there are some concerns about technical issues.\n\nThe detailed comments are as follows.\n\n(1) It seems that the input of the proposed method is limited to data aggregated at grid cells. But, in my understanding, the proposed method allows data aggregated over regions that have various shapes; is that right? If that is the case, it would be helpful to clarify it; this idea may be similar to those of [e.g., 10, R1]. \n\nActually, the spatially aggregated data (e.g., disease incidence [7]) are often associated with irregularly shaped regions such as districts. As in [10, R1], can the proposed model handle the effects from the region's size and shape? That could be critical for estimating prediction uncertainty. \n\n(2) I support you in that GP-based modeling is useful, but it would be better to mention the recently published downscaling methods based on DNN [e.g., R2] and clarify the advantages of the proposed approach.\n\n(3) In [6], learning from aggregated data is more generally formulated. Can the authors consider the extension of the proposed method to other machine learning tasks such as classification and ordinal regression?\n\n(4) Could you please elaborate on the reason why the predictive performance of BAGG-GP and VBAGG is worse than vanilla GP?\n\n[R1] M. T. Smith, M. A. Álvarez, and N. D. Lawrence. Gaussian process regression for binned data. In arXiv e-prints, 2018.\n\n[R2] Vandal, T.; Kodra, E.; Ganguly, S.; Michaelis, A.; Nemani, R.; and Ganguly, A. R. 2017. DeepSD: Generating high resolution climate change projections through single image super-resolution. In KDD, 1663–1672.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Please see comments (1), (2), and (3) in my main review.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "12 hours"}, "extra_text": {}}, {"review_id": "FDrpvek88g", "reviewer": "Reviewer_szhM", "summary": "The paper considers the problem of learning a function \n$$f:\\mathcal{X}\\to \\mathbb{R}$$ from the relationship \n$$g(y) = \\int_x f(x) \\mathbb{P}[x|y]dx$$\nwhere the training data consists of two datasets. The first set is  $\\mathcal{D}_1 = ( x^b_i, y_i )$ where each $x^b_i$ is a bag of points in $\\mathcal{X}$ and $y_i$ is a single point in a space $\\mathcal{Y}$, these are empirical samples from some distribution $\\mathbb{P}[x|y]\\mathbb{P}[y]$. The second dataset is of $\\mathcal{D}_2 = (\\tilde y_j, \\tilde z_i)$ with $\\tilde y_j \\in \\mathcal{Y}$ and $\\tilde z_i \\in \\mathbb{R}$ is a set of input-output pairs to from $\\tilde z_i = g(\\tilde y_i)$.\n\nLearning models for $f(x)$ from aggregate data, (input bag, output) pairs e.g. $( x^b_i, z_i )$ has been studied, similarly the use of  there are intermediate covariates $y_i \\in \\mathcal{Y}$ with a dataset of the form $( x^b_i, y_i, z_i )$ has been studied. As I undertsand, the novel constributino of this work is the extension to the case where there are two datasets, $( x^b_i, y_i)$ and $(\\tilde y_j, \\tilde z_j )$ that both contain points from the domain $\\mathcal{Y}$ yet they are not the same points, they are mismatched.\n\nThe paper proposes a \"conditional mean process\", a Gaussian process model that interpolates aggregated data as well as the marginalization operater in order to infer the un-aggregated latent function under this above problem setting. \n\n", "questions": "", "limitations": "", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I have two major concerns, the overall writing and the novelty of the proposed conditional mean process.\n\nRegarding the Conditional Mean Process\n\n  -  the integral of a Gaussian process is another Gaussian process is a standard result (eg convolution of a GP is another GP), I feel I have seen variations of equation (3) in multiple previous works [1] [2].\n\n  - I believe the following \"deconditioning\" step is standard (as rightly acknowledged to be Gaussian conditioning) and equations (5), (6) are textbook results.\n\n  - the sparse inducing point variational approach was applied in Hsu and Ramos\n\n  - I believe the main practical novelty mainly lies in using the kernel $l_y$ as a means to model $\\mathbb{P}[X|Y]$ away from the training data thereby bridging the gap between the unmatched $y_i$ and $\\tilde y_j$, I think this is really cool and the inclusion of VBAGG baseline (where $l_y$ is effectively \"switched off\") demonstrates the efficacy. However this seems rather unsurprising and rather incremental.\n\n[1] https://arxiv.org/abs/1902.07908\n\n[2] https://ieeexplore.ieee.org/abstract/document/9383978\n\n\nOverall I found the paper initially rather difficult to read, I felt the mathematical statements were very dense (mainly Section 2.2) with minimal explanation and required multiple passes to follow (I am largely familiar with the notation though not an expert). I had to read the paper multiple times to first acquire an intuitive understanding of the problem definition, I felt there were multiple subtle ambiguities that threw me off course and minor clarifications would have made the paper much easier to read.\n\nHere are some of my personal writing nitpicks.\n -  the first dataset is is introduced as $(x, y)$. This would traditionally imply that $x$ is predicting $y$, i.e. modelling $\\mathbb{P}[Y|X]$, Figure 1 seems to imply this as well, however $y$ is an \"input\" for the output $x^b$, it is used to Monte Carlo approximate $\\mathbb{P}[X|Y]$ instead. (I felt the introductory paragraph of Hsu and Ramos [19] was far easier to follow), perhaps using $y$ as function outputs (like [3]) as in  and $z$ (or anything else) as the \"middle man\" variables.\n\n - $z = \\mathbb{E}[f(X)|Y] + \\epsilon$ does not explicitly state what the expectation is over, though it slowly becomes clear that is is $X$.\n\n  - the word \"conditional\" is used very frequently and I believe it is referring to the _marginal_ over $X$ using $\\mathbb{P}[X|Y]$ (which happens to be a conditional distribution although $Y$ could be any parameter). Again upon first reading this led to a lot of ambiguity for me. (I would personally find \"de-marginalizing\" or \"de-aggregation\" to be much more intuitive than \"de-conditioning\")\n\n  - given the two datasets (x, y), (y, z) and the \"two stage regression\" statements, I was partly expecting deep Gaussian Processes, again, this misinterpretation took a long time to un-misinterpret\n\n  -  a \"hands on\" example application in the introduction with example $x$, $x$ and $z$ values would help a clear up much speculation from thee start avoiding many misinterpretations and blind alleys, upon first reading I didn't fully understand the problem setting until the experiments.\n\n\n\nI have given a low confidence score to reflect my difficulty in reading the paper.\n\n\n[3] Variational Learning on Aggregate Outputs with Gaussian Processes, Law et.al. NeurIPS 2019\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors do not appear to explicitly raise possible failure modes, below is one example.\n\n  - if $l_y$ is poorly chosen, this may cause the method to perform worse than VBAGG in matched data setting and worse than naive methods (like the imputation method on line 338) in the unmatched data setting.", "review_text": "", "extra_scores": {"ethical_concerns": "I have no ethical concerns.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "10"}, "extra_text": {}}, {"review_id": "C9xCmpfBTbd", "reviewer": "Reviewer_wwEN", "summary": "This paper extends on the work of deconditional mean embeddings (DMEs) and task transformed Gaussian processes (TTGPs) from Hsu and Ramos (2019) from both an application and theoretical stand point.\n\nFrom a theoretical stand point, the authors developed the framework around deconditional mean embeddings further by three ways: (1) formulating another elegant way to arrive at deconditional posteriors from Hsu and Ramos (2019), (2) establishing deconditional mean operators (DMOs) as vector-valued regressors in a manner that mirrors Grunewalder et al (2012) which highlights its interpretation as a reconstruction operator, and (3) enhancing guarantees on the convergence rate of deconditional mean operators from what was established in Hsu and Ramos (2019) using the setup and results from  Caponnetto and De Vito (2007), Szabó et al (2016), and Singh et al (2019).\n\nFrom an application stand point, the authors apply DMEs for refining low resolution spatial fields with high resolution information. This is the special case of DMEs and TTGPs where the task dataset consists of collections of low resolution covariates ($\\tilde{y}$) and aggregated targets ($\\tilde{z}$), and the transformation dataset consists of collections of (another potentially unmatched set of) low resolution covariates ($y$) each matched with bags of high resolution covariates ($^{b}\\bf{x}$). The latter constructions regarding the transformation dataset is what makes this a non-trivial special case due to having bagged high resolution covariates, which allows for a slightly different empirical estimator for the cross-covariance operator. This leads to an alternative conditional mean operator (CMO) which they call the conditional mean shrinkage operator. Finally, the authors then apply their work to toy experiments (swiss roll) and downscaling of atmospheric temperature (CMIP6), where to scale it they also derive a variational formulation to approximate the deconditional posterior.\n", "questions": "", "limitations": "", "rating": 7, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Overall, the contributions of the paper is well motivated and well developed. The paper presents non-trivial extensions to the original DME paper on a theoretical level, and at the same time successfully showcases a particularly useful type of application for the theory developed through refining low-resolution spatial fields with high-resolution information.\n\nIn particular, of all the theoretical contributions, I find the first contribution on an alternative view to arrive at deconditional posteriors using conditional mean processes very simple and elegant. In Hsu and Ramos 2019, this posterior was derived by defining and performing inference on the graphical model of the task transformed Gaussian process. Instead, this paper arrives at the posterior by first considering the conditional mean process constructed by taking the conditional mean of an original latent Gaussian process which we would like to infer. Due to linearity of the conditional mean, the resulting process is also a Gaussian process, and the authors derive their mean and covariance functions. This naturally leads to the notion of a joint Gaussian field between noisy observations of the resulting Gaussian process after the conditional mean and the original latent Gaussian process, which means the deconditional posterior can be easily derived from Gaussian conditioning. The authors show that this recovers the same deconditional posterior, also termed the posterior of the task transformed Gaussian process, from Hsu and Ramos (2019), establishing another elegant view for interpreting deconditional posteriors. An added bonus is that the authors also explicitly kept a non-zero mean function throughout this derivation.\n\nThe second and third contributions are also helpful in the sense that it casts the deconditional problem as a reconstruction problem, which provides a different two-staged reconstruction regression view to the two-staged task transformed regression view originally discussed in Hsu and Ramos (2019). The theorems developed for capturing the convergence rate under this reconstruction loss is well appreciated and indeed non-trivial.\n\nThe main weakness of the paper predominantly lies in its presentation, which hopefully should be easy to fix. Understandably, the presentation is dense due to the amount of content and contributions the authors are trying to communicate. However this is especially where the structure and presentation of the paper can vastly improve the communication of the core contributions of the paper. Unfortunately, as it currently stands it can be difficult to appreciate the contributions and their relationships all at once especially because they operate at different levels of generality.\n\nFor instance, all three theoretical contributions (as summarised in my summary section) are in fact highly general to the theory of deconditional mean embeddings, and not limited to the core application discussed in the paper, which is refining low-resolution spatial fields with high-resolution information. To me, these contributions are of higher significance compared to the specific application of focus, which is a special case (albeit non-trivial) to the original DME and TTGP setup. For the author's benefit, I would suggest to emphasize the former more than the latter, in contrast to how it is currently communicated. At the very least, it seems strange that the title, abstract, and introduction focuses heavily on setting up this specific application, while the majority of the remainder of the paper develops or extends a theorical framework that is much more general than the application aforementioned.\n\nFurthering this point with a specific example, section 3.3 seems quite out of place in that this is the only sub-section in the major contribution sections (sections 3 and 4) where the multi-resolution application setting comes into play. It seems more appropriate to delay this to Section 5 instead when the specific application is discussed, and keep section 3 and 4 focused on the more general theoretical developments.\n\n[Update] In light of the shared concerns on the paper's clarity and presentation which I agree with as initially noted, I have revised my score to give more weight to this concern, although this does not change my overall stance on the paper.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The work is motivated by applications in statistical downscaling with unmatched multi-resolution data. However, its contributions are general in nature in terms of the technique and not limited to this application. In fact, as highlighted in my review, I believe the general theoretical contributions are more significant than the specific application discussed.\n\nThe paper does address societal impact in the submission, citing positive benefits in understanding climate phenomenon from unmatched multi-resolution data which can be quite common, as well as potential negative impacts in using this technique to recover sensitive information from low resolution to high resolution where privacy is a concern.\n\nOverall, I agree with the author's self assessment on limitations and societal impact of the paper.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "7"}, "extra_text": {}}]

Author rebuttal:

Meta-review: 
After the rebuttal and discussion, the consensus among the high-confidence reviewers is that the paper should be accepted. The meta-reviewer agrees. On the theoretical side, the paper makes valuable contributions that clearly go beyond previous work. On the practical side, the application in climate science is valuable.

For the community to fully benefit from the paper, serious efforts need to go into making the paper more accessible. First, the explanation of the problem setup and the background needs to be improved to make the paper accessible to non-experts. In addition to the advice given by the reviewers, which needs to be included, the authors could make use of the appendix to further clarify the background material and previous work. Second, as promised by the authors in the rebuttal, the presentation of the two theoretical contributions to (1) the theory of deconditional mean embeddings and (2) statistical downscaling needs to be improved.



Decision: Accept (Poster)