Local Function Complexity for Active Learning via Mixture of Gaussian Processes

TL;DR

This paper introduces a Gaussian process-based local function complexity measure to enhance active learning, demonstrating improved efficiency in reconstructing quantum chemical force fields with limited data.

Contribution

It develops a robust, scalable GPR-based local complexity measure and integrates it into an active learning framework for real-world, high-dimensional problems.

Findings

Effective in low-dimensional synthetic data

Achieves state-of-the-art performance in quantum chemistry task

Reduces training data requirements significantly

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.

Equations

$$\text{MSE}\left(x,\hat{f}|{\mathbf{X}}_n\right)=\text{sideset}{}_{{\mathbf{Y}}_n}\text{mathds}\mathrm{E}\left[({\hat{f}}_{{\mathbf{X}}_n,{\mathbf{Y}}_n}(x)-f(x){)}^2\right]=\text{sideset}{}_{{\mathbf{\epsilon }}_n}\text{mathds}\mathrm{E}\left[({\hat{f}}_{{\mathbf{X}}_n,f({\mathbf{X}}_n)+{\mathbf{\epsilon }}_n}(x)-f(x){)}^2\right].$$

$$\text{MISE}\left(q,\hat{f}|{\mathbf{X}}_n\right)={\int }_{\mathcal{X}}\text{MSE}\left(x,\hat{f}|{\mathbf{X}}_n\right)q(x)dx.$$

$${\mathbf{X}}_n^{\prime }\approx {\mathrm{argmin}}_{{\mathbf{X}}_n\in {\mathcal{X}}^n}\text{MISE}\left(q,\hat{f}|{\mathbf{X}}_n\right).$$

$${\Sigma }_{\hat{f}}^n(x)={\mathrm{argmin}}_{\Sigma \in \mathcal{S}}\text{MSE}\left(x,{\hat{f}}^{\Sigma }|{\mathbf{X}}_n\right),$$

$${\hat{f}}_{\text{LPS}}^Q:=m_Q^{{\Sigma }_Q^n(x)}(x)$$

$$p_{\text{Opt}}^{Q,n}(x)\propto {\left[{\mathfrak{C}}_Q^n(x)q(x)\right]}^{\frac{2(Q+1)+d}{4(Q+1)+d}}v(x{)}^{\frac{2(Q+1)}{4(Q+1)+d}}(1+o(1)),$$

$${\mathfrak{C}}_Q^n(x):={\left[\frac{v(x)}{p(x)n}\right]}^{\frac{d}{2(Q+1)+d}}{\left|{\Sigma }_Q^n(x)\right|}^{-1}={\left[\frac{v(x)}{p(x)n}\right]}^{\frac{d}{2(Q+1)+d}}{\sigma }_Q^n(x{)}^{-d}.$$

$${\mathfrak{C}}_1^n(x)\propto \text{trace}(D_f^2(x){)}^{\frac{2d}{2(Q+1)+d}}(1+o(1))$$

$${\Sigma }_{\text{GPR}}^n\propto n^{-\frac{1}{2\alpha +d}}.$$

$$\hat{y}({\mathbf{X}}_{\ast })\sim \mathcal{N}(\cdot ;{\mathbf{\mu }}^{\ast }({\mathbf{X}}_{\ast }),{\mathbf{C}}^{\ast }({\mathbf{X}}_{\ast })|{\theta }_e)$$

$${\mathbf{\mu }}^{\ast }({\mathbf{X}}_{\ast })={\mathbf{K}}_{\ast †}{\mathbf{K}}_{†}^{-1/2}({\tilde{\mathbf{\mu }}}_{†}-{\mathbf{K}}_{†}^{-1/2}\mathbf{\mu }({\mathbf{X}}_{†}))+\mathbf{\mu }({\mathbf{X}}_{\ast }),$$

$${\mathbf{C}}^{\ast }({\mathbf{X}}_{\ast })=\lambda \left[{\mathbf{K}}_{\ast }+{\mathbf{K}}_{\ast †}{\mathbf{K}}_{†}^{-1/2}({\tilde{\mathbf{S}}}_{†}-{\mathcal{I}}_m){\mathbf{K}}_{†}^{-1/2}{\mathbf{K}}_{\ast †}^{\top }\right]+\text{diag}(\hat{v}({\mathbf{X}}_{\ast })),$$

$$\text{sideset}{}_{u\sim Q(\cdot |{\mathbf{X}}_n,{\mathbf{Y}}_n)}\text{mathds}\mathrm{E}\log (P({\mathbf{X}}_n,{\mathbf{Y}}_n|u))-\mathcal{K}\mathcal{L}\left[Q(\cdot |{\mathbf{X}}_n,{\mathbf{Y}}_n)\Vert P\right]$$

$$\approx$$

$$P_b:=\text{sideset}{}_{u\sim Q(\cdot |{\mathbf{X}}_n,{\mathbf{Y}}_n)}\text{mathds}\mathrm{E}\log \int P(y_b|f)P(f|u,x_b)df.$$

$${\hat{f}}_{\text{MoE}}(x)=\text{sideset}{}_{l=1}^L\sum G(x{)}_l{\hat{y}}_l(x),$$

$$\displaystyle\operatorname*{\textbf{soft\hskip 1.00006ptmax}}(\bm{a})_{i}:=\exp\{\bm{a}_{i}\}\Big{/}\sideset{}{{}_{l=1}^{L}}{\sum}\exp\{\bm{a}_{l}\},$$

$$G(x)=\mathrm{softmax}({\tilde{h}}_1(x),\dots ,{\tilde{h}}_L(x)),\text{where}{\tilde{h}}_{i_l}(x)=\{\begin{array}{ll}{h}_{i_l}(x)& ,l<\kappa \\ -\infty & ,l\ge \kappa +1\end{array}$$

$$\Theta =(\{{\theta }_{e_l}{\}}_{l=1}^L,\{{\theta }_{g_l}{\}}_{l=1}^L,\kappa ,\mathfrak{s}).$$

$${\mathfrak{C}}_{\text{GPR}}^n(x):={\left[\frac{1}{p(x)n}\right]}^{\frac{d}{2\alpha +d}}{\left|{\Sigma }_{\text{GPR}}^n(x)\right|}^{-1}.$$

$$p_{\text{Sup}}^{\text{GPR},n}(x)\propto {\left[{\mathfrak{C}}_{\text{GPR}}^n(x)q(x)\right]}^{\frac{2\alpha +d}{4\alpha +d}}v(x{)}^{\frac{2\alpha }{4\alpha +d}}(1+o(1)).$$

$${\mathfrak{C}}_{\text{GPR}}^n(x)={\left|{\Sigma }_{\text{GPR}}^n(x)\right|}^{-1}\text{and}{p}_{\text{Sup}}^{\text{GPR},n}(x)\propto {\left[{\mathfrak{C}}_{\text{GPR}}^n(x)q(x)v(x)\right]}^{\frac{1}{2}}.$$

$$\text{MISE}\left(q,\hat{f}|{\mathbf{X}}_n\right)=\text{MISE}\left(q,\hat{f}|{\mathbf{X}}_{n^{\prime }}^{\prime }\right)(1+o(1)).$$

$$\displaystyle\textstyle\widehat{\mathfrak{C}}^{n_{k}}_{\text{GPR}}(x)\propto\left[1\big{/}p_{k}(x)\right]^{\frac{1}{2\alpha+d}}\left|\widehat{\Sigma}^{n_{k}}_{\text{GPR}}(x)\right|^{-1},\;\;\text{and}$$

$${\hat{p}}_{\text{Sup}}^{\text{GPR},n_k}(x)\propto {\left[{\hat{\mathfrak{C}}}_{\text{GPR}}^{n_k}(x)q(x)\right]}^{\frac{2\alpha +d}{4\alpha +d}}\hat{v}(x{)}^{\frac{2\alpha }{4\alpha +d}}.$$

$$\text{norm}(p)={\left|{\mathbf{X}}_{pool}\right|}^{-1}\text{sideset}{}_{x\in {\mathbf{X}}_{pool}}\sum p(x)/{\hat{p}}_{\mathcal{X}}(x).$$

$$\text{norm}(p)\approx {\int }_{\mathcal{X}}\frac{p(x)}{{\hat{p}}_{\mathcal{X}}(x)}{p}_{\mathcal{X}}(x)dx=\frac{1}{c}{\int }_{\mathcal{X}}p(x)dx,$$

$${\hat{f}}_{\text{MoE}}(x)=\text{sideset}{}_{l=1}^L\sum G(x{)}_l{\hat{y}}_l(x),$$

$${\theta }_{e_l}=({\mu }_E,{\lambda }_E,\hat{v},{\sigma }_l{\Sigma }_E,{\mathbf{X}}_{†}^E,{\mathbf{\mu }}_{†}^{e_l},0).$$

$${\theta }_{g_l}=({\mu }_{g_l},{\lambda }_G,0,{\sigma }_G{\mathcal{I}}_d,{\mathbf{X}}_{†}^G,{\mathbf{\mu }}_{†}^{g_l},0).$$

$${\hat{\Sigma }}_{\text{GPR}}^n(x)={\hat{\sigma }}_{\text{GPR}}^n(x){\Sigma }_E,\text{where}{\hat{\sigma }}_{\text{GPR}}^n(x)=\exp \left\{\text{sideset}{}_{l=1}^L\sum G(x{)}_l\log ({\sigma }_l)\right\}.$$

Full text

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 treat specific kinds of inhomogeneities. For example, Kersting et al. (2007); Cawley et al. (2006) deal with heteroscedasticity by reconstructing a local noise variance function that is used to adapt the regularization of the model locally. Some approaches adjust bandwidths locally with respect to the input density (Wang & Wang, 2007; Mackenzie & Tieu, 2004; Moody & Darken, 1989; Benoudjit et al., 2002). Inhomogeneous complexity can also be captured using a combination of several 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 regularization according to the inhomogeneities in noise, complexity, and data density. While this is the path we will pursue, the focus in 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 on the informativeness of certain locations of the input space, which subsequently can be used to guide the sampling process during training—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 forecasting (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 fields 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 (Wu, 2019) over sophisticated, model-based AL approaches with strong assumptions as the former are inherently robust and model-agnostic 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 \raisebox{1.0pt}{\varrho}>0, which asymptotically equates the performance of n\cdot\raisebox{1.0pt}{\varrho} active training samples to $n$ random test samples (see Definition 2). Accordingly, we call an AL approach asymptotically superior to random test sampling, if \raisebox{1.0pt}{\varrho}<1. Unfortunately, the performance gain of model-free AL approaches that we observe at small training sizes over random test sampling eventually diminishes completely (\raisebox{1.0pt}{\varrho}=1) with growing training size.

For the described learning task, we therefore require an AL approach that is model-based but comes with mild regularity assumptions at the same time to feature robustness and model-agnosticity to a certain extent.

Recently, Panknin et al. (2021) addressed the outlined AL scenario, where the fundamental idea is to analyze the distribution of the optimal training set of a model in the asymptotic limit of the sample size.

Assuming that this limiting distribution exists, they then propose to sample training data in a top-down manner from this very distribution, knowing that with growing sample size the training set will eventually become optimal. By a top-down AL approach, we mean an (infinite) training data refinement process $x_{1}^{\prime},x_{2}^{\prime},\ldots$ such that—when optimizing an AL criterion with respect to $\{x_{1},\ldots,x_{n}\}$—$\{x_{1}^{\prime},\ldots,x_{n}^{\prime}\}$ asymptotically becomes a respective optimizer as $n\rightarrow\infty$. 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 complexity (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 locally optimal kernel bandwidth (LOB) of the LPS model, calibrated for the local effects of the training input density and noise level. Given a 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 a 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 input 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 on 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 is 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 Gaussian process regression (GPR) model to preserve this property. While we lose the strict asymptotic sampling optimality this way, we 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 inhomogeneities, namely heteroscedasticity and inhomogeneous complexity. Note that it is first and foremost superior for our adaptive bandwidth GPR model.

Fig. 1 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 a 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-agnosticity of the optimal training density of LPS, we replace the herein contained LPS-based LFC estimate for our GPR-based LFC estimate to obtain a superior training density for locally adaptive models. From this point, we can implement the AL framework by Panknin et 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 a scalable LOB estimate for GPR as 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., 1996). 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 an implementation111https://github.com/DPanknin/modelagnostic_superior_training of both, the AL framework and our model. Finally, we propose a 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.

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.i.d. sampling from the superior training density.

Force field reconstruction   Quantum interactions exhibit multi-scale behavior due to the complex electronic interactions that give rise to any observable property of interest, like the total energy or atomic forces of a system (Bereau et al., 2018; Yao et al., 2018; Grisafi & Ceriotti, 2019; Ko et al., 2021; Unke et al., 2021a). To demonstrate the scalability of our approach, we consider a force field reconstruction problem of a molecule with 27 dimensions, where the application of the LPS-based 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 visualizations of the scalar-valued LFC function.

We begin by discussing our work in the context of related work in Sec. 2. Next, we give a formal definition of the considered regression problem and the asymptotic AL task, and review asymptotic results for LPS and GPR in Sec. 3. In Sec. 4, we describe our MoE model and derive the GPR-based LFC and superior training density estimates. In Sec. 5, we then describe our experiments and results, which will be further discussed in Sec. 6. We finally conclude in Sec. 7.

2 Related work

Choice of MoE, experts and the gate   The common assumption of MoE approaches is that the overall problem to infer is too complex for a single, comparably simple expert. This is the case, for example in regression of nonstationary or piecewise 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 (2001); Meeds & Osindero (2006); Yuan & Neubauer (2009); Yang & Ma (2011); Chen et al. (2014).

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 inhomogeneous 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), where—locally dependent—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 a 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 to train the MoE in an expectation maximization loop, where the likelihood of an expert to have produced a training label functions as a pseudo-label to the gate.

Alternative nonstationary GPs   As opposed to a standard GP that features a stationary covariance structure, our MoE 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:

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 a tree.

(Gramacy & Apley, 2015) identify 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 a GP. Note that our work intends to elaborate the estimation of LFC and model-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. (2019). 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 model-agnosticity of our work in Sec. 5.1.

IP selection   In our work (see 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 IP 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 IP locations be the result of k-Means clustering. This choice of IP 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 IP selection approaches that select columns of the full kernel matrix according to a fixed distribution (Drineas et al., 2005) or one-by-one in a greedy, adaptive way (Smola & Schölkopf, 2000; Fine & Scheinberg, 2001; Seeger et al., 2003). An intensive overview of Nyström method based IP selection methods was given by Kumar et al. (2012), where they also analyzed ensembles of low-rank approximations. We compare our proposed IP choice (29) to the greedy fast forward IP selection approach by Seeger 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 a model, e.g., when they are based on a parametric model (Kiefer, 1959; MacKay, 1992; He, 2010; Sugiyama & Nakajima, 2009; Gubaev et al., 2018), or which refine training data bottom-up in a greedy way to maximize its information content at small sample size, where the information is either based on the inputs only (Seo et al., 2000; Teytaud et al., 2007; Yu & Kim, 2010; Wu, 2019; Liu et al., 2021) or also incorporates the labels (Burbidge et al., 2007; Cai et al., 2013). By a bottom-up AL approach, we mean a training data refinement process that is constructed by choosing the nth input $x_{n}$ as the optimizer of an AL criterion with respect to $\{x_{1},\ldots,x_{n}\}$, when keeping the previously drawn inputs $\{x_{1},\ldots,x_{n-1}\}$ (with labels $\{y_{1},\ldots,y_{n-1}\}$) fixed.

Our work is therefore complementary to the latter kind of approaches which can be better suited in another AL 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 knowledge 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., 2007; 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 on our 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 theoretically 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 (see, e.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 (RBF)-network (Moody & Darken, 1989), using their proposed training dataset.

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 boundaries, which is inherently a model-agnostic task. In addition, since the decision boundaries are a submanifold of the input space $\mathcal{X}$, a substantial part of $\mathcal{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 (Lewis & Gale, 1994; Roy & McCallum, 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 contrast, 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, 2011; 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 (pseudo/̄)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 a mixture of GPR experts, it is straightforward to derive its uncertainty as a mixture of Gaussian process uncertainties (MoGPU) by simply weighting the predictive variances of all experts with respect to the gate 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 Sec. 5.1.

3 Preliminaries

We will now give a formal definition of the regression task, the AL objective and LOB, and a 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 $\textbf{{diag}}(z)\in\operatorname*{\mathds{R}}^{d\times d}$ the diagonal matrix with the entries of the vector $z\in\operatorname*{\mathds{R}}^{d}$ on its diagonal and by $\mathcal{I}_{d}=\textbf{{diag}}(\mathbbmss{1}_{d})$ the identity matrix, where $\mathbbmss{1}_{d}$ is the vector of ones in $\operatorname*{\mathds{R}}^{d}$.

3.1 Formal definition of the regression task and AL objective

Let $f$ be the target regression function defined on an input space $\mathcal{X}\subset\operatorname*{\mathds{R}}^{d}$ that we want to infer from noisy observations $y_{i}=f(x_{i})+\varepsilon_{i}$, where $x_{i}\in\mathcal{X}$ are the training inputs and $\varepsilon_{i}$ is independently drawn noise from a distribution with mean $\operatorname*{\mathds{E}}[\varepsilon_{i}]=0$ and local noise variance $\operatorname*{\mathds{V}}[\varepsilon_{i}]=v(x_{i})$. We denote a training set by $(\bm{X}_{n},\bm{Y}_{n})$, where $\bm{X}_{n}=(x_{1},\ldots,x_{n})\in\mathcal{X}^{n}$ and $\bm{Y}_{n}=(y_{1},\ldots,y_{n})\in\operatorname*{\mathds{R}}^{n}$. For a given model class $\widehat{f}$ that returns a predictor $\widehat{f}_{\bm{X}_{n},\bm{Y}_{n}}$ for a training set $(\bm{X}_{n},\bm{Y}_{n})$, we can define the pointwise conditional mean squared error of $\widehat{f}$ in $x\in\mathcal{X}$, given $\bm{X}_{n}$, by

[TABLE]

Note that via marginalization the conditional mean squared error is no function of the training labels $\bm{Y}_{n}$. Given a test probability density q\in\mathcal{C}^{0}\left(\mathcal{X},\mathds{R}_{\raisebox{0.0pt}{\tiny+}}\right) such that $\mathop{{\int}}_{\mathcal{X}}q(x)dx=1$, the conditional mean integrated squared error of the model under the given training set is then defined as

[TABLE]

With these preparations, the AL task is to construct a training dataset $(\bm{X}_{n}^{\prime},\bm{Y}_{n}^{\prime})$ such that

[TABLE]

3.2 Locally optimal bandwidths, function complexity, and optimal training

Let $\widehat{f}^{\Sigma}$ be a family of kernel machines which is characterized by a positive definite bandwidth matrix parameter \Sigma\in\mathbb{S}^{d}_{\raisebox{0.0pt}{\tiny++}} of an RBF kernel $k^{\Sigma}(x,x^{\prime}):=\left|\Sigma\right|^{-1}\phi(\|\Sigma^{-1}(x-x^{\prime})\|_{\raisebox{-2.0pt}{\tiny$$}})$ for a monotonically decreasing function \phi\!:\mathds{R}_{\raisebox{0.0pt}{\tiny+}}\,\rightarrow\,\mathds{R}_{\raisebox{0.0pt}{\tiny+}}. The well known Gaussian kernel is for example implemented by $\phi(z)=\exp\{-\frac{1}{2}z^{2}\}$.

Given a bandwidth space \mathcal{S}\subseteq\mathbb{S}^{d}_{\raisebox{0.0pt}{\tiny++}} we define the LOB function of $\widehat{f}$ by

[TABLE]

assuming that this minimizer uniquely exists for all $x\in\mathcal{X}$.

Denote by $m^{\Sigma}_{Q}$ the predictor of the LPS model of order $Q$ under bandwidth $\Sigma$ and by $\displaystyle\Sigma^{n}_{Q}$ $\displaystyle:=\Sigma^{n}_{m_{Q}}$ the LOB function (4) of LPS, if it is well-defined. This is the case, e.g., for the isotropic bandwidths space $\mathcal{S}=\left\{\sigma\mathcal{I}_{d}\;\middle|\;\sigma>0\right\}$ under mild assumptions222For LOB being well-defined in the isotropic case, we generally require a non-vanishing bias and variance in terms of a bias-variance-decomposition of the MSE of the predictor in x, for all $x\in\mathcal{X}$. See, e.g., Eq. (36) for the LPS predictor $m_{Q}$, or Silverman (1986); Wand & Jones (1994) in more general., where we particularly can write $\Sigma^{n}_{Q}(x)=\sigma^{n}_{Q}(x)\mathcal{I}_{d}$. We refer to Appendix A for details on the LPS model and asymptotic results. For the optimal predictor

[TABLE]

of LPS, letting $\widehat{f}=\widehat{f}_{\text{\tiny{LPS}}}^{Q}$ in Eq. (3), Panknin et al. (2021) have shown that there exists an optimal training density $p^{Q,n}_{\text{\tiny{Opt}}}$ that allows the optimal training inputs in Eq. (3) to be asymptotically obtained by independently and identically sampling $\bm{X}_{n}^{\prime}\sim p^{Q,n}_{\text{\tiny{Opt}}}$. They have also shown that this density exhibits a closed-form

[TABLE]

where for an arbitrary training dataset $(\bm{X}_{n},\bm{Y}_{n})$ with $\bm{X}_{n}\sim p$, the LFC of LPS is defined by

[TABLE]

The LFC in (7) asymptotically solely depends on the behavior of $f$ as opposed to $p,v$, and $n$: It scales with the local variation of $f$ in the vicinity of $x$. For example,

[TABLE]

is a function of the trace of the Hessian of $f$ (Fan et al., 1997).

The optimal density $p^{Q,n}_{\text{\tiny{Opt}}}$ in (6) implies that we require more training data where the problem is locally more complex (large $\mathfrak{C}^{n}_{Q}$) or noisy (large $v$), or where test instances are more likely (large $q$). 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\in\mathcal{C}^{\alpha}\left(\mathcal{X},\operatorname*{\mathds{R}}\right)$, there is a canonical choice $Q=\lceil\alpha\rceil-1$ 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 on $Q$ negligible, as its choice is driven by the problem intrinsic regularity.

In practice, we obtain $\bm{X}_{n}^{\prime}\sim p^{Q,n}_{\text{\tiny{Opt}}}$ by estimating Eq. (6) and (7) from $(\bm{X}_{n^{\prime}},\bm{Y}_{n^{\prime}})$ with $\bm{X}_{n^{\prime}}\sim p$ for an arbitrary training density $p$, where $n^{\prime}<n$, followed by adding the remaining $n-n^{\prime}$ inputs appropriately (see Sec. 4.2).

The construction of $p^{Q,n}_{\text{\tiny{Opt}}}$ 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 dimensionality regarding robustness and computational feasibility. The goal of this work is to implement the above AL framework but based on a functional LOB estimate in the domain of GPR instead of LPS, since the GPR model class can naturally deal with high input space dimensions (Williams & Rasmussen, 1996). Relying on the model-agnosticity, we expect that LOB estimates based on LPS can be exchanged for LOB estimates based on GPR in the formulation of LFC and $p^{Q,n}_{\text{\tiny{Opt}}}$ 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 we keep a fixed model complexity—in the sense of the number of basis functions—in the former while there is varying model complexity in the latter as 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 asymptotic limit. In particular, given enough samples, there is no need for local bandwidth adaption.

However, there is a mismatch if $f\in\mathcal{C}^{\alpha}\left(\mathcal{X},\operatorname*{\mathds{R}}\right)$ is $\alpha$-times continuously differentiable since the Gaussian kernel is infinitely often continuously differentiable. As shown by Van der Vaart et al. (2007; 2009), in order to obtain optimal minimax-convergence of the predictor (except for logarithmic factors), the associated (global) bandwidth has to follow the asymptotic law

[TABLE]

Note that for $f\in\mathcal{C}^{\alpha}\left(\mathcal{X},\operatorname*{\mathds{R}}\right)$, where the theoretical results of LPS apply, the scaling factor $n^{-\frac{1}{2\alpha+d}}$ of LOB in sample size matches exactly for both classes, LPS and GPR. In our work, we will use (9) to deduce a GPR-based LFC estimate in analogy to the LPS-based LFC estimate (7) by Panknin et al. (2021).

3.4 Preliminaries on 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\in\mathcal{X}^{n}$ and $X^{\prime}\in\mathcal{X}^{m}$ as $\bm{K}^{\Sigma}(X,X^{\prime})=\left[k^{\Sigma}(x,x^{\prime})\right]_{x\in X,x^{\prime}\in X^{\prime}}$. As a shorthand notation we furthermore define $\bm{K}^{\Sigma}(X):=\bm{K}^{\Sigma}(X,X)$.

3.4.1 Sparse variational Gaussian processes

We define the sparse GPR model $\widehat{y}\sim\mathcal{SVGP}(\theta)$ (see, e.g. Williams & Rasmussen (1996); Hensman et al. (2015)) as follows: The sparse GP is described by the (hyper-) parameters $\theta=(\mu,\lambda,\widehat{v},\Sigma,\bm{X}_{\dagger},\bm{\mu}_{\dagger},\bm{S}_{\dagger})$, which are the global constant prior mean $\mu$, the regularization parameter $\lambda$, the label noise variance function $\widehat{v}$, the bandwidth matrix $\Sigma$ of the kernel and the prior distribution, given by the IP locations $\bm{X}_{\dagger}\in\mathcal{X}^{m}$ as well as their inducing value distribution, characterized by the moments $\bm{\mu}_{\dagger}$ and $\bm{S}_{\dagger}$. That is, for the inducing values $\bm{Y}_{\dagger}$ of $\bm{X}_{\dagger}$ we assume $\bm{Y}_{\dagger}=\widehat{y}(\bm{X}_{\dagger})\sim\mathcal{N}(\cdot;\bm{\mu}_{\dagger},\bm{S}_{\dagger})$. 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 $m_{n}=o[n]$ is typically much slower than $n$. If we can assume homoscedastic noise, we let $\widehat{v}(x)\equiv\sigma_{\varepsilon}^{2}$.

The sparse GP then outputs

[TABLE]

for the mean function

[TABLE]

and the covariance function

[TABLE]

where $\widetilde{\bm{\mu}}_{\dagger}=\bm{K}_{\dagger}^{-1/2}\bm{\mu}_{\dagger}$ and $\widetilde{\bm{S}}_{\dagger}=\bm{K}_{\dagger}^{-1/2}\bm{S}_{\dagger}\bm{K}_{\dagger}^{-1/2}$ are the whitened moments of the inducing value distribution (Pleiss et al. (2020), Sec. 5.1), and we have defined $\bm{K}_{*}=\bm{K}^{\Sigma}(\bm{X}_{*})$, $\bm{K}_{\dagger}=\bm{K}^{\Sigma}(\bm{X}_{\dagger})$ and $\bm{K}_{*\dagger}=\bm{K}^{\Sigma}(\bm{X}_{*},\bm{X}_{\dagger})$.

We choose $\bm{\mu}$ to be the constant mean function, i.e., $\bm{\mu}(X)=\mu\mathbbmss{1}_{n}$ for $X\in\mathcal{X}^{n}$, noting that other mean functions are possible. Note that test predictions $\widehat{f}_{\text{\tiny{GP}}}(x)=\bm{\mu}^{*}(x)$ are given by Eq. (11).

The training objective   Let $P$ denote the prior distribution of the inducing function values $\bm{Y}_{\dagger}$ of the IPs $\bm{X}_{\dagger}$ and let $Q$ denote a tractable variational distribution intended to approximate $P(\cdot)\approx Q(\cdot|\bm{X}_{n},\bm{Y}_{n})$. In variational inference, we want to minimize the Kullback-Leibler divergence $\mathcal{KL}\left[Q\|P\right]$ between $Q$ and $P$, which is equivalent to maximizing the data log-evidence $\log(P(\bm{X}_{n},\bm{Y}_{n}))$. As a tractable approximation, we maximize the evidence lower bound (ELBO), given by

[TABLE]

where $P_{b}$ is the predictive log-likelihood in $x_{b}$, marginalized over the variational distribution $Q$, that is,

[TABLE]

3.4.2 Sparse mixture of experts

Given a finite set of expert models $\widehat{y}_{l}$ that are parameterized by $\theta_{e_{l}}$, the MoE model is given by

[TABLE]

where the gate $G\!:\mathcal{X}\,\rightarrow\,[0,1]^{L}$ is a probability assignment of an input $x$ to the experts. In particular, it holds $\sideset{}{{}_{l=1}^{L}}{\sum}G(x)_{l}\equiv 1$ and $G(x)_{i}\geq 0,\forall x\in\mathcal{X}$ and $1\leq i\leq L$.

We implement the approach of Shazeer et al. (2017) to model the gate $G$ as follows: For the softmax function

[TABLE]

where $\bm{a}\in\operatorname*{\mathds{R}}^{L}$, Shazeer et al. (2017) propose to set

[TABLE]

for an adequate permutation $(i_{1},\ldots,i_{L})$ of $\{1,\dots,L\}$ such that $h_{i_{l}}(x)>h_{i_{l+1}}(x)$ are ordered decreasingly. Here, $h_{l}(x)=g_{l}(x)+\mathcal{N}(0,\mathfrak{s}_{l}^{2})$ is a noisy version of single-channel gating models $g_{l}$ with parameters $\theta_{g_{l}}$. Note that these models can be chosen freely and may also deviate from the choice of expert models $\widehat{y}_{l}$.

The cutoff value $1\leq\kappa\leq 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 each expert a chance to become an element of the top-$\kappa$ components by adding independent Gaussian noise $\mathcal{N}(0,\mathfrak{s}_{l}^{2})$ before thresholding, where \mathfrak{s}\in\mathds{R}^{L}_{\raisebox{0.0pt}{\tiny++}} is another hyperparameter to set or learn. This noisy gating prevents a premature discarding of initially underperforming experts.

The overall MoE hyperparameter set is thus given by

[TABLE]

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 a GPR-based, model-agnostic superior sampling scheme in Sec. 4.2 (Fig. 2, C). Next, we 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 of our superior sampling scheme to problems of high input space dimensions. We then give details on the training of the MoE in Sec. 4.4 and finally propose an LFC-based IP selection method in Sec. 4.5. We summarize the pseudo-code of our superior sampling scheme in Algorithm 1.

4.1 GPR-based LFC and the superior training density

Let $\Sigma^{n}_{\text{GPR}}(x)$ denote the LOB function (4) of GPR. Inspired by the results to LFC and the superior training density of LPS in Eq. (7) and (6), we are able to deduce their GPR-based analog. Here, we need to take into account that GPR adapts universally333That is, the MISE decays at the minimax-rate $n^{-\frac{2\alpha}{2\alpha+d}}$ of nonparametric models. to functions $f\in\mathcal{C}^{\alpha}\left(\mathcal{X},\operatorname*{\mathds{R}}\right)$, as 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(x)n$ with the scaling result of the global $\Sigma^{n}_{\text{GPR}}$ in (9) from Sec. 3.3, we propose an LFC estimate for GPR as follows (see Appendix C for proof details).

Theorem 1** (LFC of GPR).**

For $f\in\mathcal{C}^{\alpha}\left(\mathcal{X},\operatorname*{\mathds{R}}\right)$, $\bm{X}_{n}\sim p$ and homoscedastic noise, the GPR-based LFC estimate of $f$ in $x\in\mathcal{X}$ is asymptotically given by

[TABLE]

In analogy to Eq. (7), $\mathfrak{C}^{n}_{\text{GPR}}$ measures the structural complexity of $f$, as it asymptotically does not depend on $p$, $v$ and $n$. Note that the LPS model provides no explicit way to adapt to the local noise variance $v(x)$, such that the LOB of LPS scales with respect to $v$ to address heteroscedasticity (see (36) in Appendix A). For GPR, we have made the restriction of homoscedastic noise in the definition of $\mathfrak{C}^{n}_{\text{GPR}}$ in Theorem 1, since we are not aware of a theory on the scaling of GPR-based LOB with respect to heteroscedasticity. However, as opposed to LPS, a heteroscedastic GPR model provides an explicit way to adapt to the local noise variance $v(x)$ 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, $\mathfrak{C}^{n}_{\text{GPR}}$ will be sufficiently calibrated in a heteroscedastic scenario, making it a reasonable estimate of LFC in practice without further restrictions.

Now, when putting $\mathfrak{C}^{n}_{\text{GPR}}$ into Eq. (6) with $Q=\alpha-1$, we obtain the superior training density

[TABLE]

For even $\alpha\in\operatorname*{\mathds{N}}$ with $Q=\alpha-1$ and $\mathfrak{C}^{n}_{\text{GPR}}\equiv\mathfrak{C}^{n}_{Q}$, $p^{\text{GPR},n}_{\text{\tiny{Sup}}}$ and $p^{Q,n}_{\text{\tiny{Opt}}}$ coincide, which proved to be optimal for LPS. In this sense, (20) generalizes (6) to the general case of \alpha\in\mathds{R}_{\raisebox{0.0pt}{\tiny+}}, where we expect that the true optimal training density for $f\in\mathcal{C}^{\alpha}\left(\mathcal{X}\right)$ will not deviate by a lot from $p^{\text{GPR},n}_{\text{\tiny{Sup}}}$. Since LPS and GPR are related models, we furthermore expect $\mathfrak{C}^{n}_{\text{GPR}}$ to be similar to $\mathfrak{C}^{n}_{Q}$ for the appropriate order $Q$.

Note that for $f\in\mathcal{C}^{\infty}\left(\mathcal{X}\right)$, we let $\alpha\rightarrow\infty$ in Eq. (19) and (20) to obtain

[TABLE]

While $\bm{X}_{n}\sim p^{\text{GPR},n}_{\text{\tiny{Sup}}}$ will not be optimal for our model, we expect it to be asymptotically superior to the naive random test sampling, i.e., $\bm{X}_{n}\sim q$, due to the model-agnosticity of the LPS-based result. To assess the asymptotic performance of a training density $p$ (such as $p^{\text{GPR},n}_{\text{\tiny{Sup}}}$), let us first observe the following:

For regression problems and under weak assumptions the law of the MISE does not change with respect to $p$, except for a constant multiple (Györfi et al., 2002; Willett et al., 2005). Accordingly, the number of actively selected training samples ($\sim p$) that are required to achieve the same level of accuracy of random test sampling is given by a constant \raisebox{1.0pt}{\varrho}>0. Formally, we can define $\varrho$ as follows.

Definition 2**.**

Over the space of square-integrable functions $f\in\mathcal{L}^{2}\left(\mathcal{X}\right)$, for a nonparametric regression model $\widehat{f}$ and a training density $p$, we define by \raisebox{1.0pt}{\varrho}(\widehat{f},p)>0 the relative required sample size such that for n^{\prime}=\raisebox{1.0pt}{\varrho}(\widehat{f},p)n, $\bm{X}_{n^{\prime}}^{\prime}\sim p$ and $\bm{X}_{n}\sim q$ it holds that

[TABLE]

Thus, a training density $p$ is asymptotically superior to random test sampling, if \raisebox{1.0pt}{\varrho}(\widehat{f},p)<1, since we achieve the same performance as random test sampling with only a fraction of the number of training samples. In Sec. 5 we will demonstrate the superiority of the training density $p^{\text{GPR},n}_{\text{\tiny{Sup}}}$ for our GPR-based MoE model.

Respecting the intrinsic dimension in high-dimensional input spaces   In Eq. (19), (20) and (21) we assume the input space $\mathcal{X}$ to have full degrees of freedom $d$, which in practice is particularly not the case in high-dimensional feature spaces. For an intrinsic dimension $\delta<d$ of $\mathcal{X}$, we adjust as follows: For the space $\mathcal{S}=\left\{\sigma\Sigma\;\middle|\;\sigma>0\right\}$ of bandwidth candidates that are essentially isotropic up to a fixed, shared positive definite factor \Sigma\in\mathbb{S}^{d}_{\raisebox{0.0pt}{\tiny++}} that is, e.g., calculated in a pre-processing step, let $\Sigma^{n}_{\text{GPR}}(x)=\sigma^{n}_{\text{GPR}}(x)\Sigma$ be the LOB function of GPR with respect to $\mathcal{S}$. Then we replace all occurrences of $d$ for $\delta$ and $\left|\Sigma^{n}_{\text{GPR}}(x)\right|$ for $\sigma^{n}_{\text{GPR}}(x)^{\delta}$ in Eq. (19), (20) and (21).

Besides being an ingredient to AL, LFC can also be used to reduce the required model complexity. For example, in an RBF-network or a sparse GPR model, we can coarsen or refine the model resolution by placing an adequate amount of basis functions or IPs, respecting LFC. We will discuss this choice in Sec. 4.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 high-dimensional data (see Sec. 5.2).

4.2 The active learning framework

Starting with an initial training set $\bm{X}_{n_{0}},\bm{Y}_{n_{0}}$ of size $n_{0}$ with $\bm{X}_{n_{0}}\sim p_{0}$ for some initial training distribution such as $p_{0}\equiv q$, we implement the online sampling procedure as described in Panknin et al. (2021), such that $\bm{X}_{n}\sim p^{\text{GPR},n}_{\text{\tiny{Sup}}}$ as $n\rightarrow\infty$. We grow the training set as follows:

Given the current training set $\bm{X}_{n_{k}},\bm{Y}_{n_{k}}$ we estimate $\widehat{\Sigma}^{n_{k}}_{\text{GPR}}$ as described in Sec. 4.3. Using (19), (20), it is

[TABLE]

Letting the next sample size be $n_{k+1}=2n_{k}$, we have already drawn half the samples of $n_{k+1}$ according to a potentially different distribution $p_{k}$ than the new proposed $\widehat{p}^{\text{GPR},n_{k}}_{\text{\tiny{Sup}}}$. The closest we can get in distribution to $\widehat{p}^{\text{GPR},n_{k}}_{\text{\tiny{Sup}}}$ is given by $\displaystyle\bm{X}_{n_{k+1}}\!\!\!\sim p_{k+1}$, where $p_{k+1}:=\gamma_{2}p_{k}+(1-\gamma_{2})\widehat{p}^{\text{GPR},n_{k}}_{\text{\tiny{Sup}}}$, for $\displaystyle\gamma_{2}=\max\left\{0,\frac{0.5-\gamma_{1}^{-1}}{1-\gamma_{1}^{-1}}\right\}\in[0,0.5)$ and $\displaystyle\gamma_{1}=\max_{x\in\mathcal{X}}\frac{p_{k}(x)}{\widehat{p}^{\text{GPR},n_{k}}_{\text{\tiny{Sup}}}(x)}$. This is achieved by sampling $x_{n_{k}+1},\ldots,x_{n_{k+1}}\sim\widetilde{p}_{k+1}$ for $\widetilde{p}_{k+1}=2p_{k+1}-p_{k}$, 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 $\mathcal{X}$ is restricted to a large, unlabeled pool of samples $\bm{X}_{pool}\in\mathcal{X}^{N}$. This pool follows a distribution $\bm{X}_{pool}\sim p_{\mathcal{X}}$, for which it is common to assume an (unnormalized) density estimate $\widehat{p}_{\mathcal{X}}$ to be given: Unlabeled inputs are considered cheaply accessible, whereas querying labels is expensive.

For our AL framework to be applicable, it suffices to keep all considered densities such as $\widehat{p}^{\text{GPR},n_{k}}_{\text{\tiny{Sup}}}$ at equal norm, which we can enforce via normalizing a density $p$ by $\bar{p}=p/\text{norm}(p)$, where

[TABLE]

To see this, note that first of all $\widehat{p}_{\mathcal{X}}$ is an unnormalized estimate of $p_{\mathcal{X}}$ such that we can write $\widehat{p}_{\mathcal{X}}\approx c\cdot p_{\mathcal{X}}$ for some unknown constant $c>0$. On the one hand, it is $\displaystyle\int_{\mathcal{X}}\widehat{p}_{\mathcal{X}}(x)dx=c$ by definition. On the other hand, it is

[TABLE]

such that also $\displaystyle\int_{\mathcal{X}}\bar{p}(x)dx=\frac{1}{\text{norm}(p)}\int_{\mathcal{X}}p(x)dx\approx c$ holds for any unnormalized density $p$.

Subsequently, the required samples $x_{n_{k}+1},\ldots,x_{n_{k+1}}\sim\widetilde{p}_{k+1}$ are obtained via importance sampling from the pool with importance weights $\displaystyle\operatorname*{\mathds{P}}(x_{i}=x)\propto\widetilde{p}_{k+1}(x)/[\text{norm}(\widetilde{p}_{k+1})\widehat{p}_{\mathcal{X}}(x)]$ for $x\in\bm{X}_{pool}$ and $n_{k}+1\leq i\leq n_{k+1}$.

4.3 Sparse mixture of Gaussian processes

Recall the sparse MoE model (15) from Sec. 3.4.2, given by

[TABLE]

where the gate $G\!:\mathcal{X}\,\rightarrow\,[0,1]^{L}$ is a probability assignment of an input $x$ to the expert models $\widehat{y}_{l}$. In particular, it holds $\sideset{}{{}_{l=1}^{L}}{\sum}G(x)_{l}\equiv 1$ and $G(x)_{i}\geq 0,\forall x\in\mathcal{X}$ and $1\leq i\leq L$. According to (18), besides the expert and gate model parameters $\{\theta_{e_{l}}\}_{l=1}^{L}$ and $\{\theta_{g_{l}}\}_{l=1}^{L}$, this MoE approach has two hyperparameters, $\kappa$ and $\mathfrak{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 as well as the single channel gating models to be sparse variational GPs (Williams & Rasmussen, 1996; Hensman et al., 2015), that is, $\widehat{y}_{l}\sim\mathcal{SVGP}(\theta_{e_{l}})$ and $g_{l}\sim\mathcal{SVGP}(\theta_{g_{l}})$, which are parameterized by $\theta_{e_{l}}$ and $\theta_{g_{l}}$, as described in Sec. 3.4.1. The overall MoE hyperparameter set is thus given by $\Theta=(\{\theta_{e_{l}}\}_{l=1}^{L},\{\theta_{g_{l}}\}_{l=1}^{L},\kappa,\mathfrak{s})$.

We will keep certain hyperparameters of $\Theta$ constant after initialization, and share some hyperparameters across experts and the channels of the gate: While the covariances of the inducing value distributions $\bm{S}_{\dagger}\in\theta,\theta\in\Theta$ could be full positive definite matrices, we apply $\bm{S}_{\dagger}=0$ throughout, giving favorable stability and computational efficiency. For the same reasons, we fix the inducing point (IP) locations $\bm{X}_{\dagger}\in\theta,\theta\in\Theta$ after initialization. Furthermore, we share the IP locations among the experts, respectively the gate channels, such that for $\bm{X}_{\dagger}\in\theta_{e_{l}}$ we apply $\bm{X}_{\dagger}=\bm{X}^{E}_{\dagger}$ and for $\bm{X}_{\dagger}\in\theta_{g_{l}}$ we apply $\bm{X}_{\dagger}=\bm{X}^{G}_{\dagger}$, for all $1\leq l\leq L$.

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 $\mu_{E}$, the regularization parameter $\lambda_{E}$, and the noise variance function $\widehat{v}$, respectively the global noise variance $\sigma_{\varepsilon}^{2}$ with $\widehat{v}(x)\equiv\sigma_{\varepsilon}^{2}$ in case of homoscedasticity. Furthermore, we apply a fixed, logarithmically spaced set of individual expert bandwidth scaling factors $\sigma_{1}<\ldots<\sigma_{L}$ that are multiplied by a fixed, shared bandwidth matrix $\Sigma_{E}$. Our expert parameters therefore reduce to

[TABLE]

Remark 3**.**

Recall from Sec. 2 that it is possible to replace the variational GPR expert models for full as well as sparse analytic GPR formulations (see Appendix B). With slight abuse of notation, these cases are subsumed by setting $\bm{\mu}^{e_{l}}_{\dagger}=\emptyset$ or $\bm{X}^{E}_{\dagger}=\bm{\mu}^{e_{l}}_{\dagger}=\emptyset$ 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 $\widehat{v}\equiv 0$ for $\widehat{v}\in\theta_{g_{l}}$ and all $1\leq l\leq L$. Each output channel of the gate poses its own classification 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

[TABLE]

After training as described in Sec. 4.4, this MoE can cope with a varying structural complexity through the individual bandwidth scaling factors $\sigma_{l}$ of the experts and heteroscedastic 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 $\Sigma^{n}_{\text{GPR}}(x)$ as

[TABLE]

Due to the finite candidate set $\sigma_{1},\ldots,\sigma_{L}$ we are limited to measure a quantization of $\Sigma^{n}_{\text{GPR}}(x)$ through $G(x)_{l}=\operatorname*{\mathds{P}}(\Sigma^{n}_{\text{GPR}}(x)=\sigma_{l}\Sigma_{E})$. If, in fact, $\Sigma^{n}_{\text{GPR}}\in\{\sigma_{1}\Sigma_{E},\dots,\sigma_{L}\Sigma_{E}\}$ holds true, then there exists an index function $j(x)\in\{1,\dots,L\}$ such that $\Sigma^{n}_{\text{GPR}}(x)=\sigma_{j(x)}\Sigma_{E}$. In this case, we are able to exactly recover LOB with $G(x)_{l}=\begin{cases}1,&l=j(x)\\ 0,&\text{else}\end{cases}$. In any other case, the estimate (24) of LOB is a reasonable interpolation, which deviation from $\Sigma^{n}_{\text{GPR}}$ 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.4.1 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 systematically on the initial training dataset in Sec. 4.4.3.

4.4.1 The training objective

First, we will set up the objective function for training our MoE model in batch mode.

The main objective   Denote by $\varnothing\subsetneq\mathcal{B}\subseteq\{1,\dots,n\}$ the indices of a batch, and let $w_{\mathcal{B}}=\sideset{}{{}_{b\in\mathcal{B}}}{\sum}w(x_{b})$ for the training importance weight function $w\propto q/p$. Let $P_{l}$ be the prior distribution of the inducing function values of the l-th expert and $Q_{l}$ the corresponding variational distribution as defined in Sec. 3.4.1. We choose the (through the gate G) weighted sum of the individual expert negative ELBO objectives (13), denoted by

[TABLE]

as our main objective, where $n_{l}=nw_{\mathcal{B}}/\nu_{\mathcal{B},l}$ for $\nu_{\mathcal{B},l}=\sideset{}{{}_{b\in\mathcal{B}}}{\sum}v_{l}(x_{b})$ with $v_{l}(x)=G(x)_{l}w(x)$ and

[TABLE]

is the predictive log-likelihood (14) of the l-th expert in $x_{b}$, marginalized over its variational distribution $Q_{l}$.

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 to enforce this, we penalize smaller bandwidth choices by adding the following term:

[TABLE]

Note that $\text{pen}_{\sigma}(\bm{X}_{n},\bm{Y}_{n},\mathcal{B},w,\Theta)=1$ if $\nu_{\mathcal{B},1}=\ldots=\nu_{\mathcal{B},L}$. Our total objective then amounts to

[TABLE]

Remark 4**.**

If we assume our problem to be (almost) noise-free, we replace the Obj in our objective (26) for the mean squared error

[TABLE]

4.4.2 Training procedure

We implement our model in PyTorch (Paszke et al., 2019), using the GPyTorch-package (Gardner et al., 2018). Given the training set $(\bm{X}_{n},\bm{Y}_{n})$, we minimize the objective described in Sec. 4.4.1 via ADAM-optimization (Kingma & 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 hyperparameters, $\kappa$ for enforcing sparse gate responses and a noise term on the gate responses during the training, which is controlled by $\mathfrak{s}$. Instead of learning $\mathfrak{s}$ as a parameter during the training—like proposed by Shazeer et al. (2017)—we propose to shrink $\mathfrak{s}\leftarrow\mathfrak{s}\eta_{\mathfrak{s}}$ after each training epoch, for a multiplicative factor $\eta_{\mathfrak{s}}<1$ and an initial value $\mathfrak{s}:=\mathfrak{s}_{0}$ 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 $\eta$, where we shrink $\eta\leftarrow\eta_{i}:=\frac{1}{2}\eta_{i-1}$ for an initial base learning rate $\eta_{0}$ during the training as soon as the validation performance gets stuck until $\eta_{i}$ crosses a lower threshold, e.g., $\eta_{i}<\eta_{0}/1000$. Note that, relative to the base learning rate, good learning rates for the individual types of tunable parameters should be deployed: Within a GP component $\mathcal{SVGP}(\theta)$ with $\theta=(\mu,\lambda,\widehat{v},\Sigma,\bm{X}_{\dagger},\bm{\mu}_{\dagger},\bm{S}_{\dagger})$, the hyperparameters $(\mu,\lambda,\widehat{v})$ must be updated on a smaller scale than the inducing value distribution, given by $(\bm{\mu}_{\dagger},\bm{S}_{\dagger})$. In this regard, let $\eta_{H}\leq 1$ be the factor such that, if we update $\bm{\mu}_{\dagger}$ at rate $\eta$, then we update $(\mu,\lambda,\widehat{v})$ at rate $\eta_{H}\eta$.

Similarly, we need to update the gate parameters $\theta_{g_{l}}$ on a smaller scale than the expert parameters $\theta_{e_{l}}$. In this regard, let $\eta_{G}\leq 1$ be the factor such that, if we update $\theta_{e_{l}}$ at rate $\eta$, then we update $\theta_{g_{l}}$ at rate $\eta_{G}\eta$.

The set of hyperparameters that require off-training selection (e.g., via cross-validation) is thus given by

[TABLE]

whereas the overall set of parameters that get tuned while training is given by

[TABLE]

We provide further details on 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 et al., 2017) we can train our model using well-established software libraries (Kingma & Ba, 2015; Paszke et al., 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.

Note that some hyperparameters impact the computational complexity rather than the model performance. Thus, as long as they are not underestimated, their tuning is optional and will therefore be postponed:

• •

Since our MoE is robust concerning unnecessarily large choices of the gate output sparsity $\kappa$, we initialize $\kappa\equiv L$ while choosing the remaining hyperparameters, followed by tuning $\kappa$ as the last hyperparameter, where we successively reduce $\kappa$ until we observe a significant loss of performance of the MoE.

• •

The numbers $m_{E}=\left|\bm{X}^{E}_{\dagger}\right|$, $m_{G}=\left|\bm{X}^{G}_{\dagger}\right|$ 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 to the smallest value that leads to no significant loss of performance to keep the computational complexity of the model moderate at larger training sizes. In the initial iteration, we use $m_{E}=n_{0}$ for the experts and $m_{G}=\frac{n_{0}}{4}$, where the locations of the IPs $\bm{X}^{E}_{\dagger},\bm{X}^{G}_{\dagger}\sim p$ are chosen diverse as described in Appendix E. In the second iteration, where we have first estimates of LFC, $\bm{X}^{E}_{\dagger}$ and $\bm{X}^{G}_{\dagger}$ are drawn as discussed in Sec. 4.5. Here, we gradually decrease $m_{E}$ and $m_{G}$ until we observe a significant loss of validation performance of the MoE. We hold $m_{E}$ and $m_{G}$ 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 ($\eta_{0},\eta_{H}$) and the batch size $B$ that are related to a single GPR expert rather than the whole MoE model:

• •

First, we hold $\eta_{H}=0.2$, $B=n_{0}$ fixed at reasonable initial values and perform line search over $\eta_{0}$ according to the resulting validation performance of a single, isotropic, sparse GPR expert. Here, too small values of $\eta_{0}$ result in slow convergence of the objective, in which case we interrupt the training immediately and increase $\eta_{0}$ as long as the first objective updates are consistently decreasing.

• •

Next, we choose $\eta_{H}$ according to the resulting validation performance of a single, isotropic, sparse GPR expert, where we gradually decrease $\eta_{H}$, starting from $\eta_{H}=1$. Again, we interrupt the training for too large choices of $\eta_{H}$, where the training objective will diverge.

• •

Finally, we choose $B$ according to the resulting validation performance of a single, isotropic, sparse GPR expert, where we gradually decrease $B$, starting from $B=n_{0}$.

Next, we have to deal with the remaining hyperparameters that are related to the MoE model. For this, we first initialize the less crucial hyperparameters at reasonable values, tuning them afterward: We apply a set of experts with $\{\sigma_{l}\}_{l=1}^{L}$, where $L=7$, and $\sigma_{l}=2^{\frac{l-4}{\delta}}$, which are logarithmically spaced around the global bandwidth estimate $\Sigma_{E}$ 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 ($\mathfrak{s}_{0},\eta_{\mathfrak{s}}$) as well as the regularization $\vartheta_{\sigma}$ of the bandwidth function are about fine-tuning of the model. We set them to $\mathfrak{s}_{0}=0.1,\eta_{\mathfrak{s}}=1/\sqrt{2}$ and $\vartheta_{\sigma}=0.01$. We suggest to keep $\eta_{\mathfrak{s}}=1/\sqrt{2}$ throughout without further tuning.

• •

Via grid search, we choose $\sigma_{G},\lambda_{G}$ according to the resulting validation performance of the MoE model, where we gradually decrease the gate learning rate from $\eta_{G}=1$.

As a last step, we choose the hyperparameters for fine-tuning of the MoE model:

• •

First, we perform line search over $\vartheta_{\sigma}$ according to the resulting validation performance of the MoE.

• •

Next, we tune $\{\sigma_{l}\}_{l=1}^{L}$: We observe that unreasonable bandwidths will be automatically dropped during the training. Therefore, if the minimal or maximal candidate associated with $\sigma_{1},\sigma_{L}$ is not chosen during the training, we remove the respective expert and retrain the MoE. Vice versa, we expand the bandwidth candidate range beyond $\sigma_{1},\sigma_{L}$ with a factor of $2^{\mp\frac{1}{\delta}}$ as long as the boundary candidates are not dropped during the training.

• •

Finally, we perform line search over $\mathfrak{s}_{0}$ according to the resulting validation performance of the MoE.

4.5 Initializing the IP locations

Since the numbers of IPs of the gate $\bm{X}^{G}_{\dagger}$, as well as the experts $\bm{X}^{E}_{\dagger}$ 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-dispersion) so 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:

Definition 5**.**

The dispersion, given by $\sup_{x\in\mathcal{X}}\min_{1\leq i\leq n}\|x-x_{i}\|_{\raisebox{-2.0pt}{\tiny$$}}$ (Niederreiter, 1988) is a measure of how well spread out the training sample is. We say that a sequence has low-dispersion if its dispersion is lower than the dispersion of random uniform sampling.

Indeed, by sampling the IPs in this manner, we can reduce the distance of an evaluation point $x$ to its closest neighbor in $\bm{X}^{E}_{\dagger}$, which is known to reduce the reconstruction error of a full kernel matrix by a 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 summary, we therefore propose to choose the IPs

[TABLE]

as the geometric mean of LFC and the training density $\bm{X}_{n}\sim p$ in a diverse way. Here, we ensure diversity by implementing distribution preserving clustering or particle repulsion as described in Appendix E.

5 Experiments

In this section, we will first analyze our approach on toy-data, regarding the MoE model, LFC, and the superior training density. Then, we apply our approach to a high-dimensional MD simulation dataset from quantum chemistry, by which we can deduce deeper insights into this regression problem.

We denote the root mean squared error (RMSE) and the maximum absolute error (max​ AE) of a model $\widehat{f}$ for a test set $\bm{X}_{\text{T}}\in\mathcal{X}^{N}$ with $\bm{X}_{\text{T}}\sim q$ by

[TABLE]

As already discussed in Sec. 4.1, the learning rate is invariant under change of the training density $\bm{X}_{n}\sim p$ in the considered scenario. For our MoE model and $f\in\mathcal{C}^{\alpha}\left(\mathcal{X},\operatorname*{\mathds{R}}\right)$ we, thus, can write

[TABLE]

where $C_{p}>0$ is a constant depending on the training density $p$. Note that we can theoretically bound the asymptotic RMSE from below by $C^{*}n^{-\tau}$, where we have defined $C^{*}:=C_{p^{*}}$ with $\bm{X}_{n}^{\prime}\sim p^{*}$ being the optimal training set from (3). Unfortunately, since $p^{*}$ is unknown—even when given the ground truth—we are not able to estimate $C^{*}$ and, thus, provide a lower bound of the RMSE beyond the known learning rate $n^{-\tau}$.

As an AL performance measure, we use the relative required sample size from Definition 2 which can be estimated for a GPR-based model such as our MoE and $f\in\mathcal{C}^{\alpha}\left(\mathcal{X},\operatorname*{\mathds{R}}\right)$ according to

[TABLE]

where it is $\bm{X}_{n}^{\prime}\sim p$ and $\bm{X}_{n}\sim q$ with respective labels $\bm{Y}_{n}^{\prime}$ and $\bm{Y}_{n}$. Using (31), we can compare the asymptotic AL performance of different AL sampling schemes in the following experiments. For example, we can quantify the AL performance of our proposed AL framework by sampling $\bm{X}_{n}^{\prime}\sim\widehat{p}^{\text{GPR},n}_{\text{\tiny{Sup}}}$.

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 et al., 2013). For $x\in\mathcal{X}=[0,1]$, let

[TABLE]

where $\epsilon=0.05$, $C$ is chosen such that \|f\|_{\raisebox{-2.0pt}{\tiny2}}=7 and $\mathcal{N}(\cdot;\mu,\sigma^{2})$ denotes the Gaussian distribution with mean $\mu$ and variance $\sigma^{2}$. We assume a uniform test distribution $q\sim\mathcal{U}(\mathcal{X})$ in all Doppler function experiments.

This one-dimensional, homoscedastic toy-example allows for an easy and intuitive visualization. Fig. 3 shows an example dataset as blue dots and the true function $f$ to infer in black. Due to the strong variation of structural complexity, a single-scale GPR model does not cope well with the Doppler function (see Appendix G.1 for a comparison of single-scale to multi-scale GPR).

We implement our proposed MoE model as 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 $\sigma_{j}=10^{(j-10)/3},1\leq j\leq 7$, as the expert bandwidths, $\lambda_{E}=20$ as the initial expert regularization, and $\sigma_{G}=0.05$ and $\lambda_{G}=10$ for the gate. For the training, we apply a batch size of $B=512$, a terminal expert sparsity $\kappa=2$, a penalty factor of $\vartheta_{\sigma}=0.5$ for small bandwidth choices, gate noise parameters $\mathfrak{s}_{0}=0.1$ and $\eta_{\mathfrak{s}}=1/\sqrt{2}$, and learning rate parameters $\eta=0.01$, $\eta_{H}=0.2$, $\eta_{G}=1$.

Fig. 4 shows the gate function after training of the MoE model, as described in Sec. 4.4.1, 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\in\mathcal{C}^{\infty}\left(\mathcal{X},\operatorname*{\mathds{R}}\right)$, our deduced superior training density estimate is given by Eq. (21). In Fig. 5 we plot our estimates of LOB and the superior training density in comparison to the LPS-based results for polynomial degrees of order $Q=1,3$, and with implementation and hyperparameters as described in Panknin et al. (2021). 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 in Sec. 4.2, we additionally observe quantitative benefits in Fig. 6 over random test sampling—quite similar to the LPS-based result for $Q=3$: When estimating the relative sample size (31) we require to achieve the same RMSE via active sampling compared to random test sampling, we obtain \textstyle\raisebox{1.0pt}{\varrho}(\widehat{f}_{\text{\tiny{MoE}}},\widehat{p}^{\text{GPR},n}_{\text{\tiny{Sup}}})=0.58\pm 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 $\mathcal{X}=[0,1]$ we mean the deterministic construction of the training inputs, where for $n$ given training samples the subsequent $n$ training inputs get placed halfway between all nearest neighbors of the former $n$ 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, indeed, equidistant sampling is superior to random test sampling at small training sizes. As claimed, however, with growing training size, this advantage gradually diminishes until it has vanished completely at $n=2^{15}$ 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 was obtained for $2^{15}$ training samples and which we now keep fixed—we define the uncertainty estimate of our model MoGPU as a straightforward 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 $x$ with respect to the gate values $G(x)$ from (15), we derive

[TABLE]

where $\bm{C}_{\theta_{l}}^{*}$ is the predictive variance of the l-th expert (see (12)). Note that we consider MoGPU as a baseline competitor to our superior sampling scheme.

Intuitively, the uncertainty estimate in $x\in\mathcal{X}$ increases as the applied bandwidth $\sigma^{n}_{\text{GPR}}(x)$ decreases. Now, in order to equalize uncertainty over the input space, MoGPU will sample more in regions where $\sigma^{n}_{\text{GPR}}$ is smaller. For $\bm{X}_{n}$ drawn according to MoGPU, we expect $\bm{X}_{n}\sim[\sigma^{n}_{\text{GPR}}]^{-d}$. This expectation holds as can be seen at the top in Fig. 8.

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). In all error measures the beneficial effect of the low-dispersion property of $\bm{X}_{n}$ drawn according to MoGPU has already vanished for about $512$ training samples, from where the asymptotic law dominates. As expected, our approach is superior to MoGPU when comparing RMSE. Interestingly, MoGPU is superior to our approach regarding the max​ AE, suggesting that $\bm{X}_{n}\sim[\sigma^{n}_{\text{GPR}}]^{-d}$ is the preferable training distribution under the supremum-norm.

AL** Performance on deep Gaussian processes**   To demonstrate the model-agnosticity of our AL approach, we deploy the DGP model of Sauer et al. (2023b) using the CRAN package deepgp444See https://cran.r-project.org/web/packages/deepgp/index.html. This package also implements an aggregate variance-based AL criterion (Cohn, 1994), which they named active learning Cohn (ALC) after the originator. We deploy a 3-layer DGP model, using the Gaussian kernel. For test evaluation, We train the model using Vecchia-approximation (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 equidistant 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. 9. 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, we 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 $\bm{X}_{n}\sim p^{\text{GPR},n}_{\text{\tiny{Sup}}}$ can be performed in batch mode. This result emphasizes the complementary nature of our asymptotic work to the classic bottom-up AL literature.

Necessity of the small bandwidth penalty   We impose a penalty on small bandwidth choices through the factor $\vartheta_{\sigma}=0.5$ to regularize the bandwidth function and prevent overfitting, as described in Appendix F. We demonstrate this overfitting issue in Appendix G.2 that results from applying no regularization ($\vartheta_{\sigma}=0$).

Parsimonious modeling using LFC   In Sec. 4.1 we mentioned that LFC can also be used to coarsen or refine the model resolution adequately to 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 $m=1024$ IPs under active, respectively random test sampling, we here investigate the influence of the number of IPs and their distribution on 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\ll 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.

In Fig. 10, we compare the RMSE for the fixed training size $n=2^{15}$ for both settings, active and passive, when sampling the IPs according to the training density $p$, the LFC and their geometric mean (29). First of all, we observe that we generally require less IPs with active sampling compared to random test sampling, which originates from the fact that the superior training density $\widehat{p}^{\text{GPR},n}_{\text{\tiny{Sup}}}$ already 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, non-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 a factor of $8$ for the Doppler function without a significant loss of performance, when respecting LFC in the model design.

Comparing our proposed IP selection method to a greedy fast forward selection   In Sec. 2, we discussed other IP selection approaches. For comparison, we have implemented the greedy fast forward (GFF) IP selection method of Seeger 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 $\mathcal{GP}(\theta)$ distribution. Here, the information of an IP candidate $x_{i}\in\bm{X}_{n}\setminus\bm{X}_{\dagger}$ is measures by

[TABLE]

which is the Kullback-Leibler divergence between the posterior distributions based on the IPs $\bm{X}_{\dagger}\cup\{x_{i}\}$ and $\bm{X}_{\dagger}$. Accordingly, the updated set of IPs is given by $\bm{X}_{\dagger}\leftarrow\bm{X}_{\dagger}\cup\{x_{i}^{*}\}$, where

[TABLE]

The procedure converges, when the remaining IP candidates carry no further information, that is, $J(x_{i}^{*})<\varepsilon_{\mathcal{J}}$, up to a specified threshold $\varepsilon_{\mathcal{J}}\geq 0$.

At a given threshold $\varepsilon_{\mathcal{J}}$, we observe that the number of selected IPs is very small for the experts with a large bandwidth, while it increases drastically ($\propto\sigma_{i}^{-1}$) for experts with a small bandwidth. Now that the overall complexity of the MoE is dominated by the expert with the most IPs, it is fair to compare the number of IPs of the expert at bandwidth $\sigma_{1}$ with our statically specified number $m$ of IPs in Fig. 10. Here, we will vary the threshold $\varepsilon_{\mathcal{J}}$ to obtain a 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, as opposed to our proposed IP selection method.

In any case, even under training according to $p^{\text{GPR},n}_{\text{\tiny{Sup}}}$, the selected IPs by Seeger et al. (2003) are uniformly distributed. In particular, the need for less IPs of the experts at smaller bandwidths to the right of $\mathcal{X}$ is not recognized and, thus, we observe no IP savings at an acceptable performance over a random IP selection at all for this inhomogeneously complex problem.

Heteroscedastic noise treatment   While the treatment of heteroscedastic 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(x)=(3-4\left|x-0.5\right|))^{2}\in[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. 4.1, the influence of $v$ on the LOB estimates obtained from heteroscedastic GPR experts is relatively small. Likewise, we proceed with the evaluation of our proposed AL scheme under heteroscedasticity. Here, we compare to random test sampling but also to $p^{\text{GPR},n}_{\text{\tiny{Sup}}}$ under 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 to materialize later from $n=2^{13}$ 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.

5.2 Force field reconstruction

We now turn our attention to a real-world example in which we predict the potential energy surface (PES) and corresponding force field (FF) of a molecule from first-principles calculations. The PES function links the geometry $x=\left[R_{1},\ldots,R_{\bm{a}}\right]\in\operatorname*{\mathds{R}}^{3\times\bm{a}}$ of a molecule to its potential energy $E\in\operatorname*{\mathds{R}}$, where $R_{i}$ are the Cartesian positions of the $\bm{a}$ atoms of the molecule. In ab initio computations, this mapping is achieved by solving the time-independent 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 configurations. 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 equations of motion. To this end, energy-conserving forces acting on each atom are required. These forces are the negative derivative of the PES with respect to the atomic positions $F\in\operatorname*{\mathds{R}}^{3\times\bm{a}}$.

This type of proxy for the prohibitively expensive ab initio quantum mechanical calculations is commonly used to enable 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 structure-function relationships of physical systems at atomic scale. They are widely used in molecular biology research and play a 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 in the way of scaling to larger system sizes, both, from a data acquisition cost and training perspective. Here, we demonstrate how our method can be used to construct smaller, yet more effective training datasets.

In this experiment, we reconstruct a FF for the molecule malonaldehyde, which has $\bm{a}=9$ atoms and the chemical formula $C_{3}H_{4}O_{2}$ (see Fig. 13 (A)). Formally, we try to infer the high-dimensional target function $f\!:\mathcal{X}\,\rightarrow\,\mathcal{Y},R\mapsto[E,F]$, where $\mathcal{X}=\operatorname*{\mathds{R}}^{3\bm{a}}$ and $\mathcal{Y}=\operatorname*{\mathds{R}}^{1+3\bm{a}}$. 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 (aldehyde groups) (Chmiela et al., 2018; Sauceda 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 degrees of freedom to obtain a minimal energy configuration. We will refer to these geometries as the relaxed configurations in the following.

To reconstruct the FF, we consider the broadly adopted symmetric gradient-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 to 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 on this, we derive the superior training density, which we finally input into our AL framework to refine the training dataset in a superior way.

Experimental setup   All experiments use an extensive pre-computed reference trajectory (almost a million data points ($\bm{X}_{pool},\bm{Y}_{pool}$)) as ground truth, as opposed to generating new data points on demand. This test setup allows a 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 $\bm{X}_{pool}\sim p_{\mathcal{X}}$ since we are dealing with a pool-based AL scenario. We estimate $\widehat{p}_{\mathcal{X}}$ by standard kernel density estimation, based on the energy-to-energy entry of the sGDML kernel $\widetilde{\bm{k}}$ from (46) at $\sigma=0.03$. Fig. 13 (C) shows the density estimate of the relaxed configurations, where we observe that $p_{\mathcal{X}}$ is very unbalanced, with a strong concentration of mass near the stable configurations.

We implement our MoE approach, using the sGDML kernel $\widetilde{\bm{k}}$ 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 as ground truth), we tune the experts (and MoE model) with respect to MSE rather than the Obj objective. For stability, we will apply $\widehat{v}(x)=10^{-9}$ even though we assume no noise.

Anisotropic bandwidths   sGDML operates on $\bm{d}=\bm{a}(\bm{a}-1)/2=36$ features that are based on the interatomic distances of the molecule. In contrast to the work of Chmiela et al. who restrict themselves to an isotropic bandwidth $\Sigma_{E}=\sigma_{E}\mathcal{I}_{\bm{d}}$, our implementation of sGDML in GPyTorch naturally enables us to tune an anisotropic bandwidth $\Sigma_{E}=\textbf{{diag}}(\sigma_{1},\ldots,\sigma_{\bm{d}})$ in the preprocessing step.

We partially offset the increased memory footprint of the model due to the tunable $\Sigma_{E}$ by implementing the sparse GPR model from Sec. 3.4.1 under the sGDML kernel $\widetilde{\bm{k}}$ from (46) and limiting the number of IPs to $m=128$ configurations. Since all our features are of the same type—pairwise interatomic distances—they are inherently calibrated in terms of scale. Hence, the reciprocal entries of $\Sigma_{E}$ directly translate into the importance of the features, which we display in Fig. 14.

We observe, that the importance assigned to some pairs 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 we 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 $\Sigma_{E}$, we apply dense sGDML experts with $\Sigma_{j}=\sigma_{j}\Sigma_{E}$, where $\sigma_{j}=2^{-5/4+j/2},1\leq j\leq 8$ as the individual expert bandwidths, $\lambda_{E}=1$ as the initial expert regularization, and $\sigma_{G}=0.1$ and $\lambda_{G}=10^{4}$ for the sparse gate with $1024$ IPs. For the training, we apply a batch size of $B=1024$, a terminal expert sparsity $\kappa=8$, a penalty factor of $\vartheta_{\sigma}=0.01$ for small bandwidth choices, gate noise parameters $\mathfrak{s}_{0}=0.01$ and $\eta_{\mathfrak{s}}=1/\sqrt{2}$, and learning rate parameters $\eta=0.005$, $\eta_{H}=0.05$, $\eta_{G}=0.1$. As we discuss in Appendix F, for tuning the MoE with dense (sGDML) experts, we either require an additional gate training set, which is independent of the training set for the experts, or we could provide leave-one-out (LOO) responses of the experts for the training of the gate. In our experiment, we use an additional gate training set $\bm{X}_{n_{G}}^{G}$ of fixed size $n_{G}=2^{14}$. The anisotropic 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.

Active learning   We assume an intrinsic dimension of $\delta=2$ (the two aldehyde rotor angles, the most salient features of malonaldehyde) and a smooth target function $f\in\mathcal{C}^{\infty}\left(\mathcal{X},\mathcal{Y}\right)$. The test distribution is given by the MD trajectory such that $q=p_{\mathcal{X}}$. Prior to the AL procedure, we separate the validation samples $\bm{X}_{\text{val}}$ and test samples $\bm{X}_{\text{T}}$ at random from the pool $\bm{X}_{pool}$. We apply an initial expert training size of $n_{0}=2^{9}$, doubling the sample size with each iteration of the AL procedure. The initial expert training set $\bm{X}_{n_{0}}$ and the gate training set $\bm{X}_{n_{G}}^{G}$ are drawn via importance sampling from the remaining pool with weights $\widehat{p}_{\mathcal{X}}^{-1/2}(\bm{X}_{pool}\setminus(\bm{X}_{\text{val}}\cup\bm{X}_{\text{T}}))$. By this it is $\bm{X}_{n_{0}}\sim q^{1/2}$, which is more in alignment with the superior training density (21) than sampling $\bm{X}_{n_{0}}\sim q$.

In Fig. 13 (D, E) we 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 \textstyle\raisebox{1.0pt}{\varrho}(\widehat{f}_{\text{\tiny{MoE}}},\widehat{p}^{\text{GPR},n}_{\text{\tiny{Sup}}})=0.920\pm 0.013. This means that we save about 8% of samples under the MoE model with our active sampling scheme compared to 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 et 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 practical importance.

6 Discussion

Active learning   Recall that in this work we have restricted ourselves to the scenario of model-agnostic AL with persistent performance at large training size. This scenario is complementary to the more common small sample size regime. And while both cases are important, a lot of AL 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 asymptotically come with no benefit over random test sampling. Since our proposed model is GPR-based, 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 state-of-the-art performance to recent, sophisticated, model-agnostic AL approaches. This was done by training different models on the actively constructed training sets of their and other AL approaches and assessing their performance. In particular, they compared favorably to Goetz et al. (2018)—a random tree-based AL approach—in a heteroscedastic setting, using a regression forest model and to Bull et al. (2013)—a wavelet-based AL approach—in a setting of inhomogeneous complexity, using an RBF-network. This demonstrates the flexibility 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.

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, a 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 a low-dimensional projection.

Parsimonious modeling   We proposed a novel, model-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 $v$ on the LOB of GPR and, thus, the adequate adjustment with respect to $v$ in the derivation of the LFC. As we argued, GPR already treats heteroscedasticity through local adaptions of the regularization. Hence, we assume the influence of $v$ on 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 $\delta\leq d$ and the smoothness $\alpha\in(0,\infty]$ of $f\in\mathcal{C}^{\alpha}\left(\mathcal{X},\operatorname*{\mathds{R}}\right)$ to be given through domain knowledge.

If $\delta$ is unknown, we can estimate it from unlabeled input instances, such as $\bm{X}_{pool}$ in a pre-processing step. However, this is beyond the scope of our work and we refer to the approach of Facco et al. (2017)555For an implementation in Python see https://scikit-dimension.readthedocs.io/en/latest/index.html for an estimate of the $\delta$ for unbalanced input distributions in high-dimensional input spaces $\mathcal{X}$.

If we have no ground truth knowledge about the smoothness $\alpha$ of the target function $f\in\mathcal{C}^{\alpha}\left(\mathcal{X},\operatorname*{\mathds{R}}\right)$, we resort to $\widehat{\alpha}=\infty$ as default in practice. This assumption is justified, as long as $f$ happens to be rougher in at most finitely many locations of the input space. Since, asymptotically, the violation of $\alpha=\infty$ 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 $\alpha$ is to deploy the Matérn kernel

[TABLE]

where $\Gamma$ is the gamma function and $\widetilde{k}_{\nu}$ is the modified Bessel function of the second kind or order $\nu$. After finding the best fitting $\nu^{*}$, we obtain by $\widehat{\alpha}:=\lceil\nu^{*}\rceil-1$ a reasonable estimate to $\alpha$. We defer this idea to future work as it is beyond the scope of this work.

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 a fixed number $L$ 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 representations scale poorly in $d$ as their IPs must be space-filling (Binois & Wycoff, 2022). Likewise, our density-based AL approach encounters decaying power for large $d$. A low intrinsic dimension ($\delta<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 is characterized by features $x\in\mathcal{X}$ with labels $f(x)\in\operatorname*{\mathds{R}}$, a hypothesis space of locally adaptive models, and the MISE as the loss function (or rather the pointwise MSE from (1) due to the localization of LFC). In this sense, LFC is formally a property of the combination of model and loss according to the characterization of (Jung (2022), 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 $\bm{X}_{pool}$ to be given that already follows the true molecular distribution. This will not be given in practice, since building the input trajectory 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 a 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, e.g., in multiscale problems from the sciences such as turbulence (Brunton et al., 2020) or quantum chemistry (Noé et al., 2020; von Lilienfeld et al., 2020; Unke 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 datasets 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 a model that naturally can cope with high input space dimensions on the other hand. Furthermore, we have shown how respecting LFC in the selection of IPs contributes to parsimonious modeling.

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 a 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 a 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(x)$. This claim should be supported by further theoretical investigation. While we deployed our estimates of LFC and the superior training density, using $\widehat{\alpha}=\infty$, if the smoothness of the target function $f\in\mathcal{C}^{\alpha}\left(\mathcal{X},\operatorname*{\mathds{R}}\right)$ is unknown, it is possible to (re-)estimate $\widehat{\alpha}$, e.g., by tuning the regularity of the Matérn kernel of a 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 idea is also to combine our LFC estimate with the IP selection approach by Moss et al. (2023) to 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 AI (e.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.

Acknowledgments

D. Panknin, S. Chmiela, 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 III, Autonomous Learning in Complex Environments (01IS18049B). K.-R. Müller was also supported by the BMBF Grants 01GQ1115 and 01GQ0850, under the Grants 01IS14013A-E, 031L0207A-D; DFG under Grant Math+, EXC 2046/1, Project ID 390685689 and by the Institute of Information & Communications Technology Planning & Evaluation (IITP) grants funded by the Korea Government (No. 2017-0-00451, Development of BCI 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.

Appendix 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 \Sigma\in\mathbb{S}^{d}_{\raisebox{0.0pt}{\tiny++}} in $x\in\mathcal{X}$ can be understood as follows: First, the regression problem is localized around $x$ according to weights $k^{\Sigma}(\cdot,x)$ that decrease with growing distance to $x$. Then we search for the polynomial up to order $Q$ that fits the localized regression problem best. Finally, the evaluation of this polynomial in $x$ is returned as the prediction. Formally, it is

[TABLE]

and $\mathcal{P}_{Q}(\operatorname*{\mathds{R}}^{d})$ is the space of the real polynomial mappings $\mathfrak{p}\!:\operatorname*{\mathds{R}}^{d}\,\rightarrow\,\operatorname*{\mathds{R}}$ up to order $Q$.

The localization is controlled by $\Sigma$ through the kernel weights $\bm{K}^{\Sigma}(x,x_{i})$ for $x_{i}\in\bm{X}_{n}$: For an RBF-kernel, $k^{\Sigma}(x,x^{\prime})$ decays monotonically with growing distance of $x^{\prime}$ to $x$. This decay is dampened or amplified as $\Sigma$ increases or decreases, respectively (in the sense of the Loewner order).

For readability, since $\Sigma$ will be replaced by terms with more involved notation, we redefine (1) by

[TABLE]

For a bandwidth space \mathcal{S}\subseteq\mathbb{S}^{d}_{\raisebox{0.0pt}{\tiny++}}, Panknin et al. (2021) proposed to minimize the AL objective

[TABLE]

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\in\mathcal{X}$ there exists a unique $\Sigma^{\prime}\in\mathcal{S}$ such that

[TABLE]

we are able to define the LOB function

[TABLE]

This function exists, for example, in the isotropic case $\mathcal{S}=\left\{\sigma\mathcal{I}_{d}\;\middle|\;\sigma>0\right\}$ for LPS under mild conditions, where we denote $\Sigma^{n}(x)=\sigma^{n}(x)\mathcal{I}_{d}$ (see, e.g., Masry (1996; 1997); Fan et al. (1997) or Panknin et al. (2021) for an overview).

Assuming the isotropic bandwidths candidate space $\mathcal{S}=\left\{\sigma\mathcal{I}_{d}\;\middle|\;\sigma>0\right\}$, the LOB as in Eq. (4) is an asymptotically well-defined function under mild assumptions666We require non-vanishing leading bias- and variance-terms of $m_{Q}(x)$, which is guaranteed if $\forall x\in\mathcal{X}$ it holds that $\text{b}_{Q}\left[x,\mathcal{I}_{d}\right]\neq 0$ from Eq. (38) and $v(x)>0$.: Denoting the LOB of LPS of order $Q$ by $\Sigma^{n}_{Q}(x)=\sigma^{n}_{Q}(x)\mathcal{I}_{d}$ such that

[TABLE]

asymptotically it holds

[TABLE]

where $C_{Q}$ is a constant, and $\text{b}_{Q}\left[x,\mathcal{I}_{d}\right]$ is a function of $x$ taken from the asymptotic conditional bias $f(x)-\operatorname*{\mathds{E}}\left[m^{h_{n}\mathcal{I}_{d}}_{Q}(x)\middle|\bm{X}_{n}\right]$ of LPS (Masry, 1996; 1997). That is, for a sequence $h_{n}\rightarrow 0$ as $n\rightarrow\infty$ we can write the conditional bias, which is of order $Q+1$, as

[TABLE]

Eq. (37) shows how LOB scales asymptotically with respect to the training size $n$, the local noise level function $v(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 properties and LOB itself, we are able to formulate LFC in a closed form.

Definition 6** (Panknin et al. (2021)).**

For LPS of order $Q$, the LFC of $f$ in $x\in\mathcal{X}$ is asymptotically given by

[TABLE]

As already mentioned in Eq. (35), given a test density $q$, the AL task is to minimize $\text{MISE}\left(q,m_{Q}|\bm{X}_{n}\right)$. Now, if LOB is well-defined, we can rewrite

[TABLE]

Finally, when solving for the optimal training dataset

[TABLE]

as in Eq. (3), the optimal training inputs $\bm{X}_{n}^{\prime}$ can be written asymptotically as an independent and identically distributed sample from the optimal training distribution, whose density $p^{Q,n}_{\text{\tiny{Opt}}}$ possesses an asymptotic closed form.

Theorem 7** (Panknin et al. (2021)).**

Let v,q\in\mathcal{C}^{0}\left(\mathcal{X},\mathds{R}_{\raisebox{0.0pt}{\tiny+}}\right) for a compact input space $\mathcal{X}$, where $q$ is a test probability density. Additionally, assume that $v$ and $q$ are bounded away from zero. I.e., $v,q\geq\epsilon$ for some $\epsilon>0$. Let $k$ be a RBF-kernel with bandwidth parameter space $\mathcal{S}=\left\{\sigma\mathcal{I}_{d}\;\middle|\;\sigma>0\right\}$. Let $Q\in\operatorname*{\mathds{N}}$ be odd and $f\in\mathcal{C}^{Q+1}\left(\mathcal{X}\right)$ such that the bias of order $Q+1$ does not vanish almost everywhere. Then the optimal training density for LPS of order $Q$ is asymptotically given by

[TABLE]

We will use this optimal distribution to sample $\bm{X}_{n}^{\prime}\sim p^{Q,n}_{\text{\tiny{Opt}}}$ with a proposed estimator for $\mathfrak{C}^{n}_{Q}$ that is scalable with respect to the input space dimension.

For LPS with $\bm{X}_{n}^{\prime}\sim p$ and $\bm{X}_{n}\sim q$, we can asymptotically calculate the relative required sample size from Definition 2 in Sec. 4.1 by

[TABLE]

Appendix B Analytic GPR formulations

B.1 Classical Gaussian process regression

The GPR model $\widehat{y}\sim\mathcal{GP}(\theta)$ (see, e.g. Williams & Rasmussen (1996)) is defined as follows: The GP is described by the hyperparameters $\theta=(\mu,\lambda,\widehat{v},\Sigma)$, which are the global constant prior mean $\mu$, the regularization parameter $\lambda$, the label noise variance function $\widehat{v}$ and the bandwidth matrix $\Sigma$ of the kernel. If we can assume homoscedastic noise, we let $\widehat{v}(x)\equiv\sigma_{\varepsilon}^{2}$.

The GP prior then assumes the labels $\bm{Y}_{n}$ of $\bm{X}_{n}$ to be distributed according to $\bm{Y}_{n}=\widehat{y}(\bm{X}_{n})\sim\mathcal{N}(\cdot;\bm{\mu}(\bm{X}_{n}),\bm{C}(\bm{X}_{n})|\theta)$, for the constant mean function $\bm{\mu}(\bm{X}_{n})=\mu\mathbbmss{1}_{n}$, and the covariance function

[TABLE]

where $\bm{K}_{n}=\bm{K}^{\Sigma}(\bm{X}_{n})$ is the kernel matrix of $\bm{X}_{n}$.

For test inputs $\bm{X}_{*}$, the posterior predictive distribution of $\bm{Y}_{*}$ is then given by

[TABLE]

where the predictive mean and covariance are given by

[TABLE]

and we have defined

[TABLE]

B.2 Analytic sparse Gaussian processes

We define the sparse GPR model $\widehat{y}\sim\mathcal{SGP}(\theta)$ as follows, following Snelson & Ghahramani (2005): The sparse GP is described by the (hyper-) parameters $\theta=(\mu,\lambda,\widehat{v},\Sigma,\bm{X}_{\dagger})$, which are the global constant prior mean $\mu$, the regularization parameter $\lambda$, the label noise variance function $\widehat{v}$, the bandwidth matrix $\Sigma$ of the kernel and the prior distribution, given by the IP locations $\bm{X}_{\dagger}\in\mathcal{X}^{m}$.

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 $m_{n}=o[n]$ is typically much slower than $n$. If we can assume homoscedastic noise, we let $\widehat{v}(x)\equiv\sigma_{\varepsilon}^{2}$.

The sparse GP then outputs

[TABLE]

for the mean function

[TABLE]

and the covariance function

[TABLE]

where we have defined $\bm{K}_{\dagger}=\bm{K}^{\Sigma}(\bm{X}_{\dagger})$, $\bm{K}_{n}=\bm{K}^{\Sigma}(\bm{X}_{n})$, $\bm{K}_{*\dagger}=\bm{K}^{\Sigma}(\bm{X}_{*},\bm{X}_{\dagger})$, $\bm{K}_{n\dagger}=\bm{K}^{\Sigma}(\bm{X}_{n},\bm{X}_{\dagger})$, $\bm{Q}_{\dagger}=\bm{K}_{\dagger}+\bm{K}_{n\dagger}^{\top}(\Lambda+\textbf{{diag}}(\widehat{v}(\bm{X}_{n})))^{-1}\bm{K}_{n\dagger},$ and $\Lambda=\textbf{{diag}}(\bm{\lambda})$ with $\bm{\lambda}=\textbf{{diag}}(\bm{K}_{n}+\bm{K}_{n\dagger}\bm{K}_{\dagger}^{-1}\bm{K}_{n\dagger}^{\top}).$

We choose $\bm{\mu}$ to be the constant mean function, i.e., $\bm{\mu}(X)=\mu\mathbbmss{1}_{n}$ for $X\in\mathcal{X}^{n}$, noting that other mean functions are possible.

Appendix C LFC of GPR

See 1

Proof.

Let $\mathcal{X}=\sideset{}{{}_{i=1}^{k}}{\biguplus}\mathcal{X}^{k}_{i}$ be a segmentation of the input space with non-empty interiors $(\mathcal{X}^{k}_{1})^{\circ},\ldots,(\mathcal{X}^{k}_{k})^{\circ}\neq\emptyset$, over which we can define the restricted bandwidth function search space

[TABLE]

Here, $\sideset{}{{}_{A}}{\operatorname*{\scalerel*{\mathbbmss{1}}{\textstyle\sum}}}(z)$ is the indicator function, returning $1$ for $z\in A$ and [math], else. Furthermore let $\Sigma^{k,n}\in\bm{\mathcal{S}}_{k}$ be the minimizer of the MISE over $\bm{\mathcal{S}}_{k}$ with $\Sigma^{k,n}(x)=\sideset{}{{}_{i=1}^{k}}{\sum}\sideset{}{{}_{\mathcal{X}^{k}_{i}}}{\operatorname*{\scalerel*{\mathbbmss{1}}{\textstyle\sum}}}(x)\Sigma^{k,n}_{i}$ such that

[TABLE]

Recall from (9) that $\Sigma^{n}_{\text{GPR}}\propto n^{-\frac{1}{2\alpha+d}}$ generally holds for arbitrary input spaces. Due to this, asymptotically, $\widehat{f}_{|_{\mathcal{X}^{k}_{i}}}$ does not depend on training samples outside $\mathcal{X}^{k}_{i}$. Hence, letting $\bm{X}_{i,n}:=\left\{x\in\bm{X}_{n}\;\middle|\;x_{i}\in\mathcal{X}^{k}_{i}\right\}$, the individual $\Sigma^{k,n}_{i}$ are asymptotically found by solving the isolated segments of the objective

[TABLE]

First of all, it is $\operatorname*{\mathds{E}}\bm{X}_{i,n}=p(\mathcal{X}^{k}_{i})n$, where $p(A):=\sideset{}{{}_{A}}{\operatorname*{{\int}}}p(x)dx$ is the probability for a training sample to fall into $A\subset\mathcal{X}$. In addition, we need to account for the expanse of $\mathcal{X}^{k_{n}}_{i}$, which we measure by $\textbf{{Vol}}(\mathcal{X}^{k_{n}}_{i})$. Here, $\textbf{{Vol}}(A):=\textstyle\sideset{}{{}_{A}}{\operatorname*{{\int}}}dx$ is the volume of $A\subset\mathcal{X}$. Again with (9), it is therefore

[TABLE]

Subsequently, we can slowly refine the segmentation $\mathcal{X}=\sideset{}{{}_{i=1}^{k_{n}}}{\biguplus}\mathcal{X}^{k_{n}}_{i}$, where $\max_{1\leq i\leq k_{n}}\textbf{{Vol}}(\mathcal{X}^{k_{n}}_{i})\rightarrow 0$ for $k_{n}\rightarrow\infty$ slow enough (with $k_{n}=o(n)$). Then, for almost every $x\in\mathcal{X}$, there exists a sequence $(i_{k,x})_{k\in\operatorname*{\mathds{N}}}$ with $x\in\mathcal{X}^{k}_{i_{k,x}}$ for all $k\in\operatorname*{\mathds{N}}$ such that

[TABLE]

By construction, it is $\Sigma^{k_{n},n}_{i_{k_{n},x}}=\Sigma^{n}_{\text{GPR}}(x)(1+o_{p}\left[1\right])$. It follows $\Sigma^{n}_{\text{GPR}}(x)=p(x)n(1+o_{p}\left[1\right])$ such that $\left|\Sigma^{n}_{\text{GPR}}(x)\right|=\left[p(x)n\right]^{\frac{d}{2\alpha+d}}(1+o_{p}\left[1\right])$. Therefore, asymptotically, $\mathfrak{C}^{n}_{\text{GPR}}(x):=\left[\frac{1}{p(x)n}\right]^{\frac{d}{2\alpha+d}}\left|\Sigma^{n}_{\text{GPR}}(x)\right|^{-1}$ does not depend on $n$ and $p$. Under homoscedasticity, asymptotically, $\mathfrak{C}^{n}_{\text{GPR}}$ is necessarily a function that only depends on $f$, which justifies its use as a measure of LFC. ∎

Appendix D Algorithmic summary of the proposed AL framework

Appendix E Finding diverse IP locations

In order to obtain diverse IP locations with a certain distribution, we consider two approaches, Stein variational gradient descent (SVGD) (Liu & 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, of which we require the density, as well as its derivative (Liu & 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 to 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 a considerable number of times, we 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 is similar to the known k-means clustering (Gan et al., 2020) but solves a different inertia objective, that is modified such that asymptotically, as the number of cluster centers $\left|\bm{X}_{\dagger}\right|\rightarrow\infty$, the distribution of the training data is preserved (Krishna et al., 2019). Under the standard k-means clustering objective, we would observe $\bm{X}_{\dagger}\sim p^{\frac{d}{2+d}}$ (Graf & Luschgy, 2007), where it was $\bm{X}_{n}\sim 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 to obtain a representative set of IPs. Due to very mild assumptions on the problem, DC is specifically easy to perform in higher dimensions.

Dealing with local optima of DC   The inertia objective of DC is given by

[TABLE]

where

[TABLE]

are those elements in $\bm{X}_{n}$ 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 is given in the case of DC by

[TABLE]

It is a known problem that k-means-related 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 running multiple repetitions of k-means, we will extend the state-of-the-art method k-means $++$ for choosing the initial set of clusters in a more sophisticated way, where we additionally account for the training distribution. Given the inertia objective

[TABLE]

of the cluster centers $\bm{c}$, the k-means $++$ procedure builds the set of initial cluster centers as follows: Draw the first center $c_{1}$ randomly from $\bm{X}_{n}$. Then keep track of the current closest squared distance

[TABLE]

of each element $x_{i}\in\bm{X}_{n}$ to the so far drawn centers $c_{1},\ldots,c_{m}$ and sample the next center $c_{m+1}$ with probability $\propto(d_{i}^{m})_{i=1}^{n}$ from $\bm{X}_{n}$. 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 $\bm{X}_{n}$. However, in its standard form, the centers initialized by k-means $++$ are themselves distributed flatter than $\bm{X}_{n}$. And so, in the case of DC, we propose the following adjustment for a distributional k-means $++$:

We sample with probability $\propto\left(d_{i}^{m}p(x_{i})^{2/d}\right)_{i=1}^{n}$ from $\bm{X}_{n}$, where $\bm{X}_{n}\sim p$.

Symmetrized DC for molecules   Since any symmetric molecule has multiple equivalent representations, care must be taken 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.1, let

[TABLE]

be the symmetrized distance between two molecule representations. The symmetrized DC algorithm is then obtained by replacing all occurrences of $\|z-z^{\prime}\|_{\raisebox{-2.0pt}{\tiny$$}}$ with $d(z,z^{\prime})$ in the cluster assignments $I_{c}$, the objective $\text{inertia}_{\text{DC}}[\bm{c}|\bm{X}_{n}]$, the cluster updates $c^{*}$ and closest distances $d_{i}^{m}$ from Equations 43, 42, 44 and 45.

Appendix F Design choices of the sparse MoE model

In Sec. 4.3 we have made several design choices with computational feasibility in mind. We will discuss these summarized in this section.

The gate model   While in Sec. 4.3 we have chosen the gate $g_{l}\sim\mathcal{SVGP}(\theta_{g_{l}})$ to be a 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 MoE. 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 $\bm{X}^{G}_{\dagger}$ across all gate channels. While this is not necessary, it simplifies our method without costs as the MoE is rather insensitive concerning the gate IP 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. 4.3, whether these experts should be sparse or dense. Here, the deciding factor is 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 $\bm{X}^{E}_{\dagger}$ across all experts, which simplifies our model. In contrast to the gate situation, the MoE is sensitive to the choice of expert IP 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 gate, or we use leave-one-out expert responses on a shared training set.

In the sparse expert case, it is necessary to learn reasonable inducing values $\bm{\mu}_{\dagger}$ prior to the actual learning procedure of the MoE to not get stuck in a spurious solution. Therefore, there should be a short pre-training phase for each individual expert.

In addition—whether or not the experts are sparse—the shared expert parameters $\mu_{E},\lambda_{E},\widehat{v},\Sigma_{E}$ should be initialized reasonably. In this regard, we suggest training a single, global expert model before the (pre-)training 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 $\Sigma_{E}=\sigma_{E}\mathcal{I}_{d}$ and learn the scalar $\sigma_{E}>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 $\Sigma_{E}$ 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 $\bm{S}_{\dagger}=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 no 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 IP 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 one that is capable of modeling the function. Now that we are able to model a comparably flat function by small bandwidths, as long as 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 in Sec. 5.1.

The gate noise $\mathfrak{s}$   Like already mentioned in Sec. 4.4.1, it is possible to tune $\mathfrak{s}$ in the training process:

Remark 8**.**

Shazeer et al. (2017)** proposed to learn the $\mathfrak{s}$ parameter by adding a penalty term to the main objective that penalizes the imbalance of how likely training inputs are assigned to each expert: Let $\pi_{b}\in[0,1]^{L}$ be the expert assignment probabilities of $x_{b}$ and define $\pi_{\mathcal{B}}=\sideset{}{{}_{b\in\mathcal{B}}}{\sum}\pi_{b}$. Then they add a penalty \operatorname*{\mathds{V}}[\pi_{\mathcal{B}}]\big{/}[\operatorname*{\mathds{E}}\pi_{\mathcal{B}}]^{2} to the objective, which is the squared coefficient of variation—a coefficient that accounts for the non-uniformity of a set of positive variables.

We justify our simple heuristic to shrink $\mathfrak{s}$ in a static way as follows: Recall from Sec. 4.3 that $\mathfrak{s}$ prevents premature commitment to a spurious solution. When treating $\mathfrak{s}$ as a 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 $\mathfrak{s}$ behaves best when decaying towards zero as the training progresses.

Appendix G Supplemental results on the Doppler experiment

G.1 The single-scale GPR model

When training a single-scale GPR model on the Doppler dataset, the tuned bandwidth parameter will typically take an intermediate value, trying to compromise between more complex and simpler regions. This is reflected in the predictions in Fig. 16, 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 F 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 $\vartheta_{\sigma}\text{pen}_{\sigma}(\bm{X}_{n},\bm{Y}_{n},\mathcal{B},w,\Theta)$ with the penalty term $\text{pen}_{\sigma}$ from (25) and a scaling factor $\vartheta_{\sigma}\geq 0$.

Now, while the LOB estimate with $\vartheta_{\sigma}=0.5$ (see Fig. 4) consistently behaves as expected, we show for comparison a typical LOB estimate in Fig. 18 that results from applying no regularization ($\vartheta_{\sigma}=0$). By chance—here, the flat region of the Doppler function to the right—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 a detrimental oversampling of these locations by the AL procedure.

Appendix H Supplemental on the malonaldehyde MD simulation experiment

H.1 The sGDML model

The GDML model by Chmiela et al. (2017) represents the geometry $x=\left[R_{1},\ldots,R_{\bm{a}}\right]\in\operatorname*{\mathds{R}}^{3\times\bm{a}}$ of each molecule in terms of the reciprocal distances $\Phi(x)_{kl}=\|R_{k}-R_{l}\|_{\raisebox{-2.0pt}{\tiny$$}}^{-1}$ of all atom-pairings to achieve roto-translational invariance of the input. This representation gives us a total $\bm{d}=\bm{a}(\bm{a}-1)/2$ input features. The similarity of a pair of configurations $(z,E,F)$ and $(z^{\prime},E^{\prime},F^{\prime})$ is then given by the extended covariance function

[TABLE]

Hence, we denote the overall kernel function of two configurations by

[TABLE]

Atoms of the same type are physically identical and therefore exchangeable, albeit only a small subset of such symmetries is exercised at a 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.

The symmetric extension sGDML (Chmiela et al., 2018; 2019) automatically identifies all accessed atom permutations from the training set and adds this symmetric prior to the covariance function. Formally, let $(\pi_{s})_{s=1}^{\bm{s}}$ be atomic permutations that lead to an equivalent molecular representation. Then, the extended symmetric kernel of sGDML is given by

[TABLE]

Malonaldehyde possesses $\bm{s}=4$ such permutations.

Remark 9**.**

The identified set of permutations is transitively closed to form a group. Under isotropy, it suffices to permute only one of the two configurations given to the kernel: Permuting both entries (as in (46)) equals permuting one entry and multiplying by the constant $\bm{s}$. However, if the applied bandwidth is not of the form $\Sigma=\sigma\mathcal{I}_{\bm{d}}$, this property does not hold.

H.2 The malonaldehyde MD simulation experiment under a uniform test distribution

In this scenario, we assume a uniform test density $q=\mathcal{U}(\mathcal{X})$. Accordingly, we weight the validation and test MSE by the importance weights $1/\widehat{p}_{\mathcal{X}}(\bm{X}_{\text{val}})$ and $1/\widehat{p}_{\mathcal{X}}(\bm{X}_{\text{T}})$. We draw the initial expert training set $\bm{X}_{n}$ of size $n=2^{9}$ and the gate training set $\bm{X}_{n_{G}}^{G}$ via importance sampling from the remaining pool with weights $1/\widehat{p}_{\mathcal{X}}(\bm{X}_{pool}\setminus(\bm{X}_{\text{val}}\cup\bm{X}_{\text{T}}))$. By this it is $\bm{X}_{n}\sim\mathcal{U}(\mathcal{X})$.

In Fig. 20 we show the estimates of LOB, LFC, and the superior training density under the pool test distribution, evaluated on the relaxed configurations of malonaldehyde. The LFC estimates in Fig. 20 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. 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 \textstyle\raisebox{1.0pt}{\varrho}(\widehat{f}_{\text{\tiny{MoE}}},\widehat{p}^{\text{GPR},n}_{\text{\tiny{Sup}}})=0.965\pm 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 ($\bm{X}_{n_{0}},\bm{Y}_{n_{0}}$) as well as the validation set ($\bm{X}_{\text{val}},\bm{Y}_{\text{val}}$) by simulating the true MD trajectory solving the computationally expensive Schrödinger equation. Estimate the initial MD density $p_{\mathcal{X},0}$ based on ($\bm{X}_{n_{0}}\cup\bm{X}_{\text{val}}$).

For k\in\operatorname*{\mathds{N}_{\raisebox{0.0pt}{\tiny0}}}:

• •

Set $q^{k}\leftarrow p_{\mathcal{X},k}^{1-\frac{1}{k+1}}$ to encourage exploration in early iterations and exploitation in later iterations

• •

Estimate the model $\widehat{f}_{k}:=\widehat{f}_{\text{\tiny{MoE}}}$ based on ($\bm{X}_{n_{k}},\bm{Y}_{n_{k}}$)

• •

Estimate $p^{\text{GPR},n_{k}}_{\text{\tiny{Sup}}}$ based on $q^{k}$ and $\widehat{f}_{k}$

• •

Sample a large pool ($\bm{X}_{pool,k+1},\widehat{\bm{Y}}_{pool,k+1}$) of, e.g., 100,000 candidates by simulating the approximate MD trajectory using the computationally cheap model $\widehat{f}_{k}$. While simulation, avoid unreliable out-of-distribution predictions, e.g., by resetting the trajectory in some $x$ whenever $p_{\mathcal{X},k}(x)<\epsilon$ drops below a reasonable threshold.

• •

Estimate the trajectory density $p_{\mathcal{X},k+1}$ of $\bm{X}_{pool,k+1}$

• •

Update the training set ($\bm{X}_{n_{k+1}},\bm{Y}_{n_{k+1}}$) by selecting input candidates from the pool $\bm{X}_{pool,k+1}$ with distribution $p^{\text{GPR},n_{k}}_{\text{\tiny{Sup}}}$ and estimating the respective labels solving the Schrödinger equation