summaries: This paper presents GraphRetro, a single-step retrosynthesis model that represents an advance from last years' NeurIPS state of the art model on semi-template-based modeling of the USPTO-50k dataset. GraphRetro uses a graph neural network that predicts edits to transform a target into synthons and then expands the synthons into full molecules by attaching leaving groups (with the synthons and leaving groups learned from the training set). The authors compare the predictions on the USPTO-50k to existing baselines, and show examples of predictions, including wrong ones. Because of the particular construction of the model, the authors can decompose the sources of errors as wrong edit predictions or wrong synthon completions.

summaries: ## Summary.  The authors propose a new method for single-step retrosynthesis prediction. Their method first predicts graph edits using a message passing network yielding synthons which are then completed by adding leaving groups. The leaving groups are chosen from a fixed set using another MPN.

summaries: This paper deals with molecular graph generation via synthon, incomplete molecular graphs.  Unlike existing works, this paper formulates the completion of incomplete molecular graphs not as a generation problem but as a classification problem to select one of a set of pre-computed leaving groups fragments. The training and evaluation are performed using the standard dataset USPTO-50k. The proposed model has high retrosynthesis performance, covering a large number of responses with a small set of leaving groups.

summaries: This paper proposes a new retrosynthetic prediction model using neural networks. The algorithm proposed to make the prediction based on breaking the target product into synthons and adding leaving groups to the synthons. The main difference of the proposed algorithm to the existing synthon-based retrosynthetic model is on introducing the vocabulary of leaving groups to complete the synthons in an efficient way.  Extensive evaluation verifies that the proposed method is competitive with the existing algorithms, and each algorithmic component plays a meaningful role.

summaries: This paper proposes a new strategy for reactants generation for the retrosynthesis prediction. Motivated by the fact that the graph topology of precursor molecules is largely unchanged during the chemical reaction, the authors propose to expand intermediate synthons into reactants by selecting and attaching leaving groups (molecule sub-graphs extracted from training data) to the synthons. The product molecule is transformed into synthons through graph edits. The proposed method achieves 53.7% top-1 accuracy.

abstractText: Retrosynthesis prediction is a fundamental problem in organic synthesis, where the task is to identify precursor molecules that can be used to synthesize a target molecule. A key consideration in building neural models for this task is aligning model design with strategies adopted by chemists. Building on this viewpoint, this paper introduces a graph-based approach that capitalizes on the idea that the graph topology of precursor molecules is largely unaltered during a chemical reaction. The model first predicts the set of graph edits transforming the target into incomplete molecules called synthons. Next, the model learns to expand synthons into complete molecules by attaching relevant leaving groups. This decomposition simplifies the architecture, making its predictions more interpretable, and also amenable to manual correction. Our model achieves a top-1 accuracy of 53.7%, outperforming previous template-free and semi-template-based methods.

authors:
B. Chen
T. Shen
T.S. Jaakkola
R. Barzilay

title: Learning to Make Generalizable and Diverse Predictions for Retrosynthesis

venue: In Submission,

authors:
C.W. Coley
R. Barzilay
T.S. Jaakkola
W.H. Green
K.F. Jensen

title: Prediction of Organic Reaction Outcomes Using Machine Learning

venue: In ACS Central Science. ACS Publications,

authors:
C.W. Coley
L. Rogers
W.H. Green
K.F. Jensen

title: Computer-Assisted Retrosynthesis Based on Molecular Similarity

venue: ACS Central Science,

authors:
C.W. Coley
W. Jin
L. Rogers
T.F. Jamison
T.S. Jaakkola
W.H. Green
R. Barzilay
K.F. Jensen

title: A graph-convolutional neural network model for the prediction of chemical reactivity

venue: Chemical Science,

authors:
E. Corey
W.T. Wipke

title: Computer-assisted design of complex organic syntheses

venue: Science,

title: The Logic of Chemical Synthesis: Multistep Synthesis of Complex Carbogenic Molecules (Nobel Lecture)

venue: Angewandte Chemie International Edition,

authors:
H. Dai
C. Li
C. Coley
B. Dai
L. Song

title: Retrosynthesis Prediction with Conditional Graph Logic Network

venue: In Advances in Neural Information Processing Systems (NeurIPS),

authors:
S. Genheden
A. Thakkar
V. Chadimová
J.-L. Reymond
O. Engkvist
E. Bjerrum

title: Aizynthfinder: a fast, robust and flexible open-source software for retrosynthetic planning

venue: Journal of cheminformatics,

authors:
J. Gilmer
S.S. Schoenholz
P.F. Riley
O. Vinyals
G.E. Dahl

title: Neural Message Passing for Quantum Chemistry

venue: In International Conference on Machine Learning (ICML),

authors:
M. Hartenfeller
M. Eberle
P. Meier
C. Nieto-Oberhuber
K.-H. Altmann
G. Schneider
E. Jacoby
S. Renner

title: A Collection of Robust Organic Synthesis Reactions for In Silico Molecule Design

venue: In Journal of Chemical Information and Modeling,

authors:
W. Jin
C. Coley
R. Barzilay
T. Jaakkola

title: Predicting Organic Reaction Outcomes with WeisfeilerLehman Network

authors:
W. Jin
R. Barzilay
T. Jaakkola

title: Junction Tree Variational Autoencoder for Molecular Graph Generation

title: Composing Molecules with Multiple Property Constraints

authors:
J. Law
Z. Zsoldos
A. Simon
D. Reid
Y. Liu
S.Y. Khew
A.P. Johnson
S. Major
R.A. Wade
H.Y. Ando

title: Route Designer: A Retrosynthetic Analysis Tool Utilizing Automated Retrosynthetic Rule Generation

authors:
B. Liu
B. Ramsundar
P. Kawthekar
J. Shi
J. Gomes
Q. Luu Nguyen
S. Ho
J. Sloane
P. Wender
V. Pande

title: Retrosynthetic Reaction Prediction Using Neural Sequence-to-Sequence Models

venue: In ACS Central Science,

authors:
N. Schneider
N. Stiefl
G.A. Landrum

title: What’s What: The (Nearly) Definitive Guide to Reaction Role Assignment

authors:
M.H. Segler
M.P. Waller

title: Neural-Symbolic Machine Learning for Retrosynthesis and Reaction Prediction

venue: Chemistry–A European Journal,

authors:
C. Shi
M. Xu
H. Guo
M. Zhang
J. Tang

title: A graph to graphs framework for retrosynthesis prediction, 2020

authors:
R. Sun
H. Dai
L. Li
S. Kearnes
B. Dai

title: Energy-based view of retrosynthesis, 2021

authors:
S. Szymkuć
E.P. Gajewska
T. Klucznik
K. Molga
P. Dittwald
M. Startek
M. Bajczyk
B.A. Grzybowski

title: Computer-assisted synthetic planning: The end of the beginning

authors:
A. Vaswani
N. Shazeer
N. Parmar
J. Uszkoreit
L. Jones
A.N. Gomez
Ł. Kaiser
I. Polosukhin

title: Attention is All You Need

title: SMILES, a Chemical Language and Information System

venue: Journal of Chemical Information and Computer Sciences,

authors:
R.J. Williams
D. Zipser

title: A Learning Algorithm for Continually Running Fully Recurrent Neural Networks

venue: In Neural Computation,

authors:
C. Yan
Q. Ding
P. Zhao
S. Zheng
J. YANG
Y. Yu
J. Huang

title: Retroxpert: Decompose retrosynthesis prediction like a chemist

venue: Advances in Neural Information Processing Systems,

authors:
S. Zheng
J. Rao
Z. Zhang
J. Xu
Y. Yang

title: Predicting Retrosynthetic Reactions using SelfCorrected Transformer Neural Networks

venue: In Journal of Chemical Information and Modeling. ACS Publications,

heading: 1 Introduction

text:
Retrosynthesis prediction, first formalized by E. J. Corey [Corey, 1991] is a fundamental problem in organic synthesis that attempts to identify a series of chemical transformations for synthesizing a target molecule. In the single-step formulation, the task is to identify a set of reactant molecules given a target. Beyond simple reactions, many practical tasks involving complex organic molecules are difficult even for expert chemists. As a result, substantial experimental exploration is needed to cover for deficiencies of analytical approaches. This has motivated interest in computer-assisted retrosynthesis [Corey and Wipke, 1969], with a recent surge in machine learning methods [Chen et al., 2019, Coley et al., 2017b, Dai et al., 2019, Zheng et al., 2019, Genheden et al., 2020].
Computationally, the main challenge is how to explore the combinatorial space of reactions that can yield the target molecule. Largely, previous methods for retrosynthesis prediction can be divided into template-based [Coley et al., 2017b, Dai et al., 2019, Segler and Waller, 2017] and template-free [Chen et al., 2019, Zheng et al., 2019] approaches. Template-based methods match a target molecule against a large set of templates, which are molecular subgraph patterns that highlight changes during a chemical reaction. Despite their interpretability, these methods fail to generalize to new reactions. Template-free methods bypass templates by learning a direct mapping from the SMILES [Weininger, 1988] representations of the product to reactants. Despite their greater generalization potential, these methods generate reactant SMILES character by character, increasing generation complexity.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Another important consideration in building retrosynthesis models is aligning model design with strategies adopted by expert chemists. These strategies are influenced by fundamental properties of chemical reactions, independent of complexity level: (i.) the product atoms are always a subset of the reactant atoms1, and (ii.) the molecular graph topology is largely unaltered from products to reactants. For example, in the standard retrosynthesis dataset, only 6.3% of the atoms in the product undergo any change in connectivity.
This consideration has received more attention in recent semi-template-based methods [Shi et al., 2020, Yan et al., 2020], that generate reactants from a product in two stages: (i.) first identify intermediate molecules called synthons, (ii.) and then complete synthons into reactants by sequential generation of atoms or SMILES characters.. Our model GRAPHRETRO also uses a similar workflow. However, we avoid sequential generation for completing synthons by instead selecting subgraphs called leaving groups from a precomputed vocabulary. This vocabulary is constructed during preprocessing by extracting subgraphs that differ between a synthon and the corresponding reactant. The vocabulary has a small size (170 for USPTO-50k) indicating remarkable redundancy, while covering 99.7% of the test set. Operating at the level of these subgraphs greatly reduces the complexity of reactant generation, with improved empirical performance. This formulation also simplifies our architecture, and makes our predictions more transparent, interpretable and amenable to manual correction.
The benchmark dataset for evaluating retrosynthesis models is USPTO-50k [Schneider et al., 2016], which consists of 50000 reactions across 10 reaction classes. The dataset contains an unexpected shortcut towards predicting the edit, in that the product atom with atom-mapping 1 is part of the edit in 75% of the cases, allowing predictions that depend on the position of the atom to overestimate performance. We canonicalize the product SMILES and remap the existing dataset, thereby removing the shortcut. On this remapped dataset, GRAPHRETRO achieves a top-1 accuracy of 53.7% when the reaction class is not known, outperforming both template-free and semi-template-based methods.

heading: 2 Related Work

text:
Retrosynthesis Prediction Existing machine learning methods for retrosynthesis prediction can be divided into template-based, template-free and recent semi-template-based approaches.
Template-Based: Templates are either hand-crafted by experts [Hartenfeller et al., 2011, Szymkuć et al., 2016], or extracted algorithmically from large databases Coley et al. [2017a], Law et al. [2009]. Exhaustively applying large template sets is expensive due to the involved subgraph
1ignoring impurities
matching procedure. Template-based methods therefore utilize different ways of prioritizing templates, by either learning a conditional distribution over the template set [Segler and Waller, 2017], ranking templates based on molecular similarities to precedent reactions [Coley et al., 2017b] or directly modelling the joint distribution of templates and reactants using logic variables [Dai et al., 2019]. Despite their interpretability, these methods fail to generalize outside their rule set.
Template-Free: Template-free methods [Liu et al., 2017, Zheng et al., 2019, Chen et al., 2019] learn a direct transformation from products to reactants using architectures from neural machine translation and a string based representation of molecules called SMILES [Weininger, 1988]. Linearizing molecules as strings does not utilize the inherently rich chemical structure. In addition, the reactant SMILES are generated from scratch, character by character. Attempts have been made to improve validity by adding a syntax correcter [Zheng et al., 2019] and a mixture model to improve diversity of suggestions [Chen et al., 2019], but the performance remains worse than [Dai et al., 2019] on the standard retrosynthesis dataset. Sun et al. [2021] formulate retrosynthesis using energy-based models, with additional parameterizations and loss terms to enforce the duality between forward (reaction prediction) and backward (retrosynthesis) prediction.
Semi-Template-Based: Our work is closely related to recently proposed semi-template-based methods [Shi et al., 2020, Yan et al., 2020], which first identify synthons and then expand synthons into reactants through sequential generation using either a graph generative model [Shi et al., 2020] or a Transformer [Yan et al., 2020]. To reduce the complexity of reactant generation, we instead complete synthons using subgraphs called leaving groups selected from a precomputed vocabulary. This allows us to view synthon completion as a classification problem instead of a generative one. We also utilize the dependency graph between possible edits, and update edit predictions using a message passing network (MPN) [Gilmer et al., 2017] on this graph. Both innovations together yield a 4.8% and 3.3% performance improvement respectively over previous semi-template-based methods.
Reaction Center Identification The reaction center covers a small number of participating atoms involved in the reaction. Our work is also related to models that predict reaction outcomes by learning to rank atom pairs based on their likelihood to be in the reaction center [Coley et al., 2019, Jin et al., 2017]. The task of identifying the reaction center is related to the step of deriving the synthons in our formulation. Our work departs from [Coley et al., 2019, Jin et al., 2017] as we utilize the property that new bond formations occur rarely (~0.1%) from products to synthons, allowing us to predict a score only for existing bonds and atoms and reduce prediction complexity from O(N2) to O(N). We also utilize the dependency graph between possible edits, and update edit predictions using a MPN on this graph.
Utilizing Substructures Substructures have been utilized in various tasks from sentence generation by fusing phrases to molecule generation and optimization [Jin et al., 2018, 2020]. Our work is closely related to [Jin et al., 2020] which uses precomputed substructures as building blocks for property-conditioned molecule generation. However, instead of precomputing, synthons —analogous building blocks for reactants— are indirectly learnt during training.

heading: 3 Model Design

text:
Our approach leverages the property that graph topology is largely unaltered from products to reactants. To achieve this, we first derive suitable building blocks from the product called synthons, and then complete them into valid reactants by adding specific functionalities called leaving groups. These derivations, called edits, are characterized by modifications to bonds or hydrogen counts on atoms. We first train a neural network to predict a score for possible edits (Section 3.1). The edit with the highest score is then applied to the product to obtain synthons. Since the number of unique leaving groups are small, we model leaving group selection as a classification problem over a precomputed vocabulary (Section 3.2). To produce candidate reactants, we attach the predicted leaving group to the corresponding synthon through chemically constrained rules. The overall process is outlined in Figure 1. Before describing the two modules, we introduce relevant preliminaries that set the background for the remainder of the paper.
Retrosynthesis Prediction A retrosynthesis pair R is described by a pair of molecular graphs (Gp,Gr), where Gp are the products and Gr the reactants. A molecular graph is described as G = (V, E) with atoms V as nodes and bonds E as edges. Prior work has focused on the single product case, while reactants can have multiple connected components, i.e. Gr = {Grc}Cc=1. Retrosynthesis pairs are atom-mapped so that each product atom has a unique corresponding reactant atom. The retrosynthesis task then, is to infer {Grc}Cc=1 given Gp.
Edits Edits consist of (i.) atom pairs {(ai, aj)} where the bond type changes from products to reactants, and (ii.) atoms {ai} where the number of hydrogens attached to the atom change from products to reactants. We denote the set of edits by E. Since retrosynthesis pairs in the training set are atom-mapped, edits can be automatically identified by comparing the atoms and atom pairs in the product to their corresponding reactant counterparts.
Synthons and Leaving Groups Applying editsE to the product Gp results in incomplete molecules called synthons. Synthons are analogous to rationales or building blocks, which are expanded into valid reactants by adding specific functionalities called leaving groups that are responsible for its reactivity. We denote synthons by Gs and leaving groups by Gl. We further assume that synthons and leaving groups have the same number of connected components as the reactants, i.e Gs = {Gsc}Cc=1 and Gl = {Glc}Cc=1. This assumption holds for 99.97% reactions in the training set. Formally, our model generates reactants by first predicting the set of edits E that transform Gp into Gs, followed by predicting a leaving group Glc to attach to each synthon Gsc . The model is defined as
P (Gr|Gp) = ∑ E,Gl P (E|Gp)P (Gl|Gp,Gs), (1)
where Gs,Gr are deterministic given E,Gl, and Gp.

heading: 3.1 Edit Prediction

text:
For a given retrosynthesis pair R = (Gp,Gr), we predict an edit score only for existing bonds and atoms, instead of every atom pair as in [Coley et al., 2019, Jin et al., 2017]. This choice is motivated by the low frequency (~0.1%) of new bond formations in the training set examples. Coupled with the sparsity of molecular graphs, this reduces the prediction complexity from O(N2) to O(N) for a product with N atoms. Our edit prediction model has variants tailored to single and multiple edit prediction. Since 95% of the training set consists of single edit examples, the remainder of this section describes the setup for single edit prediction. A detailed description of our multiple edit prediction model can be found in Appendix ??.
Each bond (u, v) in Gp is associated with a label yuvk ∈ {0, 1} indicating whether its bond type k has changed from the products to reactants. Each atom u is associated with a label yu ∈ {0, 1} indicating a change in hydrogen count. We predict edit scores using representations that are learnt using a graph encoder.
Graph Encoder To obtain atom representations, we use a variant of the message passing network (MPN) described in [Gilmer et al., 2017]. Each atom u has a feature vector xu indicating its atom type, degree and other properties. Each bond (u, v) has a feature vector xuv indicating its aromaticity, bond type and ring membership. For simplicity, we denote the encoding process by MPN(·) and describe architectural details in Appendix ??. The MPN computes atom representations {cu|u ∈ G} via
{cu} = MPN(G, {xu}, {xuv}v∈N (u)), (2)
where N (u) denotes the neighbors of atom u. The graph representation cG is an aggregation of atom representations, i.e. cG = ∑ u∈V cu. When G has connected components {Gi}, we get a set of graph representations {cGi}. For a bond (u, v), we define its representation cuv = (ABS(cu, cv)||cu+cv), where ABS denotes absolute difference and || refers to concatenation. This ensures our representations are permutation invariant. These representations are then used to predict atom and bond edit scores using corresponding neural networks,
su = ua T τ(Wacu + b) (3)
suvk = uk T τ(Wkcuv + bk), (4)
where τ(·) is the ReLU activation function.
Updating Bond Edit Scores Unlike a typical classification problem where the labels are independent, edits can have possible dependencies between each other. For example, bonds part of a stable system such as an aromatic ring have a greater tendency to remain unchanged (label 0). We attempt to leverage such dependencies to update initial edit scores. To this end, we build a graph with bonds (u, v) as nodes, and introduce an edge between bonds sharing an atom. We use another MPN(·) on this graph to learn aggregated neighborhood messages muv, and update the edit scores suvk in a manner similar to how LSTMs update representations,
fuvk = σ(W f kxxuv +W f kmmuv) (5)
iuvk = σ(W i kxxuv +W i kmmuv) (6)
m̃uvk = umτ(W m kxxuv +W m kmmuv) (7)
s̃uvk = fuvk · suvk + iuvk · m̃uvk. (8)
Training We train by minimizing the cross-entropy loss over possible bond and atom edits
Le = − ∑
(Gp,E)  ∑ ((u,v),k)∈E yuvklog(s̃uvk) + ∑ u∈E yulog(su)  . (9) The cross-entropy loss enforces the model to learn a distribution over possible edits instead of reasoning about each edit independently, as with the binary cross entropy loss used in [Jin et al., 2017, Coley et al., 2019].

heading: 3.2 Synthon Completion

text:
Synthons are completed into valid reactants by adding specific functionalities called leaving groups. This involves two complementary tasks: (i.) selecting the appropriate leaving group, and (ii.) attaching the leaving group to the synthon. As ground truth leaving groups are not directly provided, we extract the leaving groups and construct a vocabulary X of unique leaving groups during preprocessing. The vocabulary has a limited size (|X | = 170 for a standard dataset with 50, 000 examples, and 72000 synthons) indicating the redundancy of leaving groups used in accomplishing retrosynthetic transformations. This redundancy also allows us to formulate leaving group selection as a classification problem over X , while retaining the ability to generate diverse reactants using different combinations of leaving groups.
Vocabulary Construction Before constructing the vocabulary, we align connected components of synthon and reactant graphs by comparing atom mapping overlaps. Using aligned pairs Gsc = (Vsc , Esc) and Grc = (Vrc , Erc) as input, the leaving group vocabulary X is constructed by extracting subgraphs Glc = (Vlc , Elc) such that Vlc = Vrc \ Vsc . Atoms {ai} in the leaving groups that attach to synthons are marked with a special symbol. We also add three tokens to X namely START, which indicates the start of synthon completion, END, which indicates that there is no leaving group to add and PAD, which is used to handle variable numbers of synthon components in a minibatch.
Leaving Group Selection For synthon component c ≤ C, where C is the number of connected components in the synthon graph, we use three inputs for leaving group selection – the product representation cGp , the synthon component representation cGsc , and the leaving group representation for the previous synthon component, elc−1 . The product and synthon representations are learnt using the MPN(·). For each xi ∈ X , representations can be learnt by either training independent embedding vectors (ind) or by treating each xi as a subgraph and using the MPN(·) (shared). In the shared setting, we use the same MPN(·) as the product and synthons.
The leaving group probabilities are then computed by combining cGp , cGsc and elc−1 via a single layer neural network and softmax function
q̂lc = softmax ( Uτ ( W1cGp +W2cGsc +W3el(c−1) )) , (10)
where q̂lc is distribution learnt over X . Using the representation of the previous leaving group elc−1 allows the model to understand combinations of leaving groups that generate the desired product from the reactants. We also include the product representation cGp as the synthon graphs are derived from the product graph.
Training For step c, given the one hot encoding of the true leaving group qlc , we minimize the cross-entropy loss
Ls = C∑
c=1
L(q̂lc , qlc). (11)
Training utilizes teacher-forcing [Williams and Zipser, 1989] so that the model makes predictions given correct histories. During inference, at every step, we use the representation of leaving group from the previous step with the highest predicted probability.
Leaving Group Attachment Attaching leaving groups to synthons is a deterministic process and not learnt during training. The task involves identification of the type of bonds to add between attaching atoms in the leaving group (marked during vocabulary construction), and the atom(s) participating in the edit. These bonds can be inferred by applying the valency constraint, which determines the maximum number of neighbors for each atom. The attachment process does not modify any stereochemistry. Given synthons and leaving groups, the attachment process has a 100% accuracy. The detailed procedure is described in Appendix ??.

heading: 3.3 Inference

text: Inference is performed using beam search with a log-likelihood scoring function. For a beam width n, we select n edits with highest scores and apply them to the product to obtain n synthons, where each synthon can consist of multiple connected components. The synthons form the nodes for beam search. Each node maintains a cumulative score by aggregating the log-likelihoods of the edit and predicted leaving groups. Leaving group inference starts with a connected component for each synthon, and selects n leaving groups with highest log-likelihoods. From the n2 possibilities, we select n nodes with the highest cumulative scores. This process is repeated until all nodes have a leaving group predicted for each synthon component.

heading: 4 Evaluation

text:
Evaluating retrosynthesis models is challenging as multiple sets of reactants can be generated from the same product. To deal with this, previous works [Coley et al., 2017b, Dai et al., 2019] evaluate the ability of the model to recover retrosynthetic strategies recorded in the dataset.
Data We use the benchmark dataset USPTO-50k [Schneider et al., 2016] for all our experiments. The dataset contains 50, 000 atom-mapped reactions across 10 reaction classes. We use the same dataset version and splits as provided by [Dai et al., 2019]. The USPTO-50k dataset contains a shortcut in that the product atom with atom-mapping 1 is part of the edit in ~75% of the cases. If the product SMILES is not canonicalized, predictions utilizing operations that depend on the position of the atom or bond will be able to use the shortcut, and overestimate performance. We canonicalize the product SMILES, and reassign atom-mappings to the reactant atoms based on the canonical ordering, which removes the shortcut. Details on the remapping procedure can be found in Appendix ??.
Evaluation We use the top-n accuracy (n = 1, 3, 5, 10) as our evaluation metric, defined as the fraction of examples where the recorded reactants are suggested by the model with rank ≤ n. Following prior work [Coley et al., 2017b, Zheng et al., 2019, Dai et al., 2019], we compute the accuracy
by comparing the canonical SMILES of predicted reactants to the ground truth. Atom-mapping is excluded from this comparison, but stereochemistry, which describes the relative orientation of atoms in the molecule, is retained. The evaluation is carried out for two settings, with the reaction class being known or unknown.
Baselines For evaluating overall performance, we compare GRAPHRETRO to nine baselines — four template-based, three template-free, and two semi-template-based methods. These include:
Template-Based: RETROSIM Coley et al. [2017b] ranks templates for a given target molecule by computing molecular similarities to precedent reactions. NEURALSYM [Segler and Waller, 2017] trains a model to rank templates given a target molecule. GLN [Dai et al., 2019] models the joint distribution of templates and reactants in a hierarchical fashion using logic variables. DUALTB [Sun et al., 2021] uses an energy-based model formulation for retrosynthesis, with additional parameterizations and loss terms to enforce the duality between forward (reaction prediction) and backward (retrosynthesis prediction). Inference is carried out using reactant candidates obtained by applying an extracted template set to the products.
Template-Free: SCROP [Zheng et al., 2019], LV-TRANSFORMER [Chen et al., 2019] and DUALTF [Sun et al., 2021] use the Transformer architecture [Vaswani et al., 2017] to output reactant SMILES given a product SMILES. To improve the validity of their suggestions, SCROP include a second Transformer that functions as a syntax correcter. LV-TRANSFORMER uses a latent variable mixture model to improve diversity of suggestions. DUALTF utilizes additional parameterizations and loss terms to enforce the duality between forward (reaction prediction) and backward (retrosynthesis prediction).
Semi-Template-Based: G2GS [Shi et al., 2020] and RETROXPERT [Yan et al., 2020] first identify synthons, and then expand the synthons into reactants by either sequential generation of atoms and bonds (G2Gs), or using the Transformer architecture (RETROXPERT). The training dataset for the Transformer in [Yan et al., 2020] is augmented with incorrectly predicted synthons with the goal of learning a correction mechanism.
Results for NEURALSYM are taken from [Dai et al., 2019]. The authors in [Yan et al., 2020] report their performance being affected by the dataset leakage2. Thus, we use the most recent results from their website on the canonicalized dataset. For remaining baselines, we directly use the values reported in their paper. For the synthon completion module, we use the ind configuration given its better empirical performance.
2https://github.com/uta-smile/RetroXpert

heading: 4.1 Overall Performance

text:
Reaction class unknown As shown in Table 1, when the reaction class is unknown, GRAPHRETRO outperforms G2GS by 4.8% and and RETROXPERT by 3.3% in top-1 accuracy. Performance improvements are also seen for larger n, except for n = 5. Barring DUALTB, the top-1 accuracy is also better than other template-free and template-based methods. For larger n, one reason for lower top-n accuracies than most template-based methods is that templates already contain combinations of leaving group patterns. In contrast, our model learns to discover these during training. A second hypothesis to this end is that simply adding log-likelihood scores from edit prediction and synthon completion models may be suboptimal and bias the beam search in the direction of the more dominating term. We leave it to future work to investigate scoring functions that rank the attachment.
Reaction class known When the reaction class is known, GRAPHRETRO outperforms G2GS and RETROXPERT by a margin of 3% and 2% respectively in top-1 accuracy. GRAPHRETRO also outperforms all the template-free methods in top-n accuracy. for GRAPHRETRO are also better than most template-based and template-free methods. When the reaction class is known, RETROSIM and GLN restrict template sets corresponding to the reaction class, thus improving performance. The increased edit prediction performance (Section 4.2) for GRAPHRETRO helps outweigh this factor, achieving comparable or better performance till n = 5.

heading: 4.2 Individual Module Performance

text:
To gain more insight into the working of GRAPHRETRO, we evaluate the top-n accuracy (n = 1, 2, 3, 5) of edit prediction and synthon completion modules, along with corresponding ablation studies, with results shown in Table 2.
Edit Prediction For the edit prediction module, we compare the true edit(s) to top-n edits predicted by the model. We also consider two ablation studies, one where we directly use the initial edit scores without updating them, and the other where we predict edits using atom-pairs instead of existing bonds and atoms. Both design choices lead to improvements in performance, as shown in Table 2. We hypothesize that the larger improvement compared to edit prediction using atom-pairs is due to the easier optimization procedure, with lesser imbalance between labels 1 and 0.
Synthon Completion For evaluating the synthon completion module, we first apply the true edits to obtain synthons, and compare the true leaving groups to top-n leaving groups predicted by the model. We test the performance of both the ind and shared configurations. Both configurations perform similarly, and are able to identify ~ 97% (close to its upper bound of 99.7%) of the true leaving groups in its top-5 choices, when the reaction class is known. The top-1, 3 and 5 accuracies of the synthon completion for unknown reaction classes for G2Gs are 61.1%, 81.5% and 86.7% respectively, while ours are 75.6%, 92.5% and 96.1%, indicating a 10-14% performance improvement using a classification formulation over the generative one adopted by G2Gs.

heading: 4.3 Example Predictions

text: In Figure 2, we visualize the model predictions and the ground truth for three cases. Figure 2a shows an example where the model identifies both the edits and leaving groups correctly. In Figure 2b, the correct edit is identified but the predicted leaving groups are incorrect. We hypothesize this is due to the fact that in the training set, leaving groups attaching to the carbonyl carbon (C=O) are small (e.g. -OH, -NH2, halides). The true leaving group in this example, however, is large. The model is unable to reason about this and predicts the small leaving group -I. In Figure 2c, the model identifies the edit and consequently the leaving group incorrectly. This highlights a limitation of our model. If the edit is predicted incorrectly, the model cannot suggest the true precursors.

heading: 4.4 Limitations

text: The simplified and interpretable construction of GRAPHRETRO comes with certain limitations. First, the overall performance of the model is limited by the performance of the edit prediction step. If the predicted edit is incorrect, the true reactants cannot be salvaged. This limitation is partly remedied by our model design, that allows for user intervention to correct the edit. Second, our method is reliant on atom-mapping for extracting edits and leaving groups. Extracting edits directly based on substructure matching currently suffer from false positives, and heuristics to correct for these result in correct edits in only ~90% of the cases. Third, our formulation assumes that we have as many synthons as reactants, which is violated in some reactions. We leave it to future work to extend the model to realize a single reactant from multiple synthons, and introduce more chemically meaningful edit correction mechanisms.

heading: 5 Conclusion

text:
Previous methods for single-step retrosynthesis either restrict prediction to a template set, are insensitive to molecular graph structure or generate molecules from scratch. We address these shortcomings by introducing a graph-based semi-template-based model inspired by a chemist’s
workflow, enhancing the interpretability of retrosynthesis models. Given a target molecule, we first identify synthetic building blocks (synthons) which are then realized into valid reactants, thus avoiding molecule generation from scratch. Our model outperforms previous semi-template-methods by significant margins on the benchmark dataset. Future work aims to extend the model to realize a single reactant from multiple synthons, and introduce more chemically meaningful components to improve the synergy between such tools for retrosynthesis prediction and a practitioner’s expertise.

heading: Acknowledgements

text: This research was supported by the Machine Learning for Pharmaceutical Discovery and Synthesis Consortium at MIT. V.R.S. was also supported by the Zeno Karl Schindler Foundation. C.B. was supported by the Swiss National Science Foundation under the National Center of Competence in Research (NCCR) Catalysis under grant agreement 51NF40 180544. We thank the Leonhard scientific computing cluster at ETH Zürich for providing computational resources.

title: Learning Graph Models for Retrosynthesis Prediction


summaries: The authors consider the problem of refining low-res (LR) spatial observations using high-res (HR) covariates.   Basically, 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.  In 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).  The 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.

summaries: 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.  From 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).  From 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.

summaries: The paper considers the problem of learning a function  $$f:\mathcal{X}\to \mathbb{R}$$ from the relationship  $$g(y) = \int_x f(x) \mathbb{P}[x|y]dx$$ where 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)$.  Learning 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.  The 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.

summaries: 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.

summaries: 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.

abstractText: 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 [1]. 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.

authors:
Kelvin Hsu
Fabio Ramos

title: Bayesian Deconditional Kernel Mean Embeddings

venue: Proceedings of Machine Learning Research. PMLR,

authors:
L.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
B.N. Holben

title: The MODIS Aerosol Algorithm, Products, and Validation

venue: Journal of the Atmospheric Sciences,

authors:
S. Platnick
M.D. King
S.A. Ackerman
W.P. Menzel
B.A. Baum
J.C. Riedi
R.A. Frey

title: The MODIS cloud products: algorithms and examples from Terra

venue: IEEE Transactions on Geoscience and Remote Sensing,

authors:
Graeme L Stephens
Deborah G Vane
Ronald J Boain
Gerald G Mace
Kenneth Sassen
Zhien Wang
Anthony J Illingworth
Ewan J O’connor
William B Rossow
Stephen L Durden

title: The CloudSat mission and the A-train: A new dimension of space-based observations of clouds and precipitation

venue: Bulletin of the American Meteorological Society,

authors:
William L. Barnes
Thomas S. Pagano
Vincent V. Salomonson

title: Prelaunch characteristics of the Moderate Resolution Imaging Spectroradiometer (MODIS) on EOS-AM1

authors:
Zbyněk Malenovský
Helmut Rott
Josef Cihlar
Michael E. Schaepman
Glenda García-Santos
Richard Fernandes
Michael Berger

title: Sentinels for science: Potential of Sentinel-1, -2, and -3 missions for scientific observations of ocean, cryosphere, and land

venue: Remote Sensing of Environment,

authors:
Edward J. Knight
Geir Kvaran

title: Landsat-8 operational land imager design

authors:
Yivan Zhang
Nontawat Charoenphakdee
Zhenguo Wu
Masashi Sugiyama

title: Learning from aggregate observations

venue: In Advances in Neural Information Processing Systems,

authors:
Leon Ho Chung Law
Dino Sejdinovic
Ewan Cameron
Tim C.D. Lucas
Seth Flaxman
Katherine Battle
Kenji Fukumizu

title: Variational learning on aggregate outputs with Gaussian processes

authors:
Oliver Hamelijnck
Theodoros Damoulas
Kangrui Wang
Mark A. Girolami

title: Multiresolution multi-task Gaussian processes

authors:
Fariba Yousefi
Michael Thomas Smith
Mauricio A. Álvarez

title: Multi-task learning for aggregated data using Gaussian processes

authors:
Yusuke Tanaka
Toshiyuki Tanaka
Tomoharu Iwata
Takeshi Kurashima
Maya Okawa
Yasunori Akagi
Hiroyuki Toda

title: Spatially aggregated Gaussian processes with multivariate areal outputs

authors:
Arto Klam Ville Tanskanen
Krista

title: Longi. Non-Linearities in Gaussian Processes with Integral Observations

venue: IEEE international Workshop on Machine Learning for Signal,

title: Variational learning of inducing variables in sparse gaussian processes

venue: In AISTATS,

authors:
Anna Vaughan
Will Tebbutt
J Scott Hosking
Richard E Turner

title: Convolutional conditional neural processes for local climate downscaling

venue: Geoscientific Model Development Discussions,

authors:
Brian Groenke
Luke Madaus
Claire Monteleoni

title: Climalign: Unsupervised statistical downscaling of climate variables via normalizing flows

venue: In Proceedings of the 10th International Conference on Climate Informatics,

authors:
Thomas Vandal
Evan Kodra
Sangram Ganguly
Andrew Michaelis
Ramakrishna Nemani
Auroop R Ganguly

title: Deepsd: Generating high resolution climate change projections through single image super-resolution

venue: In Proceedings of the 23rd acm sigkdd international conference on knowledge discovery and data mining,

authors:
Veronika Eyring
Sandrine Bony
Gerald A Meehl
Catherine A Senior
Bjorn Stevens
Ronald J Stouffer
Karl E Taylor

title: Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization

venue: Geoscientific Model Development,

title: Earth system models: an overview

venue: WIREs Climate Change,

authors:
Marko Scholze
J. Icarus Allen
William J. Collins
Sarah E. Cornell
Chris Huntingford
Manoj M. Joshi
Jason A. Lowe
Robin S. Smith
Oliver Wild

title: Earth system models

authors:
Krikamol Muandet
Kenji Fukumizu
Bharath Sriperumbudur
Bernhard Schölkopf

title: Kernel mean embedding of distributions: A review and beyond

venue: arXiv preprint arXiv:1605.09522,

authors:
Steffen Grünewälder
Guy Lever
Luca Baldassarre
Sam Patterson
Arthur Gretton
Massimilano Pontil

title: Conditional Mean Embeddings as Regressors

venue: In Proceedings of the 29th International Coference on International Conference on Machine Learning,

authors:
Zoltán Szabó
Bharath K. Sriperumbudur
Barnabás Póczos
Arthur Gretton

title: Learning theory for distribution regression

venue: Journal of Machine Learning Research,

authors:
Rahul Singh
Maneesh Sahani
Arthur Gretton

title: Kernel instrumental variable regression

authors:
Le Song
Kenji Fukumizu
Arthur Gretton

title: Kernel embeddings of conditional distributions: A unified kernel framework for nonparametric inference in graphical models

venue: IEEE Signal Processing Magazine,

authors:
Kenji Fukumizu
Francis R Bach
Michael I Jordan

title: Dimensionality reduction for supervised learning with reproducing kernel hilbert spaces

authors:
Kenji Fukumizu
Arthur Gretton
Xiaohai Sun
Bernhard Schölkopf

title: Kernel measures of conditional dependence

authors:
Bharath K. Sriperumbudur
Kenji Fukumizu
Arthur Gretton
Bernhard Schölkopf
Gert R.G. Lanckriet

title: On the empirical estimation of integral probability metrics

venue: Electronic Journal of Statistics,

authors:
Le Song
Jonathan Huang
Alex Smola
Kenji Fukumizu

title: Hilbert space embeddings of conditional distributions with applications to dynamical systems

venue: In Proceedings of the 26th Annual International Conference on Machine Learning,

title: Gaussian measure in Hilbert space and applications in numerical analysis

venue: The Rocky Mountain Journal of Mathematics,

authors:
François-Xavier Briol
Chris J Oates
Mark Girolami
Michael A Osborne
Dino Sejdinovic

title: Probabilistic integration: A role in statistical computation

authors:
Tim GJ Rudner
Dino Sejdinovic
Yarin Gal

title: Inter-domain deep Gaussian processes

venue: In International Conference on Machine Learning,

authors:
Krikamol Muandet
Bharath Sriperumbudur
Kenji Fukumizu
Arthur Gretton
Bernhard Schölkopf

title: Kernel mean shrinkage estimators

authors:
Vern Paulsen
Mrinal Raghupathi

title: An introduction to the theory of reproducing kernel Hilbert spaces

authors:
Andrea Caponnetto
Ernesto De Vito

title: Optimal rates for the regularized least-squares algorithm

venue: Foundations of Computational Mathematics,

authors:
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 DeVito
Martin Raison
Alykhan Tejani
Sasank Chilamkurthy
12 Benoit Steiner
Lu Fang
Junjie Bai
Soumith Chintala

title: PyTorch: An Imperative Style, High-Performance Deep Learning Library

venue: Advances in Neural Information Processing Systems

authors:
Jacob Gardner
Geoff Pleiss
Kilian Q Weinberger
David Bindel
Andrew G Wilson

title: Gpytorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration

authors:
F. Pedregosa
G. Varoquaux
A. Gramfort
V. Michel
B. Thirion
O. Grisel
M. Blondel
P. Prettenhofer
R. Weiss
V. Dubourg
J. Vanderplas
A. Passos
D. Cournapeau
M. Brucher
M. Perrot
E. Duchesnay

title: Scikit-learn: Machine learning in Python

authors:
Diederik P. Kingma
Jimmy Ba

title: Adam: A Method for Stochastic Optimization

venue: In 3rd International Conference on Learning Representations,

authors:
D. Watson-Parris
N. Schutgens
N. Cook
Z. Kipling
P. Kershaw
E. Gryspeerdt
B. Lawrence
P. Stier

title: Community Intercomparison Suite (CIS) v1.4.0: a tool for intercomparing models and observations

title: MOHC HadGEM3-GC31-HM model output prepared for CMIP6 High- ResMIP hist-1950

venue: Earth System Grid Federation,

title: CNRM-CERFACS CNRM-CM6-1-HR model output prepared for CMIP6 HighResMIP hist-1950

authors:
Ali Rahimi
Benjamin Recht

title: Random features for large-scale kernel machines

authors:
Charles A Micchelli
Massimiliano Pontil

title: On learning vector-valued functions

venue: Neural computation,

authors:
Ho Chung Law
Peilin Zhao
Leung Sing Chan
Junzhou Huang
Dino Sejdinovic

title: Hyperparameter learning via distributional transfer

authors:
Ingo Steinwart
Andreas Christmann

title: Support Vector Machines

heading: 1 Introduction

text:
Spatial observations often operate at limited resolution due to practical constraints. For example, remote sensing atmosphere products [2, 3, 4, 5] provide measurement of atmospheric properties such as cloud top temperatures and optical thickness, but only at a low resolution. Devising methods to refine low-resolution (LR) variables for local-scale analysis thus plays a crucial part in our understanding of the anthropogenic impact on climate
When high-resolution (HR) observations of different covariates are available, details can be instilled into the LR field for refinement. This task is referred to as statistical downscaling or spatial disaggregation and models LR observations as the aggregation of an unobserved underlying HR field. For example, multispectral optical satellite imagery [6, 7] typically comes at higher resolution than atmospheric products and can be used to refine the latter.
Statistical downscaling has been studied in various forms, notably giving it a probabilistic treatment [8, 9, 10, 11, 12, 13], in which Gaussian processes (GP) [14] are typically used in conjunction with a sparse variational formulation [15] 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 [16, 17, 18] study this problem from a different setting, where they often assume that both HR and LR matched observations are available for training. Then, their approaches follow a standard supervised learning setting in learning a mapping between different resolutions. ∗Indicates equal contribution †Department of Statistics, Oxford, UK, OX1 3LB. <siu.chau@stats.ox.ac.uk, shahine.bouabid@stats.ox.ac.uk, dino.sejdinovic@stats.ox.ac.uk>
35th Conference on Neural Information Processing Systems (NeurIPS 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 a 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 [19, 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 bx = {x(1), . . . , x(n)} ⊂ X be a general notation for bags of HR covariates, f : X → R the field of interest we wish to recover and z̃ the LR aggregate observations from the field f . We suppose that bx and z̃ are unmatched, but that there exists mediating covariates y, ỹ ∈ Y , such that (bx, y) are jointly observed and likewise for (ỹ, z̃) as illustrated in Figure 1. We assume the following aggregation observation model z̃ = EX [f(X)|Y = ỹ] + ε with some noise ε. Our goal in mediated statistical downscaling is then to estimate f given (bx, y) and (ỹ, z̃), which corresponds to the deconditioning problem introduced in [1].
Motivated by applications in likelihood-free inference and tasktransfer regression, Hsu and Ramos [1] first studied the deconditioning 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 deconditioning that results into a much simpler and elegant way to arrive to the DME-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 a 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 a posterior distribution of the downscaled field as a principled Bayesian solution to the downscaling task on indirectly matched data. For scalability, we provide a tractable variational inference approximation and an alternative approximation to the conditional mean operator (CMO) [22] to speed up computation for large multi-resolution datasets.
From a theoretical stand point, we further develop the framework of DMEs by establishing it as a two-staged vector-valued regressor with a natural reconstruction loss that mirrors Grünewälder et al. [23]’s work on conditional mean embeddings. This perspective allows us to leverage distribution regression theory from [24, 25] 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 statistical-computational efficiency.
Our contributions are summarized as follows:
• We propose a Bayesian formulation of the deconditioning problem of Hsu and Ramos [1]. We establish its posterior mean as a DME-based estimate and its posterior covariance as a 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 a 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.

heading: 2 Background Materials

heading: 2.1 Notations

text:
Let X,Y be a pair of random variables taking values in non-empty sets X and Y , respectively. Let k : X × X → R and ` : Y × Y → R be positive definite kernels. The closure of the span of their canonical feature maps kx := k(x, ·) and `y := `(y, ·) for x ∈ X and y ∈ Y respectively induces RKHSHk ⊆ RX andH` ⊆ RY endowed with inner products 〈·, ·〉k and 〈·, ·〉`.
We observe realizations bD1 = {bxj , yj}Nj=1 of bags bxj = {x (i) j } nj i=1 from conditional distribution PX|Y=yj , with bag-level covariates yj sampled from PY . We concatenate them into vectors x :=[ bx1 . . . bxN ]> and y := [y1 . . . yN ] >.
For simplicity, our exposition will use the notation without bagging – i.e. where D1 = {xj , yj}Nj=1 – when the generality of our contributions will be relevant to the theory of deconditioning from Hsu and Ramos [1], in Sections 2, 3 and 4. We will come back to a bagged dataset formalism in Sections 3.3 and 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 Φx := [kx1 . . . kxN ] and Ψy := [`y1 . . . `yN ] and we denote Gram matrices as Kxx = Φ>x Φx = [k(xi, xj)]1≤i,j≤N and Lyy = Ψ > y Ψy = [`(yi, yj)]1≤i,j≤N . The notation abuse (·)>(·) is a shorthand for the elementwise RKHS inner products when it is clear from the context.
Let Z denote the real-valued random variable stemming from the noisy conditional expectation of some unknown latent function f : X → R, as Z = EX [f(X)|Y ] + ε. We suppose one observes another set of realizations D2 = {ỹj , z̃j}Mj=1 from PY Z , which is sampled independently from D1. Likewise, we stack observations into vectors ỹ := [ỹ1 . . . ỹM ] >, z̃ := [z̃1 . . . z̃M ] > and define feature map Ψỹ := [`ỹ1 . . . `ỹM ].

heading: 2.2 Conditional and Deconditional Kernel Mean Embeddings

text:
Marginal and Joint Mean Embeddings Kernel mean embeddings of distributions provide a powerful framework for representing and manipulating distributions without specifying their parametric form [22, 26]. The marginal mean embedding of measure PX is defined as µX := EX [kX ] ∈ Hk and corresponds to the Riesz representer of expectation functional f 7→ EX [f(X)]. It can hence be used to evaluate expectations EX [f(X)] = 〈f, µX〉k. If the mapping PX 7→ µX is injective, the kernel k is said to be characteristic [27], a property satisfied for the Gaussian and Matérn kernels on Rd [28]. In practice, Monte Carlo estimator µ̂X := 1N ∑N i=1 kxi provides an unbiased estimate of µX [29].
Extending this rationale to embeddings of joint distributions, we define CY Y := EY [`Y ⊗ `Y ] ∈ H` ⊗ H` and CXY := EX,Y [kX ⊗ `Y ] ∈ Hk ⊗ H`, which can be identified with the crosscovariance operators between Hilbert spaces CY Y : H` → H` and CXY : H` → Hk. They correspond to the Riesz representers of bilinear forms (g, g′) 7→ Cov(g(Y ), g′(Y )) = 〈g, CY Y g′〉` and (f, g) 7→ Cov(f(X), g(Y )) = 〈f, CXY g〉k. As above, empirical estimates are obtained as ĈY Y = 1 NΨyΨ > y = 1 N ∑N i=1 `yi ⊗ `yi and ĈXY = 1 NΦxΨ > y = 1 N ∑N i=1 kxi ⊗ `yi . Again, notation abuse (·)(·)> is a shorthand for element-wise tensor products when clear from context.
Conditional Mean Embeddings Similarly, one can introduce RKHS embeddings for conditional distributions referred to as Conditional Mean Embeddings (CME). The CME of conditional probability measure PX|Y=y is defined as µX|Y=y := EX [kX |Y = y] ∈ Hk. As introduced by Fukumizu et al. [27], it is common to formulate conditioning in terms of a Hilbert space operator CX|Y : H` → Hk called the Conditional Mean Operator (CMO). CX|Y satisfies by definition CX|Y `y = µX|Y=y and C>X|Y f = EX [f(X)|Y = ·] ,∀f ∈ Hk, where C > X|Y de-
notes the adjoint of CX|Y . Plus, the CMO admits expression CX|Y = CXY C −1 Y Y , provided `y ∈ Range(CY Y ), ∀y ∈ Y [26, 27]. Song et al. [30] show that a nonparametric empirical form of the CMO writes
ĈX|Y = Φx(Lyy +NλIN ) −1Ψ>y , (1)
where λ > 0 is some regularisation ensuring the empirical operator is globally defined and bounded.
As observed by Grünewälder et al. [23], since CX|Y defines a mapping from H` to Hk, 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 [1] as a new class of embeddings, Deconditional Mean Embeddings (DME) are natural counterparts of CMEs. While CME µX|Y=y ∈ Hk allows to take the conditional expectation of any f ∈ Hk through inner product EX [f(X)|Y = y] = 〈f, µX|Y=y〉k, the DME denoted µX=x|Y ∈ H` solves the inverse problem3 and allows to recover the initial function of which the conditional expectation was taken, through inner product 〈EX [f(X)|Y = ·], µX=x|Y 〉` = f(x). The associated Hilbert space operator, the Deconditional Mean Operator (DMO), is thus defined as the operator DX|Y : Hk → H` such that D>X|Y EX [f(X)|Y = ·] = f , ∀f ∈ Hk. It admits an expression in terms of CMO and cross-covariance operators DX|Y = (CX|Y CY Y )>(CX|Y CY Y (CX|Y )>)−1 provided `y ∈ Range(CY Y ) and kx ∈ Range(CX|Y CY Y C>X|Y ) ,∀y ∈ Y and ∀x ∈ X . A nonparametric empirical estimate of the DMO using datasets D1 and D2 as described above, is given by D̂X|Y = Ψỹ ( A>KxxA +M IM )−1 A>Φ>x where > 0 is a regularisation term and A := (Lyy + NλI)−1Lyỹ can be interpreted as a mediation operator. Applying the DMO to expected responses z̃, Hsu and Ramos [1] are able to recover an estimate of f as
f̂(x) = k(x,x)A ( A>KxxA +M IM )−1 z̃. (2)
Note that since separate samples ỹ can be used to estimate CY Y , this naturally fits a mediating variables setup where x and the conditional means z̃ are not jointly observed.

heading: 3 Deconditioning with Gaussian processes

text:
In this section, we introduce Conditional Mean Process (CMP), a stochastic process stemming from the conditional expectation of a GP. We provide a 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 5.
For simplicity, we put aside observations bagging in Sections 3.1 and 3.2, our contributions being relevant to 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 a Borel space and feature maps kx and `y are Borel-measurable functions for any x ∈ X , y ∈ Y . All proofs and derivations of the paper are included in the appendix.

heading: 3.1 Conditional Mean Process

text:
Bayesian quadrature [14, 31, 32] is based on the observation that the integral of a GP with respect to some marginal measure is a Gaussian random variable. We start by probing the implications of integrating with respect to conditional distribution PX|Y=y 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.e.∫ X |f |dPX < ∞ a.s. The CMP induced by f with respect to PX|Y is defined as the stochastic
3the slightly unusual notation µX=x|Y is taken from Hsu and Ramos [1] and is meant to contrast the usual conditioning µX|Y =y
process {g(y) : y ∈ Y} given by g(y) = EX [f(X)|Y = y] = ∫ X f(x) dPX|Y=y(x). (3)
By linearity of the integral, it is clear that g(y) is a Gaussian random variable for each y ∈ 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 [|m(X)|] <∞ and EX [‖kX‖k] <∞ and let (X ′, Y ′) ∼ PXY . Then g is a Gaussian process g ∼ GP(ν, q) a.s. , specified by
ν(y) = EX [m(X)|Y = y] q(y, y′) = EX,X′ [k(X,X ′)|Y = y, Y ′ = y′] (4) ∀y, y′ ∈ Y . Furthermore, q(y, y′) = 〈µX|Y=y, µX|Y=y′〉k a.s.
Intuitively, the CMP can be understood as a GP on the conditional means where its covariance q(y, y′) is induced by the similarity between the CMEs at y and y′. Resorting to the kernel ` defined on Y , we can reexpress the covariance using Hilbert space operators as q(y, y′) = 〈CX|Y `y, CX|Y `y′〉k. A natural nonparametric estimate of the CMP covariance thus comes using the CMO estimator from (1) as q̂(y, y′) = `(y,y) (Lyy +NλIN ) −1 Kxx (Lyy +NλIN ) −1 `(y, y′). When m ∈ Hk, the mean function can be written as ν(y) = 〈µX|Y=y,m〉k for which we can also use empirical estimate ν̂(y) = `(y,y) (Lyy +NλIN ) −1 Φ>xm. Finally, one can also quantify the covariance between the CMP g and its integrand f , i.e. Cov(f(x), g(y)) = EX [k(x,X)|Y = y]. Under the same assumptions as Proposition 3.2, this covariance can be expressed using mean embeddings, i.e. Cov(f(x), g(y)) = 〈kx, µX|Y=y〉k and admits empirical estimate k(x,x) (Lyy +NλIN ) −1 `(y, y).

heading: 3.2 Deconditional Posterior

text:
Given independent observations introduced above, D1 = {x,y} and D2 = {ỹ, z̃}, we may now consider an additive noise model with CMP prior on aggregate observations z̃|ỹ ∼ N (g(ỹ), σ2IM ). Let Qỹỹ := q(ỹ, ỹ) be the kernel matrix induced by q on ỹ and let Υ := Cov(f(x), z̃) = Φ>xCX|Y Ψỹ be the cross-covariance between f(x) and z̃. The joint normality of f(x) and z̃ gives[
f(x) z̃
] | y, ỹ ∼ N ([ m(x) ν(ỹ) ] , [ Kxx Υ
Υ> Qỹỹ + σ 2IM
]) . (5)
Using Gaussian conditioning, we can then readily derive the posterior distribution of the underlying GP field f given the aggregate observations z̃ corresponding to ỹ. The latter can naturally be degenerated if observations are paired, i.e. y = ỹ. This formulation can be seen as an example of the inter-domain GP [33], where we utilise the observed conditional means z̃ as inter-domain inducing features for inference of f . Proposition 3.3 (Deconditional Posterior). Given aggregate observations z̃ with homoscedastic noise σ2, the deconditional posterior of f is defined as the Gaussian process f |z̃ ∼ GP(md, kd) where
md(x) = m(x) + k > x CX|Y Ψỹ(Qỹỹ + σ 2IM ) −1(z̃− ν(ỹ)), (6) kd(x, x ′) = k(x, x′)− k>x CX|Y Ψỹ(Qỹỹ + σ2IM )−1Ψ>ỹC>X|Y kx′ . (7)
Substituting terms by their empirical forms, we can define a nonparametric estimate of the md as
m̂d(x) := m(x) + k(x,x)A(Q̂ỹỹ + σ 2IM ) −1(z̃− ν̂(ỹ))) (8)
which, when m = 0, reduces to the DME-based estimator in (2) by taking the noise variance σ 2
N as the inverse regularization parameter. Hsu and Ramos [1] recover a similar posterior mean expression in their Bayesian interpretation of DME. However, they do not link the distributions of f and its CMP, which leads to much more complicated chained inference derivations combining fully Bayesian and MAP estimates, while we naturally recover it using simple Gaussian conditioning.
Likewise, an empirical estimate of the deconditional covariance is given by
k̂d(x, x ′) := k(x, x′)− k(x,x)A(Q̂ỹỹ + σ2IM )−1A>k(x, x′). (9)
Interestingly, the latter can be rearranged to write as the difference between the original kernel and the kernel undergoing conditioning and deconditioning steps k̂d(x, x′) = k(x, x′) − 〈kx, D̂X|Y ĈX|Y kx′〉k. This can be interpreted as a measure of reconstruction quality, which degenerates in the case of perfect deconditioning, i.e. D̂X|Y ĈX|Y = IdHk .

heading: 3.3 Deconditioning and multiresolution data

text:
Downscaling application would typically correspond to multiresolution data, with bag dataset bD1 = {(bxj , yj)}Nj=1 where bxj = {x (i) j } nj i=1. In this setup, the mismatch in size between vector concatenations x = [x(1)1 . . . x (nN ) N ] > and y = [y1 . . . yN ] > prevents from readily applying (1) to estimate the CMO and thus infer the deconditional posterior. There is, however, a straightforward approach to alleviate this: simply replicate bag-level covariates yj to match bags sizes nj . Although simple, this method incurs a O(( ∑N j=1 nj) 3) cost due to matrix inversion in (1).
Alternatively, since bags bxj are sampled from conditional distribution PX|Y=yj , unbiased Monte Carlo estimators of CMEs are given by µ̂X|Y=yj = 1 nj ∑nj i=1 kx(i)j . Let M̂y = [µ̂X|Y=y1 . . . µ̂X|Y=yN ] > denote their concatenation along columns. We can then rewrite the crosscovariance operator as CXY = EY [EX [kX |Y ] ⊗ `Y ] and hence take 1N M̂yΨ > y as an estimate for CXY . Substituting empirical forms into CX|Y = CX|Y C −1 Y Y and applying Woodbury identity, we obtain an alternative CMO estimator that only requires inversion of a N × N matrix. We call it Conditional Mean Shrinkage Operator and define it as
SĈX|Y := M̂y(Lyy + λNIN ) −1Ψ>y . (10)
This estimator can be seen as a generalisation of the Kernel Mean Shrinkage Estimator [34] to the conditional case. We provide in Appendix C modifications of (8) and (9) including this estimator for the computation of the deconditional posterior.

heading: 4 Deconditioning as regression

text:
In Section 3.2, we obtain a DMO-based estimate for the posterior mean of f |z̃. 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 a two-staged vector-valued regression. This frequentist turn enables us to leverage rich theory of vector-valued regression and establish under mild assumptions a convergence rate on the DMO estimator, with conditions to fulfill minimax optimality in terms of statistical-computational 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 to [35] for a comprehensive overview of vector-valued RKHS theory. As for Sections 3.1 and 3.2, this section contributes to the general theory of DMEs and we hence put aside the bag notations.
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>X|Y is equivalent to solving a vector-valued kernel ridge regression problem in the hypothesis space of Hilbert-Schmidt operators fromHk toH`, denoted as HS(Hk,H`). Specifically, we may consider the operator-valued kernel defined overHk as Γ(f, f ′) := 〈f, f〉k IdH` . We denoteHΓ theH`-valued RKHS spanned by Γ with norm ‖ · ‖Γ, which can be identified to HS(Hk,H`)4. Singh et al. [25] frame CMO regression as the minimisation surrogate discrepancy Ec(C) := EX,Y [ ‖kX − C>`Y ‖2k ] , to which they substitute an empirical regularised version restricted toHΓ given by
Êc(C) := 1
N N∑ i=1 ‖kxi − C>`yi‖2k + λ‖C‖2Γ C ∈ HΓ λ > 0 (11)
ThisHk-valued kernel ridge regression problem admits a closed-form minimiser which shares the same empirical form as the CMO, i.e. Ĉ>X|Y = arg minC∈HΓ Êc(C) [23, 25].
4HΓ = Span {Γfh, f ∈ Hk, h ∈ H`} = Span {f ⊗ h, f ∈ Hk, h ∈ H`} = Hk ⊗H` ∼= HS(Hk,H`)
Stage 2 : Regressing the Deconditional Mean Operator The DMO on the other hand is defined as the operator DX|Y : Hk → H` such that ∀f ∈ Hk, D>X|Y C > X|Y f = f . Since deconditioning corresponds to finding a pseudo-inverse to the CMO, it is natural to consider a reconstruction objective Ed(D) := EY [ ‖`Y −DCX|Y `Y ‖2` ] . Introducing a novel characterization of the DMO, we propose to minimise this objective in the hypothesis space of Hilbert-Schmidt operators HS(Hk,H`) which identifies to HΓ. As per above, we denote ĈX|Y the empirical CMO learnt in Stage 1, and we substitute the loss with an empirical regularised formulation onHΓ
Êd(D) := 1
M M∑ j=1 ‖`ỹj −DĈX|Y `ỹj‖2` + ‖D‖2Γ D ∈ HΓ > 0 (12)
Proposition 4.1 (Empirical DMO as vector-valued regressor). The minimiser of the empirical reconstruction risk is the empirical DMO, i.e. D̂X|Y = arg minD∈HΓ Êd(D)
Since it requires to estimate the CMO first, minimising (12) can be viewed as a two-staged vector value regression problem.
Convergence results Following the footsteps of [24, 25], this perspective enables us to state the performance of the DMO estimate in terms of asymptotic convergence of the objective Ed. 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 Hk is infinite dimensional. The family of distribution considered is a general class of priors that does not assume parametric distributions and is parametrized by two variables: b > 1 controls the effective input dimension and c ∈]1, 2] controls functional smoothness. Mild regularity assumptions are also placed on the original spaces X ,Y , their corresponding RKHS Hk,H` and the vector-valued RKHS HΓ. We discuss these assumptions in details in Appendix D. Importantly, whileHk can be infinite dimensional, we nonetheless have to assume the RKHSH` is finite dimensional. In further research, we will relax this assumption. Theorem 4.2 (Empirical DMO Convergence Rate). Denote DPY = arg minD∈HΓ Ed(D). Assume assumptions stated in Appendix D are satisfied. In particular, let (b, c) and (0, c′) be the parameters of the restricted class of distribution for PY and PXY respectively and let ι ∈]0, 1] be the Hölder
continuity exponent inHΓ. Then, if we choose λ = N− 1 c′+1 , N = M a(c′+1) ι(c′−1) where a > 0, we have the following result,
• If a ≤ b(c+1)bc+1 , then Ed(D̂X|Y )− Ed(DPY ) = O(M −ac c+1 ) with = M −a c+1
• If a ≥ b(c+1)bc+1 , then Ed(D̂X|Y )− Ed(DPY ) = O(M −bc bc+1 ) with = M −b bc+1
This theorem underlines a trade-off between the computational and statistical efficiency with respect to the datasets cardinalities N = |D1|, M = |D2| and the problem difficulty b, c, c′. For a ≤ b(c+1)bc+1 , smaller ameans less samples fromD1 at fixedM and thus computational savings. But it also hampers convergence, resulting in reduced statistical efficiency. At a = b(c+1)bc+1 < 2, convergence rate is a minimax computational-statistical 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 ≥ b(c+1)bc+1 does not improve the convergence rate but only increases the size of D1 and hence the computational cost it bears.

heading: 5 Deconditional Downscaling Experiments

text:
We demonstrate and evaluate our CMP-based downscaling approaches on both synthetic experiments 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 a GP regression [14] baseline (GPR) modified to take bags centroids as the input. Experiments are implemented in PyTorch [37, 38], all code and datasets are made available5 and computational details are provided in Appendix E.
5https://github.com/shahineb/deconditional-downscaling
Variational Aggregate Learning VBAGG is introduced by Law et al. [9] as a variational aggregate learning framework to 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 zj and bag covariates {x(i)j }i by bag-wise averaging of a GP prior on the function of interest. In fact, the Gaussian VBAGG corresponds exactly to a special case of CMP on matched data, where the bag covariates are simply one hot encoded indices with kernel `(j, j′) = δ(j, j′) where δ is 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 to as BAGG-GP.

heading: 5.1 Swiss Roll

text:
The scikit-learn [39] swiss roll manifold sampling function allows to generate a 3D manifold of points x ∈ R3 mapped with their position along the manifold t ∈ R. Our objective will be to recover t for each point x by only observing t at an aggregate level. In the first experiment, we compare our model to existing weakly supervised spatial disaggregation methods when all high-resolution covariates are 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. [9], we regularly split space along height B − 1 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 {(bxj , tj)}Bj=1 where the jth bag contains nj points bxj = {x(i)j } nj i=1 and their corresponding targets tj = {t(i)j } nj i=1. We then construct bag aggregate targets by taking noisy bag targets average zj := 1nj ∑nj i=1 t (i) j + εj , where εj ∼ N (0, σ2). We thus obtain matched weakly supervised bag learning dataset D◦ = {(bxj , zj)}Bj=1. Since each bag corresponds to a height-level, the center altitude of each height split yj ∈ R is a natural candidate bag-level covariate that informs on relative positions of the bags. We can
augment the above dataset as D = {(bxj , yj , zj)}Bj=1. Using these bag datasets, we now wish to downscale aggregate targets zj to recover the unobserved manifold locations {tj}Bj=1 and be able to query the target at any previously unseen input x.
Models We use a zero-mean prior on f and choose a Gaussian kernel for k and `. 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 σ2 are learnt on D by optimising the marginal likelihood. For VARCMP, they are learnt jointly with variational distribution parameters by maximising an evidence lower bound objective. While CMP-based methods can leverage bag-level covariates yj , baselines are restricted to learn from D◦. Adam optimiser [40] is used in all experiments.
Results We test models against unobserved groundtruth {tj}Bj=1 by evaluating the root mean square error (RMSE) to the posterior mean. Table 1 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 ` is reflected even in a setting where each bag is matched with an aggregate target.

heading: 5.1.2 Indirect matching

text:
Experimental Setup We now impose indirect matching through mediating variable yj . We randomly select N = bB2 c bags which we consider to be the N first ones without loss of generality and split D into D1 = {(bxj , yj)}Nj=1 and D2 = {(ỹj , z̃j)} B−N j=1 = {(yN+j , zN+j)} B−N j=1 , such that no pair of covariates bag bxj and aggregate target z̃j are jointly observed.
Models CMP-based methods are naturally able to learn from this setting and are trained by independently 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 g ∼ GP(0, `) and fit regression model z̃j = g(ỹj) + εj over D2. We then use the predictive posterior mean to augment the first dataset as D′1 = {( bxj ,E[g(yj)|D2] )}N j=1
. This dataset can then be used to train BAGG-GP, VBAGG and GPR.
Results For comparison, we use the same evaluation as in the direct matching experiment. Table 1 underlines an 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 a substantial margin. This highlights how using a separate prior on g to mediate D1 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 a relative drop in performance for S-CMP.

heading: 5.2 Mediated downscaling of atmospheric temperature

text:
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, 4, 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 [42, 43] for the following variables: air temperature at cloud top (T), mean cloud top pressure (P), total cloud cover (TCC) and cloud albedo (α). First, we collocate TCC and α onto a HR latitudelongitude grid of size 360×720 to obtain fine grained fields (latitudeHR, longitudeHR, altitudeHR, TCCHR, αHR) augmented with a static HR surface altitude field. Then we collocate P and T onto a LR grid of size 21×42 to obtain coarse grained fields (latitudeLR, longitudeLR, PLR, TLR). 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 to an aggregate target zj := TLRj corresponding to bag j. Each bag includes HR covariates bxj = {x(i)j } nj i=1 := {(latitude HR(i) j , longitude HR(i) j , altitude HR(i) j , TCC HR(i) j , α HR(i) j )} nj i=1. To emulate unmatched observations, we randomly select N = bB2 c of the bags { bxj}Nj=1 and keep the opposite half of LR observations {zN+j}B−Nj=1 , 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 yj := (latitudeLRj , longitude LR j , P LR j ) which can separately be matched with both above sets to obtain datasets D1 = {(xj , yj)}Nj=1 and D2 = {(ỹj , z̃j)} B−N j=1 = {(yN+j , zN+j)} B−N j=1 .
Models Setup We only consider variational methods to scale to large number of pixels. VARCMP is naturally able to learn from indirectly matched data. We use a Matérn-1.5 kernel for rough spatial covariates (latitude, longitude) and a Gaussian kernel for atmospheric covariates (P, TCC, α) and surface altitude. k and ` are both taken as sums of Matérn and Gaussian kernels, and their hyperparameters are learnt along with noise variance during training. A high-resolution noise term is 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 [15]. Since baseline methods require a matched dataset, we proceed as with the unmatched swiss roll experiment and fit a GP regression model g with kernel ` on D2 and then use its predictive posterior mean to obtain pseudo-targets for the bags of HR covariates from D1.
Results Performance is evaluated by comparing downscaling deconditional posterior mean against original high resolution field THR available in CMIP6 data [43], which we emphasise is never observed. We use random Fourier features [44] approximation of kernel k to scale kernel evaluation to the HR covariates grid during testing. As reported in Table 2, VARCMP substantially outperforms both baselines with statistical significance provided in Appendix E. Figure 3 shows the reconstructed image with VARCMP, plots for other methods are included in the Appendix E. The model resolves statistical patterns from HR covariates into coarse resolution temperature pixels, henceforth reconstructing a faithful HR version of the temperature field.

heading: 6 Discussion

text:
We introduced a scalable Bayesian solution to 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, we demonstrated the advantages over spatial disaggregation baselines in synthetic data and in a challenging atmospheric fields downscaling problem.
In future work, exploring theoretical guarantees of the computationally efficient shrinkage formulation in a multi-resolution setting and relaxing finite dimensionality assumptions for the convergence rate will have fruitful practical and theoretical implications. Further directions also open up to quantify uncertainty over the deconditional posterior since it is computed using empirical estimates of the CMP covariance. This may be problematic if the mediating variable undergoes covariate shift between the two datasets.

heading: Acknowledgments

text: SLC is supported by the EPSRC and MRC through the OxWaSP CDT programme EP/L016710/1. SB is supported by the EU’s Horizon 2020 research and innovation programme under Marie SkłodowskaCurie grant agreement No 860100. DS is supported in partly by Tencent AI lab and partly by the Alan Turing Institute (EP/N510129/1).

title: Deconditional Downscaling with Gaussian Processes