Local Function Complexity for Active Learning via Mixture of Gaussian Processes
Authors: Danny Panknin
Stefan Chmiela
Klaus Robert Muller
Shinichi Nakajima
Venue: Accepted by TMLR ()

Abstract: Inhomogeneities in real-world data, e.g., due to changes in the observation noise level or variations in the structural complexity of the source function, pose a unique set of challenges for statistical inference. Accounting for them can greatly improve predictive power when physical resources or computation time is limited. In this paper, we draw on recent theoretical results on the estimation of local function complexity (LFC), derived from the domain of local polynomial smoothing (LPS), to establish a notion of local structural complexity, which is used to develop a model-agnostic active learning (AL) framework. Due to its reliance on pointwise estimates, the LPS model class is not robust and scalable concerning large input space dimensions that typically come along with real-world problems. Here, we derive and estimate the Gaussian process regression (GPR)-based analog of the LPS-based LFC and use it as a substitute in the above framework to make it robust and scalable. We assess the effectiveness of our LFC estimate in an AL application on a prototypical low-dimensional synthetic dataset, before taking on the challenging real-world task of reconstructing a quantum chemical force field for a small organic molecule and demonstrating state-of-the-art performance with a significantly reduced training demand.

Published in Transactions on Machine Learning Research (12/2023)
Local Function Complexity for Active Learning via Mixture of Gaussian Processes
Danny Panknin
danny-penkkiing dennknin@ttteerriinnno
Uncertainty, fnverse Modeling and Machine Learning Group, Berlin Institute of Technology, 10587 Berlin, Germany Physikalisch- Techhsssche Bundesonstall, 10587 Berlin, Germany
Stefan Chmicla
Stefan@chhicia..comm com
Machine Learning Deportment, Berlin Institute of Technology, 10587 Berlin, Germany BIFOLD-Berd Berlin Institute for the Foundations of Learning and Data, Germany
Müller
Machine Learning Deportment, Berlin Institute of Technology, 10587 Berlin, Germany BIFOLDDBerlisn Institute for the Foundations of Learning and Data, Germang
Department of Artificial Intelligence, Korea University, Scoul 136-715 South Korea Max Planck Institute for Informatics, 66123 Saarbrücken, Germany
Shinichi Nakajima
nakajjmnoBuu-brriin bereinnninneee
Machine Learning Deportment, Berlin Institute of Technology, 10587 Berlin, Germany BIFOL.D-Berlin Institute for the Foundations of Learning and Data, Germany RIKEN ATP, 1-4-1 Nihonbashi, Chuo-ku, Tokyo, Japan
Reviewed on OpenReview: https: openreviee. net/ forum? td- #MMMQ
Abstract
Inhomogeneities in real-world data, e.g., due to changes in the observation noise level or variations in the structural complexity of the source function, pose unique set of challenges for statistical inference. Accounting for them can greatly improve predictive power when physical resources or computation time limited. In this paper, we draw on recent thcoretical results on the estimation of local function complexity (LFC), derived from the domain of local polynomial smoothing (LPS), to establish notion local structural complexity, which used develop model-agnostic active learning (AL) framework. Due to its reliance on pointwise estimates, the LPS model class is not robust and scalable concerning large input space dimensions that typically come along with real-world problems. Here, we derive and estimate the Gaussian process regression (GPR)-based analog of the LPS-based LFC and use it as a substitute in the above framework to make it robust and scalable. We assess the effectiveness our LFC estimate an AL application on a prototypical low-dimensional synthetic dataset, before taking on the challenging real-world task of reconstructing quantum chemical force field for small organic molecule and demonstrating state-of-the-art performance with significantly reduced training demand.
1 Introduction
Inference problems from real-world data often exhibit inhomogeneities, e.g., the noise level, the density of the data distribution, or the complexity of the target function may change over the input space. There exist different approaches from various domains that treal specific kinds of inhomogeneitiess. For example, Kersting et al. (2007); Cawley et al. (2006) deal with heteroscedasticity by reconstructing a local noise variance function that used adapt the regularization of the model locally. Some approaches adjust bandwidths locally with respect to the input density (Wang & Wang 20077 Mackenzie & Ticu [2004] Moody & Darken 1989 Benoudjit et. al.] 2002). Inhomogeneous complexity can also be captured using combination of several

---

Published in Transactions on Machine Learning Research (12/2023)
kernel-linear models with different bandwidths, either learned jointly (Zheng et al. 2006 Guigue et al. 2005 or hierarchically (Ferrari et al. 2010 Bellocchio et al.] 2012). The most widely applicable models treat all types of aforementioned inhomogeneities in a unified way (Tresp) 2001 Panknin et al. 2021). Namely, they locally adapt bandwidths or rogularization according to the inhomogencitics in noise, complexity, and data density. While this is the path we will pursue, the focus this work will be on inhomogeneous complexity under the assumption of homoscedastic noise. In addition, we will investigate our proposed estimates in a heteroscedastic setting to demonstrate negligible practical limitations.
Exposing inhomogeneities sheds light the informativeness of certain locations of the input space, which subsoquently can be used to guide the sampling process during training valoong-also known as active learning (AL).
AL (Kiefer 1959; MacKay 1992; Seung et al. 1992 Seo et al.) 2000) is a powerful tool to enhance the training process of a model when the acquisition of labeled training data is expensive. It has been successfully implemented in various regression applications like reinforcement learning (Teytaud et al. 2007). wind speed
forccasting (Douak et al. 2013), and optimal control (Wu et al. 2020).
Nowadays, machine learning (ML) methods are increasingly deployed in physical modeling applications across various disciplines. In that setting, the labels that are necessary for model training are typically expensive as they stem, e.g., from computationally expensive first-principles calculations (Chmiela et al.) 2017) or even laboratory experiments. Due to the need for effective training datasets, AL has become an integral part of ever-growing importance in real-world applications, e.g., in the domains of pharmaceutics (Warmuth
et al. 2003) and quantum chemistry (Gubaev et al. 2018 Tang & de Jong, 2019 Huang & von Lilienfeld 2020
which raises the demand for AL solutions and the importance of AL research in general.
Through the advance of ML in scientific ficlds that hold the potential for significant impact, new regression problems emerge, for which there is initially only scarce domain knowledge while they simultaneously require thousands to tens of thousands of training samples for ML models to operate at an acceptable performance level. Regarding AL, these two characteristics of regression problems are hard to reconcile:
Due to insufficient domain knowledge on the one hand, a suitable AL approach shall be robust, since unjustified assumptions may result in a training performance that is even worse than random test sampling. By random test sampling, we refer to the naive training data construction that draws samples i.i.d. according to the test distribution. Additionally, the AL approach shall be model-agnostic since the state-of-the-art. is ever-evolving for this particular kind of regression problem. For these reasons, practitioners prefer model-free AL approaches for regression Ww 2019) over sophisticated, modellbased AL approaches with strong assumptions as the former are inherently ly robust and model-ngnostic by ignoring label information.
On the other hand, it is preferable that an AL approach outperforms random test sampling even at large training sizes. In the following, we will measure the AL performance by the relative required sample size e O, which asymptotically equates the performance of n e active training samples n random test samples (sce Definition 2 Accordingly, we call an AL approach asymptotically superior to random test sampling, if e <<< 1. Unfortunately, the performance gain of model-free AL approaches that we observe at small training sizes over rendom test sampling eventually diminishes completely (e= (@=1) 1) with growing training size.
For the described lcarning task, we therefore require an AL approach that is model-based but comes with mild regularity assumptions al the same time to feature robustness and model-agnosticity to a certain extent.
Recently, Panknin et al. (2021) h addressed the outlined AL scenario, where the fundamental idea is analyze the distribution of the optimal training set of model in the asymptotic limit of the sample size.
Assuming that this limiting distribution exists, they then propose to sample training data in top-down manner from this very distribution, knowing that with growing sample size the training set will eventually become optimal. By a lop-down AL approach, we mean an (infinite) training data refinement process æ1, 26.... such that when optimizing an AL criterion with respect to (#1...., (xf..... ,n1}} asymptotically becomes a respective optimizer as n - co. They have shown for the local polynomial smoothing (LPS) model class [Cleveland & Devlin 1988) that the asymptotically optimal distribution exists, whose density furthermore factorizes into contributions of the test density, heteroscedastic noise, and local function complerity (LFC) a measure of the local structural complexity of the regression function. Intuitively, LFC scales with the local amount of variation of the regression function. It is essentially estimated as the reciprocal determinant of the
2

---

Published in Transactions on Machine Learning Research (12/2023)
locally optimal kernel bandwidth (LOB) of the LPS model, calibrated for the local effects of the training input density and noise level. Given small but sufficient training set, these factors can be estimated, allowing the construction of the optimal training density and subsequently enabling the refinement of the training data towards asymptotic optimality.
While the previous work by Panknin et al. (2021) provides theoretically sound solution to our considered AL scenario, the required pointwise estimates that are inherent to the LPS model class prevent scalability with regard to the imput space dimension d. The goal of our work is to extend the above approach in a scalable way. The key idea is to build the required estimate of LFC based the LOB of the related Gaussian process (GP) model class that can naturally deal with high input space dimensions. Subsequently, we plug our scalable LFC estimate into the AL framework of the existing method, whose functioning is justified by the method's model-agnostic nature.
It particularly the almost assumption-free nature of LPS that made the results of Panknin et al. (2021) model-agnostic. To that effect, we base our results on the nonparametric, adaptive bandwidth Gaussion process regression (GPR) model to preserve this property. While we lose the strict asymptotic sampling optimality this way, expect it to be reasonably close due to the model-agnosticity nevertheless. We refer to the resulting training density as the superior training density for locally adaptive models, by which we refer to models that adapt to the considered inhomogeneilies, namely heleroscedasticity and inhomogeneous complexity. Note that it is first and foremost superior for our adaptive bandwidth GPR model.
Fig. summarizes all steps of our contribution and shows how they are interlinked and in which sections they will be discussed. Specifically, we contribute in two ways:
Theoretical contribution Assuming homoscedastic data, we propose GPR-based LFC estimate which is inspired by the design of the LPS-based LFC estimate. Here, we need to respect the scaling behavior of LOB of GPR that differs from the LPS case, where we use asymptotic results on the scaling of optimal bandwidths for GPR, as described in Sec. 3.3 Making use of the model-agnostieiit of the optimal training density of LPS, we replace the herein contained L.PS-based LFC estimate for our GPR- basedased LFC estimate to obtain a superior training density for locally adaptive models. From this point, we can implement the AL framework by Panknin al. (2021), for which we propose a novel pool-based formulation. Both, our LFC and density estimate will inherit. the scalability of the deployed GPR-based LOB estimate.
Methodological contribution We propose scalable LOB estimate for GPR the weighted average of bandwidth candidates, where the weights are given by the gate function of a sparse mixture of GPs model-a special case of a mixture of experts (MoE) model (Jacobs et al. 1991 Jordan & Jacobs 1994 Pawelzik et al.) 1996p. Here, each expert of the MoE is a GPR model that holds an individual, fixed bandwidth candidate. We construct the MoE in PyTorch (Paszke et al. 2019) out of well-established components from the related work and design a training objective that is regularized with respect to small bandwidth choices to obtain a robust and reasonable LOB estimate in the end. In addition, we provide implementation of both, the AL framework and our model. Finally, we propose novel model-agnostic way of choosing inducing points (IPs) of sparse GPR models: Respecting LFC, we place additional basis functions of a kernel method in more complex regions while removing basis functions in simpler regions of the input space.
MoE Mixture Experts
LOB Locally Optimal Bandwidth
LFC Local Function Complexity
Superior Training Density
Active Learning Framework
Sec.
Sec.
4
Figure 1: An overview of the steps of our contribution and how they are interlinked. We will elaborate on the main steps in the specified sections.
3

---

Published in Transactions on Machine Learning Research (12/2023)
To show the capabilities of our approach, we consider two inhomogeneously complex regression problems:
The Doppler function In a controlled setting of 1-dimensional synthetic data, we will first analyze our MoE model and the proposed estimates of LFC and the superior training density, where we will demonstrate the asymptotic superior performance of our superior sampling scheme and compare to related work. By the superior sampling scheme, we refer to i.id. sampling from the superior training density.
Force field reconstruction Quantum interactions exhibit multi-scale behavior due the complex electronic interactions that give rise to any observable property of interest, like the total energy atomic forces of
system (Bereau et al. 2018 Yao et al. 2018 Grisali & Ceriotti 2019 Ko et al. 2021 Unke et al. 2021a). To
demonstrate the scalability of our approach, we consider force field reconstruction problem of a molecule with 27 dimensions, where the application of the LPseb AL framework by [Panknin et al. (2021) is
intractable. Besides the asymptotic superior AL performance, we gain insights into the local structural complexity of this high-dimensional molecular configuration space through risualizations of the scalar-valued LFC function.
We begin by discussing our work in the context of related work in Sec. 2 Next, we give formal definition of the considered regression problem and the asymptotic AL lask, and review asymptolic results for LPS and GPR in Sec. B In Sec. 4 we describe our MoE model and derive the GPR-based LFC and superior training
densily estimates. In Sec. 5 we then describe our experiments and results, which will be further discussed in
Sec. We finally conclude in Sec.
2 Related work
Choice of MoE, experts and the gate The common assumption of MoF approaches is that the overall problem to infer is too complex for single, comparably simple expert. This is the case, for example in regression of nonstationary or piecewisc continuous data, and naturally in classification where each cluster shape may follow its own pattern. In such a scenario each expert of the MoE model can specialize in modeling an individual, (through the lens of a single expert) incompatible subset of the data, where the gate learns a soft assignment of data to the experts. Under these assumptions, the hyperparameters of each expert can be tuned individually on the respective assigned data subset. In the light of this paradigm, there exist several instances of mixture of GPs, for example, Tresp 2007
Meeds & Osindero 2006); Yuan & Neubauer 2009)
Yang & Ma 2011 Chen et al. (0044)
In our work, we aim to infer a single regression problem, where there is no such segmentation as described above: Each individual (reasonably specified) expert of our mixture model is eventually capable of modeling the whole problem on its own. Yet, if the problem possesses an inhomogencous structure, the prediction performance can be increased by allowing for a local individual bandwidth choice. Therefore, we deviate from the common MoE paradigm, sharing all those parameters across the experts that describe the regression function. This less common assumption was also made by Pawelzik et al. (1996), schre--localy dependent dependent-some some experts are expected to perform superior compared to the others.
For the expert and gate components of our MoE model, we focus on the sparse, variational GPR model (see, e.g., Hensman et al. (2015)) trained by stochastic gradient descent. However, there exist other (sparse) GPR approaches that could be considered for the gate or the experts of our MoE model. Some are computationally appealing as they solve for the inducing value distribution analytically (Seeger et al., 2003 Snelson & Ghahramani 2005; Titsias 2009) or do not require inducing points in the first place in the case of full GPR model (see, e.g. Williams & Rasmussen 1996)).
Particularly the expert models of our MoE model can be exchanged for arbitrary sparse and full formulations of GPR, as long as we can access the posterior predictive distribution. We give a short summary of these model alternatives in Appendix B.1 and B.2 which are also included in our provided implementation.
For the gate, however, analytic approaches come with complications as they require labels. Such labels do not exist for the gate, and so we would need train the MoE in an ezpectation mozimization loop, where the likelihood of an expert. to have produced training label functions a pseudo-label to the gate.
4

---

Published in Transactions on Machine Learning Research (12/2023)
As opposed to standard GP that features stationary covariance structure, our MoF, model with GP experts of individual bandwidths can be interpreted as a nonstationary GP. Apart from MoE model architectures, there exist other approaches to construct a nonstationary GP:
Alternative nonstationary GPs
Closely related to the MoE model, (Rullière et al. 2018) proposed to aggregate individual expert model outputs through a nested GP. Note that this approach is not well-suited for our purpose as the aggregation weights lack the interpretation of a hidden classifier, leaving the subsequent LOB estimation (as in (24) open. (Gramacy & Lee 2008) deploy individual stationary GPs on local patches of the input space that are given by an input space partitioning of tree.
Gramacy & Apley] 2015) identily subsets of the training data that are necessary to resemble the GP covariance structure at each individual evaluation point. Careful bandwidth choice for each subset then yields a nonstationary GP as well as local bandwidths. Roininen et al. (2019) obtain local bandwidths by imposing a hyperprior on the bandwidth of GP. Note that our work intends to elaborate the estimatiio of LFC and agnel-agnostic superior training, given any estimate of LOB of GPR. We deployed an MoE approach as a simple means to obtain these estimates. The MoE component in our LFC and superior training density estimates may be readily replaced by the approaches of Gramacy & Apley
2015 or Roininen et al. 2019
In deep Gaussian process (DGP) regression (Damianou & Lawrence. 2013). the inputs are mapped through one or more hidden (stationary) GPs. This warping of the input space yields a nonstationary covariance structure of DGP. An approximate DGP model through random feature expansions was exercised by Roininen et al. (019)) Sauer et al. (2023b) discuss active learning for DGP regression. We will implement the DGP model as well as the AL scheme of Sauer et al. (2023b) and compare the AL performance of our superior AL scheme on this very model to demonstrate the modellagnostieity of our work in Sec. 5.1
IP selection In our work (sce Sec. 4.5 we choose the IP locations of the gate and the experts of our MoE model in a diverse and representative way but also in alignment with the structural complexity of the target function, interpreting this choice as a nested AL problem. There are a variety of TP selection approaches in the literature. Zhang et al. (2008) interpreted the choice of IP locations from a geometric view that is similar to ours: They derived a bound on the reconstruction error of a full kernel matrix by a Nyström low-rank approximation in terms of the sum of distances of all training points to their nearest IP. This exposes a local minimum by letting the TP locations be the result. of k-Means clustering. This choice of TP locations is representative and diverse, while it solely considers input space information. In this sense, our approach extends their work by additionally considering label information. This and our approach draw a fixed number of IPs at once. There are also a lot of Nyström method based TP selection approaches that select columns of the full kernel matrix according to a fixed distribution (Drincas ct al., 2005) or one-by-one in a greody, adaptive way (Smola & Schölkopl 2000 Fine & Scheinberg, 2001 Seeger et al. 2003). An intensive overview of Nyström method based IP selection methods was given by Kunar et al. (2012), where they also analyzed ensembles of low-rank approximations. We compare our proposed IP choice 29 to the greody fast forward IP selection approach by SSeger et al. (2003) in Sec. 5.1
Moss et al. 2023) incorporate a quality function into a diverse IP selection process that can be specified flexibly. They consider Bayesian optimization rather than regression, they exercise a quality function proportional to the label. However, other measures of informativeness that are better suited for regression could be deployed. Note that LFC would be a possible candidate for this purpose.
The AL scenario In this work, we consider model-agnostic AL with persistent performance at large (or even asymptotic) training size as opposed to the common AL paradigm that is concerned with small sample sizes. In this sense, we delimit ourselves from AL approaches that are tied to model, e.g., when they are based on a parametric model (Kiefer 1959 MacKay
1992 He 2010 Sugiyama & Nakajima 2009 Gubaev el al. 2018), ; or which refine training data bollom-up in a greedy way to maximize its information content at sample size, where the information is either based on the inputs only (Seo et al. 2000 Teytaud et al,
small
2007
Yu & Kim 2010 Wu| 2019 Liu et al., (021) or also incorporates the labels (Burbidge et al. 2007 Cai 2013). By botlom-up AL approach, we mean a training data refinement process that is constructed by choosing the nh imput 2, the optimizer of an AL criterion with respect to {I1,..., Zn), when keeping
ee al.
the previously drawn inputs In-11 (with labels (yI...., Bn-1) }) fixed.
5

---

Published in Transactions on Machine Learning Research (12/2023)
Our work is therefore complementary to the latter kind of approaches which can be better suited in another AT scenario. For example, if there is enough domain knowledge such that we can deduce a reasonable parametric model without the need for a model change in hindsight, an active sampling scheme based on this model will be best. Our category of interest is for the other case, when domain knowlodge is scarce, where we have no idea about the regularity or structure of the problem to decide on a terminal model. Here, for small training sizes (and particularly from scratch), input space geometric arguments (Teytaud et al. 20077 Yu & Kim 2010 Wu, 2019 Liu et al., 2021) are applied in practice. However, as already noted in the introduction, their benefit is limited to this small sample size regime, which we will demonstrate
synthetic dataset. They serve reasonably for the initialization of supervised AL. approaches, including ours, nevertheless.
Regarding our considered AL scenario, Panknin et al. (2021) have recently proposed an AL framework based on the LPS model class, where training samples are added so as to minimize the mean integrated squared error (MISE) in the asymptotic limit. This approach is therefore provably asymptotically superior to random test sampling. Additionally, it is robust since the LPS model is almost free of regularity assumptions. Finally, their LPS-based solution then showed to be model-agnostic: On the one hand, this is indicated tcalle by the fact that the LPS model has only indirect influence on the asymptotic form of LFC and the optimal training density since the predictor is asymptotically not involved (sec, c.g., Eq. (8) On the other hand, this is validated empirically by assessing the performance of their LPS-based training dataset construction under reasonable model change in hindsight. This model change is restricted to locally adaptive models. Here, Panknin et al. (2021) observed a consistent performance superior to random test sampling when training a random forest model and a radial basis function proposed training dataset.
(Moody & Darken 1989). using their
AL for classification Note that the outlined AL scenario can be solved more easily for classification:
Here, AL is intuitively about the identification and rendering of the decision boundarics, which is inherently a model-agnostic lask. In addition, since the decision boundaries are a submanifold of the input space X, a substantial part of X' can be spared when selecting training samples. Therefore, AL for classification leverages the decay of the generalization error from a polynomial to an exponential law (Seung et al.) 1992) over random test sampling. For the above reasons, AL for classification has been applied successfully in practice Roy & McCallum,
Lewis & Gale; 1994
2001
Goudjil et. al., 2018 Warmuth et al.] 2003 Pasolli & Melgani
2010 Saito et al. 2015 Bressan et al.) 2019 Sener & Savarese, 2018 Beluch et al., 2018 Haut et al.) 2018
Tong & Chang, 2001 He
2010). In rontrast, the performance gain of AL for regression is more limited in the sense that, under weak assumptions, we are tied to the decay law of the generalization error of random test
sampling (Györfi et al. 2002 Willett et al. 2005).
GP uncertainty sampling There exists a lot of research on AL for GPR (Seo et al. 2000 Pasolli & Melgani |011] [Schreiter et al., [2015] Yue et al., 2020), which is typically based on minimizing prediction uncertainties of the model. With our proposed AL approach being based on GPR models, this research area is the most related competitor to our work.
For a standard GPR model, the prediction uncertainty is the higher the farther away we move from training inputs. In this way, GP uncertainty sampling samples (pscudo-)uniformly from the input space which makes up for a low-dispersion sequence (Niederreiter 1988) (see Definition 5 Note that standard GP uncertainty sampling is an input space geometric argument since it does not depend on the regression function to infer. As already indicated in the introduction and as we will show in Sec. 5.1 input space geometric arguments feature no benefit regarding asymptotic AL performance.
Since our model is mixture of GPR experts, it is straightforward to derive its uncertainty as mixture of Goussian process uncertainties (MoGPU) by simply weighting the predictive variances of all experts with respect to the gale output (see 32 for a definition). As opposed to GP uncertainty sampling, MoGPU can cope with structural inhomogeneities. Therefore, we consider MoGPU as a fair baseline competitor to our superior sampling scheme and compare both in See. 5.1
6

---

Published in Transactions on Machine Learning Research (12/2023)
3
Preliminaries
We will now give formal definition of the regression task, the AL objective and LOB, and short review of the asymptotic results on LFC and the optimal training distribution of the LPS model in Sec. 3.1 and
3.2
Then we recap asymptotic results on the optimal bandwidth of GPR in Sec. 3.3 and known models in Sec. 3.4 that will serve as building blocks of our proposed adaptive bandwidth MoE model later on.
In the following, we denote by diag(z) E Rdxd the diagonal matrix with the entries of the vector z Rd on its diagonal and by La diag(1,) the identity matrix, where 1 is the vector of ones in R
3.1 Formal definition of the regression task and AL objective
Let f be the target regression function defined on an input space X' C Rd that we want to infer from noisy observations y; = fx;) + En where z; E are the training inputs and E; is independently drawn noise from distribution with mean F[&y] = and local noise variance V[s,] = v(x;). We denote a training set. by (X,, n. Yn), where X, = (II...., In) E x and Y, = (yI,..,, .n)) E R For a given model class ʃ that returns predictor fx,,x, for training set (X,,, Y,), can define the pointwise conditional mean squared error of J in x E 2, given X,, by
MSE (a.flx,)) Ey, (x)
Ec, ((x,, (2) (x)
f(x))²
(1)
Note that via marginalization the conditional mean squared error is no function of the training labels Y.. Given a test probability density g E Ct (xt R,) such that q(x))dx = 1, the conditional mean integrated squared error of the model under the given training set is then defined as
MISE (I.
/ MSE (I, f1,,) g(x)dxx
(2)
With these preparations, the AL task is to construct a training dataset (X,, Y2 such that
X ~ argminx, x,Ex" MISE (q. f(x.).
(3)
3.2 Locally optimal bandwidths, function complexity, and optimal training
Let. - be a family of kernel machines which characterized by a positive definite bandwidth matrix parameter Σ E sd of an RBF kernel (=xxx') == ΣΣ |زーー x(| E((E-1(< - z'Il)) for a monotonically decreasing
function Φ: R. + R, The well known Gaussian kernel example implemented by d(z) = expl -
Given bandwidth space S C sf we define the LOB function of f by
my(0) = argminges MSE ((x
(4)
assuming that this minimizer uniquely exists for all x E X.
Denote by m the predictor of the LPS model of order Q under bandwidth Σ and by Eg = Em the LOB function [4] of LPS, if it is well-defined. This is the case, e.g., for the isotropic bandwidths space S = {oI, σ >> 0f under mild assumptions]] where we particularly can write (x) o (z) We refer to Appendix for details on the LPS model and asymptotic results, For the optimal predictor
FQ
== ma
(g(a)
(x)
5))
"For LOB being well-defined in the isotropic case, we generally require a non-vanishing bias and variance in terms a cnla-variannedeommaittioon of the MSE of the predictor in x, for all a E X. See, e.g., Eq. 36 for the LPS predictor m or Silverman (1986): Wand & Jones (1994) in more general.
7

---

Published in Transactions on Machine Learning Research (12/2023)
of LPS, letting fups in Eq. (3). Panknin et al. (2021) have shown that there exists an optimal training density Popt that allows the optimal training imputs in Eq. (3) to be asymptotically obtained by independently
and identically sampling X ~ 2nn. They have also shown that this density exhibits a closed-form
HQ+1+4
2(++1) RI+HTt v(() TG 1T42 (1 + o(1),
Pop (x) (z)(2)((z)
(6)
where for an arbitrary training dataset (X,, 12 Yn with X, ~ ps the LFC of LPS is defined by
v(x)
aQ+1++
[a(+)| = ((()
212-1+
CQ (c) = p(r)
(7)
(p(x))
The LFC in 7 asymptotically solcly depends on the behavior of f as opposed to U, and n: It scales with the local varialion of f in the vicinity of x. For example,
(x) x trace(D)(())
अक (1 (1+ o1))
(8)
is a function of the trace of the Hessian of f (Fan et al. 1997
Q.n The optimal density Pop- in 6 implics that we require more training data where the problem is locally more complex (large (%) or noisy (large v), or where test instances are more likely (large g). As already noted in the introduction, the results to LFC and the optimal training density of LPS indicate their problem intrinsic nature, as they reflect no direct dependence on the LPS model except for the order Q. Note that for f E C (X, R), there a canonical choice Q = [a] - of the LPS model order. When deriving LFC under this canonical-order model, we consider the dependence of the associated LFC and the optimal training density Q negligible, its choice is driven by the problem intrinsic regularity.
In practice, we obtain X popt by estimating Eq. (6) and 7 from (X,,, Y.)) with X ~ p for an arbitrary training density p, where n' < n, followed by adding the remaining n - n inputs appropriately (see Sec. 4.2)
Q.n
The construction of pon crucially depends on reliable estimates of LOB as the key ingredient for the estimation of LFC. While Panknin et al. (2021) provide such an estimate based on Lepski's method (Lepski 1991 Lepski Spokoiny] 1997), it does not scale well with increasing input space dimension d. This is because pointwise estimates suffer from the curse of dimensionalitt regarding robustness and computational feasibility. The goal of this work is to implement the above AL framework but based on functional LOB estimate in the domain of GPR instead of TPS,, since the GPR model class can naturally deal with high input space expect that LOB estimates
dimensions (Williams & Rasmussen 1996). Relying the model-aggosticityy, based on LPS can be exchanged for LOB estimates based on GPR in the formulation of LFC and popt when matching the degree Q to the smoothness of the regression function appropriately.
3.3
On the scaling of GPR bandwidths
The major difference between LPS and GPR is that keep fixed model complexity in the sense of the number of basis functions in the former while there is varying model complexity in the latter we add further training instances. E.g., under the Gaussian kernel the model complexity of GPR grows infinitely. When the regularity of the kernel and the target function f match, then, as soon as the training size n becomes large enough, there is no need for further shrinkage of the bandwidth to reproduce f with GPR in the asymplotic limit. In particular, given enough samples, there no need for local bandwidth adaption.
However, there mismatch if f Ca (X, a-times continuously differentiable since the Gaussian kernel is infinitely often contimuously differentiable. As shown by Van der Vaart et al. (2007) 2009), in order to obtain optimal mininaaz-convergemn of the predictor (except for logarithmic factors), the associated (global) bandwidth has to follow the asymptotic law
E I n
Note that for f E C* (X, R), where the thcoretical results of LPS apply, the scaling factor n 2012 of LOB in sample size matches exactly for both classes, LPS and GPR. In our work, we will use ( to deduce a GPR-based LFC estimate in analogy to the L.PS-based LFC estimate (7) by Panknin et al. 2021
8

---

Published in Transactions on Machine Learning Research (12/2023)
3.4
Preliminaries the applied models
We will now introduce the models that we implement in this work. For the RBF-kernel k, we define the kernel matrix between X E X and X'' E x as K (X,, = [k* (2,x²)) xex As shorthand notation we furthermore define K (X) = K (X,
3.4.1 Sparse variational Gaussian processes
We define the sparse GPR model ì ~ SV) (sce, c.g. Williams & Rasmussen (1996); Hensman et al. (J,, X, V, E, 1 ,4, S²), which are the global constant prior mean μ, the regularization parameter A, the label noise variance function 0, the bandwidth matrix of the kernel and the prior distribution, given by the IP locations X, E Am as well their inducing value distribution, characterized by the moments μ4 and ST. That is, for the inducing values Y of X, we assume Y, (xX,) ~ M,, Here, the degree of sparsity is described by m IPs: This number can be fixed in advance gradually increased with training size n, where the increase m. o[n] is typically much slower than n. If we can assume homoscedastic noise, we let 0(x) = o.
(2015)) as follows: The sparse GP is described by the (hyper-) parameters
The sparse GP then outputs
(XX)) ~ (X
(10)
for the mean function
('())) K,,KK,1/(11 K, 1/1
u(xx)) (T-111((x,)) +M(X.,,
(11)
and the covariance function
Imm) C*(X.) K. + K.; 1/2 -m, 1/2
K, KT + diag( (X.)),
(12)
where i K 1/2 , and S K. 1- K. 1/2 are the whitened moments of the inducing value distribution Pllciss et al. (2020). Sec. 5.1), and we have defined K. = (x,), K, = (xx,,, and K,, = (X,,, X,).
We choose la to be the constant mean function, i.e. (x) plln for X Xm, noting that other mean functions are possible. Note that test predictions for(x) = 1 (x) are given by Eq. [1
The training objective Let P denote the prior distribution of the inducing function values Y, of the IPs X, and let denote tractable variational distribution intended to approximate P(.) Q(- Y,).- In variational inference, we want to minimize the Kullback-Leibler Leibler divergence K.C [Q][P] between Q and P, which is equivalent to maximizing the data log-evidence log(P(X.,, Y,)). As tractable approximation, we maximize the evidence lower bound (ELBO), given by
E Y,) (X Y,lw)) KCIQOXX Y.
B KC Y.),IP),
(13)
where Ps the predictive likelihdl in Th, marginalized over the variational distribution Q, that is,
Ph = Eu-Q(),,,..) log f P(wl))))(lu 2n)df.
(14)
3.4.2 Sparse mixture of experts
Given a finite set of expert models ÿi that are parameterized by oe,, the MoE model is given by
Pruee (x)
Σ 1= (())
(15)
where the gate G: X + 1 is a probability assignment of an input x to the experts. In particular, it holds Σ G(r); = 1 and G(z), 2 0, E X' and 1 L.

---

Published in Transactions on Machine Learning Research (12/2023)
We implement the approach of Shazeer et al. 2017 to model the gate G as follows: For the softmar function
soft max a); == expfai) 1 exp(ai):
(16)
where a E R", Shazeer et al. (0017) propose to set
his (z)
,I <
G(x) mtf (s),...,,,((()))
where (x(
(17)
-00
,12nt x
for an adequate permutation (f1,...,, (L) of {I,...: L such that h2 (x) > hics (x) are ordered decreasingly. Here, hy(x) (x) + -V(0, st) noisy version of single-channel gabing models Φi with parameters 0gr Note that these models can be chosen freely and may also deviate from the choice of expert models ğ
The cutoff value 1 S K < L controls the sparsity of the MoE, as it enforces the minor mixture weights to strictly equal zero. For stability reasons, during the training, we give cach expert a chance to become an element of the top-s components by adding independent Gaussian noise (O, sf) before thresholding, where E R4, another hyperparameter to set or learn. This noisy gating prevents a premature discarding of initially underperforming experts.
The overall MoF hyperparameter set is thus given by
= (f0c)) 10 1,,,,,))))) s),
(18)
4 Estimating locally optimal bandwidths via mixture of Gaussian processes
In this section, we derive our main contribution, namely the GPR-based AL framework, which we summarized
in Fig. 2 We first derive our GPR-based estimates of LFC and the superior training density in Sec. 4.1
(Fig. 2 B). Combining this estimate with the AL framework from Panknin et al. (2021), we obtain GPR-
based, model-agnostic superior sampling scheme in Sec. 4.2 (Fig- CC. Next, describe our GPR-based
MoE model in Sec. 4.3 (Fig. 2 A) of which we obtain the required LOB estimate of GPR for the estimation of superior training density. The scalability of this estimate enables the application
our superior sampling scheme to problems of high input space dimensions. We then give details on the training of the MoE in Sec. and finally propose an L.FC-based TP selection method in Sec. 4.5 We summarize the pseudo-code of our superior sampling scheme in Algorithm !
4.1 GPR-based LFC and the superior training density
Let EGPR (x) denote the LOB function 44) of GPR. Inspired by the results LFC and the superior training densily of LPS in Eq. 7 and ( we are able to deduce their GPR-based analog. Here, we need to take into account that GPR adapts universally[] to functions f E C (X, R.,, opposed to LPS, whose decay rate is determined by the specified polynomial order Q. The idea of the LFC estimate was to adjust LOB appropriately so that it becomes invariant under the influence of the training density, heteroscedasticity, and its global decay with respect to the training size n.
Combining the local effective sample size p(z)n with the scaling result the global EPR in @ffrom Sec. 3.3 we propose an LFC estimate for GPR as follows (see Appendix C for proof details).
Theorem 1 (LFC of GPR). For Ca (X, X, ~ p and homoscedaslic noise, the GPR-based LFC estimate of f in x E X is asymptotically given by
2o4d
CapR(()) = [())n
DGPa(())1
(19)
In analogy to Eq. 7
COPR measures the structural complexity of f, as it asymplotically does not depend on p. D and n. Note that the LPS model provides no explicit way to adapt to the local noise variance v(z),
9That the MISE decays the minimax-rate
20r4 of nonparametric models.
10

---

Published in Transactions on Machine Learning Research (12/2023)
A
MoE
Experts
CP()))
GiDi
(gP
Gate
G (a)
X Xnb, Y nkk
GPR, (
nkk
CCR (x)
Psup
C
B
GPR, Lj Psup
k ++=1
query label Yi of Di nk + +1,- 7k+11
Figure 2: The proposed AL framework.
such that the LOB of L.PS scales with respect to u address heterosccdasticity (see (36) in Appendix A For GPR, we have made the restriction of homoscedastic noise in the definition of C in Theorem since not aware of theory on the scaling of GPR-based LOB with respect to heteroscedaatiiityy However, opposed to LPS, heteroscedasticc GPR model provides explicit way to adapt to the local noise variance v(z) via regularization. As a result, we observe only very little influence of heteroscedasticity on LOB function, which we will demonstrate in Sec. 5.1 Thus, GPR will be sufficiently calibrated in a heteroscedastic scenario, making a reasonable estimate LFC in practice without further restrictions.
Now, when putting CRR into Eq. 6 with Q -a-1, we obtain the superior training density
Psup
(P).
[G (x) (z))(())
2adid lotd v(x) Bo-d (1+0(1)))
(20)
GPR,1 and Pops coincide, which proved be optimal for For even a E IN with a and CGPR Cz Psup
LPS. In this sense, (20 generalizes 66) to the general case of RS, where expect that the true optimal training density for f C (X) will not deviate by a lot from Paup
GPR.n Since LPS and GPR are related models, we furthermore expect CRR to be similar to Cg for the appropriate order Q.
Note that for C (X), let 0o in Eq.
19 and 20 to obtain
Co (x = nn (T)(
and
Psup
GPR, (x) x [C(mn())(())())]
(21)
While X ~ Psup ~ naive random test sampling, i.e,, X, ~, due to the model-agnosticity the LPS-based result. To assess the asymptotic performance of a training density p (such as Psup GPR."). let us first observe the following:
GPR, will not be optimal for our model, we expect to be asymptotically superior to the
For regression problems and under weak assumptions the law the MISE does not change with respect to p, except for constant multiple (Györfi al.) 20022 Willett al. 2005p. Accordingly, the number of actively selected training samples (~ p) that are required to achieve the same level of accuracy of random test sampling is given by constant 0 0. Formally, we define & follows.
11

---

Published in Transactions on Machine Learning Research (12/2023)
Algorithm 1: Superior training data process (X,, WNEN with labels Y, of training inputs X, 4 GPR,
Paup
Input
1: Intermediate training sizes (ME)EE, with ny nk+1, vk IN, for rocstimation
2: labeled validation set X vvll
3: An imput generating process X pool E X with pool ~ Px
4: The label oracle y:: X R
5: (optional) The test density
6: (optional) The intrinsic dimension d d of the imput space X'
7: (optional) The regularity of the target function E C (X,R)
Output
8: (Iufinite) training data process (x) with labels y, of training imputs X, PSug
GPR,
Procedure
9: & Initialization
10: Estimate pool density Px based on X pool
e.g., using kernel density estimation
q not specified then
Set q Px
is not specified then 14: Estimate 5 based on pool
Di e.go, following the work of Facco et: al. 2017
15: if is not specified then
16: Set a-
as discussed in Sec.
17: Set 4-
18: Draw initial training inputs X
19: Query labels Y, y( np) from the oracle
Set (0H,E Σx) hyper init(XXng ng Yng Po. vaal Yval) Sample Process
Algorithm
Appendix
for &t NN
> then
24:
Update TP locations X X E On, where X8, XG ~
CaPR
29
in Sec.
4.5
25:
if k == 1 then
26:
Gradually decrease
x and mG x as long as the validation performance of F does not degrade
as discussed in Sec. 143
Train the model MMUE from Sec. 4.3 with hyperparameters Oy on (Xng Yng) as described Sec. 4.4
Estimate the LOB GGR of GPR according
24
Estimate LPC CPPR
[x)
Zard
S" GGPR (1)
according 22
-CPR.ng
30:
Estimate the superior training density Psup (GPR.ME(s) (x) and 7/2 max {0, (0.5 -0----71 7/11/(11 (1-
6(x) उu according to
23
31:
Set 7 max Pt(tt)
Päw
Sec.
4.2
Set P4+1 7/Pk + P/2 GPR.
32:
Pk+1 < 2px+1 Ph
D see Soc.
4.2
34:
Draw Xngti...., Xmess ~ pk41 via importance sampling from X poot as descrihed Sec.
4.2
35:
Query labels tngtli-.., Umats =(Nnatls..... Xnuss) from the oracle
Set X Xny U (Xux +l. Xne+i) and Yuuss Yng U (Mog +1..... Mo
Then X Pk-1
Definition 2. Over the space of nquare-integrable functions f L (X), for a nonparametric regression model f and training density p, define by e(f.p) > 0 the relative required sample size such that for n! = e(f.p) X', ~ and X,, ~ q it holds that
MISE fl.x)))) MISE (1.f%%))) (1+ +(())
Thus, a training density p is asymptotically superior to random test sampling, if e(f.p) p) < 1, since we achieve the same performance as random lest. sampling with only fraction of the number of training samples. GPR.n for our GPR- based MoE model.
Sec. we will demonstrate the superiority the training density Pau
Respecting the intrinsic dimension in high-dimensional input spaces In Eq. 19 assume the input space X' to have full degrees of freedom which in practice is particularly not the case
20 and (21 we
12

---

Published in Transactions on Machine Learning Research (12/2023)
in high-dimensional feature spaces. For an intrinsic dimension 8 < d of X, we adjust follows: For the space S = {aE] >> 0} of bandwidth candidates that are essentially isotropic up to fixed, shared positive definite factor E' E S that is, e.g., calculated in a pre-processing step, let (T) OGPR (2)2 be the LOB function of GPR with respect to S. Then we replace all occureences of d for and (_GPR()) (2) for (x)
in Eq. 19
20 and 21
Besides being an ingredient AL, LFC can also be used to reduce the required model complexity. For example, in an RBF-network or a sparse GPR model, we coarsen or refine the model resolution by placing adequate amount of basis functions IPs, respecting LFC. We will discuss this choice in Sec. 1.5 and demonstrate its ability to reduce the overall model complexity in Sec. 5.1 Finally, LFC can be inspected to obtain deeper insights into the research field of the regression problem, which is particularly hard for nigh-dimensional data (see Sec. 5.2
4.2 The active learning framework
Starting with an initial training set X Y.. of size no with Xn, ~ Po for some initial training distribution such as Po = q, we implement the online sampling procedure as described in Panknin al. (2021) such that
X, ~ Psup
GPR, as n + oc. We grow the training set as follows:
Given the current training set X ny my we estimate
Sny -GPR as described in Sec. 4.3 Using
19
20
it is
Ch (x) x [1/px(x)]
1642
Sm, (x) and
(22)
2a+d (CPR(())(())
Psup (x) (x))).";())>
GPR..
n
((x) 4044
(23)
Letting the next sample size be ns+1 2m we have already drawn half the samples of 7x+1 according to a potentially different distribution Pk. than the new proposed Poup
,PRR,, The closest we can get in distribution
.PR,.2 is given by X Xnxt Pk+1, where P6+1 == 722k + (1 (1-7))11 Pau Peup
GPRMM for 72 = max {
{0, 0.5 71 E [0, 0.5)
7/
(x)
and 71 = max GPR.n (x) SEX Psup
This is achieved by sampling Inethy Ingi1 p for Px+1 2px+1 Pk.
which is a valid probability density (Panknin et al., 2021
Adaptions in the pool-based active learning scenario In the AL framework described above, we deal with properly normalized probability densities. But in the pool-based AL scenario such normalization is usually impossible since our information about the input space X' is restricted to a large, unlabeled pool of samples X poot X This pool follows a distribution X poot Px, for which it is common to assume an (unnormalized) density estimate px to be given: Unlabeled inputs are considered cheaply accessible, whereas querying labels is expensive.
GPR,,, at equal For our AL framework to be applicable, it suffices to keep all considered densities such as PeuP norm, which we can enforce via normalizing a density p by p/ p/norm(p), where
norm(p) X pool Σ REX p(x)/Dx(x).
To see this, note that first of all Px is an unnormalized estimate of Px such that we can write Px s c Px for some unknown constant e >>0 O. On the one hand, it is f (x(x))d by definition. On the other hand, it is
norm(p) ~ p(x) x Px (2)
/ p(z)dr,
1 such that also / ((r)dr= norm(p) / p(r)dr s c holds for any unnormalized density p.
Subsequently, the required samples Fnstl..... Ire-s Ps+1 are obtained via importance sampling from the pool with importance weights P(x;; æ) pk+t (x) (tcl)//oommp/ (p (x)) for E X poot and ni i & net1-
13

---

Published in Transactions on Machine Learning Research (12/2023)
4.3 Sparse mixture of Gaussian processes
Recall the sparse MoE model
15 from Sec. 3.4.2 given by
PM (2)
Σ L G(+)((())
where the gate G: X + [0, 1/- a probability assignment of an input x to the expert models ih- In particular, it holds Σ t-1 G(z); = and G(æ) 2 0, Vx E X' and 1 iiis L. According to (18) besides the expert and gale model parameters (0., and (0 H= this MoE approach has two hyperparameteess, x. and s, for controlling the sparsity of the gate and adding noise to the gate responses during the training to escape local optima.
We choose the expert models well as the single channel gating models to be sparse variational GPs Williams & Rasmussen 1996 Hensman et al. 2015), that is, Vi (V,))(e..) and ai SVGP((g:) which are parameterized by Oeg and On as described in Sec. 3.4.1 The overall MoE hyperparameter set is thus given by θ = ((ee)/i (8 k)..
We will keep certain hyperparameters of θ constant after initialization, and share some hyperparameters across experts and the channels of the gate: While the covariances of the inducing value distributions S E f, 0 E @ could be full positive definite matrices, we apply S throughout, giving favorable stability and computational efficiency. For the same reasons, we fix the inducing point (IP) locations X, E 00,0 O after initialization. Furthermore, we share the IP locations among the experts, respectively the gate channels, such that for X, E 0e, we apply X, x and for X E Og we apply X; =x for all 1 <I
In this work, our goal is to fit a single, coherent regression problem by a MoE approach. Therefore, we propose to share all the parameters across the experts that characterize the regression function rather than the expert model. That is, we share the mean (E,, the regularization parameter λE, and the noise variance function é, respectively the global noise variance of with 0(x) of in case of homoscedasticityy Furthermore, we apply fixed, logarithmically spaced set of individual expert bandwidth scaling factors σi <... < ar that are multiplied by fixed, shared bandwidth matrix Eg. Our expert parameters therefore reduce to
oz = (ME: dE. G, D, x (1.00)
Remark 3. Recall from Sec. 2 that it is possible to replace the variational GPR expert models for full well as sparse analytic GPR formulations (see Appendia D With slight abuse of notation, these cases are subsumed by setting 1 = @ or X = x = 0 for sparse, respectively full analytic GPR.
Since our objective does not incorporate any likelihood about the gate's output, there is no noise function to fit for the gate, such that we set Ô = 0 for ē E 0g. and all 1 < l < L. Each output channel of the gate poses its own classification ion problem, which is why we do not share the means. Yet, we share the regularization parameter and the bandwidth, as the individual channels should be structurally similar. Our gate parameters therefore reduce to
0g: (Ugn do, O, aoIa, X,, 14,4
After training described in Sec. 4.4 this MoE' can cope with varying structural complexity through
the individual bandwidth scaling factors oi of the experts and heteroscodastic noise through the adaptive regularization. Additionally, we can now use the gate of our MoE to propose an LOB estimate of GPR.
A GPR-based LOB estimate After training of the MoE, we use the learned gate G from 15 to predict
En EPR (x) as
Son G = GCP Ex) where
@th @G) = exp Σ 1. " G()))oo(()))
(24)
Due to the finite candidate set o G(z); P(EGPR Eg If, in fact, EaPR Eg, Egt holds true, then there exists an index function j(x)) (1,.... I.) such that In (x) = Hay In this case, we able to exactly recover LOB
......,,,, Tr we are limited to mcasure a quantization of Σn E (x) through
14

---

Published in Transactions on Machine Learning Research (12/2023)
1, 1, )=x()) In any other case, the estimate (24 of LOB a reasonable interpolation, which
with G(z))
0, clsc
deviation from 2Grn can be controlled by the number of bandwidth candidates of the MoE.
4.4 Model training
This section is devoted to the training of the model described in Sec. 4.3 We first set up the training objective
in Sec. [4..1..e and describe the training procedure of our model in Sec. 4.4.2 where we identify hyperparameters
of the approach. Then, we discuss how to choose the essential hyperparameters systemalically on the initial training dataset in Sec. 4.4.3
4.4.1 The training objective
First, will set up the objective function for training our MoE model in batch mode.
The main objective Denote by Ø C (I....,n] nu the indices of a batch, and let wg E bEB w(xs) for
the training importance weight function w x q/p. Let Pt be the prior distribution of the inducing function
values of the 1-th expert and Qi the corresponding variational distribution defined in Sec. 3.4.1 We choose
the (through the gale G) weighted sum of the individual expert negalive ELBO objectives
13 denoted by
Obj(X Y,, B, () = - "B Eveo (or)) 1 LKC[Q((IXX, Y, x.MlIP]]
as our main objective, where nu nwB/MB, for VB,1 Eseg va(os) with u((x) = (())w((x and
Po.1 Eu X, log / P((mll)) m))
is the prodictive log-likelihood 14 of the 1-th expert in 2s, marginalized over its variational distribution Qi-
A penalty on small bandwidth choices In the spirit of Lepski's method (Lepski 1991 Lepski & Spokoiny 1997], we prefer the largest choice of bandwidth out of all candidates that perform comparably well. In order enforce this, we penalize smaller bandwidth choices by adding the following term:
Σ 1=1 PB, (-1--)) Σ VB,1-
pen, (X n. Y,, B, w. 4) = (L. 1I-
(25)
Note that pen,((,, Y, B, w 0) = 1 VB,1 =... VB Our total objective then amounts to
Obj X n. Y, B, w. +) = Obj(X n. Y,, B, +) + Oxpen, Y, B, w. O).
(26)
Remark 4. If we assume our problem to be (almost) noise-free, we replace the Obj in our objective the mean squared error
26 for
MSF X Y B, w. 4 = W/' Σ hes w(ra)llyb 2
4.4.2 Training procedure
We implement our model in PyTorch (Paszke et al. 2019), using the CPyTorch-package (Gardner el al. 2018 Given the training set (X,, Y,), minimize the objective described in Sec. 4.4.1 via. ADAM-optimization (Kigmma & Ba 2015). It remains to identify those variables of the MoE that will be adapted as parameters during the training. Then, the remaining variables are hyperparameters that need to be specified or tuned through an external validation step.
Recall from Sec. 4.3 that the MoE has two further hyperparamcters, x for enforcing sparse gate responses and a noise term on the gate responses during the training, which controlled by s. Instead of learning as parameter during the training like proposed by Shazeer et al. (007 propose to shrink 91/, after
15

---

Published in Transactions on Machine Learning Research (12/2023)
each training epoch, for multiplicative factor 7/ 1 and an initial value == so as hyperparameters. We discuss this heuristic choice in Appendix F
We require appropriate learning rates for the optimization of the parameters and tunable hyperparameters of the model. Generally, we suggest applying an adaptive base learning rate 1) where we shrink <- 7/ := 2 for an initial base learning rate yo during the training as soon as the validation performance gets stuck until Th crosses a lower threshold, e.g., T < 70/0000 Note that, relative the base learning rate, good learning rates for the individual types of tunable parameters should be deployed: Within GP component SVGP(() with 0 == (j,, A, E, X,, 1, the hyperparameters (μ, A, ū) must be updated
smaller scale than the inducing value distribution, given by S,). In this regard, let лn S 1 be the factor such that, if we update 2 at rate 7,, then we update (f,, O) at rate T)и.)
Similarly, we need to update the gate parameters 0gr on a smaller scale than the expert parameters E. In this regard, let ne I be the factor such that, if we update de, at rate n,, then we update 0g. rate 167)
The set of hyperparameters that require off-training selection (e.g., via cross-validation) is thus given by
2 (B, (ar) ag, Ao, X x So, Ma, Oor lo, qH, nG);
(27)
whereas the overall set of parameters that get tuned while training is given by
θ (Mg, λE 0, ME u4 PPG. u().
(28)
We provide further details the design choices for our MoE model in Appendix F
4.4.3
Choosing the hyperparameters
Since our MoE approach is based on known building blocks (Williams & Rasmussen 1996 Hensman et al,
2015
Shazeer ct al.] 2017) we train our model using well-ostablished software libraries Kingma & Ba
2019 Gardner et al. 2018), with the hyperparameters chosen by following best practice. While the set of hyperparameters (27 appears to be large, most of them can be tuned in advance on the initial training dataset of moderate size and held fixed in the subsequent training data refinement process.
2015
Paszke et al.
Note that some hyperparameters impact the computational complexity rather than the model performance. Thus, as long as they are not underestimated, their tuning optional and will therefore be postponed:
Since our MoE is robust concerning unnecessarily large choices of the gate output sparsity K, we initialize x = L while choosing the remaining hyperparameters, followed by tuning k as the last hyperparameter, where we successively reduce λ until we observe a significant loss of performance of the MoE.
The numbers ma = x me = x [x of IPs of the expert and the gate are the main driver of the computational complexity of our MoE. While unnecessarily large numbers will not hurt the model performance, they should therefore be set the smallest value that leads to no significant loss of performance to keep the computational complexity of the model moderate at larger training sizes. Tn the initial iteration, XF, X p are chosen diverse described in Appendix E. In the second iteration, where we have
use mg Ro for the experts and mg "4, where the locations of the IPs
first estimates of LFC, X and X, are drawn as discussed in Sec. 4.5 Here, we gradually decrease mæ and my until we observe significant loss of validation performance of the MoE. We hold me and my fixed in subsequent iterations. Note that the IP locations are not subject to optimization.
The initial base learning rate and the expert's internal hyperparameters learning rate (no- nH) and the batch size B that are related to single GPR expert rather than the whole MoE model:
First, we hold ли = 0.2, = no fixed at reasonable initial values and perform line search over 7o0 according to the resulting validation performance of single, isotropic, sparse GPR export. Here, too small values of no result in slow convergence of the objective, in which case we interrupt the training immediately and increase fo long the first objective updates are consistently decreasing.
16

---

Published in Transactions on Machine Learning Research (12/2023)
Next, we choose Huи according to the resulting validation performance of a single, isotropic, sparse GPR expert, where we gradually decrease 1)8, starting from т/и = 1. Again, we interrupt the training for too large choices of I]H, where the training objective will diverge.
Finally, we choose B according to the resulting validation performance of single, isotropic, sparse GPR expert, where we gradually decrease B, starting from B = nl-
Next, we have to deal with the remaining hyperparameters that are related to the MoF model. For this, we first initialize the less crucial hyperparameters reasonable values, tuning them afterward: We apply a set of experts with {ai where L = 7, and or = 2 which are logarithmically spaced around the global bandwidth estimate Σg of the best-performing model that we obtained from the above tuning of the hyperparameters related to a single GPR expert. The noise added to the gate (5p, as well the regularization 9r of the bandwidth function are about fine-tuning of the model. We set them to 50 0.1, 1/v2 and 0,- 0.01. We suggest to keep 7 1/v2 throughout without further tuning.
Via grid scarch, we choose aa, ÀG according to the resulting validation performance of the MoE model, where we gradually decrease the gate learning rate from BG 1. 1.
As a last step, we choose the hyperparameters for fine-tuning of the MoE model:
First, we perform line search over 9, according to the resulting validation performance of the MoE.
Next, we tune {ai)/ We observe that unreasonable bandwidths will be automatically dropped during the training. Therefore, if the minimal or maximal candidate associated with 01, oz is not chosen during the training, we remove the respective expert and retrain the MoE. Vice versa, we expand the bandwidth candidate range beyond 01, TL with factor of 2F as long as the boundary candidates are not dropped during the training.
Finally, we perform line search over $o according to the resulting validation performance of the MoE.
4.5 Initializing the IP locations
Since the numbers of IPs of the gate X as well as the experts X are the computational bottleneck of our model, they should be chosen advisedly. We can interpret the choice of IPs as a nested AL task at small sample size. In the small sample size regime, input space geometric arguments have proven to be robust and superior in comparison to naive approaches like random sub-sampling from the training inputs (Teytaud et al.) 2007: Yu & Kim 2010 Wu 2019 Liu et al. 2021). They are representative, respecting the training distribution, and diverse (with low-disperes dispersiopr that they achieve an acceptable representation of the dataset at the smallest. possible number of IPs. By low-dispersion, we resort to the following definition:
The dispersion, given by supzex minigign || - 24 [Nederveiter 1988) a measure of how well spread out the training sample is. We say that sequence has low-dispersion if its dispersion is lower than the dispersion of random uniform sampling.
Definition 5.
Indeed, by sampling the IPs in this manner, we reduce the distance of evaluation point æ to its closest neighbor in X.F. which is known to roduce the reconstruction error of full kernel matrix by sparse representation (Zhang et al. 2008).
In addition, recall that our derived LFC measure of local structural complexity quantifies the local variation of the target function. Intuitively, we require more IPs to sense and reconstruct the target function where this local variation is higher. In summaary, we therefore propose to choose the IPs
x X ~ vp- &GP
(29)
as the geometric mean of LFC and the training density X,, ~ p in diverse way. Here, we ensure diversity by implementing distribution preserving clustering particle repulsion as described in Appendix P
17

---

Published in Transactions on Machine Learning Research (12/2023)
5 Experiments
In this section, we will first analyze our approach on toy-data, regarding the MoE model, LFC, and the superior training densily. Then, we apply our approach to a high-dimensional MD simulation dataset from quantum cherisstr, by which we can deduce deeper insights into this regression problem.
We denote the root mean squared error (RMSE) and the marimum absolute envor (max AE) of a model f for Lest set X E X with XT q by
RMSE(S, X x,
1 N
Σ xEXT
- -
[(x)
(x) and max AB(J. X, Y,) max zEX,
[(x) - Fx,,,, (x)
As already discussed in Sec. 4.1 the learning rate is invariant under change of the training density X,, p in the considered scenario. For our MoE model and f E Ca (K, we, thus, can write
RMSE MOES X,, Y,,) Cpn (1 + (1+)(1)) with T=
2a+d'
a <0000 00
(30)
1122
==
where Cp > 0 is a constant depending the training density p. Note that we can theoretically bound the asympiotic RMSF from below by C"n-*, where we have defined C* == Cp- with X p' being the optimal training set from 33. Unfortumately, since p" is unknown even when given the ground truth we are not able to estimate C and, thus, provide a lower bound the RMSE beyond the known learning rate n n=T.
As an AL performance measure, we use the relative required sample size from Definition 2 which can be estimated for based model such our MoE and f f Ca (X, (K, according
f MORE p) ~
RM X X..Y) RMSE MOR X X
(31)
where it is X ~ and X, ~ with respoctive labels Y and Y, Using 31 we can compare the asymptotic AL performance different AL sampling schemes in the following experiments. For example, we can quantify the AL performance of our proposed AL framework by sampling X', -Pup
~GPR,n
5.1
Doppler function
We will first demonstrate our approach on the Doppler function (see, for example, Donoho & Johnstone 1994)), which was also discussed in related work that deals with inhomogeneous complexity (Panknin et al, 2021 Bull al. 2013). For x E X' 1], let
P(yl=) = N(y: ((), 1), f(x)= C sin (2m(1+ e))
where 0.05, C' is chosen such that I1/II, 7 and Ni: o²) denotes the Gaussian distribution with mean a and variance o². We assume a uniform test distribution U(X) in all Doppler function experiments.
This one-dimensional, homoscedastic toy-example allows for an easy and intuitive visualization. Fig. D shows an example dataset as blue dots and the true function f to infer black. Due the strong variation of structural complexity, a single-scale GPR model does not cope well with the Doppler function (scc Appendix G.1 for comparison of single-scale to multi-scale GPR).
We implement our proposed MoE model described in Sec. 4.3 with sparse GPs as the expert and gate models and using the Gaussian kernel k. We apply 512, respectively 128 IPs for the experts and the gate, which are chosen via SVGD (see Appendix E Furthermore we apply oj = 100-10)/3 , S j < 7, as the expert bandwidths, Äg = 20 as the initial expert regularization, and og = 0.05 and ÅG = 10 for the gate. For the training, we apply a batch size of B 512, a terminal export sparsity к = 2, a penalty factor of 0, = 0.5 for small bandwidth choices, gate noise parameters su = 0.1 and 7, 1/v2, and learning rate parameters y = 0.01, ?/ 0.2, n0 1.
18

---

Published in Transactions on Machine Learning Research (12/2023)
10
10.
random IPs
-10-
-10-
optimized IPs
10a
X
logescale
Figure 3: The Doppler experiment: An exemplary dataset and the locations of 128 IPs, once sampled most naively-that is, random according to the test distribution-ani and once optimized regarding diversity as well as structural complexity and representativeness as described in Sec. 4.5 shown on natural x-scale (left) and on logarithmic x-scale (right).
Learned gate function
Predicted local properties
1.0
10a
HGP
0.8
6
0.6
Puu
10-1_
0.4
0.2
10-2
6.0
10
10
10
10 10a
10
10 2
10
10
log scale x
Icale
Figure 4: The Doppler experiment: The gate function (left) and the associated estimates of LOB, LFC and
superior braining densily (right) trained the dataset from Fig. g and shown logarithmic x-scale.
Fig. 4 shows the gate function after training of the MoE model, as described in Sec. 4.4. and the associated estimates of LOB, LFC and the superior training density, calculated according to 24 and (21)
Comparing the active learning framework in the LPS and GPR domain Since f C (X,R), our deduced superior training density estimate given by Eq. (21. In Fig. and the superior training densily in comparison to the LPS-based results for polynomial degrees of order
plot our estimates of LOB
Q = 1,3, and with implementation and hyperparameters as described in Panknin et al. (021). Here, we can observe the qualitative similarity of the LPS- and GPR-based estimates of LOB.
When conducting the proposed GPR-based active sampling scheme as described Sec. 4.2 we additionally observe quantitative benefits Fig. f over random lest sampling quite similar the L.PS-based result. for Q 3: When estimating the relative sample size (31) we require achieve the same RMSE via active sampling compared to random test sampling, we obtain e(Fr.rr GGPR.n = 0.58 + 0.04. This means that we save about 42% of samples via our active sampling scheme.
This provides evidence for the effectiveness of our superior sampling scheme, combining the theoretical foundation of the LPS domain with the efficient access to LOB estimates in the GPR domain.
Comparing random test sampling to equidistant sampling In the introduction, we indicated that the advantage of the robust and model-agnostic input space geometric arguments (Teytaud et al.) 2007 Yu & Kim 2010 Wu, 2019 Liu et al. 2021) diminishes as the training size grows. We can substantiate this claim by comparing random test sampling to equidistant sampling on the Doppler dataset. By equidistant sampling over X' [0, 1] we mean the deterministic construction of the training inputs, where for n given training samples the subsequent training inputs get placed halfway between all nearest neighbors of the former n.
19

---

Published in Transactions on Machine Learning Research (12/2023)
Exemplary Bandwidth Estimates
Training Density
10
10"2
GPR..
10a
10
6
OPR
10
0.0
0.2
0.4
0.6
1.0
0.4
0.6
0.8
1.0
x
Figure 5: The Doppler experiment: The LOB estimates (left) and the resulting superior training density of our proposed GPR-based approach in comparison to the LPS-based approach of order Q 1,3 (right). The results are averaged over 20 repetitions.
RMSE
max
MOE, MoE. 6GPk.n
Rips.
1 Popt
C n
10f
10-1
512
2048
8192
32768
2048
8192
32768
Training Size n
Training Size
Figure 6: The Doppler experiment: The RMSE (left) and the max AE (right) of our proposed GPR-based approach in comparison to the LPS-based approach of order Q == 3 (sce Eq. 55) and (6) once using the respective AL scheme and once, applying random test sampling. The results are averaged over 20 repetitions. The long-dashed, gray line is for illustration of the optimal learning rate т from (30), where the offset C* is imaginary. It shall therefore not be confused with a true lower bound.
samples. In this way, the training input inter-distances are halved exactly with each iteration. We regard this construction as the optimal input space geometric choice, which result will subsume all AL competitors of this type. We now observe in Fig. 7 that, indced, equidistant sampling superior to random test sampling small training sizes. As claimed, however, with growing training size, this advantage gradually diminishes until has vanished completely at 215 training samples.
On Gaussian process uncertainty In Sec. 2 we mentioned that GP uncertainty sampling is the most related approach to our superior sampling scheme since both build on GPR models. As also discussed therein, we can regard standard GP uncertainty sampling as an input space geometric argument, whose performance we can subsume by equidistant sampling in the Doppler experiment. Particularly this implies that standard GP uncertainty sampling provides no benefits regarding asymptotic AL performance.
Instead given the gate function of our MoE from the previous part of this experiment, which obtained for 215 training samples and which we now keep fixed we define the uncertainty estimate of our model MoGPU straight.foward extension of GP uncertainty sampling which takes the inhomogeneous complexity of data into account: By simply weighting the predictive variances of all experts in some input æ with respect to the gate values G(z) from (5), we derive
MoGPU(x) Σ'_ ((z))O)
(32)
20

---

Published in Transactions on Machine Learning Research (12/2023)
RMSE of MoE
max AE of MoE
10
X 9-
X ~ q
Equidistant Xn
Equidistant Xn
10-1
512
2048
8192
32768
512
2048
8192
32768
Training Size n Figure 7: The Doppler experiment: The RMSE (left) and the max AE (right) of our proposed MoE model when comparing random test sampling to equidistant sampling. The results averaged over 20 repetitions.
Training Size n
Training Data Distributions
10
1/ MoGPU @GPR.n
100
X PSup Equidistant
10-3
10-2
10-1
100
log-scale
RMSE
10
MoE, Sup MoE, MoGPU MoE, Equidistant
@ MoE, Pup
MoE, MoGPU
MoE, Equidistant
C n"n
100
10-1
128
512
128
2048
8192
Training Size n
Training Size
Figure 8: A comparison of the mixture of Gaussian process uncertainty and the equidistant sampling baselines to our proposed active sampling scheme for the Doppler experiment: (Top) The training data hist ogramss after 23 samples, contrasted with functions of oGpR, and the RMSE (bottom left) and the max AE (bottom right) at several training sizes of the compared schemes. The results are averaged over 20 repetitions. The longedashed, gray line is for illustration of the optimal learning rate т from (30) where the offset C* is imaginary. It shall therefore not be confused with a true lower bound.
where Cū is the predictive variance of the 1-th expert (sce competitor to our superior sampling scheme.
12
Note that we consider MoGPU as bascline
Intuitively, the uncertainty estimate in I X increases as the applied bandwidth (x) decreases. Now, in order equalize uncortainty over the input space, MoGPU will sample more in regions where a is smaller. For X, drawn according MoGPU, we expect X,, [oGra]-d This expectation holds as can be at the top in Fig. g
21

---

Published in Transactions on Machine Learning Research (12/2023)
Training Data Distributions
10
X ~ ALC
P n
109
X q
10
10-2
10
100
log-scale
RMSE
max AE
DGP, PSup
DGP, Sup
10a
DGP, ALC
DGP, ALC
10a
DGP, q
DGP, q
100
128
512
2048
Training Size
Training Size
Figure 9: A comparison of the DGP model performance, using the ALC criterion for DGP, the random test sampling baseline, and our proposed active sampling scheme for the Doppler experiment: (Top) The training data histograms after 512 samples, and the RMSE (bottom left) and the max AE (bottom right) at several training sizes of the compared schemes. The results are averaged over 5 repetitions.
For evaluation, we combine the fixed gate function with full GPR experts and compare our proposed sampling scheme with MoGPU (both determined through the gate). all error mcasures the beneficial effect of the low-dispersion property of X, drawn according to MoGPU has already vanished for about 512 training samples, from where the asymptotic law dominates. As expected, comparing RMSE. Interostingly, MoGPU is superior to our approach regarding the max AE, suggesting that X n [orrn-- is the preferable training distribution under the supremum- norm..
approach is superior to MoGPU when
AL Performance on deep Gaussian processes To demonstrate the model-agnosticitty of our AL approach, we deploy the DGP model Sauer et al. (2023b) using the CRAN package deepyr This package also implements an aggregate variance- brasedce-based AL criterion (Cohn 1994), which they named active learning Cohn (ALC) after the originator. We deploy 3-layer DGP model, using the Gaussian kernel. For test evaluation, We train the model using Verchia-approziimatiin (Sauer et al. 2023a with a total of 10,000 Gibbs-sampling steps, burning the initial 8000, and thinning the remaining steps to 1,000.
Beginning with 128 cquidistant samples, we refine the training data of the DGP model using the ALC criterion, random test sampling, and our proposed superior training scheme. The resulting training data distributions and the performance of the DGP model are plotted in Fig. g As expected, our superior training scheme performs superior to random test sampling. While the ALC criterion that is particularly designed for the DGP model performs best, observe only very little difference at 512 training samples. Note that sampling according to the ALC criterion becomes computationally challenging already at this point since the DGP model has to be re-trained after each new sample. In contrast, sampling X,, ~ Psup can be performed in batch mode. This result emphasizes the complementary nature of our asymptotic work to the classic bottom-up AL literature.
GPR.a
22

---

Published in Transactions on Machine Learning Research (12/2023)
RMSE at 2 training samples
IP
10a
V vp 6
GPR
10-1
p q
p=
GPA.n Sup
64
128
256
512
1024
Number of expert inducing points m
Figure 10: The Doppler experiment: The curves show the RMSE at training size n = 2 for a varying number of expert IPs m. The colors correspond to different IP distributions, whereas the line styles correspond to the underlying training distribution. The results are averaged over 20 repetitions.
Necessity of the small bandwidth penalty We impose a penalty on small bandwidth choices through
the factor 0g = 0.5 to regularize the bandwidth function and prevent overfitting, as described in Appendix We demonstrate this overfitting issuc in Appendix G.2 that results from applying no regularization (0, = 0).
modeling using LFC In Sec. 4.1 we mentioned that LFC can also be used to coarsen or refine the model resolution adequately reduce the overall complexity of the model. While we fixed the IPs to reasonable numbers in the other parts of the Doppler experiment, that is, m 512 and under active, respectively random test sampling, we here investigate the inffuence of the number of IPs and their distribution the capability to resemble the Doppler function. Recall from Sec. 4.5 that we interpret the choice of the IPs as a nested AL task at small sample size (m < n),, where it is reasonable for them to be sampled in a diverse way, respecting the training distribution but also the structural complexity of the target function. In Fig. 3 we show a naive choice and our optimized choice of IPs.
Parsimonious
1024 IPs
In Fig. 10 we compare the RMSE for the fixed training size n 2 for both settings, active and passive, when sampling the IPs according the training density p, the LFC and their geometric mean (29 First of all, we observe that generally require less IPs with active sampling compared to random test sampling, which originates from the fact that the superior training density Paup cGPR,, alrcady respects LFC to some degree. Next, we observe that the geometric mean of the training density and LFC performs best, provided that the number of IPs m is large enough. Finally, we observe that, nou-surprisingly, we can shrink m the most under the LFC distribution, namely to m = 128, before the performance of the model degrades substantially.
In summary, we are able to shrink the model complexity up to factor of 8 for the Doppler function without. significant loss of performance, when respecting LFC in the model design.
Comparing our proposed IP selection method to greedy fast forward selection In Sec. 2 discussed other IP selection approaches. For comparison, we have implemented the greedy fast forward (GFF) IP selection method of [Seger et al. (2003), in which, beginning from scratch, the most informative training inputs are gradually added to the set of IPs as a means to approximate the full GP(O) distribution. Here, the information of an IP candidate xa X, X, is measures by
J KL [(x
which is the Kullback-Leible divergence between the posterior distributions based on the IPs X, U fxz) and X,. Accordingly, the updated set of IPs is given by X. + X, U { where
= arg max, Ex, (xx)).
The procedure converges, when the remaining IP candidates carry no further information, that is, J(xt) < Eg, up to a specified threshold E 2
At a given threshold Ey, we observe that the number of selected IPs is very small for the experts with large bandwidth, while it increases drastically (xx 1) for experts with small bandwidth. Now that the overall
23

---

Published in Transactions on Machine Learning Research (12/2023)
RMSE 2 training samples. q
IP P
-Vp 6a
IP GFF
10-1.
Number of expert inducing points
Figure 11: The Doppler experiment: The curves show the RMSE at training size 2 at varying number of expert IPs m, where the IPs are either chosen at random, according to our proposed selection method, and via GFF. The results are averaged over 20 repetitions.
Heteroscedastic Doppler problem
LOB
0.10
0.08
0.06
0.04
f+ 2VV
Caaal IV
-10
/v
0.2
0.4
0.6
0.8
1.0
0.2
0.1
0.6
0.8
1.0
Training densities
X
RMSE
GAR. V
MoE/ SGAR. 1oa
Pup
4-
GGRR. V|V =
MoElV P
MoElv
3
2-
10-1
512
2048
8192
32768
Training Size
Figure 12: The Doppler experiment under heteroscedaastic noise: (Top left.) An exemplary dataset; (Top right) The LOB estimates, when comparing the homoscedastic to the heteroscedastic Doppler experiment; (Bottom left) The training densities of random test sampling, Paup homoscedasticity; (Bottom right) The RMSE at several training sizes of the compared sampling schemes. The results are averaged over 20 repetitions.
GPR,a and Psup
GPR,n when wrongly assuming
complexity of the MoE dominated by the expert with the most IPs, is fair to compare the mumber of IPs of the expert at bandwidth a with our statically specified number m of IPs in Fig. I0 Here, we will vary the threshold εg to obtain curve that maps the associated number of IPs to the achieved RMSE. The results in Fig. 11 show that with GFF almost no IPs can be saved for this inhomogeneously complex problem, opposed to our proposed IP selection method.
GPR.R the selected IPs by Seeger et al. (0033) are uniformly In any case, even under training according to Psup distributed. In particular, the need for less IPs of the experts at smaller bandwidths to the right of A' not recognized and, thus, we observe no IP savings at an acceptable performance over a random IP selection all for this inhomogeneously complex problem.
24

---

Published in Transactions on Machine Learning Research (12/2023)
Heteroscedastic noise treatment While the treatment of heteroscodastic noise is not the main focus of this work, we will now demonstrate our approach on a heteroscedastic version of the Doppler experiment. For this, we let v(z) = (3 - | - 0.55)) E [1.9], which we plot in Fig. 12 (top left) together with the resulting dataset. Here, we assume the local noise variance (or an estimate of it) to be provided externally, again, since its estimation is out of the scope of this work. However; note that the estimation of v is well-studied in the literature, especially for GPs (Kersting et al., 2007 Cawley et al. 2006 Tresp, 2001
In Fig. 12 (top right) we compare the LOB estimates obtained from the homoscedastic dataset and the heteroscedastic version. As we have suggested in Sec. 1.1 the influence of v the LOB estimates obtained from heteroscedastic GPR experts is relatively small. Likewise, we proceed with the evaluation of our proposed AL scheme under heteroscedasticityy Here, we compare to random lest sampling but also to Psup the wrong assumption of homoscedastic noise. Note that we use the heteroscedastic MoE in all cases since the wrong assumption of homoscedastic noise in the experts makes the MoE. very volatile. The respective training densities and RMSE learning curves can be seen in the bottom row in Fig. 12 Due to the stronger noise (compared to the homoscedastic experiment); the asymptotic behavior begins materialize later from = 2 training samples. Until this point, sampling only with respect to the structural complexity looks also promising. However, as soon as the training size becomes large enough to roughly resemble the target function, respecting the inhomogeneity in the noise level becomes crucial to achieve a homogeneous pointwise convergence and, thus, maintaining asymptotic superiority.
CPR,a under
5.2 Force field reconstruction
We now turn our attention to a real-world example in which we predict the potential eneryy surface (PES) and corresponding force field (FF) of molecule from first-principles calculations. The PES function links the geometry z = [R...... Ra] E R3xa of a molecule to its potential energy E E R. where R₂ are the Cartesian positions of the atoms of the molecule. In ab initio computations, this mapping is achieved by solving the indeendneee Schrödinger equation. The PES encodes essential information on the properties of a molecule. Due to thermal and quantum effects, molecules are never perfectly rigid but assume different configurat ions. The distribution of these configurations is determined by the shape of the PES. For example, the minima of the PES will be sampled more frequently than other regions and correspond to stable structures. This has practical implications since many experimental techniques measure an expectation value over molecular distributions. In order to achieve a meaningful comparison, sampling needs to be taken into account in theoretical simulations as well. One of the most successful approaches to sample molecular distributions is MD simulation. They model the time evolution of the atomic positions, sampling the PES by integrating Newton's 's equations of motion. To this end, energy- energy-conserving forces acting on each alom are required. These forces are the negative derivative of the PES with respect to the atomic positions F E 3xa
This type of proxy for the prohibitively expensive ab initio quantum mechanical calculations is commonly used to enable long- long-timescale MD simulations that consist of millions of steps, each requiring the evaluation of the PES and FF for a new geometry. Converged MD trajectories give unique insights into the dynamic behavior and structuro-function relationships of physical systems at atomic scale. They are widely used in molecular biology research and play crucial role in applications such as protein folding and drug discovery. ML has the potential to profoundly advance this field, as it bears the promise of offering a unique cost-accuracy trade-off that is not achievable with traditional methods (Noé et al.) 2020 von Lilienfeld et al., 2020 Unke et al. 2021b: Keith et al. 2021). However, some commonly deployed ML- based FFs rely on rather naive exhaustive sampling schemes to gather training data, which stands the way of scaling to larger system sizes, both, from data acquisition cost and training perspective. Here, we demonstrate how our mothod can be used to construct smaller, yet more effective training datasets.
In this experiment, we reconstruct FF for the molecule malonaldehyde, which has 9 atoms and the chemical formula C3 O2 (see Fig. 13 (A)). Formally, we try to infer the high-dimensional target function ʃ: + y/ R. + F], where X === 3a and y = R1+3a For visualization purposes, we only show a two-dimensional subspace of the PES, which is characterized by the two main features of this molecule, its two rotors (aldchyde groups) (Chmicla et al.) 2018 [auuoda et al. 2020). Their relative orientation is the dominant driver of the potential energy in this case and therefore most descriptive. Each point on the surface depicted in Fig. 13 (B) is generated by fixing the rotor pair at a particular angle and relaxing all remaining
25

---

Published in Transactions on Machine Learning Research (12/2023)
A
J
egg. 500K
B
D
E
LFC Estimate
Test density
Superior training density
Geometry sampling
Energies
x
Forces
Figure 13: Reconstructing ML-based FFs using our MoF approach: (A-B) The inputs and outputs of the regression task are the geometries and energies (including forces, i.e., energy gradients) of malonaldehyde. As an example, we highlight the geometries of the two energetically stable states found in the local minima of the energy surface. (C) The density estimate of the true MD geometry distribution.
(E) The superior training density estimate (21) based on our approach. All properties are evaluated at the relaxed malonaldchyde configurations and plotted with respect to the angles of the two aldehyde rotors of malonaldehydde (see Chmiela et al. (2018) Sauceda al. (2020)).
degrees of freedom to obtain a minimal energy configuration. We will refer to these geometries the relazed configurations in the following.
To reconstruct the FF, we consider the broadly adopted symmelric grudient-domain machine learning (sGDML.) method (Chmiela et al., 2018 2019), which is a GPR model that takes energy and force labels and also roto-translational and permutational invariances of the geometries into account (see Appendix H.1 for details). We anticipate that sGDML will benefit from our MoE approach, where we deploy sGDML as the expert model, since the transition paths along the PES vary in complexity, due the interplay between distinct atom types with different characteristic interaction length scales. Our AL approach can only improve training efficiency if there are inhomogeneities in the data. Using our LFC estimate, we therefore first verify our intuition that the PES of malonaldehyde varies in complexity. Based this, we derive the superior training density, which we finally input into our AL framework to refine the training dataset in superior way.
Experimental setup All experiments use an extensive pre-computed reference trajectory (almost a million data points (X poot poot)) as ground truth, as opposed to generating new data points on demand. This test setup allows post-hoc verification of the training distribution generated by our AL approach, while still providing ample redundancy and therefore sampling freedom.
Recall from Sec. 4.2 that we require an unnormalized density estimate of the trajectory X pool Px since we are dealing with a base- AL scenario. We estimate Px by standard kernel density estémation,
based on the enorgy-to-energy entry of the sGDML kernel 1 from (46) at 0.03. Fig. 13 (C) shows the density estimate of the relaxed configurations, where observe that Px is very unbalanced, with strong concentration of mass near the stable configurations.
We implement our MoE approach, using the sGDML kernel ie from 46 with a Gaussian base kernel function k. While we sample the training data randomly (with appropriate weights) from the pool, we will draw sub-samples (i.e., for choosing the IPs of sparse expert and gate models) via symmetrized distributional clustering (DC) with distributional k-means++ initialization (see Appendix E..
Since this dataset comes with practically noise-free labels (we consider the first principle calculations ground truth), we tune the experts (and MoE model) with respect to MSE rather than the Obj objective. For stability, we will apply 6(x)- = 10-9 even though we assume no noise.
Anisotropic bandwidths sGDML operates on = a(a- - 12/2 = 36 features that are based on the interatomic distances of the molecule. In contrast to the work of Chmicla al. who restrict themselves to an isotropic bandwidth Σg = TgIa, our implementation of sGDML in GPyTorch naturally enables us to tune anisotropic bandwidth Σg diag(01,.....u) in the preprocessing step.
26

---

Published in Transactions on Machine Learning Research (12/2023)
H
Figure 14: A visualization of the individual feature importance of malonaldehyde the anisotropic case: The structural formula of the molecule is plotted in black. The im- portance of the individual interatomic dis- tances is reciprocal to Σx, which is the bandwidth estimate obtained by training the anisotropic sGDML model. Hence, we ex- press the importance of cach interatomic dis- lance of the molecule in red, where the im- portance corresponds the line saturation
H
H
H
(in-E)))mi((-oo(2-A.. 12 E [0. max log( -min(- log())g))
We partially offset the increased memory footprint of the model due the tunable Σg by implementing the sparse GPR model from Sec. 3.4.1 under the sGDML kernel k from (46) and limiting the number of IPs to 128 configurations. Since all our features are inherently calibrated in terms of scale. Hence, the reciprocal entries of Σ5 directly translate into the importance of the features, which we display in Fig. 14
of the same type-pairwise interatomic distances-ther
We observe, that the importance assigned to some pains of atoms agrees with chemical intuition, e.g., interactions with light hydrogen atoms are generally weaker. Furthermore, the important role of the opposing aldehyde groups in malonaldehyde emerges in the form of a heavily weighted path that connects the O-C-C-C-O backbone of the molecule.
In Fig. 15 see that our anisotropic variant of sGDML, performs consistently better than the original isotropic sGDML model. Similar to the calculation of the relative sample size in (31) we can compare two models of equal asymptotic MSE law. When comparing anisotropic to isotropic sGDML, both under random test sampling, we can save about 10% of samples.
Setting up the MoE model After having trained Eg, we apply dense sGDML experts with o where oj = 2 3/
j < as the individual expert bandwidths, Äf 1 as the initial expert regularization, and og 0.1 and λg 10 for the sparse gate with 1024 IPs. For the training, batch size of B = 1024, terminal expert sparsity x = 8, a penalty factor of 0, = 0.01 for small bandwidth choices, gate noise parameters su = 0.01 and = 1/v2, and learning rate paramcters 17 0.005, 7a = 0.05, JiG = 0.1. As we discuss in Appendix F for tuning the MoF. with dense (sGDMT.) experts, we either require an additional gate training set, which independent of the training set for the experts, or we could provide leavo-ono-out (LOO) responses of the exports for the training of the gate. In our experiment, we use an additional gate training set XG of fixed size ni 214 The rnpso MoE model performs consistently better than anisotropic sGDML, as can be seen in Fig. 15 When comparing the anisotropic MoE model to isotropic sGDML, both under random test sampling, we can save about 21% of samples.
apply
Active learning We assume an intrinsic dimension of d = 2 (the two aldehyde rotor angles, the most salient features of malonaldehyde) and a smooth target function f E C (X, The test distribution is given by the MD trajectory such that q = Px- Prior to the AL procedure, we separate the validation samples X and test samples XT at random from the pool X poot- We apply an initial expert training size of no = 22, doubling the sample size with each iteration of the AL procodure. The initial expert training set no are drawn via importance sampling from the remaining pool with weights Bool (X,a U X7)). By this it is X,,, g1/2, which more in alignment with the superior training
No and the gate training sel XG
1x (x
Px
density (2 than sampling X,, X,,, ~ q-
27

---

Published in Transactions on Machine Learning Research (12/2023)
RMSE, Q Pponl
0.6
Isotropic SGDML,
Ppsol
0.5
2x10
Anisotropic SGDML, p Ppoot
0.4
Anisotropic MoË. p p Ppeol Anisotropic MoE, p @EPR.n
0.3
30
1024
4096
Training Size n
Figure 15: The RMSE under the true MD trajectory test distribution for different variants of sGDML and training distribution at varying training size: The performance is given for passive sampling, using the original isotropic sGDML (dotted), anisotropic sGDML (dash-dotted) and our MoE model with anisotropic sGDML experts (dashed), and for the MoE model, applying the proposed superior sampling scheme (solid). The results are averaged over repetitions.
In Fig. 13 (D, E) show the estimates of LFC and the superior training density under the pool test distribution, evaluated on the relaxed configurations of malonaldehyde. The LFC estimates confirm our expectation that the transition areas are more complex to model than the regions near the stable configurations. Subsequently, our active sampling scheme shifts sample mass away from the stable regimes in favor of the transition areas.
We have plotted the error curves of passive and active sampling schemes in Fig. 15 When estimating the relative sample size (31 that we require to achieve the same RMSE via active sampling compared to random test sampling, we obtain e(S e(Snue:PEup Paup -GPR.a) = 0.920 = 0.013. This means that we save about 8% of samples under the MoE model with our active sampling scheme compared random test sampling. In total, when comparing our actively trained MoE approach to the passively trained, original sGDML model, we can save about 31% of samples. Notably, DFT level calculations (Perdew et al. 1996 Blum al., 2009 Tkatchenko & Scheffler 2009) for the studied system require minutes to hours of computation per sample, CCSD(T) level computations even require days of computation per sample. So in the field of quantum chemistry saving roughly a third of computing power is of icaltical importance.
6 Discussion
Active learning Recall that in this work we have restricted ourselves the scenario of model-agnostic AL with persistent performance large training size. This scenario is complementary to the more common small sample size regime. And while both cases are important, lot of AT related work (as discussed in Sec. 2 does not apply to our AL scenario. We also discussed and demonstrated in our experiments that input space geometric arguments which are model-free asymplotically come with no benefit over random test sompling. Since our proposed model is GPRRbasedd we also analyzed uncertainty sampling (MoGPU) for our MoE model to show that our proposed superior sampling scheme differs from uncertainty sampling. Moreover, our superior sampling scheme showed to be superior to uncertainty sampling. Finally, in the regime of our AL scenario, the approach by Panknin et al. (2021) recently has demonstrated maate-of-the-at performance to recent, sophisticated, model-agnostic AL approaches. This was done by training different models the actively constructed training sets of their and other AL approaches and assessing their performance. In particular, they compared lavorably to Goetz el al. (2018) random tree-based AL approach in heteroscedastic setting, using a regression forest model and to Bull et al.) (2013) a wavelet-based AL approach in a setting of inhomogencous complexity, using an RBF-network. This demonstrates the ficxibility of this AL approach in terms of learning problem specifications as well as model choices. Now that our work builds on the previous work of Panknin et al. (2021). state-of-the-art performance of our work is implied.
28

---

Published in Transactions on Machine Learning Research (12/2023)
Interpretability of LFC Due to model-agnosticity, we consider LFC to be an intrinsic, interpretable property of the regression problem, which can be used as an analysis tool by domain experts:
When looking at a single point in a high-dimensional input space, visual assessment of the local structural complexity (e.g., a human can visually detect more complexity to the left of the Doppler function) is challenging. Here, the scalar LFC value gives a human assessable, quantitative description. The LFC (as a scalar-valued function) then even allows for an easy visualization of the local structural complexity in high-dimensional input spaces, if the input space features a reasonable low-dimensional projection. This benefit was demonstrated in the high-dimensional FF reconstruction experiment, where the two-dimensional visualization of LFC provides new insights into the regression problem. Note that the LFC function cannot. be visualized in the absence of low-dimensional projection.
Parsimonious modcling We proposed a novel, modol-agnostic approach to select the IPs of GPR, sampling them in a diverse way from a distribution that is representative for the training data and respects the LFC. In the experiments, we have seen that for problems of inhomogeneous complexity, our approach sustains the expressive power of the model at a considerably smaller number of IPs, compared to the GFF IP selection method of (Seeger et al. (2003).
Heteroscedasticity While inhomogeneities in noise are not the focus of our work, note that both, our model as well as the original AL framework upon which we built our approach naturally deal with heteroscedasticity It, therefore, suffices to complement our work with an estimate of the noise variance function v as, e.g., given in Kersting et al.] (2007); Cawley et al. (2006). The only aspect left open is to elaborate on the impact of e on the LOB of GPR and, thus, the adequate adjustment with respect to in the derivation of the LFC. As we argued, GPR already treats heteroscedast through local adaptions of the regularization. Hence, we assume the influence of v the LOB to be negligible to not existent. This is opposed to the LOB of LPS whose only way to deal with heteroscedasticity is through adaption of its LOB.
Intrinsic dimension and smoothness of the problem In our derivation of LFC and the superior training density, we assumed the intrinsic dimension 8 I d and the smoothness a E (0, oo] of f E C RR to be given through domain knowledge.
If 8 is unknown, we can estimate it from unlabeled input instances, such X poot in pre-processing step. However, this is beyond the scope of our work and we refer to the approach of Facco et al. (2017) for an estimate of the 8 for unbalanced input distributions in nigh-dimensional input spaces X.
If we have ground truth knowledge about the smoothness of the target function f Ca (X, R), we resort to â = oc as default in practice. This assumption is justified, as long as f happens to be rougher in at most finitely many locations of the imput space. Since, asymptotically, the violation of a co affects only a set of measure zero, the influence on our AL setting that addresses large training sizes is marginal. We consider the restriction to target functions that are at most rough on a finite set of input space locations a weak assumption which, thus, comes with no practical limitations.
One way to deal with an unknown smoothness is deploy the Matérn kernel
21-v (z, = [(v) (1 x' 211//)) To (1x-2/11/o).
where I is the gomma function and k is the modified Bessel function of the second kind or order v. After finding the best fitting v", we obtain by â = [w*] future work as it is beyond the scope of this work.
reasonable estimate to a. We defer this idea to
Dimensional scaling As opposed to nonstationary GP approaches (e.g., the tree-based or local GP by Gramacy & Lee (2008); Gramacy & Apley (2015) that suffer from the curse of dimensionality through input space localization, we segment the input space into fixed number of patches, given by the L experts of our MoE. Thus, if we were to instantiate our MoE with dense GPR models, our approach scales well concerning the input space dimension d. However, in real-world applications, we typically deal with training sizes that are too large for dense modeling. In this regime, sparse GPR represcntations scale poorly in d as their IPs
For an implementation in Python see http::///eiiikiiiiiiiieaaeeeeeee
29

---

Published in Transactions on Machine Learning Research (12/2023)
must be space-filing [Bioois & Wycoff 2022). Likewise, our density-based AL approach encounters decaying power for large d. A low intrinsic dimension (6 < d) of the regression problem is therefore crucial for our work to apply.
LFC and the superior sampling scheme as a ML concept Recall that we consider LFC to be a problem intrinsic property. Here, the problem characterized by features æ E X with labels f(x)) E R, a hypothesis space of locally adaptive models, and the MISE as the loss function (or rather the pointwise MSE from ( due to the localization of LFC). In this sense, LFC is formally property of the combination of model and loss according to the characterization of (Jung (202)), Chapter 2). Similarly, the superior training density is a property of the combination of model and loss with respect to the same hypothesis space and MISE instead of the pointwise MSE as loss function.
On a realistic implementation in ab initio FF reconstruction In our FF reconstruction experiment, we assumed a. large unlabeled reference trajectory poot to be given that already follows the true molccular distribution. This will not be given in practice, since building the input trajeclory already requires the computationally expensive estimation of the respective labels. At this point, the actual task behind the regression problem would already be solved. We outline realistic ab initio FF reconstruction scenario in Appendix H.3.
7
Conclusion
Standard ML tasks implicitly assume a certain homogeneity in the data scales. However, in practice this structural property of the learning problem may not be fulfilled, c.g., in multiscale problems from the sciences such lurbulence (Brunton et al. 2020 or quantum chemistry (Noé et al. 2020 von Lilienfeld et al. 2020
Cnke et al. 2021b Keith et al.) 2021).
In this work, we aimed to identify local inhomogeneities in regression tasks, which can be used to construct better models and training datasels and for domain interpretation. To this end, we combined recent results on model-agnostic LFC estimates and asymptotically optimal sampling, which are founded in the domain of LPS, with estimates of LOB, which are derived in the GPR domain. By this, we benefit from both sides, having a theoretically sound superior sampling scheme on the one hand, and having access to the required estimates from model that naturally can cope with high input space dimensions the other hand. Furthermore, we have shown how respecting LFC in the selection of IPs contributes to parsimonious modcling.
On synthetic data, we showcased and validated our approach, where we analyzed similarities with the LPS-based analog but also compared to the most related GP uncertainty sampling concepts for AL. To show the full potential of our approach, we studied a real-world, high-dimensional force field reconstruction task. Our approach not only gave access to an interpretable visualization of the inhomogeneous structural complexity but also guided the sampling process in way that takes the structural changes into account, enhancing the quality of the training data. Here, we additionally identified the multi-scale structure of the individual atomic interactions, whose treatment also results in substantial performance gain of the broadly adopted method sGDML.
Future work In Sec. 4.1 we conjecture that the LOB of heteroscedastic GPR is invariant or scales at most weakly with respect to the local noise level v(z). This claim should be supported by further theoretical investigation. While we deployed our estimates of LFC and the superior training density, using = co, if the smoothness of the target function f E C (X, ()) is unknown, it is possible to (re-)estimate a, e.g., by tuning the regularity of the Matérn kernel of GPR. model after the acquisition of each new training data batch. While we have compared to baseline IP selection methods, a thorough comparison to more sophisticated approaches remains open. A promising idca is also to combine our LFC estimate with the IP selection approach by Moss et al.] (2023) lo obtain informative and diverse IPs in GPR. Finally, we will focus on the application of our approach to real-world problems from chemistry, physics, and further domains also applying techniques from explainable AT (c.g. [Samek et al. (2021); Letzgus et al. (2022)) In particular, recent advances on sGDML regarding the scalability by Chmiela et al. (2023) will enable the application of our approach to large molecular systems.
30

---

Published in Transactions on Machine Learning Research (12/2023)
Acknowledgments
D. Panknin, S. Chmicla, S. Nakajima, and K.-R. Müller were funded by the German Ministry for Education and Research as BIFOLD Berlin Institute for the Foundations of Learning and Data (ref. BIFOLD23B). D. Panknin was also supported by the BMBF project ALICE I11, Autonomous Learning in Complex Environ- ments (01IS18049B). K-R. Müller was also supported by the BMBF Grants 01GQ1115 and O1GQ0850, under the Grants 011S14013A-E 031L0207A-D; DFG under Grant Math+, EXC 2046/1, Project ID 390685689 and by the Institute of Information & Communications Technology Planning & Evaluation (ITPP) grants funded by the Korea Government (No. 2017-0-00451, Development of BC1 based Brain and Cognitive Computing Technology for Recognizing User's Intentions using Deep Learning) and funded by the Korea Government (No. 2019-0-00079, Artificial Intelligence Graduate School Program, Korea University).
All funding sources were not involved in the process of writing and submitting this work.
References
David Arthur and Sergei Vassilvitskii. K-means++:
The advantages of careful seeding. In Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA '07, pp. 1027-103: Society for Industrial and Applicd Mathematics, 2007.
Francesco Bellocchio, Stefano Ferrari, Vincezoo Piuri, and Nunzio Alberto Borghese. Hicrarchical approach for multiscale support vector regression. IEEE Transactions on Neural Networks and Learning Systems, 23 (9), 2012.
William H Beluch, Tim Genewein, Andreas Nürnberger, and Jan M Köhler. The power of ensembles for active learning in image classification. In Proceedings of the IEEE Conference on Computer Vision and Pottern Recognition, pp. 9368 9378,9377,77 2018.
Nabil Benoudjit, Cédric Archambeau, Amaury Lendasse, John Aldo Lee, Michel Verleysen, et al. Width optimization of the gaussian kernels in radial basis function net In ESANN, volume 2, pp. 425-432, 2002.
Tristan Bereau, Robert A DiStasio Jr, Alexandre Tkatchenko, and Anatole Von Lilienfeld. Non-covalent interactions across organic and biological subsets of chemical space: Physics-based potentials paramctrizod from machine learning. J. Chem. Phys., 148 (48(2))2241706, 2018.
Mickael Binois and Nathan Wycoff. A survey high-dimensional gaussian process modeling with application to bayesian optimization. ACM Transactions on Evolutionary Learning and Optimization, 2(2):1-26, 2022.
Volker Blum, Ralf Gehrke, Felix Hanke, Paula Havu, Ville Havu, Xinguo Ren, Karsten Reuter, and Matthias Scheffler. Ab initio molecular simulations with nuneric atom-centered orbitals. Computer Physics Communications, 180(11):2175 22966, 2009.
Rafael S Bressan, Pedro H Bugatti, and Priscila TM Saito. Breast cancer diagnosis through active learning in content-based image retrieval. Neurocomputing. 357:1-10, 2019.
Steven L. Brunton, Bernd R Noack, and Petros Koumoutsakos, Machine learning for fluid mechanics. Annual review of fluid mechanics, 52:477 508, 2020.
Adam D Bull et al. Spatially-adaptive sensing in nonparametrie regression. The Annals of Statistics, 41(1): 41 41-62, 2013.
Robert Burbidge, Jem J Rowland, and Ross D King. Active learning for regression based on query by committee. In International conference on intelligent data engineering and automated learning, pp. 209 218. Springer, 2007.
Wenbin Cai, Ya Zhang, and Jun Zhou. Maximizing expected model change for active learning in regression. In 2013 IEFE 13th international conference on data mining, pp. 51-60. TEEE, 2013.
31

---

Published in Transactions on Machine Learning Research (12/2023)
Gavin C. Cawley, Nicola L. C. Talbot, and Olivier Chapelle. Estimating predictive variances with kernel ridge regression. In Machine Learning Challenges. Enaluating Predictive Uncerlainly, Visual Object Classification, and Recognising Tectual Entailment, pp. 56 77, Berlin, Heidelberg, 2006. Springer Berlin Heidelberg. ISBN 978-8-540-33428--.
Ziyi Chen, Jiuwen Ma, and Yalong Zhou. precise hard-cut em algorithu for mixtures of gaussian processes. In International Conference on Intelligent Computing, pp. 68-75 75. Springer, 2014.
Stefan Chmiela, Alexandre Tkatchenko, Huziel F.. Sauceda, Tgor Poltavsky, Kristof T. Schütt, and Klaus- Robert Müller. Machine learning of accurate energy-conserving molecular force fields. Science Advances, 3 (5):e600011, 2017.
Stefan Chmiiea, Huziel E. Sauceda, Klaus-Robert Müller, and Alexandre Tkatchenko. Towards exact molecular dynamics simulations with machine-learned force fields. Nature Communications, 9(1):8887, 2018. doi: 10.1038/s1446--01-0066666
Stefan Chmicla, Huziel E. Saucoda, Igor Poltavsky, Klaus-Robert Müller, and Alexandre Tkatchenko. sgdml: Constructing accurate and data efficient molecular force fields using machine learning. Computer Physics Communications, 240:38 45, 2019. doi: 10.1016/..pc:2019.02.000
Stefan Chmicla, Valentin Vassilev-Galindo, Oliver T. Unke, Adil Kabylda, Huziel E. Sauceda, Alexandre Tkatchenko, and Klaus-Robert Müller. Accurate global machine learning force fields for molecules with hundreds of atoms. Science Advances, 9readf0873, 2023. doi: /111////siidv.... adf0873.
William S Cleveland and Susan Devlin. Locally weighted regression: an approach to regression analysis by local fitting. Journal of the Americon statistical association, 83(403):596-66) 610, 1988.
David Cohn. Neural network exploration using optimal experiment design. Advances in Neural Information Processing Systems, 6:679-655 (669, 1994.
Andreas Damianou and Neil D Lawrence. Deep gaussian processes. In Artificial intelligence and stalislics, pp. 207-201 215. PMLR, 2013,
David Donoho and Jain M Johnstone. Ideal spatial adaptation by wavelet shrinkage. biometrika, 81 (1)): 425 455, 1994.
Fouzi Douak, Farid Melgani, and Nabil Benoudjit. Kernel ridge regression with active learning for wind speod prediction. Applied Energy, 103(0):328 340, 2013.
Petros Drineas, Michael W Mahoney, and Nello Cristianini. On the nyström method for approximating gram matrix for improved kernel-based learning. journal of machine learning research, 6(12), 2005.
Elena Facco, Maria d' d'Errico, Alex Rodriguez, and Alessandro Laio. Estimating the intrinsic dimension of datasets by a minimal neighborhood information. Scientific reports, 7(1):12140, 2017.
Jianqing Fan, Theo Gasser, Irène Gijbels, Michael Brockmann, and Joachim Engel. Local polynomial regression: optimal kernels and asymptotic minimax efficiency. Annols of the Institute of Statistical Mathematics, 49(1):79-99, 1997.
Stefano Ferrari, Francesco Bellocchio, Vincenzo Piuri, and N Alberto Borghese. Multi-scale support vector regression. In Neural Networks (IJCNN), The 2010 International Joinl Conference on, pp. 1-7. IEEE, 2010.
Shai Fine and Katya Scheinberg. Efficient svm training using low-rank kernel representations. Journal of Machine Learning Research, 2(Dec):243 264, 2001.
Guojun Gan, Chaoqun Ma, and Jianhong Wu. Data clustering: theory, algorithms, and applications. SIAM, 2020.
32

---

Published in Transactions on Machine Learning Research (12/2023)
Jacob R Gardner, Geoff Pleiss, David Bindel, Kilian Q Weinberger, and Andrew Gordon Wilson. Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration. In Advances in neural informalion processing systems, pp. 7576 7586, 2018.
Jack Goetz, Ambuj Tewari, and Paul Zimmerman. Active learning for uon-parametric regression using purely random trees. In Advances Neural Information Processing Systems, pp. 2537-25-46, 2018.
Mohamed Goudjil, Mouloud Koudil, Mouldi Bedda, and Noureddine Ghoggali. A novel active learning method using svm for text classification. International Journal of Aulomation and Compuling, 15(3): 290 290,299,-00 2018.
Siegfried Graf and Harald Luschgy. Foundations of quantization for probability distributions. Springer, 2007.
Robert Gramacy and Daniel W Apley, Local gaussian process approximation for large computer experiments. Journal of Computational and Graphical Statistics, 24(2)5561-577 2015.
Robert B Gramacy and Herbert K H Lee. Bayesian treed gaussian process models with an application to computer modeling. Journal of the American Statistical Association, 103(483): 111--110,, 2008.
Andrea Grisafi and Michele Ceriotti. Incorporating long-range physics atomic-scale machine learning. J.
Chem. Phys., 151((0)2204105, 2019.
Konstantin Gubaev, Evgeny V Podryabinkin, and Alexander V Shapeev. Machine learning of molecular properties: Locality and active learning. The Journal of chemical physics, 148(24):241727, 2018.
Vincent Guigue, Alain Rakotomamonjy, and Stéphane Canu. Kernel basis pursuit. In European Conference on Machine Learning, pp. 146 157. Springer, 2005.
László Györfi, Michael Kohler, Adam Krzyčak, and Harro Walk. A distribution-fre theory of nonparametric regression, volume 1. Springer, 2002.
Jun Han and Qiang Liu. Stein variational gradient descent without gradient. In International Conference on Machine Learning, pp. 1900 1908. PMLR, 2018.
Juan Mario Haut, Mercedes E Paoletti, Javier Plaza, Jun Li, and Antonio Plaza. Active learning with convolutional neural networks for hyperspectral image classification using a new bayesian approach. IEEE Transactions Geoscience and Remote Sensing, 56(11):6440-6440 64461, 2018.
Xiaofci Hc. Laplacian regularized d-optimal design for active learning and its application to image retrieval. Image Processing, IEEE Transactions on, 19(1):254 263, 2010.
James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational gaussian process classification. In Artificial Intelligence and Statistics, pp. 351-360. PMLR, 2015.
Bing Huang and O Anatole von Lilienfeld. Quantum machine learning using atom-in www.iii-oooleeeeeeeeeee based fragments selected on the fly. Nature Chemistry, 12(10))945-951, 2020.
Robert A Jacobs, Michael I Jordan, Steven J Nowlan, Geoffrey E Hinton, et al. Adaptive mixtures of local experts. Neural computation, 3(1):79 87, 1991.
Michael I Jordan and Robert A Jacobs. Hierarchical mixtures of experts and the algorithm. Neural compudation, 6(2):181-21) 1994.
Alexander Jung, Machine learning: The basics. In Machine Learning: Foundations, Methodologies, and Applications. Springer, 2022.
John A Keith, Valentin Vassilev-Galinen Galindo, Bingqing Cheng, Stefan Chmicla, Michael Gastegger, Klaus-Robert Müller, and Alexandre Tkatchenko. Combining machine learning and computational chemistry for predictive insights into chemical systems. Chem. Rev., 121 (16))999816 6-9872, 2021.
33

---

Published in Transactions on Machine Learning Research (12/2023)
Kristian Kersting, Christian Plagemann, Patrick Pfaff, and Wolfram Burgard. Most likely heteroscedastic
gaussian process regression. In International Conference on Machine Learning, 2007.
Jack Kiefer. Optimum experimental designs. Journal of the Royal Statistical Society. Series B (Methodological), pp. 272-30 272-319, 1959.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA,, USA, May 7-9, 2015, Conference Track Proceedings, 2015. URL
Tsz Wai Ko, Jonas A Finkler, Stefan Goedecker, and Jörg Behler. A fourth-gencration high-dimensional neural network potential with accurate rostaloss including non-local charge transfer. Nat. Comman. 12(1))398, 2021.
Arvind Krishna, Simon Mak, and Roshan Joseph. Distributional clustering: A distribuution-preeerving clustering method. arXie preprint arXiv:1911.0500 1.55999410 2019.
Sanjiv Kumar, Mehryar Mohri, and Ameet Talwalkar. Sampling methods for the nyström method. The Journal of Machine Learning Research, 13(1):981-1066, 2012.
Oleg V Lepski. On a problem of adaptive estimation in gaussian white noise. Theory of Probability & Its Applications, 35(3):454 466, 1991.
Oleg V Lepski and Vladimir G Spokoiny. Optimal pointwisc adaptive methods in nonparamctric estimation.
The Annals of Statistics, pp. 2512-2546, 1997.
Simon Letzgus, Patrick Wagner, Jonas Lederer, Wojciech Samek, Klaus-Robert Müller, and Gregoire Montavon. Toward explainable artificial intelligence for regression models: A methodological perspective. IEEE Signal Processing Magazine, 39(4):40-58 2022.
David D. Lewis and William A. Gale. A sequential algorithm for training text classifiers. In Proceedings of the 17th annual international ACM SIGIR conference on Research and development in information retrieval, pp. 3-12. Springer- Verlag New York, Inc., 1994.
Qiang Liu and Dilin Wang. Stein variational gradient descent: A general purpose bayesian inference algorithm. In D. Lee, M. Sugiyama, U. Luxburg, Guyon, and R. Garnetl (eds.), Advances in Newral Information Processing Systems, volume 29. Curran Associates, Inc., 2016. URL httpgss//prrceeddngg..euu www.inng..e.urs...om http:////.......0.....2222e222..
.pdf
Ziang Liu, Xue Jiang, Hanbin Luo, Woili Fang, Jiajing Liu, and Dongrui Wu. Pool-basod unsupervised active learning for regression using iterative representafivenes-diveeiis maximization (irdm). Pallern Recognition Letters, 142:11-19, 2021.
David J. C. MacKay. Information- basead objective functions for active data selection. Neural computation, 4 (4):590-00) 604, 1992.
Mark Mackenzie and A Kiet Tieu. Asymmetric kernel regression. IEEE transactions neural networks, 15 (2):276-282, 2004.
Elias Masry. Multivariate local polynomial regression for time series: uniform strong consistency and rates. Journal of Time Series Analysis, 17(6):571 599, 1996.
Elias Masry. Multivariate regression estimation: local polynomial fitting for time scries. Nonlinear Analysis: Theory, Methods & Applications, 30(6):3575-3558 3581, 1997.
Edward Meeds and Simon Osindero. An alternative infinite mixture of gaussian process experts. In Advances in Neural Information Processing Systems, pp. 883-890, 2006.
John Moody and Christian J Darken. Fast learning networks of locally-tuned processing units. Neural computation, 1(2):281 1294, 1989.
1
34

---

Published in Transactions on Machine Learning Research (12/2023)
Henry B Moss, Sebastian W Ober, and Victor Picheny. Inducing point allocation for sparse gaussian processes in high-throughput bayesian optimisation. In International Conference on Artificial Intelligence and Statistics, pp. 5213 5230. PMLR, 2023.
Harald Niederreiter. Low-discrepancy and low-dispersion sequences. Journal of namber theory, 30(1):51-33 1988.
Frank Noé, Alexandre Tkatchenko, Klaus-Robert Müller, and Cocilia Clementi. Machine learning for molecular simulation. Annual review of plyysccal chemistry, 71:361 390, 2020.
Danny Panknin, Klaus-Robert Müller, and Shinichi Nakajima. Optimal sampling density for nonparametric regression. Xivi preprint NNw:2105.11990 2021.
E. Pasolli and F. Melgani. Active learning methods for clectrocardiograpiic signal classification. Information
Technology in Biomedicine, IEEE Transactions on, 14(6):1405-14 1416, 11 2010. ISSN 1089-7771. doi:
(LL109//TTB.200 201.229892..
Edoardo Pasolli and Farid Melgani. Gaussian process regression within active learning scheme. In 2011 IEEE International Geoscience and Remote Sensing Symposium, pp. 3574 3577. LEEE, 2011.
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, Benoit Steiner, Lu Fang, Junjic Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. In Advonces in 2019.
neural information processing systems, pp.
Klaus Pawelzik, Jens Kohlmorgen, and Klaus-Robert Müller. Annealed competition of experts for a segmentation and classification of switching dynamics. Newral Computation, 8(2):340 356, 1996.
John P Perdew, Kieron Burke, and Matthias Ernzerhof. Generalized gradient approximation made simple.
Physical review letters, 77(18)33865 3868, 1996.
Geoff Pleiss, Martin Jankowiak, David Eriksson, Anil Damle, and Jacob Gardner. Fast matrix square roots with applications to gaussian processes and bayesian optimization. arXiv preprint arxvv:0006. 1000006.11200 11267, 2020.
Lassi Roininen, Mark Girolami, Sari Lasanen, and Markku Markkanen. Hyperpriors for matérn fields with applications bayesian inversion. Inverse Problems and Imaging, 13(1):1 29, 2019. ISSN 1930- 8337. doi: 10.3934/bpp.2019001 438d 138d-af0a-b712cf45355
URL https://wwa aimsciences. http:////......e////////////////
Nicholas Roy and Andrew McCallum. Toward optimal active learning through sampling estimation of error reduction. In Proceedings of the 18th International Conference on Machine Learning, pp. 441 448. Morgan
Kaufmann Publishers Inc., 2001.
Didier Rullière, Nicolas Durrande, François Bachoc, and Clément Chevalier. Nested kriging predictions for datasets with a large number of observations. Statistics and Computing, 28:849 867, 2018
Priscila. TM Saito, Celso TN Suzuki, Jancarlo F Gomes, Pedro J de Rezende, and Alexandre X Falcao. Robust active learning for the diagnosis of parasites. Pattern Recognition, 48(11):3572-35 111133523333,,,,, 2015.
Wojeiech Samek, Grégoire Montavon, Sebastian Lapuschkin, Christopher J Anders, and Klaus-Robert Müller. Explaining deep neural networks and beyond: A review of methods and applications. Proceedings of the IEEE, 109(3))247-278, 2021.
Huziel E Sauceda, Stefan Chmiela, Igor Poltavsky, Klaus-Robert Müller, and Alexandre Tkatchenko. Con- struction of machine learned force fields with quantum chemical accuracy: Applications and chemical insights. In Machine Learning Meels Quantum Physics, pp. 277-307 3077 Springer, 2020.
Annie Sauer, Andrew Cooper, and Robert Gramacy. approia-approximaaee deep gaussian processes for computer experiments. Journal of Compulational and Graphical Statistics, 32(3):824 837, 2023a.
35

---

Published in Transactions on Machine Learning Research (12/2023)
Annie Sauer, Robert B Gramacy, and David Higdon. Active learning for deep gaussian process surrogates. Technomelrics, 65(1):4-18, 2023b.
Jens Schreiter, Duy Nguyen-Tuong, Mona Eberts, Bastian Bischoff, Heiner Markert, and Mare Toussaint. Safe exploration for active learning with gaussian processes. In Joint European conference on machine learning and knowledge discovery in databases, pp. 183-149. Springer, 2015.
Matthias W Seeger, Christopher KT Williams, and Neil Lawrence. Fast forward selection speed up sparse gaussian process regression. In International Workshop on Artificial Intelligence and Statistics, pp- 261-261. PMLR, 2003.
Ozan Sener and Silvio Savarese. Active learning for convolutional neural networks: A core-set approach. In International Conference on Learning Represenlalions, 2018.
Sambu Seo, Marko Wallat, Thore Graepel, and Klaus Obermayer. Gaussian process regression: Active data selection and test point rejection. In Mustererkenmungg 2000, pp. 27-34. Springer, 2000.
H. Scbastian Soung, Manfred Opper, and Haim Sompolinsky. Query by committec. In Proceedings of the fifth anmual workshop on Compulational learning theory, pp. 287 294. ACM, 1992.
Noam Shazeer, Azalia Mirhoseini, Krzysztof Maziarz, Andy Davis, Quoc Le, Geoffrey Hinton, and Jeff Dean. Outrageously large neural networks: The sparsely-gated mixturo-of-oxeert layer. orXiv preprint WrXiv:100.006338 2017.
Bernard W Silverman. Densily estimation for statistics and dala analysis, volume 26 of Monographs Statistics & Applied Probability. Chapman & Hall/CRC Press, 1986.
AJ. Smola and B. Schölkopf. Sparse greedy matrix approximation for machine learning. In Proceedings of the 1 17th International Conference on Machine Learning, pp. 911-918, San Fransisco, CA, USA, 2000. Morgan Kaufman.
Edward Snelson and Zoubin Ghahramani. Sparse gaussian processes using pseudo-inputs. Advances in neural information processing systems, 18, 2005.
Masashi Sugiyama and Shinichi Nakajima. Pool-based active learning in approximate linear regression. Machine Learning, 75(3):249 274, 2009.
Yu-Hang Tang and Wibe A de Jong. Prediction of atomization energy using graph kernel and active learning. The Journal of chemical physics, 150(4):044107, 2019.
Olivier Teytaud, Sylvain Gelly, Jérémie Mary, et al. Active learning in regression, with an application to stochastic dynamic programming. ICINCO 2007, 2007.
Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In Artificial intelligence and statistics, pp. 567-574. PMLR, 2009.
Alexandre Tkatchenko and Matthias Scheffler. Accurate molecular van der waals interactions from ground-state electron density and free-atom reference data, Physscca review letters, 102(7):073005, 2009.
Simon Tong and Edward Chang. Support vector machine active lcarning for image retrieval. In Proceedings of the ninth ACM international conference on Multimedia, pp. 107-118. ACM, 2001.
Volker Tresp Mixtures of gaussian processes. In Aduances in neural information processing systems, pp- 654 10066600222 2001.
Oliver T Unke, Stefan Chmicla, Michael Gastegger, Kristof T Schütt, Huziel Saucoda, and Klaus-Robert Müller. SpookyNet: Learning force fields with electronic degrees of freedom and nonlocal effects. Nat. Commun., 12:7273, 2021a.
36

---

Published in Transactions on Machine Learning Research (12/2023)
Oliver T. Unke, Stefan Chmicla, Huziel E. Sauceda, Michael Gastegger, Igor Poltavsky, Kristof T. Schüitt, Alexandre Tkatchenko, and Klaus-Roberr Müller. Machine learning force fields. Chem. Rev., (26( 10142 10186, 2021b.
Aad W Van der Vaart, J Harry Van Zanten, et al. Bayesian inference with rescaled gaussian process priors. Electronic Journal of Statistics, 1:433-405 448, 2007.
Aad W Van der Vaart, J Harry Van Zanten, et al. Adaptive bayesian estimation using a gaussian random field with inverse gamma bandwidth. The Annols of Statistics, 37(5B):2655 2675, 2009.
O Anatole von Lilienfeld, Klaus-Robert Müller, and Alexandre Tkatchenko. Exploring chemical compound space with quantum- based machine learning. Nat. Rev. Chem., 4:347-358. 2020.
Matthew P. Wand and M Chris Jones. Kernel Smoothing, volume 60 of Monographs Statistics & Applied Probability. Chapman & Hal//CRC Press, 1994.
Bin Wang and Xiaofeng Wang. Bandwidth Selection for Weighted Kernel Density Estimation. arXiv e-prints, art. arXiv:0709.1616 09 2007.
Manfred K Warmuth, Jun Liao, Gunnar Rätsch, Michael Mathieson, Santosh Putta, and Christian Lemmen. Active learning with support vector machines in the drug discovery process. Journal of chemical information and computer sciences, 43(2):667 673, 2003.
Rebecca Willett, Robert Nowak, and Rui M. Castro. Faster rates in regression via active learning. Advances in Neural Information Processing Systems, 18, 2005.
Christopher KT Williams and Carl Edward Rasmussen. Gaussian processes for regression. Tn Advances in neural information processing systems, pp. 514 520, 1996.
Dongrui Wu. Pool-based sequential active learning for regression. IEEE transactions newral networks and learning systems, (W5):1348-1350) 2019.
Yadong Wu, Zengming Meng, Kai Wen, Chengdong Mi, Jing Zhang, and Hui Zhai. Active learning approach to optimization of experimental control. Chinese Physics Letters, 37110)103301, 2020.
Yan Yang and Jinwen Ma. An elficient approach to parameter learning of the mixture of gaussian processes. In International Symposium on Neural Networks, pp. 165-174. Springer, 2011.
Kun Yao, John E Herr, David W Toth, Ryker Mckintyre, and John Parkhill. The Tol.orMol-0.1 model chemistry: a neural network augmented with long-range physics. Chem. Sci., 98)):2261- 2018.
Hwanjo Yu and Sungehul Kim. Passive sampling for regression. In 2010 IEEE International Conference Data Mining, pp. 1151-1156. IEEE, 2010.
Chao Yuan and Claus Neubauer. Variational mixture of gaussian process experts. In Advances in Newral Information Processing Systems, pp. 1897-1904, 2009.
Xiaowei Yue, Yuchen Wen, Jeffrey H Hunt, and Jianjun Shi. Active learning for gaussian process considering uncertaintics with application to shape control of composite fuselage. IEEE Transactions on Automation Science and Engineering, 18(1):36-46) 46, 2020.
Kai Zhang, Tvor W Tsang, and James T Kwok. Improved nyström low-rank approximation and error analysis. In Proceedings of the 25th international conference on Machine learning, pp. 1232 1239, 2008.
Danian Zheng, Jiaxin Wang, and Yannan Zhao. Non-flat function estimation with multi-scale support vector regression. Neurocomputing, 70(1):420-425 -29, 2006.
37

---

Published in Transactions on Machine Learning Research (12/2023)
A Asymptotic results for local polynomial smoothing
In this section, we will review the theory of Panknin et al. (2021)
The prediction of the LPS model of order Q under the bandwidth Σ S in x E X can be understood as follows: First, the regression problem localized around according to weights (.,x) that decrease with growing distance z. Then we search for the polynomial up order Q that fits the localized regression problem best. Finally, the evaluation of this polynomial in is returned as the prediction. Formally, it is
m²(²) =2 where
(33)
PQ arg min PRP)(R")
Σ (x,,,)) (v p(a, x))²
and PQ(RR) is the space of the real polynomial mappings Rd - R up to order Q.
The localization is controlled by Σ through the kernel weights K xz) for I; E X,; For an RBF-kernel, (2 x²) decays monotonically with growing distance of x' to a. This docay is dampened or amplified as E increases or decreases, respectively (in the sense of the Loewner order).
For readability, since will be replaced by terms with more involved notation, we redefine ! by
MSE (a, (x,.)) = MSE (x, (1..))
(34)
For a bandwidth space S C s Panknin et al. [2021] proposed to minimize the AL objective
MISE (1.
infxes MSE (,, g(x)dxx
(35)
which is the optimal MISE, obtained by predictions that are based on locally optimal chosen bandwidths. If these locally optimal bandwidth choices are well-defined, that is, if for all x E there exists a unique Y' E S
such that
MSE (x, [.Ei,,))))))) infses MSE (s, (.
we are able to define the LOB function
md Ex(x) (x) = argminges MSE Jx [1
This function exists, for example, in the isotropic case S = faIa a > 0} for LPS under mild conditions,
where we denote [xxx)= = o" (x) (see, e.g., Masry 1996 for an overview).
1997
Fan et al. (1977 or Panknin et 2021
Assuming the isotropic bandwidths candidate space S Io a 0}, the LOB in Eq. is asymptotically well-defined function under mild assumptions Denoting the LOB of LPS of order Q by Eg = of (z) such that
(x()) = argmin,,,. MSE (z, mq: oTal X
(36)
asymptotically it holds
(ā(x) Co
o(x) p(zja
200+144
by [x, Iz] 200-14 + +0p अلांग
(37)
where Co is constant, and bo xx, I,l ((x) E m Aa (x) x of LPS Masry 1996 1997). That is, for sequence hn 0 - ao we
is function of taken from the asymptotic conditional bias
write the conditional bias, which is of order Q I, as
f(x) - E
Am Tu (x) m
X, , = h2 bq [z, Ial + op ( [h?+1].
(38)
We require non-vanishing leading bias- and varianoe-terms of m (a), which is guaranteed if Viz E X' it holds that by [z, Id] * 0 from Eq. and n(xx) > 0.
38

---

Published in Transactions on Machine Learning Research (12/2023)
shows how LOB scales asymptotically with respect to the training size the local noise level function e(x) and the training density p(x)). The remaining bias component depends on the local structural complexity, which can be characterized by the derivatives of f in a non-trivial way. Therefore it encodes the local structural complexity of f. Given all other propertics and LOB itself, we are able to formulate LFC in a closed form.
Eq. 37
Definition (Panknin el al. 2021 ). For LPS of order Q, the LFC of f in x E is asymplotically given by
=(x))
v(x)
C
oQ(x)
p(x)m
[(()))
35 given test density q, the AL task is to minimize MISE (q. m x,)) Now, if LOB is well-defined, we can rewrite
As already mentioned in Eq.
- infxes MSE (1,, mq x(X,) q())dz
MISE (g. (f,molX,) x,,)
_x
| MSE (z, (x)x)) (x)) (())dr.
Finally, when solving for the optimal training dataset
X = argminx, ex" MISE (g. m
as in Eq. (3, the optimal training inputs X can be written asymptotically as an independent and identically distributed sample from the optimal training distribution, whose density Popt possesses an asymptotic closed form.
Theorem Panknin al. (2021)). Let u, Co (x, for a compact input space X, where q is a test probability densily. Additionally, assume that and q are bounded away from zero. Lee, v, 2.q 2 e for some e> O. Let k be RBF-kernel with bandwidth parameter space S = {oIs o > >0. Let Q IN be odd and E CQC111 (X)) such that the bias of order Q Q+1 1 does not vanish almost everywhere. Then the optimal training density for LPS of order Q is asymptotically given by
20+11 2(++1) #Q+1+4 u(T) पद पर o(1)).
prr (z) x [()())(())
Q. We will use this optimal distribution to sample X popt with proposed estimator for eg that is scalable with respect to the input space dimension.
For LPS with X, ~pp and X, y, we can asymplotically calculate the relative required sample size from Definition in Sec. 4.1 by
2
MISE
9. X²)
e p) p) =
(39)
MISE
g. X.)
B Analytic GPR formulations
B.1 Classical Gaussian process regression
The GPR model ~~ (P( (see, e.g. Williams & Rasmussen (1996)) defined follows: The GP is described by the hyperparameters
(n, X, E), which are the global constant prior mean μ., the regularization parameter X, the label noise variance function T and the bandwidth matrix Σ of the kernel. If we can assume homoscedastic noise, we let 0(x)= a.
The GP prior then assumes the labels Y, of X, to be distributed according to Y, = ((x,)
(XX,), ((X,)))) for the constant mean function μ(X,)) = ph,s and the covariance function
C(X,) AK, diag (0( (X,)),
39

---

Published in Transactions on Machine Learning Research (12/2023)
where K2 (X,,, is the kernel matrix of x, For test inputs X., the posterior predictive distribution of Y. is then given by
(X.) (x.),
where the predietive mean and covariance are given by
(X.) = (xX.) Can C (,,),,
(40)
('(X,) = C. - -Cun Cn C.T
(41)
and we have defined
C
C
C(X, UX.)
Con
C.
B.2 Analytic sparse Gaussian processes
We define the sparse GPR model 7 SGP(() follows, following Snellon & Ghahramani 20055:
The sparse GP is described by the (hyper-) parameters 0 (n, X, X,), which are the global constant prior mean μ4, the regularization parameter X, the label noise variance function 0 the bandwidth matrix E of the kernel and the prior distribution, given by the IP locations X, E xx"
Here, the degree of sparsity is described by m IPs: This number can be fixed in advance or gradually increased with training size n, where the increase mn = o/n] is typically much slower than n. If we can assume homoscedastic noise, we let ((=)=
The sparse GP then outputs
(XX)) ~ ((x11
for the mean function
1 (X.) - KK (iag(((XX))))) ((x,))
and the covariance function
(X,) = K. - K.t(K, -
) + diag ((..))
where we have defined K₁ (Xt), K, (x((),), Kr: K* X,), Knt K² X,), Q, = K + K af (A + diag(0l (i(X)))) Knt, and A = diag(A) ()) with X = diag Kn + Kn K; kTt)
We choose f to be the constant mean function, i.e., m(X) = μll, for X E X", noting that other mean functions are possible.
C LFC of GPR
Theorem 1 (LFC of GPR). For f CO (X, X, ~ and homoscedastic noise, the GPR-based LFC estimate of f in x E X is asymptotically given by
2a
(x)R(T) = [(())n
PGPR(())
E
(19)
Proof. Let x= _ i-1
n be a segmentation of the input space with non-empty interiors (2.)))...., (2))) which we can define the restricted bandwidth function search space
#0,
[(x) Σx= (x) Σr...., Eg EES..
S
40

---

Published in Transactions on Machine Learning Research (12/2023)
Here, (z) is the indicator function, returning 1 for A and else. Furthermore let D4." E S₁ be the minimizer of the MISE over Sx with [x²((x) E.., 1 xx) such that
/ MSE (2, (0 (x,) q(x)dr minges. / MSE z, | x,, g(x)dx.
Recall from 9 that PPP 20t4 generally holds for arbitrary input spaces. Due to this, asymptotically, Fige does not depend on training samples outside X.- Hence, letting X,, == { E x,, xa E x the
individual " are asymptotically found by solving the isolated segments of the objective
MSE
x fz. " (X g(x)dr = minses
MSE (x, ( ((x)ddr.
x
First of all, it. F. E.X,,, = =((2))), where p(A) = J. p(ardde is the probability for training sample to fall into C X. addition, we need to account for the expanse of x which we measure by Vol(X."*). Here, Vol(A) (A) = S, dr is the volume of A X. Again with 99, it is therefore
Ih.n [p(xptn) Vol(X. )ml 2042
Subsequently, we can slowly refine the segmentation x= -1+ i=1 Xxn,, where maxigisk. Vol(x. + 0 for
kn - xo slow enough (with k o(n)). Then, for almost every I E X, there exists sequence (ik,)keN with x E i4, for all & E IN such that
" =p(x)n(1 + [1]).
By construction, it is kn.m Zin. (x) + op [1]. It follows Eg (x) p(z)n(1 + Op [1] such that
Eern = [p(x)m] =o (1 o [1)) Therefore, asymptotically Cürn (x) = not depend on and p. Under homoscedasticity, asymptotically, CăPR is necessarily a function that only depends on f, which justifies its measure of I.F..
Fa
[ (x) does
p(zjn
41

---

Published in Transactions on Machine Learning Research (12/2023)
D Algorithmic summary of the proposed AL framework
A Algorithm 2: (0 Ex) < hyper init (X, no? no Po, vaal Yval)
Input
1: Initial training data (X Xng))
2: Training data density pu
3: labeled validation set (xval, Yval
Output
4: Initial hyperparameters 0µ K, (0 Go AG. >G, X 50, Me Uo
5: Global (anisotropic) expert bandwidth Eg e
Procedure
6: b Initialize secondary hyperparameters related to computational complexity
7: Identify x ==
8: Set mg ny and mg
Recall my
X me
x
the number of IPs
9: Draw TP locations x X
described in Appendix F
Tuue expert-related hyperparameters
11: Choose (B, to, nir) as described in Sec. 4.4.3 according the validation performance SVGP(0), where μ. N.A,0,EE E are learned with respect (x ng: unp and (B, nH)
12: Set Eg E Dr, where choose Fg from the best performing SVGP(() of the previous step
Initialize secondary hyperparameters related fine-tuning
Set 2- and dt - for 15 S described in Sec. 44.3
15: Set + 0.1, 1/v2 and d, 0.01 described in Sec. 1.1.3
Tune MoE related hyperparameters
17: Choose (ar,, Ac, MG) as described in Sec. 4.4.3 according to the validation performance of JML from 15
where the model
parameters 0r [xEx] from (88 are learned with respect to (X ngg Ymn) and the model hyperparameters Oh (oa AG.ic)
from 27 are fixed
18: Choose (fo, 40, L, {=1)
) as described in Sec. 4.4.3 according to the validation performance of fM where the model
parameters θ7 \{Ex) from 22 learned with respect to (X Yng) and the model hyperparameters @ (io, 5o, [or)
from
27 fixed
19: Decrease @ (beginning from L) as long as the validation performance of fRNE does not degrade
E Finding diverse IP locations
In order to obtain diverse IP locations with certain distribution, consider two approaches, Stein variational gradient descent (SVGD) (Lu & Wang 2016 Han & Liu 2018) and distributional clustering (DC) (Krishna et al. 2019
Stein Variational Gradient Descent SVGD takes a particle swarm and tries to align the empirical distribution of the particles with a target distribution, which we require the density, as well as its derivative LLu Wang 2016). In addition, the individual particles repel each other, such that we have both diversity and representativeness. In our scenario we have no access to this derivative, such that we resort the work of Han & Liu (2018) that is solely based on the density. Since the particles move freely in the input space and we have to evaluate the target density considerable number of times, suggest applying SVGD, when we deal with well-behaved input spaces and target densities that are easy to evaluate. If the input space is only given through high-dimensional features from a finite set of samples, SVGD might move particles into regions far apart from the data manifold.
Distributional Clustering DC similar to the known k-means clustering (Gan et al. 2020) but solves different ánertía objective, that is modified such that asymptotically, the number of cluster centers x xc, the distribution of the training data is preserved (Krshhna et al. 2019). Under the standard
k-means clustering objective, we would observe X, ~ p=4a GGraf & Luschgy 2007), where it was X,, ~ p. Since we intend to use clustering for sub-sampling rather than identifying a fixed number of true cluster centers, we deal with a comparably large number of cluster centers, here. Thus, we will use DC so as obtain representative set of IPs. Due to very mild assumptions on the problem, DC specifically easy perform in higher dimensions.
42

---

Published in Transactions on Machine Learning Research (12/2023)
Dealing with local optima of DC
The inertia objective of DC is given by
inertiapc/c(Xx]= Ecee Exer. aa log || c|l,
(42)
where
== f= X, ||0 dll s | Vel, Ec]
(43)
are those elements in X, that are closest to the center c.
In the classical Lloyd-step the centers are updated so as to minimize the intra-cluster inertia, which given in the case of DC by
= argmin Exer. 1mye log | - =11--
(44)
It is a known problem that m-menns-rllated inertia objectives suffer from local optima (Arthur & Vassilvitskii 2007]: The converged solution of cluster centers will typically lie close to their initialization. One way to tackle this issue in practice is to run multiple repetitions of the procedure, followed by choosing the solution with minimal inertia. Unfortunately, the amount of local optima increases with the number of cluster centers. In our case, where we use DC for sub-sampling rather than clustering in its usual sense, we deal with a large number of clusters such that this strategy becomes computationally tedious.
Complementary to rumning multiple repetitions of k-means, we will extend the state-of-the-art method k-means++ for choosing the initial set of clusters in more sophisticated way, where we additionally account for the training distribution. Given the inertia objective
[ Σn min.ge ||ا; c|l
of the cluster centers c, the k-means++ procedure builds the set of initial cluster centers as follows: Draw the first center e₁ randomly from X,, Then keep track of the current closest squared distance
d
min || -cll jl.......
(45)
of cach clement Hj Xx, to the far drawn centers C1y..., Cm and sample the next center Cm+11 with probability x (d from X,, This procedure is repeated until the desired number of cluster centers is reached.
The advantage of k-means + is that the initial centers are more diverse than if they were sampled at random from X,, However, its standard form, the centers ialized by k-means ++ are themselves distributed flatter than X,. And so, the case of DC, we propose the following adjustment for a distributional k-means++:
We sample with probability x (d.))2z)))1) from Xn,, where X, p.
Symmetrized DC for molecules Since any symmetric molecule has multiple equivalent representations, care must. be laken when measuring distances in DC. The key idea is to always compare the two configurations in its closest representation. Using the notation from Appendix H.I let
d(z,2') 2') minig.ss [1) ((T,22)]|
be the symmetrized distance between two molecule representations. The symmetrized DC algorithm is then obtained by replacing all occurrences of || - z'|l with d(z, z') in the cluster assignments Les the objective inertiapple/X,,, the cluster updates e* and closest distances d,n from Equations 43
42
44 and 45
F Design choices of the sparse MoE model
In Sec. we have made several design choices with computational feasibility in mind. We will discuss these summarized in this section.
43

---

Published in Transactions on Machine Learning Research (12/2023)
The gate model While in Sec. 4.3 we have chosen the gate gi gV) to be GPR model, note that any choice of model with sufficient flexibility would have been possible. GPR features universal approximation properties, which makes it a favorable choice.
Furthermore, the gate should come with a small degree of freedom to prevent compared to the experts to prevent those from overfitting during the training of the MoF. For this reason, and the fact that we have no ground truth labels for the training of the gate anyhow, we choose our GPR-based gate to be sparse.
Finally, note that we share the set of gate IP locations X across all gate channels. While this is not necessary, it simplifies our method without costs as the MoE is rather insensitive concerning the gate TP locations, as long as these are well-spread.
The expert models While we made clear why we use GPR experts in our work, we left open in Sec. 1.3 whether these experts should be sparse or dense. Here, the deciding factor the amount of n training samples that we have to deal with: When n goes beyond a few thousand, we suggest switching to sparse GPR experts for computational reasons. Note that after training of the MoE, it is also possible to switch back to full GPR experts, if one aims for a high accuracy predictor. For the purpose of AL, this is not necessary.
Similar to the gate, we share the IP locations X across all experts, which simplifies our model. In contrast to the gate situation, the MoE is sensitive to the choice of expert TP locations. Now, if we were to allow individual IP locations for each expert, an elsewise locally underperforming expert might work better than the remaining experts due to a lucky choice of its individual IPs. Subsequently, this would result in a sub-optimal gate and, hence, ultimately in a wrong superior training density estimate.
For better generalization, if our MoE model comprises dense GPR experts, we will either have to rely on individual training sets for the experts and the gale, or use leave-one-oul expert responses on a shared training set.
In the sparse expert case, it is necessary to learn reasonable inducing values b4 prior to the actual learning procedure of the MoE to not gel stuck in a spurious solution. Therefore, there should be a short pre-training phase for each individual expert.
In addition whether not the experts are sparso-the the shared expert parameters μg, Ag 0, Eg should be initialized reasonably. In this regard, we suggest training a single, global expert model before the traininn of the actual experts to obtain those initial parameter estimates for which we have no prior knowledge. If we assume isotropic bandwidths to be sufficient, we can simply set ΣE = IgIa and learn the scalar og >> 0 instead. Note that, from practice, the training of the MoE suffers tremendously from online changes of the expert bandwidths. Thus, we suggest to keep Eg fixed after initialization.
Finally, note that, if we stick with sparse experts after training of the MoE, it can be beneficial for the prediction accuracy to re-train the MoE, where we keep the gate fixed. In this post-processing step, we would like to apply larger learning rates on the experts to escape local optima. However, larger learning rates also lead to underperforming intermediate steps, in which an actively trained gate might reject the best fitting expert at random therefore pushing the gate towards a local optimum. Keeping the pre-trained gate fixed at this point prevents this undesired behavior.
The IPs Recall that we have set the covariance S = 0 of the inducing value distribution to zero, whereas it could have also been a diagonal or positive definite matrix. Playing around with this parameter, we have seen significant improvement that would justify the considerable amount of additional model parameters from a computational point-of-view.
In our approach we suggest. keeping the TP locations fixed, which is also for reasons of computational feasibility, but, more importantly, adaptive IP locations come along with heavy prediction instabilities during the training.
We found it necessary and sufficient to initialize the IP locations by state-of-the-art methods, as described in Appendix E
The MoE objective For the training of our MoE in Sec. 4.4.1
we added a penalty on small bandwidth choices. As described in Lepski 1991 Lepski & Spokoiny 1997 the optimal bandwidth choice is the largest
44

---

Published in Transactions on Machine Learning Research (12/2023)
one that is capable of modeling the function. Now that are able to model comparably flat function by small bandwidths, as long we have got enough training support, it can occur that, with no regularization, we choose a too-small bandwidth for such a flat region. A too-small bandwidth choice might cause overfitting. But even worse, in the subsequent AL loop the flat region is falsely identified as complex, leading to more training queries in this location, which then allow for even smaller bandwidths to model this flat region. We will demonstrate this pathological behavior for the unregularized case on toy-data Sec. 5.1
The gate noise Like already mentioned in Sec. 4.4.1 it is possible to tune in the training process:
Remark 8. Shazeer et ol. (2017) proposed to learn the s parameter by odding a penally term to the main objective that penalizes the imbalance of how likely training inputs are assigned to each expert: Let Ti E [0, 1] be the expert assignment probabilitiess of Ib and define лB E Tb. Then they add penalty V [ to the objective, which is the squared coefficient of variation a coefficient that accounls for the non-uniformily of a set of positive variables.
We justify our simple heuristic to shrink in a static way as follows: Recall from Sec. 4.3 that prevents premature commitment to spurious solution. When treating s as trainable parameter, it does not decay towards zero. Maintaining the noise then prevents the locally best-performing experts from converging by randomly withholding training samples. For this reason, we find that s behaves best when docaying towards zero as the training progresses.
G
Supplemental results on the Doppler experiment
G.1 The single-scale GPR model
When training single-scale GPR model on the Doppler dataset, the tuned bandwidth parameter will typically take intermediate value, trying to compromise between more complex and simpler regions. This is reflected in the predictions in Fig. I6 where the single-scale GPR model suffers from the inhomogeneous structure, underfitting the complex region to the left while simultaneously overfitting the simple region to the right.
In Fig. 17 we compare the performance of our multi-scale MoE approach to the single-scale GPR model. The consistently inferior performance of the single-scale GPR model shows that the issue above persists even for large training sizes.
G.2 Necessity of the small bandwidth penalty
In Appendix we discuss overfitting issues with too small local bandwidth estimates as a consequence of inadequate regularization of LOB. To address this issue, we have proposed to penalize such small bandwidth choices by Ozpen, ((,,,,(,,,, X,, Y,, B, w. θ) with the penalty term pen, from (25 and scaling factor 0, 0.
-5
-5
-10
-10
0.2
10a1
10-7
10-
10a
log-scale
Figure 16: The Doppler experiment: An exemplary dataset and the predictions a global GPR model, shown on natural x-scale (left) and on logarithmic x-scale (right).
45

---

Published in Transactions on Machine Learning Research (12/2023)
max AE
RMSE
10a
GPR, D
MOE, Q
MoE, PSup
10-1
10
10a
10
Training Size Figure 17: The Doppler experiment: The max AE (left) and the RMSE (right) of our proposed MoE model in comparison to a single-scale GPR model. The results are averaged over 20 repetitions.
Training Size n
Predicted local bandwidth estimates
10
10
10
0.6
Figure 18: The Doppler experiment: An actively sampled dataset (top) with our MoF fit at = 212 training samples without small bandwidth penalty (9, = 0), and the associated LOB estimate (bottom).
Now, while the LOB estimate with 0, 0.5 (see Fig. consistently behaves as expected, we show for comparison typical LOB estimate in Fig. 18 that results from applying no regularization (0,, = 0). By chance- here, the flat region of the Doppler function to the right-te the trained model suffers from massive overfitting by too small LOB estimates. These falsely obtained small LOB estimates then lead to overestimation of LFC, which subsequently results in detrimental oversampling of these locations by the AL procedure.
Supplemental on the malonaldehyde MD simulation experiment
H.1 The sGDML model
The GDML model by Chmicla al. (2017) represents the goometry x [ Ral R of each molccule
in terms of the reciprocal distances ((z) R Rill- of all com--parings achieve roto-translational invariance of the input. This representation gives total d a(a 1)/2 input features. The similarity of
46

---

Published in Transactions on Machine Learning Research (12/2023)
RMSE, Q 1/A1
0.7
Isotropic SGDML, p p= U1A)
0.6
Anisotropic SGDML, p= 1-(1)
Anisotropic MoE, p= 1/(2) Anisotropic MoE, p pera.n
0.4
0.3
512
1024
2048 Training Size n
4096
Figure 19: The RMSE under the uniform test distribution for different variants of sGDML and training distribution at varying training size: The performance is given for passive sampling, using the original isorropcc sGDML (dottod), anisotropic sGDML (dash-dotted) and our MoE model with anisotropic sGDML experts (dashed), and for the MoE model, applying the proposed superior sampling scheme (solid). The results are averaged over repetitions.
pair of configurations (z, F) and (z', F") ) is then given by the extended covariance function
Cov(E, E*) = k()(2), ((='),
Con(E, F*) dk(((2), ((2')) d#(
db dr
Cov(F, [=)(()]
T db (l( do(2') dx ddd dx
Hence, we denote the overall kernel function of two configurations by
k<2, 2') = P). E I ((3++1)x(3a+1)
Atoms of the same type are physically identical and therefore exchangeable, albeit only small subset of such symmetrics is excrcised at given (low) MD simulation temperature. Full permutational invariance is only needed when enough energy is put into the system for all bonds to break and all atoms to disassociate.
automatically identifies all accessed atom permutations from the training set and adds this symmetric prior to the covariance function. Formally, lot (T.).=1 be alomic permutations that lead to an equivalent molecular representation. Then, the extended symmetric kernel of sGDML is given by
The symmetric extension sGDML (Chmiela et al. 2018 2019
((2,²')= Σ
Σ (1222
(46)
Malonaldehyde possesses s such permutations.
Remark 9. The identified set of permutations is transitively closed to form group. Under isotropy, it suffices to permute only one of the two configurations given to the kernel: Permuting both entries (as in (46) equals permuling one enlry and mulliphying by the conslant 8. However, if the applied bandwidth is not of the form Σ- = o.Id this property does not hold.
H.2 The malonaldehyde MD simulation experiment under uniform test distribution
In this sccnario, we assume a uniform test density q q=)(x)) Accordingly, we weight the validation and test MSE by the importance weights 1/px X X and (x,x We draw the initial expert training set X, of size n = 29 and the gate training set X via importance sampling from the remaining pool with weights ng 1/px (X poot \ (X U X))). By this is X, (Xx).
In Fig. 20 we show the estimates of LOB, LFC, and the superior training density under the pool test distribution, evaluated the relaxed configurations of malonaldehyde. The LFC estimates in Fig. 20 confirm
47

---

Published in Transactions on Machine Learning Research (12/2023)
LOB
Optimal Dersily Estimate 20Rk
1275
063
122
1175
132
090
120
108
036
051
048
0
-150
100
50
120
-100
150
-150
-300
Figure 20: Estimates of LOB (left), LFC (middle) and the superior training density (right) under the pool test distribution ((), evaluated at the relaxed malonaldehyde configurations, plotted with respect to the angles of the two aldehyde rotors of malonaldehyde.
our expectation that the transition areas are more complex model than the regions near the stable configurations. Subsequently, active sampling scheme shifts sample mass away from the stable regimes in lavor of the transition areas.
We have plotted the error curves of passive and active sampling schemes in Fig. 19 When estimating the relative sample size (31) that we require to achieve the same RMSE via active sampling compared to random test sampling, we obtain e(Fsoon Psup = 0.965 + 0.009. This means that we save about 3.5% of samples with our active sampling scheme. With similar calculations, we save about 27%, when comparing the original sGDML approach with passive sampling to our MoE model with active sampling.
H.3 A realistic MD simulation AL scenario
In the realistic ab initio FF reconstruction AL scenario, we begin by sampling the initial training set (X no X,,)) as well as the validation set (X Y.a) by simulating the true MD trajectory solving the computationally expensive Schrödinger equation. Estimate the initial MD density Px,0 based on (X,,, U vaa))))
For :N:
Set g Px. x+1
to encourage exploration in early iterations and exploitation in later iterations
Estimate the model f = Froe based on (X nu Yay))
GPR,a based on g and Jr. Estimate Psop
Sample a large pool (X poot, Y p00), of, e.g. 100,000 candidates by simulating the approximate MD trajectory using the computationally cheap model Лr. While simulation, avoid unreliable out- of-distribution predictions, e.g., by resetting the trajectory in some whenever Px, drops below reasonable threshold.
Estimate the trajectory density Px.++1 pool
Update the training set (X X by selecting input candidates from the pool X pool with GPR, and estimating the respective labels solving the Schrödinger equation
distribution Psup
48


Reviews:
[{"review_id": "fAQ0yluBZt", "reviewer": "Reviewer_kvMm", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- Related works discussion appears satisfactory.\n- Presentation is relatively clear, although some reptition could be avoided to improve readability. \n- A toy experiments is treated in-depth to illustrate the method.\n- A real world experiment is also show-cased. \n\nWeaknesses:\n- The method is rather involved. I am too far from the field to asses if all the choices are necessary: (cf question below, is a MoE of GP for the gate functions also necessary?)", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The present work introduces an active learning framework with a mixture of expert Gaussian Processes (GP) trained variational. The proposed AL strategy follows the steps of Pankin et al . 2021, formulated for local polynomial smoothing (LPS).\nFirst, the authors adapt the former work’s framework to GPR. They give an expression for the asymptotic local function complexity under the GPR framework. They deduce from it an expression for a superior training density to be used in active learning. \n\nTo be implemented in practice the scheme requires knowing the test distribution, which is not usually the case, instead they rely on an unnormalised estimate and using an importance sampling estimate of the density. The assumption is that one has access to a large amount of unlabelled points that form a pool from which the dataset will be constructed. One the optimal training density has been derived, self-normalised importance sampling is used to estimate the probability of each point in the pool.  \n\nThe authors specify their parametrisation of the sparse mixture of Gaussian process as a mixture of experts of fixed logarithmically space bandwidths. After training they can also propose an estimation of the locally optimal bandwidth as a weighted geometric average of the experts’ bandwidths, which in turn can be used to fix the sparsity inducing parameter of the model.\n\nThe authors offer a discussion on their hyperparameter choices, and their way of choosing how to initialise the IPs in a diverse way to minimise the distance to an evaluation point x.\n\nExperiments:\n\nIn a 1d Doppler effect experiment, the authors show marginal improvement over the LPS, even with optimal training density and clear improvement over the random test sampling for all competing methods. They also use the simple example to test their model under heteroscedastic noise.\nNext the authors turn to a force field reconstruction problem for a 9 atoms molecules and find improvements brought by their methods attributed both to the mixture of experts model and the superior training density over common symmetric gradient-domain machine learning.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I have the following questions:\n\n1) I do not understand the formulations “top-down\" and “bottom up\" in the following sentences:\n“they then propose to sample training data in a top-down manner from this very distribution, “middle of page 2\n“which refine training data bottom-up in a greedy way” bottom of page 4\n\n2) p7 can read: “As already noted in the introduction, the results to LFC and the optimal training density of LPS indicate their problem intrinsic nature, as they reflect no direct dependence on the LPS model.’’  — Except for the order Q correct?\n\n3) Adding more precision on how the training objective is derived would be helpful. The $q_l$ is the marginal distribution of the $l$-th component of the MoE ? (Section 4.4.1) It is not straightforward to understand also the definition of the $P_{b,l}$. How should one interpret the marginalised $f$? In the definition of $w \\propto q/p$ the appearing $q$ is the test distribution and the appearing $p$?\n\n4) Could a simpler model for the gate functions be employed? With which loss on performance?\n\n5) Could the authors clarify how to adapt their work to “on the fly\" sampling, that is when there is not a $X_{\\rm pool}$ available from the start? That is crucial for a realistic implementation in ab initio force field as producing the unlabelled data (the configurations) as a pool a priori is requires to compute the forces to run the MD simulation and is therefore not cheaper than getting also the labels.\n\n6) It is arguably not fair to say that \"current ML-based FFs typically rely on rather naive exhaustive sampling schemes to gather training data\" page 22,  cf for instance \"Section 3.3.3. Adaptive Sampling\" of the cited review by Keith et al 2022.\n\n7) The interpretability of LFC (discussion section) for high-dimensional input spaces should be relativized, as one cannot visualise this function expect in low-dimension, or using relevant, previously known, low-dimensional projections. Or maybe the authors want to point out that one can locate high LFC points from the pool? In which case the precision would be useful. \n\nMinor comments: \n- The acronym MISE is not introduced p2.\n- From an editing point of view, keeping figures to the top of pages would improve readability by avoiding lost sentences/paragraph title between figures.", "broader_impact_concerns": "none", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "we3I5ks3Au", "reviewer": "Reviewer_ZcVs", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: \n\n* The authors put forward a principled approach for the important problem of choosing optimal training sets. Active learning techniques are needed in many important fields such as multi-task learning or federated learning. \n\n* Extensive numerical experiments illustrate the proposed methods. \n\nWeaknesses: \n\n* It would be helpful to illustrate the key concepts in the most simple ML setting. In particular it should be made more explicit how the different concepts such as LFC correspond to the key components of ML applications, i.e., is LFC a property of a data (distribution) or a specific hypothesis (target function) or a  model (hypothesis space) or a combination of model/loss function ( see Ch. 2 of [Ref1] for a discussion of these three main components). \n\n* avoid the excessive use of the word \"framework\". try to use more specific words if you mean, e.g., an algorithm or a probabilistic model or a \"complexity measure\" for a model. \n\n* The main contribution should be formulated as an algorithm pseudo-code and not only as Figure 2 . This algorithm description must clearly indicate the required input data, hyper-parameter and delivered output. You must then also clearly explain how you applied that algorithm in your numerical experiments. \n\n* There should be more discussion about the (sub-) optimality of your active learning methods. How close is your method to fundamental limits (see e.g. Sec. 6.4 of [Ref1] for fundamental lower bounds on the sample size required by plain old linear regression) \n\n* The use of language could be sharpened (avoiding jargon), e.g., \"The scalability of this estimate makes...\", \"such that the optimal training inputs in Eq. (3) can asymptotically be obtained..\", \"Let us consider a family ...of kernel machines..\", \"...of the hypothetically optimal training set of a model in the asymptotic limit of the sample size.\"\n\n[Ref1] A. Jung, \"Machine Learning: The Basics,\" Springer, 2022.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper present as its main contribution an active learning framework for GPR-based methods.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "address the weaknesses", "broader_impact_concerns": "n/a", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "JzpfXLse77", "reviewer": "Reviewer_WPPN", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- the expression for the local function complexity of a GP is derived.\n- the paper combines ideas from several different fields (local polynomial smoothing, Gaussian processes, active learning, etc)\n\nWeaknesses:\n- the paper is long and sometimes hard to follow.\n- the main theoretical paper inspiring this work is a preprint.\n- the state of the art on GP models is relatively outdated.\n- GP are well known for scaling poorly with the dimension, while this is the reason to use them here.\n- the baseline for the experimental part, uncertainty sampling, is weak.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper deals with selecting training data in a more efficient way than random sampling. It is based on a sequential active learning strategy where a surrogate model is used to approximate the underlying unknown function while new points are added using a criterion based on local function complexity. The surrogate is a mixture of sparse Gaussian process (GP) regression models. As for the criterion, it relies on the derivation of proxy for the local function complexity. A 1d toy example is presented, then a 32d real world one.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "(critical) *Clarity*. To help enhance the clarity of the paper, it may be worth separating the sparse GP aspect from the main subject. Most derivations are presented for isotropic kernels, but then are used on the anisotropic sGDML model.\nThen, the considered problem is only formally presented on page six.\n\n(critical) On the selection of inducing points selection, the state of the art to compare with can be found in more recent articles, say starting from Moss, H. B., Ober, S. W., & Picheny, V. (2023, April). Inducing point allocation for sparse gaussian processes in high-throughput bayesian optimisation. In International Conference on Artificial Intelligence and Statistics (pp. 5213-5230). PMLR.\n\n(critical) The modeling of non-stationary functions is a recurrent topic with GPs, hence the use of mixtures and sparse inputs is only one option among several, see e.g., \n- Gramacy, R. B., & Lee, H. K. H. (2008). Bayesian treed Gaussian process models with an application to computer modeling. Journal of the American Statistical Association, 103(483), 1119-1130.\n- Gramacy, R. B., & Apley, D. W. (2015). Local Gaussian process approximation for large computer experiments. Journal of Computational and Graphical Statistics, 24(2), 561-578.\n- Roininen, L., Girolami, M., Lasanen, S., & Markkanen, M. (2018). Hyperpriors for Matérn fields with applications in Bayesian inversion. Inverse Problems and Imaging, 13(1), 1-29.\n- Rullière, D., Durrande, N., Bachoc, F., & Chevalier, C. (2018). Nested Kriging predictions for datasets with a large number of observations. Statistics and Computing, 28, 849-867.\n- Cutajar, K., Bonilla, E. V., Michiardi, P., & Filippone, M. (2017, July). Random feature expansions for deep Gaussian processes. In International Conference on Machine Learning (pp. 884-893). PMLR.\n\nWhen using these methods, the typical criteria would not be the Gaussian process uncertainty but the integrated mean squared prediction error, or related criteria. See for instance Sauer, A., Gramacy, R. B., & Higdon, D. (2023). Active learning for deep Gaussian process surrogates. Technometrics, 65(1), 4-18.\nP19: why not sampling where the predictive variance is the largest?\n\n(critical) A discussion on the poor scaling of GPs with the dimension is needed.\n\n\nSome additional related works may include:\n- Lederer, A., Umlauft, J., & Hirche, S. (2019). Uniform error bounds for gaussian process regression with application to safe control. Advances in Neural Information Processing Systems, 32.\n- Kanagawa, M., Hennig, P., Sejdinovic, D., & Sriperumbudur, B. K. (2018). Gaussian processes and kernel methods: A review on connections and equivalences. arXiv preprint arXiv:1807.02582.\n- Mak, S., & Joseph, V. R. (2018). Support points. The Annals of Statistics, 46(6A), 2562-2592.\n\nThe regularity of $f$ is discussed several times, which for GPs is typically handled with Matérn kernels. \n\nThe link between (4) and (5) with the introduction of Q is, in my opinion, not straightforward.\n\n(17) is said to not depend on p, v and n while it appears in (17).\n\nA full pseudo code of the approach would be appreciated, rather than Figures.\n\nFigs. 6 and 7: perhaps use consistent x-axis\n\nFig. 12: Perhaps use letters rather than positions\n\nA moderate toy dimensional example would be welcomed, to see the scaling in 2 or 3d when strong baselines can still be used.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper presents a novel approach to address inhomogeneities in real-world data, particularly focusing on variations in observation noise levels and structural complexities. The authors propose a model-agnostic active learning framework, leveraging local function complexity (LFC) estimates derived from Gaussian Process Regression (GPR). This approach aims to enhance robustness and scalability in handling large input space dimensions.

The paper makes a substantial contribution to the field of active learning by addressing key challenges in handling inhomogeneous data complexities. Its blend of theoretical innovation and practical application positions it as a valuable resource for researchers and practitioners in machine learning and related fields.

However, the reviewers still have some concerns about the writing (clarifying in the text the meaning of bottom-up/top-down, adding an algorithm in the main part), and the experiments (discussing the lower bound of RMSE). To further enhance the quality of the paper, the authors should address these issues in the revised version.
Decision: Accept (poster)