# Estimating Local Function Complexity
via Mixture of Gaussian Processes

Danny Panknin, Shinichi Nakajima, Thanh Binh Bui, and Klaus-Robert Müller All
authors are with the Machine Learning Group, Technische Universität Berlin,
Marchstraße 23, MAR 4-1, 10587 Berlin, Germany.K.-R. Müller is also with the
Department of Brain and Cognitive Engineering, Korea University, Anam-dong,
Seongbuk-gu, Seoul 02841, Korea.K.-R. Müller is also with the Max Planck
Institute for Informatics, Stuhlsatzenhausweg, 66123 Saarbrücken,
Germany.Correspondence to<EMAIL_ADDRESS>
###### Abstract

Real world data often exhibit inhomogeneity, e.g., the noise level, the
sampling distribution or the complexity of the target function may change over
the input space. In this paper, we try to isolate _local function complexity_
in a practical, robust way. This is achieved by first estimating the _locally
optimal kernel bandwidth_ as a functional relationship. Specifically, we
propose _Spatially Adaptive Bandwidth Estimation in Regression_ (SABER), which
employs the _mixture of experts_ consisting of _multinomial kernel logistic
regression_ as a _gate_ and _Gaussian process regression_ models as _experts_.
Using the locally optimal kernel bandwidths, we deduce an estimate to the
local function complexity by drawing parallels to the theory of _locally
linear smoothing_. We demonstrate the usefulness of local function complexity
for model interpretation and active learning in quantum chemistry experiments
and fluid dynamics simulations.

###### Index Terms:

Adaptive kernel bandwidth, mixture of Gaussian processes, local function
complexity, active learning

##  I Introduction

Most inference problems intrinsically exhibit inhomogeneities of different
kinds. They might arise, for example, by non-uniform sampling of inputs,
heteroscedasticity of observation noise or differences in complexity of the
function to infer. Previous works show the benefit of treating inhomogeneities
_explicitly_. They comprise approaches from _multiple kernel learning_ (MKL)
[1], _mixture of experts_ [2, 3, 4], etc.

There are approaches that treat one kind of the aforementioned types of
inhomogeneity, while assuming the remaining components to behave
homogeneously: In [5, 6] they deal with heteroscedasticity by adapting the
regularization locally. Some approaches adjust bandwidths locally with respect
to the input density [7, 8, 9, 10]. Finally, there are approaches that tackle
inhomogeneous complexity by stacking several kernel-linear models at different
bandwidths, either learnt jointly [11, 12] or hierarchically [13, 14].
Although those approaches have shown their usefulness, they do not cope with
more realistic cases featuring multiple types of inhomogeneity.

We focus on regression problems with a fixed, given input feature
representation. More formally, we try to estimate a function $f$ from noisy
observations $y_{i}=f(x_{i})+\varepsilon_{i},i\in\\{1,\ldots,n\\}$ where the
$\varepsilon_{i}$ are independent, with mean
$\operatorname*{\mathds{E}}[\varepsilon_{i}]=0$ and local noise variance
$\operatorname*{\mathds{V}}[\varepsilon_{i}]=v(x_{i})$ and $x_{i}\sim p$ are
i. i. d. samples according to a probability density $p$ defined on
$\mathcal{X}\subset\operatorname*{\mathds{R}}^{d}$.

Figure 1: The oscillating toy data, inhomogeneous in three regards—function
complexity, noise level and input density. See Sec. V-A for further details.

Fig. 1 shows an example dataset that exhibits all inhomogeneities
simultaneously. Intuitively, the complexity of the true function shrinks
monotonically as the frequency of $f$ decreases.

Some approaches on dealing with inhomogeneous noise and input density model
the functions $v$, respectively $p$ explicitly. Treating function complexity
in that way has received less attention. Although not necessary in terms of
prediction performance, it complies with the increasing demand on model
interpretability [15, 16, 17]. Recently, [18] derived a theoretically sound,
scalar-valued, continuous _local function complexity_ (LFC) measure
$\mathfrak{C}_{n}^{\text{\tiny{LLS}}}$ for general d-dimensional data.
Basically, they define LFC as the reciprocal _locally optimal kernel
bandwidth_ (LOKB) that is adjusted for $v$ and $p$ in advance. Beyond that,
they show how $\mathfrak{C}_{n}^{\text{\tiny{LLS}}}$ can be used in active
learning.

By bandwidth one refers to the positive definite matrix $\Sigma$ of a _radial
basis function_ (RBF)-kernel such as the Gaussian kernel

$\displaystyle
k^{\Sigma}(x,y)=\exp\left\\{-\frac{1}{2}\|\Sigma^{-1}(x-y)\|_{\raisebox{-2.0pt}{\tiny$2$}}^{2}\right\\}.$

Here, $\Sigma$ controls the resolution of the kernel. For (d = 1) we will use
lower case $\sigma$ instead of $\Sigma$ throughout. We call
$\bm{\Sigma}_{\text{\tiny{Opt}}}(x)$ locally optimal in $x\in\mathcal{X}$, if
it optimizes the prediction of $f(x)$ with the kernel given above.

With the goal to improve prediction accuracy, there exist approaches on
estimating LOKB such as adaptive bandwidth RBF-networks [19, 20, 21, 22] that
tune the bandwidths individually for each basis function. However, these give
estimates of LOKB only in the basis function centers and no continuous
relation of LOKB is enforced.

The analysis of _locally linear smoothing_ (LLS) (see Sec. II-A for details)
provides an asymptotic closed-form solution to LOKB. In the one-dimensional
case it is given by

$\displaystyle\sigma_{\mathrm{theo}}(x)=\left[\frac{v(x)}{f^{{}^{\prime\prime}}(x)^{2}p(x)}\right]^{\frac{1}{5}}.$
(1)

Here, a functional relation is automatically fulfilled and we might obtain a
measure of LFC similar to [18] by

$\displaystyle\mathfrak{C}_{\infty}^{\text{\tiny{LLS}}}(x)=\left|f^{{}^{\prime\prime}}(x)\right|^{\frac{2}{5}}.$
(2)

We plotted both functions in black in Fig. 2 for the toy example. The poles
correspond to the roots $f^{{}^{\prime\prime}}=0$, which are asymptotically
justified, but misleading for finite data settings.

In the d-dimensional case $\mathfrak{C}_{\infty}^{\text{\tiny{LLS}}}(x)$
depends on the Hessian $H_{f}(x)$. There are attempts [23] to estimate LOKB
for LLS from the factorization $\sigma_{\mathrm{theo}}$ for (d = 1, 2),
providing non-parametric estimates for $p$ and $v$, and for each component of
$H_{f}$. However, such an approach is subject to the counterintuitive
modelling near the poles, that arise with respect to $H_{f}$. Additionally,
the number of components to estimate grows quadratically in $d$, making it
computationally challenging for large $d$. Beyond that, the coordinate
estimates of $H_{f}(x)$ have to be summarized, which suffers from uncertainty
propagation [24]. Both drawbacks can be circumvented by directly modelling the
locally optimal bandwidth $\bm{\Sigma}_{\text{\tiny{Opt}}}(x)$ as a functional
relationship.

The LLS model is interesting for its well elaborated analysis. However, it
lacks a robust estimate of LOKB. With respect to kernel-linear models, we
could adapt the _localized multiple kernel regression_ (LMKR) approach [25] to
directly model LOKB with functional relation. We discuss in Sec. V-A2, why
LMKR is however not well suited for this task.

Figure 2: LOKB (solid) and the associated LFC (dashed) for the oscillating
dataset (cf. Fig. 1). Once estimated by our proposed method SABER (green), and
once given by (1) and (2) (black), whose components are shown in Fig. 3.

The main contribution of our work is to estimate LOKB accurately and stably
with our novel method, called _spatially adaptive bandwidth estimation in
regressions_ (SABER). SABER is a mixture of experts (ME) model consisting of
_multinomial kernel logistic regression_ (MKLR) as a gate function, and
_Gaussian process regression_ (GPR) models as _experts_. In SABER, each GPR
expert has a preset bandwidth, and the gate function and the regression
parameters are simultaneously trained. Our estimate
$\Sigma_{\text{\tiny{SABER}}}$ to LOKB is a weighted sum of the bandwidths,
where the weights are the output of the gate. Inspired by the theory of LLS,
we then propose our estimate $\mathfrak{C}$ to LFC based on
$\Sigma_{\text{\tiny{SABER}}}$, compensating for noise and input density,
which is scalar-valued, continuous and therefore interpretable.

In Fig. 2 we show (in green) our estimate $\sigma_{\text{\tiny{SABER}}}$ for
LOKB and the resulting LFC estimate $\mathfrak{C}$ for the example dataset
from Fig. 1. Here, we can verify the intuition that the complexity of the true
function shrinks.

In our experiments, the advantages of SABER are demonstrated in terms of
prediction accuracy, problem interpretation, and the use for active learning:
SABER performs comparably to deep neural networks on a dataset of molecule
configurations. Additionally our estimate $\mathfrak{C}$ provides new insights
into physical systems. Finally, we feed our estimate $\mathfrak{C}$ as an
input to active sample selection in a fluid dynamics application.

The paper is organized as follows: The mathematical notation and some basic
machine learning concepts are introduced in Sec. II. Our method SABER and the
complexity estimate $\mathfrak{C}$ are proposed in Sec. III. Afterwards we
will discuss related work on adaptive bandwidths in Sec. IV. We will then
demonstrate the capabilities of our novel approach not only to improve
prediction performance, but also to give interpretable insights in
applications of quantum-chemistry and to guide adaptive mesh selection in
numerical simulations in Sec. V. Finally, we conclude in Sec. VI.

##  II Preliminaries

In this section, we introduce some techniques our proposed method is based on.
To start with, we define some mathematical notation: Let
$A,B\in\operatorname*{\mathds{R}}^{m\times n}$ and
$v\in\operatorname*{\mathds{R}}^{n}$.

  * •

$A\odot B=[A_{ij}B_{ij}]_{ij}$ is the Hadamard (pointwise) product

  * •

In analogy $A\oslash B=[A_{ij}/B_{ij}]_{ij}$ is the pointwise division

  * •

$\operatorname*{\bf{diag}}(v)_{ij}=\begin{cases}v_{i},i=j\\\
0,\text{else}\end{cases}$ is a matrix with $v$ on its diagonal

  * •

$\mathbbmss{1}_{n}=\begin{bmatrix}1&\cdots&1\end{bmatrix}^{\top}\in\operatorname*{\mathds{R}}^{n}$
the length-n vector of ones

  * •

$\mathcal{I}_{n}=\operatorname*{\bf{diag}}(\mathbbmss{1}_{n})$ the identity
matrix in $\operatorname*{\mathds{R}}^{n\times n}$

The local linear smoothing model is well analyzed with respect to local
properties. In order to make use of them we will now recap some theoretical
results. After that, we introduce GPR and MKLR, which are the basic components
of our proposed ME framework.

###  II-A Local Linear Smoothing

For a fix bandwidth matrix $\Sigma$ consider the linear expansion

$\widehat{f}(z)=\alpha_{\Sigma,x}+\beta_{\Sigma,x}^{\top}(z-x)$

with focus in $x$. The hyperparameters
$\alpha_{\Sigma,x}\in\operatorname*{\mathds{R}}$,
$\beta_{\Sigma,x}\in\operatorname*{\mathds{R}}^{d}$ are found by minimizing
the _local least squares problem_

$\sideset{}{{}_{i=1}^{n}}{\sum}(y_{i}-\alpha_{\Sigma,x}-\beta_{\Sigma,x}^{\top}(x_{i}-x))^{2}k^{\Sigma}(x_{i},x).$

Then the local linear smoother is given by

$m^{\Sigma}(x)=\alpha_{\Sigma,x}.$

Let $\bm{X}_{n}=\left\\{x_{i}\right\\}_{i=1}^{n}\sim p$ be a training set of
size $n$ and define the _conditional mean squared error_ in $x$ by

$\text{MSE}\left(x,\Sigma|\bm{X}_{n}\right)=\operatorname*{\mathds{E}}\left[(m^{\Sigma}(x)-f(x))^{2}|\bm{X}_{n}\right].$

Let $D_{f}$ be the derivative of $f$. [18] have shown the following:

###### Theorem 1

Let $D_{f}$ be analytic, $p$ continuously differentiable, and $v$ continuous
and bounded above $0$. On $\left\\{x\;\middle|\;p(x)>0\right\\}$,

$\Sigma_{\text{\tiny{Opt}}}^{\text{\tiny{LLS}}}(x,n)=\operatorname*{\bf{argmin}}_{\Sigma}\text{MSE}\left(x,\Sigma|\bm{X}_{n}\right)$

is well defined and continuous. The continuous function

$\displaystyle\mathfrak{S}^{n}(x)=\left[\frac{n\cdot
p(x)}{v(x)}\right]^{\frac{1}{4+d}}\Sigma_{\text{\tiny{Opt}}}^{\text{\tiny{LLS}}}(x,n)$
(3)

converges to a function $\mathfrak{S}^{\infty}$ that is independent of $p$,
$v$ and $n$, almost everywhere.

Therefore a measure of LFC is naturally defined by

$\displaystyle\mathfrak{C}_{\infty}^{\text{\tiny{LLS}}}(x)=\left|\mathfrak{S}^{\infty}(x)\right|^{-\frac{1}{d}}.$

As originally proven by [26] for one dimension, and extended to the isotropic
d-dimensional case ($\Sigma=\sigma\mathcal{I}_{d}$) by [23], it is

$\mathfrak{C}_{\infty}^{\text{\tiny{LLS}}}(x)=[\mathbbmss{1}_{d}^{\top}H_{f}(x)\mathbbmss{1}_{d}]^{\frac{2}{4+d}}.$

In the following we denote the corresponding asymptotically optimal bandwidth
function (up to $n$) by

$\displaystyle\sigma_{\mathrm{theo}}(x)=\left[\frac{v(x)}{[\mathbbmss{1}_{d}^{\top}H_{f}(x)\mathbbmss{1}_{d}]^{2}p(x)}\right]^{\frac{1}{4+d}}.$
(4)

It nicely reflects the influence of all discussed inhomogeneities.
Unfortunately, this asymptotic behavior is not always representative in case
of finite data: In the introduction, we have seen the strong impact of the
poles, where $\mathbbmss{1}_{d}^{\top}H_{f}(x)\mathbbmss{1}_{d}=0$, which do
neither match intuition, nor feature optimality [18].

A continuous approximation to LFC is given by

$\displaystyle\mathfrak{C}_{n}^{\text{\tiny{LLS}}}(x)=\left|\mathfrak{S}^{n}(x)\right|^{-\frac{1}{d}}$
(5)

Consider the _conditional mean integrated squared error_

$\text{MISE}\left(q|\bm{X}_{n}\right)=\sideset{}{{}_{\mathcal{X}}}{\int}\text{MSE}\left(x,\Sigma_{\text{\tiny{Opt}}}^{\text{\tiny{LLS}}}(x,n)|\bm{X}_{n}\right)q(x)dx,$

for a test density $q$. When optimizing $\text{MISE}\left(q|\bm{X}_{n}\right)$
with respect to the training samples $\bm{X}_{n}\sim p$, [18] have shown the
following training density to be asymptotically optimal:

###### Theorem 2

Let $q(x)$ be a test density. Under the assumptions of Theorem 1, minimizing
$\text{MISE}\left(q|\bm{X}_{n}\right)$ w.r.t. $\bm{X}_{n}\sim p$, gives

$\displaystyle
p_{\text{\tiny{Opt}}}(x)\propto\left[\mathfrak{C}_{n}^{\text{\tiny{LLS}}}(x)^{d}q(x)\right]^{\frac{4+d}{8+d}}v(x)^{\frac{4}{8+d}}.$
(6)

###  II-B Gaussian Process Regression

Let us define the GPR model (see, e.g. [27]) with the following
parameterization:
$\bm{Y}\sim\mathcal{GP}(m,\lambda_{f},\sigma_{\varepsilon}^{2},v,\Sigma)$ with
global mean $m$ and covariance function

$\bm{C}(\bm{X})=\sigma_{\varepsilon}^{2}(\exp(\lambda_{f})\bm{K}^{\Sigma}(\bm{X})+\Sigma_{\varepsilon}^{2}),$

where
$\bm{K}^{\Sigma}(\bm{X})=\left[k^{\Sigma}(x_{i},x_{j})\right]_{i,j=1}^{n}$ for
a RBF-kernel $k$, and

$\Sigma_{\varepsilon}^{2}=\operatorname*{\bf{diag}}\left(v(\bm{X})\right).$

For a test input $x^{*}$, the predictive distribution is given by

$\operatorname*{\mathds{P}}(y^{*}|x^{*},\bm{X},\bm{Y})\sim\mathcal{N}\left[\mu^{*}(x^{*}),\sigma^{*}(x^{*})\right],$

where

$\displaystyle\mu^{*}(x^{*})=m+c^{*}\bm{C}^{-1}(\bm{Y}-m\mathbbmss{1}),$ (7)
$\displaystyle\sigma^{*}(x^{*})=c^{**}-c^{*}\bm{C}^{-1}(c^{*})^{\top},$ (8)

for

$\bm{C}(\bm{X}\cup\\{x^{*}\\})=\begin{bmatrix}\bm{C}&(c^{*})^{\top}\\\
c^{*}&c^{**}\end{bmatrix}.$

The test predictions
$\widehat{y}^{*}=\operatorname*{\mathds{E}}[y^{*}|x^{*},\bm{X},\bm{Y}]$ are
given by (7), which we summarized in Alg. 4 in Appendix C.

Let $\theta_{e}=(m,\lambda_{f},\sigma_{\varepsilon}^{2},\Sigma)$ be the set of
the GPR hyperparameters. Typically one tunes $\theta_{e}$ by minimizing the
observation negative log-likelihood

$\displaystyle-\log$
$\displaystyle\operatorname*{\mathds{P}}(\bm{Y}|\bm{X},v,\theta_{e})$
$\displaystyle\propto\frac{1}{2}\log\left|\det(\bm{C})\right|+\frac{1}{2}(\bm{Y}-m\mathbbmss{1})^{\top}\bm{C}^{-1}(\bm{Y}-m\mathbbmss{1}).$

However, later on we want to implement instances of GPR as the experts of an
ME model, which simplifies for an objective that factorizes with respect to
$\bm{Y}$, given $\bm{X}$. This is why we prefer to optimize the _weighted
leave-one-out negative log predictive likelihood_ instead:

$\displaystyle\ell_{\text{e}}$
$\displaystyle(\bm{X},\bm{Y},v,\theta_{e},P)=-\sum\limits_{i=1}^{n}P_{i}\log\operatorname*{\mathds{P}}(y_{i}|x_{i},\bm{X}_{-i},\bm{Y}_{-i})$
(9)
$\displaystyle\propto\sum\limits_{i=1}^{n}\frac{P_{i}}{2}\left[\frac{looR_{i}^{2}}{\sigma^{*}_{-i}(x_{i})}+\log(\sigma^{*}_{-i}(x_{i}))\right],$

where $looR_{i}=\mu^{*}_{-i}(x_{i})-y_{i}$ are the leave-one-out residuals
with the predictive mean $\mu^{*}_{-i}$ and the predictive variance
$\sigma^{*}_{-i}$ from (7) and (8), based on the restricted training set
$(\bm{X}_{-i},\bm{Y}_{-i})$. We then minimize
$\ell_{\text{e}}(\bm{X},\bm{Y},\theta_{e},P)$ with respect to $\theta_{e}$ via
gradient descent. Details on the optimization of $\theta_{e}$ can be found in
Appendix D, where we give explicit updates of $\Sigma$ for the Gaussian kernel
in the isotropic case $\Sigma=\exp(\sigma^{x})\mathcal{I}_{d}$ or diagonal
case
$\Sigma=\operatorname*{\bf{diag}}(\exp(\sigma^{x}_{1}),\ldots,\exp(\sigma^{x}_{d}))$.

###  II-C Multinomial Kernel Logistic Regression

Let $\bm{X}$ be the training data with class membership probabilities
$\bm{P}=[\bm{P}_{-L},P_{L}]=[P_{1},\ldots,P_{L}]$ with
$\bm{P}_{-L}\in\left[0,1\right]^{n\times(L-1)}$ and
$P_{L}=\mathbbmss{1}_{n}-\bm{P}_{-L}\mathbbmss{1}_{L-1}\in\left[0,1\right]^{n}$.
Note that this is a generalization of the common case where we have class
labels rather than probabilities such that $\bm{P}\in\\{0,1\\}^{n\times L}$.
Let $\theta_{g}=(\Sigma_{g},\lambda_{g})$ be the bandwidth and regularization
parameter of the MKLR model. MKLR estimates class membership probabilities via

$\displaystyle
Q(x)_{j}=\operatorname*{\mathds{P}}(z_{j}=1|x,\bm{\omega},\theta_{g})=\frac{q(x)_{j}}{1+\sideset{}{{}_{k=1}^{L-1}}{\sum}q(x)_{k}},$
(10)

where

$\bm{\omega}=\begin{bmatrix}\bm{\alpha}\\\
\bm{b}\end{bmatrix}=\begin{bmatrix}\alpha_{1}&\cdots&\alpha_{L-1}\\\
b_{1}&\cdots&b_{L-1}\end{bmatrix}$

are the feature weights and biases and where for $j<L$ it is
$q(x)_{j}=\exp(\bm{K}^{\Sigma_{g}}(x,\bm{X})\alpha_{j}+b_{j})$ and
$q(x)_{L}=1$ to avoid over-parameterization. Define $\bm{q}$ with
$\bm{q}_{i,j}=q(x_{i})_{j}$ for $1\leq i\leq n$, $1\leq j<L$. The optimal
$\bm{\omega}$ is defined as the minimizer of the regularized negentropy

$\displaystyle\ell_{\text{g}}(\bm{\omega},$
$\displaystyle\bm{X},\bm{P},\theta_{g})=\lambda_{g}\sum\limits_{j=1}^{L-1}\alpha_{j}^{\top}\bm{K}^{\Sigma_{g}}(\bm{X})\alpha_{j}$
(11)
$\displaystyle+\sum\limits_{i=1}^{n}\left[\log\left(1+\sum\limits_{j=1}^{L-1}\bm{q}_{ij}\right)-\sum\limits_{j=1}^{L-1}\bm{P}_{ij}\log(\bm{q}_{ij})\right].$

We specify an efficient training procedure for $\bm{\omega}$ in Alg. 6 with
derivation in Appendix E and summarized the test predictions according to (10)
in Alg. 5 in Appendix C.

##  III Proposed method

In this section, we set up our prediction model $\widehat{y}_{\text{SABER}}$
from which we deduce our estimate $\Sigma_{\text{\tiny{SABER}}}$ to
$\bm{\Sigma}_{\text{\tiny{Opt}}}$. Based on this, we propose our estimate
$\mathfrak{C}$ of LFC. To begin with, we would like to discuss possible model
choices and why we committed ourselves to the ME framework.

###  III-A Preconsiderations

Let
$\mathcal{E}\\!:\mathcal{X}\times\mathcal{S}_{+}^{d}\,\rightarrow\,\operatorname*{\mathds{R}},(x,\Sigma)\mapsto
E(x,\Sigma)$ be a performance measure for predicting the label of $x$ using
the positive definite bandwidth matrix $\Sigma\in\mathcal{S}_{+}^{d}$.
$\mathcal{E}$ could be a squared error or a negative log-likelihood estimate.
Then we would like to find

$\widehat{\bm{\Sigma}}_{\text{\tiny{Opt}}}=\operatorname*{\bf{argmin}}_{\bm{\Sigma}\in\mathcal{C}^{0}(\mathcal{X},\mathcal{S}_{+}^{d})}\sideset{}{{}_{\mathcal{X}}}{\int}\mathcal{E}(x,\bm{\Sigma}(x))p(x)dx.$

For example we could try to fit a kernel model

$\log(\bm{\sigma})(x)=\sideset{}{{}_{i=1}^{n}}{\sum}\beta_{i}k(x_{i},x),$

minimizing

$\widehat{\beta}=\operatorname*{\bf{argmin}}_{\beta}\sideset{}{{}_{i=1}^{n}}{\sum}\mathcal{E}(x_{i},\bm{\sigma}(x_{i})\mathcal{I})+\gamma\beta^{\top}K\beta.$

Unfortunately, estimating $\mathcal{E}(x_{i},\bm{\sigma}(x_{i})\mathcal{I})$
for a general $\bm{\sigma}(x_{i})$-input is computationally infeasible: We
would have to estimate an individual GPR model for each training point after
each update of $\beta$, resulting in a method of at least $\mathcal{O}(n^{4})$
in computation time. We could avoid that by also fitting
$\widehat{\mathcal{E}}$ from the augmented dataset
$\left(\begin{pmatrix}x_{i}\\\
\sigma_{j}\end{pmatrix},\mathcal{E}(x_{i},\sigma_{j}\mathcal{I})\right)_{ij}$
for a finite bandwidth grid $0<\sigma_{1}<\ldots<\sigma_{L}$. However, the
augmented dataset is quite large and when it comes to adaption of
hyperparameters, $\widehat{\mathcal{E}}$ must be refitted over and over again.

Consider the more general objective

$\widehat{\mu}=\operatorname*{\bf{argmin}}_{\mu\in\mathcal{C}^{0}(\mathcal{X},\mathcal{P}(\mathcal{S}_{+}^{d}))}\sideset{}{{}_{\mathcal{X}}}{\int}\sideset{}{{}_{\mathcal{S}_{+}^{d}}}{\int}\mathcal{E}(x,\Sigma)\mu_{x}(\Sigma)d\Sigma
p(x)dx,$

where we assume $x\mapsto\mu_{x}$ to be a continuous mapping and
$\mu_{x}\in\mathcal{P}(\mathcal{S}_{+}^{d})$ is a probability measure for
$\bm{\Sigma}_{\text{\tiny{Opt}}}(x)$, that is,
$\mu_{x}(S)=\operatorname*{\mathds{P}}(\bm{\Sigma}_{\text{\tiny{Opt}}}(x)\in
S)$ for $S\subset\mathcal{S}_{+}^{d}$.

This objective is a generalization of the previous one, considering the Dirac
measures $\mu_{x}(\\{\bm{\Sigma}(x)\\})=1$. While this objective is more
powerful it is even harder to fit in its general form. However, if we restrict
$\mu\in\mathcal{C}^{0}(\mathcal{X},\mathcal{P}(\bm{S}))$ for the finite set
$\bm{S}=\\{\Sigma_{1},\ldots,\Sigma_{L}\\}$, then the objective simplifies to
a latent bandwidth classification problem, fitting the ME framework. This
classifier—also called _gate_ —is given by

$Q(x)_{j}\leftarrow\operatorname*{\mathds{P}}(\bm{\Sigma}_{\text{\tiny{Opt}}}(x)=\Sigma_{j}).$

The advantage is that we can limit the number of GPR models to be maintained
simultaneously, which additionally need no expensive refitting after each
update of hyperparameters, if preprocessed appropriately.

###  III-B Spatially Adaptive Bandwidth Estimation in Regression

In the ME paradigm we assume an observation to be generated by one of several
processes: Let $(\bm{X},\bm{Y})$ be the training set of size $n$ and let
$\bm{Z}$ be hidden multinomial variables

$\bm{Z}_{ij}=\begin{cases}1,\text{process }j\text{ generated }y_{i}\\\
0,\text{else}\end{cases}.$

indicating the generating one of $L$ potential processes. Finally, let
$\Theta=\\{\theta_{g},\theta_{e_{1}},\ldots,\theta_{e_{L}}\\}$ be the set of
all model parameters for describing both, the prediction likelihood
$\operatorname*{\mathds{P}}(\bm{Y}|\bm{X},\theta_{e_{j}})$ when assuming
process $j$ to drive the observations, and the probability
$\operatorname*{\mathds{P}}(\bm{Z}|\bm{X},\theta_{g})$ that specific processes
drive the observations. The goal is to maximize the observation likelihood
$\operatorname*{\mathds{P}}(\bm{Y}|\bm{X},\Theta)$. We refer the reader to
Appendix B for details on general tuning of the ME model.

In contrast to classical ME interpretation, we learn each expert on the full
training set regardless of the belief in affiliation to the respective
generating process. Moreover we share a single set of hyperparameters for all
experts up to an individual bandwidth, i.e.
$\theta_{e_{j}}=(m,\lambda_{f},\sigma_{\varepsilon}^{2},\Sigma_{j})$. We
discuss this choice in Sec. IV, where we delimit our method from existing
mixture of Gaussian processes approaches.

As discussed in Appendix B, the training of the ME-model simplifies
considerably, if the observations $\bm{Y}$ are mutually independent with
respect to the predictive distribution of the experts. Unfortunately, this is
not the case for the GPR model. Inspired by [28] we circumvent this problem by
optimizing the leave-one-out predictive likelihood instead:

$\displaystyle\prod\limits_{i=1}^{n}\operatorname*{\mathds{P}}(y_{i}|x_{i},\bm{X}_{-i},\bm{Y}_{-i},\Theta)=$
$\displaystyle\prod\limits_{i=1}^{n}\sum\limits_{j=1}^{L}\operatorname*{\mathds{P}}(y_{i}|x_{i},\bm{X}_{-i},\bm{Y}_{-i},\theta_{e_{j}})\operatorname*{\mathds{P}}(\bm{Z}_{ij}=1|x_{i},\bm{X}_{-i},\bm{Y}_{-i},\theta_{g})$

Given our belief about the hidden variables $\bm{Z}$ we can augment this
likelihood via

$\displaystyle\operatorname*{\mathds{P}}(\bm{Y},\bm{Z}|\bm{X},\Theta)=\prod\limits_{i=1}^{n}\prod\limits_{j=1}^{L}\big{[}$
$\displaystyle\operatorname*{\mathds{P}}(y_{i}|x_{i},\bm{X}_{-i},\bm{Y}_{-i},\theta_{e_{j}})$
$\displaystyle\cdot\,$
$\displaystyle\operatorname*{\mathds{P}}(\bm{Z}_{ij}=1|x_{i},\bm{X}_{-i},\bm{Y}_{-i},\theta_{g})\big{]}^{\bm{Z}_{ij}}.$

From that point we can formulate an _expectation maximization_ (EM) procedure,
similar to the standard framework of ME: In the E-step we calculate

$\displaystyle\bm{P}^{(t)}_{ij}$
$\displaystyle=\operatorname*{\mathds{P}}(\bm{Z}_{ij}=1|\bm{X},\bm{Y},\Theta^{(t-1)})$
$\displaystyle=\frac{\bm{Q}^{(t-1)}_{ij}\operatorname*{\mathds{P}}(y_{i}|x_{i},\bm{X}_{-i},\bm{Y}_{-i},\theta^{(t-1)}_{e_{j}})}{\sideset{}{{}_{k=1}^{L}}{\sum}\bm{Q}^{(t-1)}_{ik}\operatorname*{\mathds{P}}(y_{i}|x_{i},\bm{X}_{-i},\bm{Y}_{-i},\theta^{(t-1)}_{e_{k}})},$

where we model the leave-one-out gate responses

$\bm{Q}^{(t)}_{ij}=\operatorname*{\mathds{P}}(\bm{Z}_{ij}=1|x_{i},\bm{X}_{-i},\bm{Y}_{-i},\bm{P}^{(t)}_{-i},\theta^{(t)}_{g})$

via MKLR according to Sec. II-C. Our final objective function is given by

$\displaystyle\ell^{(t)}$
$\displaystyle(\bm{X},\bm{Y},\Theta)=-\sideset{}{{}_{\bm{Z}^{(t)}}}{\operatorname*{\mathds{E}}}\log\operatorname*{\mathds{P}}(\bm{Y},\bm{Z}^{(t)}|\bm{X},\Theta)$
(12) $\displaystyle=$
$\displaystyle\underbrace{\left[-\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{L}\bm{P}^{(t)}_{ij}\log\operatorname*{\mathds{P}}(\bm{Z}_{ij}=1|x_{i},\bm{X}_{-i},\bm{Y}_{-i},\theta_{g})\right]}_{R^{(t)}_{g}(\theta_{g})}$
$\displaystyle+\underbrace{\left[-\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{L}\bm{P}^{(t)}_{ij}\log\operatorname*{\mathds{P}}(y_{i}|x_{i},\bm{X}_{-i},\bm{Y}_{-i},\theta_{e_{j}})\right]}_{R^{(t)}_{e}(\theta_{e_{1}},\ldots,\theta_{e_{L}})}$

Note that $R^{(t)}_{g}(\theta_{g})$ and
$R^{(t)}_{e}(\theta_{e_{1}},\ldots,\theta_{e_{L}})$ can be optimized
independently.

We keep the individual bandwidths $\Sigma_{j}$ fixed for now. We therefore
minimize $R^{(t)}_{e}$ solely with respect to the shared expert
hyperparameters in the M-step via

$(m,\lambda_{f},\sigma_{\varepsilon}^{2})^{(t)}\leftarrow\operatorname*{\bf{argmin}}_{m,\lambda_{f},\sigma_{\varepsilon}^{2}}R^{(t)}_{e}(\theta_{e_{1}},\ldots,\theta_{e_{L}}).$

Recalling (9) from Sec. II-B, we can rewrite

$R^{(t)}_{e}(\theta_{e_{1}},\ldots,\theta_{e_{L}})=\sum\limits_{j=1}^{L}\ell_{\text{e}_{j}}\left(\bm{X},\bm{Y},\theta_{e_{j}},\left[\bm{P}^{(t)}_{ij}\right]_{i}\right).$

Due to linearity of the gradient the optimization of the shared parameters
$(m,\lambda_{f},\sigma_{\varepsilon}^{2})$ is therefore a straightforward
extension to the derivations in Appendix D:

For the sake of an efficient implementation, we estimate the
eigendecompositions $[U_{j},\Lambda_{j}]$ of the kernel matrices associated to
the bandwidths $\Sigma_{j}$ in advance. Note that we can reduce the
heteroscedastic case to the math of the homoscedastic case by building the
pseudo eigendecompositions $[U_{j},\Lambda_{j}]$ as follows: We estimate the
eigendecomposition
$\Sigma_{\varepsilon}^{-1}\bm{K}^{\Sigma_{j}}\Sigma_{\varepsilon}^{-1}=\widetilde{U}_{j}\Lambda_{j}\widetilde{U}_{j}^{\top}$
and set $U_{j}=\Sigma_{\varepsilon}^{-1}\widetilde{U}_{j}$. This pseudo
eigendecomposition is to be used with caution, as it is only proper with
respect to quantities that use the inverse of the covariance function. We then
define a set of help variables in Alg. 1.

Algorithm 1
$\text{helpVariables}(\lambda_{f},\bm{Y},[U_{j},\Lambda_{j}]_{j=1}^{L})$

1:[-0.4ex]0.490.4pt4pt 3pt

2:Input

3:0.490.4pt4pt 3pt

4:GPR regularization $\lambda_{f}$

5:Training labels $\bm{Y}$

6:Eigendecompositions $[U_{j},\Lambda_{j}]_{j=1}^{L}$

7:[-0.4ex]0.490.4pt4pt 3pt

8:Output

9:0.490.4pt4pt 3pt

10:Matrices
$W,loo\mathbbmss{1},looY,A^{\mathbbmss{1}},A^{Y},d^{A}\in\operatorname*{\mathds{R}}^{n\times
L}$

11:[-0.4ex]0.490.4pt4pt 3pt

12:Procedure

13:0.490.4pt4pt 3pt

14:for $1\leq j\leq L$ do

15: Set
$\overline{\bm{D}}^{j}=U_{j}[\exp(\lambda_{f})\Lambda_{j}+\mathcal{I}_{n}]^{-1}U_{j}^{\top}$

16: $W_{j}=\operatorname*{\bf{diag}}(\overline{\bm{D}}^{j})$,
$loo\mathbbmss{1}_{j}=\overline{\bm{D}}^{j}\mathbbmss{1}\oslash W_{j}$,
$looY_{j}=\overline{\bm{D}}^{j}\bm{Y}\oslash W_{j}$,

17:
$A^{\mathbbmss{1}}_{j}=\left(\mathcal{I}-\overline{\bm{D}}^{j}\right)\overline{\bm{D}}^{j}\mathbbmss{1}$,
$A^{Y}_{j}=\left(\mathcal{I}-\overline{\bm{D}}^{j}\right)\overline{\bm{D}}^{j}\bm{Y}$,

18:
$d^{A}_{j}=\operatorname*{\bf{diag}}\left(\left(\mathcal{I}-\overline{\bm{D}}^{j}\right)\overline{\bm{D}}^{j}\right)$

Let $looR=looY-m\cdot loo\mathbbmss{1}$ and $A^{R}=A^{Y}-m\cdot
A^{\mathbbmss{1}}$ in terms of these. As before, the mean and noise variance
can be updated in closed form:

$\displaystyle\frac{dR^{(t)}_{e}}{dm}=0\Leftrightarrow
m\leftarrow\frac{\mathbbmss{1}_{n}^{\top}[loo\mathbbmss{1}\odot\bm{P}^{(t)}\odot
W\odot
looY]\mathbbmss{1}_{L}}{\mathbbmss{1}_{n}^{\top}[loo\mathbbmss{1}\odot\bm{P}^{(t)}\odot
W\odot loo\mathbbmss{1}]\mathbbmss{1}_{L}}$ (13)
$\displaystyle\frac{dR^{(t)}_{e}}{d{\sigma_{\varepsilon}^{2}}}=0\Leftrightarrow\sigma_{\varepsilon}^{2}\leftarrow\frac{1}{n}\mathbbmss{1}_{n}^{\top}[looR\odot\bm{P}^{(t)}\odot
W\odot looR]\mathbbmss{1}_{L}$ (14)

It remains to perform a gradient descent for $\lambda_{f}$, using

$\displaystyle\frac{dR^{(t)}_{e}}{d{\lambda_{f}}}=$
$\displaystyle\frac{1}{2\sigma_{\varepsilon}^{2}}\mathbbmss{1}_{n}^{\top}\Big{[}\bm{P}^{(t)}\odot
d^{A}\odot\big{[}\sigma_{\varepsilon}^{2}\oslash W+looR\odot looR\big{]}$
$\displaystyle-2\big{[}looR\odot\bm{P}^{(t)}\odot
A^{R}\big{]}\Big{]}\mathbbmss{1}_{L}.$ (15)

We summarized the SABER training procedure in Alg. 2 and predictions in Alg.
3.

Algorithm 2
$\text{train}_{\text{SABER}}(\bm{X},\bm{Y},\bm{X}_{\text{v}},\bm{Y}_{\text{v}},\bm{P},v,\Theta)$

1:[-0.4ex]0.490.4pt4pt 3pt

2:Input

3:0.490.4pt4pt 3pt

4:Training data
$\bm{X}=[x_{1},\ldots,x_{n}]^{\top}\in\operatorname*{\mathds{R}}^{n\times
d},\bm{Y}=[y_{1},\ldots,y_{n}]^{\top}\in\operatorname*{\mathds{R}}^{n}$

5:Validation data $\bm{X}_{\text{v}},\bm{Y}_{\text{v}}$

6:Initial bandwidths assignment probabilities $\bm{P}\in[0,1]^{n\times L}$

7:Local noise variance $v$ and ME hyperparameters
$\Theta=\\{\theta_{g},\theta_{e_{1}},\ldots,\theta_{e_{L}}\\}$ of the experts
$\theta_{e_{j}}=(m,\lambda_{f},\sigma_{\varepsilon}^{2},\Sigma_{j})$ and the
gate $\theta_{g}=(\Sigma_{g},\lambda_{g})$

8:[-0.4ex]0.490.4pt4pt 3pt

9:Output

10:0.490.4pt4pt 3pt

11:SABER model parameters $(\bm{P},\Theta)$

12:[-0.4ex]0.490.4pt4pt 3pt

13:Procedure

14:0.490.4pt4pt 3pt

15:Estimate
$\Sigma_{\varepsilon}^{2}=\operatorname*{\bf{diag}}\left(v(\bm{X})\right)$ and
kernel matrices $\bm{K}^{\Sigma_{j}}$ for $1\leq j\leq L$

16:Estimate eigendecompositions
$\Sigma_{\varepsilon}^{-1}\bm{K}^{\Sigma_{j}}\Sigma_{\varepsilon}^{-1}=\widetilde{U}_{j}\Lambda_{j}\widetilde{U}_{j}^{\top}$

17:Set $U_{j}=\Sigma_{\varepsilon}^{-1}\widetilde{U}_{j}$

18:while $\text{Err}_{\text{v}}$ decreases do

19: $\triangleright$ E-step

20: for $1\leq i\leq n$ do$\triangleright$ Fit $\bm{Q}$ to $\bm{P}$

21:
$\bm{\omega}_{-i}\leftarrow\text{train}_{\text{MKLR}}(\bm{X}_{-i},\bm{P}_{-i},\theta_{g})$
$\triangleright$ see Alg. 6

22:
$\bm{Q}[i,:]\leftarrow\text{predict}_{\text{MKLR}}(x_{i}|\bm{X}_{-i},\bm{\omega}_{-i},\theta_{g})$
$\triangleright$ see Alg. 5

23: $\begin{bmatrix}W,loo\mathbbmss{1},looY\\\
A^{\mathbbmss{1}},A^{Y},d^{A}\end{bmatrix}\leftarrow\text{helpVariables}(\lambda_{f},\bm{Y},[U_{j},\Lambda_{j}]_{j=1}^{L})$$\triangleright$
Alg. 1

24: Set $looR=looY-m\cdot loo\mathbbmss{1}$

25: Update $\bm{P}\leftarrow\bm{Q}\odot\sqrt{W}\odot\exp\left(-looR\odot
looR\odot W\middle/(2\sigma_{\varepsilon}^{2})\right)$

26: Normalize
$\bm{P}\leftarrow\operatorname*{\bf{diag}}(\bm{P}\mathbbmss{1}_{L})^{-1}\bm{P}$

27: $\triangleright$ M-step

28: while $\lambda_{f}$ not converged do

29: Set $A^{R}=A^{Y}-m\cdot A^{\mathbbmss{1}}$

30: Estimate $\frac{dR^{(t)}_{e}}{d{\lambda_{f}}}$ according to (III-B)

31: Update
$\lambda_{f}\leftarrow\lambda_{f}-\eta^{(t)}\frac{dR^{(t)}_{e}}{d{\lambda_{f}}}$
(for some decaying step size $\eta^{(t)}$)

32: Update $\begin{bmatrix}W,loo\mathbbmss{1},looY\\\
A^{\mathbbmss{1}},A^{Y},d^{A}\end{bmatrix}\leftarrow\text{helpVariables}(\lambda_{f},\bm{Y},[U_{j},\Lambda_{j}]_{j=1}^{L})$

33: Update $m$ according to (13)

34: Set $looR=looY-m\cdot loo\mathbbmss{1}$

35: Update $\sigma_{\varepsilon}^{2}$ according to (14)

36: Predict validation labels $\triangleright$ see Alg. 3
$\widehat{\bm{Y}}_{\text{v}}\leftarrow\text{predict}_{\text{SABER}}(\bm{X}_{\text{v}}|\bm{X},\bm{Y},\bm{P},v,\Theta)$

37: Estimate
$\text{Err}_{\text{v}}\leftarrow\|\bm{Y}_{\text{v}}-\widehat{\bm{Y}}_{\text{v}}\|_{\raisebox{-2.0pt}{\tiny$$}}^{2}$
$\triangleright$ Validation sum squared error

Algorithm 3
$\widehat{\bm{Y}}^{*}\leftarrow\widehat{y}_{\text{SABER}}(\bm{X}^{*}|\bm{X},\bm{Y},\bm{P},v,\Theta)$

1:[-0.4ex]0.490.4pt4pt 3pt

2:Input

3:0.490.4pt4pt 3pt

4:Training data
$\bm{X}=[x_{1},\ldots,x_{n}]^{\top}\in\operatorname*{\mathds{R}}^{n\times
d},\bm{Y}=[y_{1},\ldots,y_{n}]^{\top}\in\operatorname*{\mathds{R}}^{n}$

5:Test input data
$\bm{X}^{*}=[x_{1}^{*},\ldots,x_{T}^{*}]^{\top}\in\operatorname*{\mathds{R}}^{T\times
d}$

6:Training bandwidths assignment probabilities $\bm{P}\in[0,1]^{n\times L}$

7:Local noise variance $v$ and ME hyperparameters
$\Theta=\\{\theta_{g},\theta_{e_{1}},\ldots,\theta_{e_{L}}\\}$ of the experts
$\theta_{e_{j}}=(m,\lambda_{f},\sigma_{\varepsilon}^{2},\Sigma_{j})$ and the
gate $\theta_{g}=(\Sigma_{g},\lambda_{g})$

8:[-0.4ex]0.490.4pt4pt 3pt

9:Output

10:0.490.4pt4pt 3pt

11:Predictions $\widehat{\bm{Y}}^{*}\in\operatorname*{\mathds{R}}^{T}$

12:[-0.4ex]0.490.4pt4pt 3pt

13:Procedure

14:0.490.4pt4pt 3pt

15:$\bm{\omega}\leftarrow\text{train}_{\text{MKLR}}(\bm{X},\bm{P},\theta_{g})$
$\triangleright$ see Alg. 6

16:$\bm{Q}^{*}\leftarrow\text{predict}_{\text{MKLR}}(\bm{X}^{*}|\bm{X},\bm{\omega},\theta_{g})$
$\triangleright$ see Alg. 5

17:for $1\leq j\leq L$ do

18:
$\widehat{\bm{Y}}^{j}\leftarrow\text{predict}_{\text{GPR}}(\bm{X}^{*}|\bm{X},\bm{Y},v,\theta_{e_{j}})$
$\triangleright$ see Alg. 4

19:Predict
$\widehat{\bm{Y}}^{*}\leftarrow\left[\bm{Q}^{*}\odot[\widehat{\bm{Y}}^{1},\ldots,\widehat{\bm{Y}}^{L}]\right]\mathbbmss{1}_{L}$

The learnt gate $Q$ contains the information of spatial inhomogeneity, that we
propose to summarize according to

$\displaystyle\Sigma_{\text{\tiny{SABER}}}(x)$
$\displaystyle=\exp_{M}\operatorname*{\mathds{E}}\left[\log_{M}(\bm{\Sigma}_{\text{\tiny{Opt}}})|x\right]$
(16)
$\displaystyle=\exp_{M}\left(\sideset{}{{}_{j=1}^{L}}{\sum}Q(x)_{j}\log_{M}(\Sigma_{j})\right),$

where $\log_{M}$ is the inverse of the matrix exponential $\exp_{M}$.

###  III-C Hyperparameter Initialization and Limitations of SABER

The initialization of the assignment matrix $\bm{P}^{(0)}$ and the shared
expert parameters $(m,\lambda_{f},\sigma_{\varepsilon}^{2})^{(0)}$ of SABER is
more crucial than one might expect at first glance: Given a reasonable initial
$\bm{P}^{(0)}$, we could perform a standard gradient descent on
$\ell^{(0)}(\bm{X},\bm{Y},\Theta)$ from (12) to solve for optimal, initial
$(m,\lambda_{f},\sigma_{\varepsilon}^{2})^{(0)}$. However, if we were to know
a reasonable $\bm{P}^{(0)}$ in advance, the overall problem would already be
solved.

We suggest a conservative initialization: We tune each expert with respect to
its individual set $(m,\lambda_{f}^{j},\sigma_{\varepsilon}^{2},\Sigma_{j})$
as described in Sec. II-B in advance, and define
$\lambda_{f}^{*}=\min_{j}\lambda_{f}^{j}$. We then set
$\bm{P}^{(0)}\equiv\frac{1}{L}$ and
$(m,\lambda_{f}^{*},\sigma_{\varepsilon}^{2})^{(0)}$ with the optimal $m$ and
$\sigma_{\varepsilon}^{2}$ given by (13) and (14), based on $\bm{P}^{(0)}$ and
$\lambda_{f}^{*}$.

In case of heteroscedasticity we assume the true—or an estimate of—local noise
variance to be given in advance, which we simply plug in, as this is not the
focus of our work. Such an estimate could be obtained by [5, 6, 29].

Obviously, there are lots of degrees of freedom in choosing appropriate
candidate bandwidth matrices $\Sigma_{1},\ldots,\Sigma_{L}$ and $\Sigma_{g}$
for the experts and the gate, respectively. In the following, we restrict
ourselves to bandwidth matrices of the form $\Sigma_{j}=\sigma_{j}\Sigma$ and
$\Sigma_{g}=\sigma_{g}\Sigma$ for a shared bandwidth matrix $\Sigma$, where
the scalar-valued $\sigma_{j}$ model spatial inhomogeneity.

In doing so, we drastically reduce the degrees of freedom making the approach
robust enough to cope with real-world problems. In addition, we obtain a
straightforward construction plan for the bandwidth matrices: In a
preprocessing step, we estimate $\Sigma$ via tuning a global GPR model, as
explained in Sec. II-B. Then we set up $\sigma_{1}<\ldots<\sigma_{L}$
logarithmically spaced around $1$ to construct the expert bandwidth matrices.

Note that it is possible to refit the shared bandwidth $\Sigma$ while training
of SABER, which is a computational disadvantage: We would have to repeatedly
decompose the kernels.

Basically we were able to tune the gate hyperparameters $\theta_{g}$,
optimizing the fit to the current assignment $\bm{P}$ in each iteration. This
is however a pitfall: As soon as $\bm{P}$ starts to converge, the optimal
$\theta_{g}$ tends to 0, leading to a perfect reconstruction of $\bm{P}$. This
is, because $\bm{P}$ is itself a prediction of the gate from the last
iteration. On the other hand, before $\bm{P}$ converges (and becomes
meaningful with respect to the problem) it is not reasonable to force the gate
to predict $\bm{P}$ well. Given the initial gate hyperparameters $\theta_{g}$,
we keep $\theta_{g}^{(t)}=\theta_{g}$ fixed while the training of SABER.

The initial $\lambda_{g},\sigma_{g}$ and $\sigma_{1},\ldots,\sigma_{L}$ are
the hyperparameters that remain to be chosen, minimizing the validation error.

###  III-D Local Function Complexity Estimation

Inspired by the theory of LLS, in order to obtain an estimate for LFC, our
LOKB estimate needs to be adjusted for the effects of local noise variance $v$
and input density $p$. We assume $p$ to be given or estimated by standard
methods (see e.g. [30]) in advance. We base our next statements on the
experimental results in Sec. V-A1.

In GPR we can treat heteroscedasticity in the regularization, as opposed to
LLS. More precisely, the GPR model is misspecified if we apply its
homoscedastic version to a factually heteroscedastic problem. Nevertheless, we
found empirical evidence that SABER will compensate for such misspecification
such that $\Sigma_{\text{\tiny{SABER}}}$ holds proportionality in $v$
analogous to LLS. But for obvious reasons it is more robust to apply
heteroscedastic GPR right from the start. As expected,
$\Sigma_{\text{\tiny{SABER}}}$ is then already adjusted for $v$.

We furthermore observe proportionality of $\Sigma_{\text{\tiny{SABER}}}$ in
$p$ analogously to LLS. Following (5), we propose our estimate for LFC via

$\displaystyle\mathfrak{C}(x)=\left|\Sigma_{\text{\tiny{SABER}}}(x)\right|^{-\frac{1}{d}}p(x)^{-\frac{1}{4+d}}.$
(17)

##  IV Related Work

Several methods were proposed for better prediction performance in
inhomogeneous data. In its basic form MKL prepares complementary kernels
$\\{k^{1},\dots,k^{L}\\}$, e.g., an RBF-kernel at different scales, and use
the mixed kernel $k_{\mu}=\sideset{}{{}_{j=1}^{L}}{\sum}\mu_{j}k^{j}$, where
the mixing weight vector $\bm{\mu}$ along with the dual variables
$\bm{\alpha}$ is optimized in a single framework. Although MKL copes with
inhomogeneity to some extent, this approach uses a single mixed kernel
globally over the input space and is therefore too restrictive. There exist
more sophisticated versions of MKL (see [31] for a summary). For example [25]
proposed LMKR—a localized MKL approach—where a gate $\eta(x)$ is learnt to
obtain a locally weighted kernel mixture

$\displaystyle
k_{\eta}(x,x^{\prime})=\sideset{}{{}_{j=1}^{L}}{\sum}\eta_{j}(x)k^{j}(x,x^{\prime})\eta_{j}(x^{\prime}),$

which is plugged in to the standard MKL model:

$\displaystyle\widehat{f}(x)=\sideset{}{{}_{i=1}^{n}}{\sum}\alpha_{i}k_{\eta}(x_{i},x)+b.$
(18)

Letting $k^{j}=k^{\sigma_{j}}$ be RBF-kernels at different bandwidths, we
obtain a localized bandwidth estimate similar to
$\Sigma_{\text{\tiny{SABER}}}(x)$ in (16) with $Q\equiv\eta$. In fact, we
considered this approach before SABER. Unfortunately, it heavily suffers from
local optima even on simple toy data, as we will show in Sec. V-A2. This makes
the approach infeasible for real-world problems. Learning such a local
weighting makes LMKR an instance of ME—the concept that we also pursue in this
paper (see [2, 3, 4]). [32] give an introduction and summary of ME approaches.
A fundamental assumption of the ME model 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 ME model can specialize on
modelling 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 it suggests itself to tune the hyperparameters of each
expert individually on the respective assigned data subset. In the light of
this paradigm there exist several instances of mixture of Gaussian processes,
for example [33, 34, 35, 36, 37].

In this work, we break with this paradigm assuming that each single expert is
capable of modelling the whole problem, given enough training data. That is,
we drop the incompatibility assumption. Yet we assume that—locally
dependent—some experts perform superior compared to the others. While we still
try to soft-partition the data with respect to the well performing experts, we
let all experts share one set of hyperparameters, since they all model the
same function.

For the LLS model there is an asymptotic closed form solution
$\sigma_{\mathrm{theo}}$ to LOKB, which [23] attempt to estimate. We refer to
the introduction, why we refrain from this approach.

There have been attempts to learn the model

$\displaystyle\widehat{f}(x)=\sideset{}{{}_{i=1}^{n}}{\sum}\sideset{}{{}_{j=1}^{L}}{\sum}\alpha^{j}_{i}k^{j}(x_{i},x)$

for individual coefficient vectors $\alpha^{j}$ corresponding to a single
kernel $k^{j}$, where $k^{j}$ may for example be a Gaussian kernel with
bandwidth $\sigma_{j}$: _Multiscale support vector regression_ (MS-SVR) [11]
minimizes a loss function similar to SVR with stacked constraints for each
candidate bandwidth. In _kernel basis pursuit_ (KBP) [12], features
$k^{j}(x_{i},\cdot)$ are selected iteratively updating the prediction. In an
iteration one selects the feature that is most correlated with the current
residuals. While [11] have shown that KBP is inferior to MS-SVR, we will show
that MS-SVR is inferior to SABER in Sec. V-A2. _Hierarchical SVR_ [13, 14] is
another approach, where several single-scale SVR models are trained, beginning
with a coarse bandwidth, followed by successively fitting the current
residuals with a finer bandwidth. These approaches provide no access to a
local bandwidth estimator. This is mainly because neither the amplitude of
coefficients $\alpha^{j}_{i}$ nor the positions of support vectors with
respect to a bandwidth $\sigma_{j}$ do necessarily match the region in the
input space where this bandwidth is locally optimal.

One may also tune the bandwidths and/or centers of RBF-networks individually,
as done in [19], _sparse multiscale Gaussian process regression_ [20], in [22]
and in the adaptive spherical Gaussian kernel _relevance vector machine_
approach by [21]. Here, no functional relationship is imposed on LOKB,
disqualifying these approaches for our purpose.

##  V Experiments

First we analyze the behavior of $\Sigma_{\text{\tiny{SABER}}}$ with respect
to changes in noise and input density, and show its superior performance to
related approaches on toy data. Then we demonstrate the capabilities of our
LFC estimate to gain insights into data from quantum chemistry, and to serve
as an input to active learning for the simulation of a physical system. We
assume homoscedasticity whenever not stated differently.

###  V-A Oscillating function toy data

Let us consider a regression problem to fit samples generated by the following
generative process:

$\displaystyle p(y|x)=\mathcal{N}(y;f(x),v(x)),\quad x\in[0,10],$
$\displaystyle f(x)=\sin\left(2\pi(10c_{1}+10)\middle/(x\cdot
c_{1}+10)\right),$ $\displaystyle
p(x)\propto(10\cdot(c_{2}+1)-c_{2}\cdot((x+3.5)\bmod 10)),$
$\displaystyle\log_{10}v(x)=-2+\frac{c_{3}}{2}(1-\cos(2\pi(x-5)/5)),$

where $\mathcal{N}(\cdot;\mu,\sigma^{2})$ denotes the Gaussian distribution
with mean $\mu$ and variance $\sigma^{2}$. The signal $f$ is a down-chirping
oscillating function, where $c_{1}$ specified the number of periods and
therefore controls the inhomogeneity in function complexity. The skewness of
the input distribution $p$ is given by $c_{2}$, whereas the degree of
heteroscedasticity is given by $c_{3}$.

This dataset is an extension to an experiment of [11] such that we can compare
our approach to their results, later on.

####  V-A1 Proportionalities of SABER

Fixing $c_{1}=2.5$, we show the individual inhomogeneities in Fig. 3.

Figure 3: The oscillating toy data. From top to bottom, we show the absolute
second derivative of the true function to infer, and the input density and
local noise variance in their homogeneous and heterogeneous states,
respectively. Figure 4: Example oscillating toy datasets with varying degree
of inhomogeneity. (top) All aspects are inhomogeneous. (middle) Removing
heteroscedasticity, (bottom) Additionally removing input density skew.

Starting from a specification that is inhomogeneous in all regards (a), we
successively remove heteroscedasticity (b), then input density skew (c), with
exemplary datasets given in Fig. 4. Note that the initial case corresponds to
the example in the introduction.

In this experiment we show how $\sigma_{\text{\tiny{SABER}}}$ behaves with
respect to changes in $v$ and $p$. We use the hyperparameter configuration
$\sigma_{j}=0.2\cdot 10^{\frac{j-3}{6}},1\leq j\leq 11$, $\sigma_{g}=5$,
$\lambda_{g}=5\cdot 10^{-4}$. The resulting bandwidth functions are based on
$1000$ training samples and averaged over 20 repetitions. Here we specify
explicitly the $v$ that is provided to SABER beforehand:

$\displaystyle\Sigma_{\text{\tiny{SABER}}}^{v}(x)$
$\displaystyle=\exp_{M}\operatorname*{\mathds{E}}\left[\log_{M}(\bm{\Sigma}_{\text{\tiny{Opt}}})|x,v\right]$

Therefore by $\Sigma_{\text{\tiny{SABER}}}^{\bm{1}}$ we enforce homoscedastic
treatment.

Figure 5:  Bandwidth estimates for the three oscillating toy data settings.
Bandwidths are given relative (top) and absolute (middle), when removing
heteroscedasticity. (bottom) Relative bandwidths, when removing input skew.

Since LLS and kernel-linear models appear related, we expect LOKB to behave
similar in both cases. More precisely, based on (3) we expect the relation

$\Sigma_{\text{\tiny{SABER}}}^{\bm{1}}(x)\propto\left[\frac{v(x)}{p(x)}\right]^{\frac{1}{4+d}}\mathfrak{S}^{n}(x).$

Let us first consider changes in local noise variance $v$: We compare
$\sigma_{\text{\tiny{SABER}}}^{\bm{1}}$ in setting (a) to (b)—which only
differ in being heteroscedastic or not—by building their ratio as can be seen
to the top of Fig. 5. The theory of LLS suggests
$\sigma_{\text{\tiny{Opt}}}^{\text{\tiny{LLS}}}$ between setting (a) and (b)
to differ proportional to $v^{1/5}$, which we added to this plot. We observe a
significant alignment of the ratio to the expected proportionality, which is
however not perfect being owed to the difficulty of the dataset.

In the middle of Fig. 5 we show $\sigma_{\text{\tiny{SABER}}}^{\bm{1}}$ in
setting (a) and (b), and $\sigma_{\text{\tiny{SABER}}}^{v}$ in setting (a).
First of all, bandwidth estimates in setting (b) are smaller than the
respective in setting (a) since the noise variance in setting (b) is equal or
less everywhere, compared to setting (a). As just discussed,
$\sigma_{\text{\tiny{SABER}}}^{\bm{1}}$ shows proportionality in $v$. However,
setting (b) is homoscedastic and therefore naturally adjusted for noise. More
interestingly, we observe that $\sigma_{\text{\tiny{SABER}}}^{v}$ is
proportional to $\sigma_{\text{\tiny{SABER}}}^{\bm{1}}(b)$, implying that
$\sigma_{\text{\tiny{SABER}}}^{v}$ is invariant with respect to $v$. We deduce

$\Sigma_{\text{\tiny{SABER}}}^{v}(x)\propto
p(x)^{-\frac{1}{4+d}}\mathfrak{S}^{n}(x).$

Let us now consider changes in the input density $p$: We compare
$\sigma_{\text{\tiny{SABER}}}^{\bm{1}}$ in setting (b) to (c)—which only
differ in featuring an inhomogeneous or uniform input distribution—by building
their ratio as can be seen to the bottom of Fig. 5. The theory of LLS suggests
$\sigma_{\text{\tiny{Opt}}}^{\text{\tiny{LLS}}}$ between setting (b) and (c)
to differ proportional to $(1/p)^{1/5}$, which we added to this plot. We also
plot the ratio of $\sigma_{\text{\tiny{SABER}}}^{v}$ in the harder setting (a)
to $\sigma_{\text{\tiny{SABER}}}^{\bm{1}}$ in setting (c) in green. Again, we
observe a significant alignment of the ratios to the expected proportionality.

####  V-A2 Comparison to Related Work

Figure 6: The oscillating dataset (top), LOKB estimates (middle) and the gate
responses of SABER (bottom). Predictions are based on 100 training examples.

With the specification $(c_{1},c_{2},c_{3})=(3.5,0,0)$ the dataset can be
reduced to an experiment in [11]. In this setting there is inhomogeneity in
function complexity only. To the top of Fig. 6 we show an example dataset.
Besides an improved performance over [11], SABER is able to provide a
structural analysis of the inhomogeneity: As can be seen in Fig. 6 (middle),
the LOKB estimates increase monotonically in $x$, matching the intuition.

They generate 30 datasets $(\bm{X}_{(i)},\bm{Y}_{(i)})$ of size 100 with
corresponding validation datasets given by $(\bm{X}_{(i)},\bm{f}_{(i)})$ with
$\bm{f}_{(i)}=f(\bm{X}_{(i)})$. Then they select hyperparameters by minimizing
the validation _root mean squared error_ (RMSE) over the first five datasets
and report performance as the average validation RMSE over all 30 datasets. We
evaluate the following experiments on the oscillating function toy data
accordingly, reporting the performance in Tab. I.

In MS-SVR [11], they use four bandwidths logarithmically spaced between
$\sigma_{1}=0.15$ and $\sigma_{4}=1.4$. However, they lack an analysis of a
spatial bandwidth segmentation. Now that our gate $Q(x)$ is not a static
hyperparameter which we can perform a grid search on, we fit the gate and the
expert hyperparameters only on the first dataset minimizing its validation
RMSE:

For several hyperparameter candidates we perform

$\text{train}_{\text{SABER}}(\bm{X}_{(1)},\bm{Y}_{(1)},\bm{X}_{(1)},\bm{f}_{(1)},\bm{P},v,\Theta)$

according to Alg. 2, and found $\sigma_{j}=\sigma_{\text{\tiny{GPR}}}\cdot
10^{\frac{j-3}{3}},1\leq j\leq 7$, $\sigma_{g}=5$, $\lambda_{g}=5\cdot
10^{-4}$ to work well. The global GPR model within the training procedure was
tuned to $\sigma_{\text{\tiny{GPR}}}=0.45$. As can be seen in Fig. 6 (bottom),
$\sigma_{2}$ to $\sigma_{5}$ are mainly used, effectively resulting in a
bandwidth range from $0.21$ to $2.1$. The expert hyperparameters are
$(m,\lambda_{f},\sigma_{\varepsilon}^{2})=(-2.7317,5.1369,0.0896)$.

We then build the SABER predictions similar to Alg. 3, reusing the gate
$\bm{Q}_{(1)}=\text{predict}_{\text{MKLR}}(\cdot|\bm{X}_{(1)},\bm{\omega},(\sigma_{g}\sigma_{\text{\tiny{GPR}}},\lambda_{g}))$
for
$\bm{\omega}=\text{train}_{\text{MKLR}}(\bm{X}_{(1)},\bm{P},(\sigma_{g}\sigma_{\text{\tiny{GPR}}},\lambda_{g}))$:

$\widehat{\bm{f}}_{(i)}=\left[\bm{Q}_{(1)}(\bm{X}_{(i)})\odot[\widehat{\bm{Y}}^{1}_{(i)},\ldots,\widehat{\bm{Y}}^{L}_{(i)}]\right]\mathbbmss{1}_{L}$

with
$\widehat{\bm{Y}}^{j}_{(i)}=\text{predict}_{\text{GPR}}(\bm{X}_{(i)}|\bm{X}_{(i)},\bm{Y}_{(i)},v,(m,\lambda_{f},\sigma_{\varepsilon}^{2},\sigma_{j}\sigma_{\text{\tiny{GPR}}}))$.
This gives a significantly better performance than MS-SVR.

As already mentioned in Sec. IV, we could have used LMKR based on (18) by [25]
for our purpose of both, treating inhomogeneities as well as modelling them
explicitly. We will however show with this controlled toy example that the
objective of LMKR suffers from local optima to such a degree which makes it
infeasible to use for real world problems. Let:

$\displaystyle\eta(x)_{j}=\frac{\exp\left\\{-\|s_{j}^{-1}(x-c_{j})\|_{\raisebox{-2.0pt}{\tiny$2$}}^{2}\right\\}}{\sideset{}{{}_{k=1}^{L}}{\sum}\exp\left\\{-\|s_{k}^{-1}(x-c_{k})\|_{\raisebox{-2.0pt}{\tiny$2$}}^{2}\right\\}}.$

In a first part we incorporate the prior knowledge of the optimal bandwidth
function to be increasing into the gate initialization. That is, we match the
random centers $c_{j}$ to the bandwidth candidates $\sigma_{j}$ such that for
$\sigma_{j}>\sigma_{k}$ it is $c_{j}>c_{k}$. We do so for 100 random gate
initializations. Beyond that, we also tune the hyperparameters $C$ and
$\varepsilon$ of the SVR model, the number of bandwidth candidates $L$ and the
lower and upper bound $\sigma_{1}$ and $\sigma_{L}$. Analogously to SABER we
chose the best set of hyperparameters (and learnt gate) on the first dataset,
being $C=100$, $\varepsilon=0.1$, $\sigma_{j}=10^{\frac{j-2}{2}},1\leq j\leq
7$.

The achieved performance is non-significantly better than SABER. But when
inspecting the so learnt gate, we observe that the corresponding initial gate
was already almost optimal in first place. This means that, if we know the
perfect solution in advance, we achieve superior performance, non-
surprisingly.

In the second part, we reuse the optimal hyperparameters, but start from
random gate initializations without the encoded prior knowledge about the
order. In Fig. 7 we plot the number of repetitions against the so far best
achieved performance, marking the performance of SABER and optimal LMKR.

Figure 7: LMKR RMSE plotted against number of random gate initializations.

We observe that LMKR beats SABER after about 9000 gate initializations. Note
that SABER does not use the degree of freedom to start from a random
assignment $P$—even though possible. Therefore the performance of SABER is
obtained in just one initialization. We claim that this problem of LMKR will
be amplified in less well-behaved, high dimensional settings.

Finally, we compare SABER to the LLS model, once using the theoretically
optimal bandwidth function $\sigma_{\mathrm{theo}}$ from (4) as an oracle,
once using $\widehat{\sigma}_{\mathrm{theo}}$ via estimation of
$f^{{}^{\prime\prime}}$ from the data and once plugging in
$\sigma_{\text{\tiny{SABER}}}$. We have plotted the applied bandwidth
functions in Fig. 6 (middle) and report the resulting performance in Tab. I.

TABLE I: Results on the oscillating function toy data experiment. The performance of MS-SVR is taken from [11]. Significantly best performances (with 95% confidence) are denoted in boldface. Method | RMSE
---|---
MS-SVR | $0.0515\pm 0.0072$
SABER | $\bm{0.0464\pm 0.0025}$
LMKR | $\bm{0.0437\pm 0.0030}$
LLS \+ $\sigma_{\mathrm{theo}}$ | $0.0570\pm 0.0029$
LLS \+ $\widehat{\sigma}_{\mathrm{theo}}$ | $0.0606\pm 0.0030$
LLS \+ $\sigma_{\text{\tiny{SABER}}}$ | $0.0530\pm 0.0025$

In accordance to [18] we observe suboptimality of $\sigma_{\mathrm{theo}}$ to
a bandwidth function estimate on finite data such as
$\sigma_{\text{\tiny{SABER}}}$. Next, we observe that approximating
$\sigma_{\mathrm{theo}}$ from the data is already instable in this rather
simple experiment on the one hand, whereas it suffers from the (true) poles of
$\sigma_{\mathrm{theo}}$ on the other hand.

###  V-B Malonaldehyde Molecular Dynamics

Figure 8: The two stable states of the malonaldehyde molecule.

The dataset by [38], [39] is a molecular dynamics (MD) simulation of
malonaldehyde—a molecule with chemical formula $C_{3}O_{2}H_{4}$, which
structure can be seen in Fig. 8. This dataset consists of Cartesian geometries
$(R_{k})_{k=1}^{9}$ with corresponding forces and potential energies, that we
use as our labels. The inputs are formed by the reciprocal distances
$\|R_{k}-R_{l}\|_{\raisebox{-2.0pt}{\tiny$$}}^{-1}$ of all atom pairings—in
total 36 features—as suggested by [38]. We visualize both, featurewise and
spatial adaption of bandwidth, which reproduce physically meaningful behavior.
Moreover, SABER raises prediction performance compared to the global energy
based model.

We conduct the experiment on the unrelaxed configurations of the true MD
simulation and illustrate our results, using a physically plausible projection
of the 27-dimensional potential energy surface (PES): When interpreting the
two aldehyde groups as rotators around the axes given by the respective carbon
bonds, their corresponding angles are the dominant drivers of the potential
energy. When relaxing the molecule while keeping these two angles fixed, we
obtain an interpretable projection of the PES to the left of Fig. 10.

Figure 9: Exemplaric malonaldehyde molecule conformation. (Left): Atomic
Bonds. (Right): Atomic pair importance depicted by line saturation. Figure 10:
Potential energies (left), estimated input density (middle) and local function
complexity (right) of the relaxed malonaldehyde configurations. Surfaces are
plotted with respect to the angles of the two aldehyde rotators of
malonaldehyde.

We apply ARD to determine the importance of specific atom pairs, which can be
seen in Fig. 9. Note that we take the symmetry of the molecule into account
here, such that we consider 17 feature groups instead of 36 features. We
observe that the features describing the geometry within each local, carbon
centered environment (that is, the aldehyde groups and the central $CH_{2}$)
have negligible influence on the energy prediction. In contrast ARD affirms
the geometry between these local environments to be most important which
matches insights from the physical point-of-view. Via tuning on a validation
set we chose the following SABER hyperparameters: $\sigma_{g}=0.5$,
$\lambda_{g}=5\cdot 10^{-3}$ and $\sigma_{j}=10^{\frac{j-3}{6}},1\leq j\leq
7$.

TABLE II: Results on the atomization energy prediction for the true malonaldehyde MD simulation with reference values taken from [38] (denoted by ref). Estimates are based on a computationally comparable number of configurations. | Performance |
---|---|---
Method | RMSE | MAE | max AE | #Configurations
GPR | 0.1738 | 0.1201 | 2.4735 | 12k
GPR \+ ARD | 0.1423 | 0.0965 | 2.1928 | 12k
SABER \+ ARD | 0.1357 | 0.0934 | 1.8658 | 12k
$\text{GDML}^{\textit{ref}}$ | 0.25 | 0.16 | — | 1k
$\text{GPR}^{\textit{ref}}$ | 0.16 | 0.11 | — | 27k

Tab. II shows the prediction performance for the energies in the malonaldehyde
MD dataset by [38] for a computationally comparable amount of configurations.

[38] proposed a gradient-domain machine learning (GDML) which is basically a
vector-valued GPR model on the derivatives of the PES rather than the energy
itself.

Since there are 27 (= 9 atoms in 3D) such derivatives, their model complexity
grows 27 times faster. This is why they suggested to compare their model using
1k configurations to energy based GPR using 27k configurations. We determined
the computationally comparable sample size for SABER to be 14k, in the sense
that $(27k)^{3}>7\cdot(14k)^{3}$—comparing the leading computationally complex
kernel matrix inversion operation which is required 7 times more often in
SABER than in GPR. Due to RAM limitations we restricted our experiment further
to 12k configurations. Non-surprisingly the reference GPR at 27k performs
better than our GPR at 12k. However, we observe an overall boost in
performance when taking featurewise importance into account, dropping
significantly below the reference errors—even using far less samples.

When inspecting LFC of the relaxed malonaldehyde configurations to the right
in Fig. 10, we observe local minima at the stable configurations of
malonaldehyde, whereas complexity grows in the transition areas of the PES.
This is also what we intuitively expect, as the energy surface in the
transition areas is more volatile compared to the stable valleys.

Treating spatial inhomogeneity has, on first glance, less impact on the
averaging measures RMSE and MAE. This is however due to the skew of the input
distribution: Most of the samples reside in the stable energy sinks whereas
transitions between those happen rarely. While the average error of a
predictor might be low, running an MD simulation with it still can
fail—meaning that the true physics are not resembled—when the performance in
the rare transition areas is poor. SABER mainly enhances prediction in these
transition areas, which is not reflected in a significant drop of RMSE and
MAE, but in the maximum AE.

###  V-C Atomization Energies of GDB-9

Now we conduct experiments on the Many-Body Descriptors for prediction of
atomization energies, discussed by [40]. These describe interactions between
sub-groups of two and three atoms, totalling 3270 features. The GDB-9 dataset
constitute the database of all stable and synthetically accessible organic
molecules composed of (C,N,O,H) with up to 9 heavy atoms (C,N,O). GDB-9 is a
subset of the GDB-17 [41] database, totalling 131722 molecules. Atomization
energies of the molecules have been calculated using hybrid density functional
theory at the B3LYP/6-31G(2df,p) level of quantum chemistry [42]. Fig. 11
(top) shows the distribution of energies.

(a) Distribution of the atomization energy. (b) Distribution of the local
function complexity.

Figure 11: Distribution of properties of the GDB-9 dataset.

[40] combined these features with a global GPR model. It is therefore
straightforward to relate (and extend) their approach using SABER, which
results in a strong performance boost. When inspecting the LFC estimates, we
identify a significant correlation to high-level features from a chemical
point-of-view.

We base our further analysis on 10k training molecules, implementing isotropic
bandwidths due to the large amount of features. Note that building feature
subsets for ARD might be considered, for example, by segmenting the features
into two groups of two-body, respectively three-body interactions or 36 groups
of specific atom-combination interactions.

The following SABER hyperparameters performed well on validation data:
$\sigma_{j}=10^{\frac{j-3}{3}},1\leq j\leq 7$, $\sigma_{g}=50$,
$\lambda_{g}=10^{-3}$. We observe a large spread of LFC estimates (see Fig.
11, bottom). Their meaningfulness is underpinned by an enormous performance
boost, as can be seen in Tab. III.

TABLE III: Average achieved errors for the atomization energy prediction of the GDB-9 dataset, given in kcal/mol. Estimates are based on 10k training samples. The reference value for SchNet is taken from ([43], Figure 3). | Performance
---|---
Method | RMSE | MAE | max AE
MB + GPR | 4.7987 | 1.4474 | 230.83
MB + SABER | 2.3809 | 1.2759 | 72.567
SchNet | — | 1.2791 | —

When comparing with the state-of-the-art method SchNet [43], we are on par.
This emphasizes the reasonability of our result. In addition, we can analyze
the bandwidth function to gain deeper insights into the problem:

(a) Molecule at high LFC.

(b) Molecule at low LFC.

Figure 12: Exemplaric molecules of the GDB-9 dataset with extremal LFC.

Inspecting molecules with high and low LFC, the molecules shown in Fig. 12 are
prototypical: The high LFC molecules form clusters of carbon which we try to
quantify by the total number of unique cycles contained in the molecule. On
the other hand, the low LFC molecules contain specific chemical functional
groups such as the nitro, azo and oxime group.

Figure 13: Correlation between LFC and the number of functional groups (left)
and the number of cycles (middle) for the GDB-9 dataset. Right: A log-linear
fit of the two features (number of functional groups, number of cycles) to
LFC.

In Fig. 13 we scatter the number of cycles, the number of functional groups
and a log-linear fit of these two features to LFC against the latter. The
correlation between LFC and its fit is strong, indicating that we identified
the correct sophisticated molecular properties which influence the complexity
of energy prediction. Note that the coefficients of the two features in the
fit have different signs. Therefore the two properties have an antagonizing
effect.

We interpret this as follows: Functional groups have itself a strong impact on
the energy of a molecule. Hence, there are several degrees of
freedom—performing changes on everything else but the functional groups—that
don’t affect the energy fundamentally. On the other hand, performing any
change in a heavy cyclic carbon configuration has strong impact on its energy
due to the high degree of interconnectivity.

###  V-D Flow in a square cavity with differentially heated side walls

In this final experiment, we will show how $\Sigma_{\text{\tiny{SABER}}}$ can
be integrated in an active learning scenario on the simulation of a physical
system that is described by partial differential equations. Since there is no
closed form solution to the system, we need to solve the problem numerically
over a mesh. A naive, uniform mesh requires an enormously large set of
simulation points in order to resemble the true physical system. Given a fix
budget, we will present a redistribution procedure of the mesh points, letting
the simulation approach the physical system with each iteration.

We feature a two dimensional test case from fluid mechanics, where variation
of spatial resolution is often required, because of different scales present
in the solution. Although defined by by a trivial domain and simple boundary
conditions, the setup produces a rather complex flow field. The test case is
known as square cavity with differentially heated side walls, which was first
proposed for incompressible flow by [44] and was extended to compressible flow
at large temperature differences by [45].

The domain is a square with four no-slip walls and subject to a gravitational
potential. The top and bottom walls (w.r.t. the gravitational potential) are
insulated. Left and right walls feature Dirichlet boundary conditions for the
temperature. The left wall is hot (i.e. $T_{h}$) and the right wall is cold
(i.e. $T_{c}$). The expansion of the heated fluid and the contraction of the
cooled fluid on left and right walls, respectively, will yield natural
convection, i.e. rise of hot fluid and descend of cold fluid. The simulation
is performed with a large temperature difference of $T_{h}/T_{c}=4$. The most
influential dimensionless number for convection flows is the Rayleigh number

$\displaystyle Ra=\dfrac{Pr~{}g~{}H^{3}~{}(T_{h}-T_{c})}{\nu^{2}}=10^{7},$

where $Pr$ is the Prandtl number, $g$ the gravitational acceleration, $H$ the
side length of the square domain and $\nu$ the kinematic viscosity. The
Rayleigh number is chosen as $10^{7}$, which is the highest value for which a
stationary solution is reported in the reference [45]. A visualization of the
reference flow field at this Rayleigh number based on [46] is shown in Fig.
14. Further, dimensionless numbers are the Barometric number $Ba$ as $0.01$,
$Pr=1.0$ and the number of internal degrees of freedom $K=3$.

Figure 14: Reference flow field of the square cavity with differentially
heated walls, based on [46] (with permission). The reference is based on the
high resolution simulation from [45].

For the simulation a two-dimensional gas kinetic scheme (GKS) [47, 48] is
applied, which was previously shown to yield high quality results for thermal
compressible flows and especially the square cavity with differentially heated
side walls [46]. The GKS is implemented on Cartesian meshes, which are
beneficial for automated mesh generation. For mesh refinement the quad tree
type method proposed by [49] is used.

In this scenario we are not interested in prediction performance, but the
resemblance of the true physical system, which we measure qualitatively by
appearance and quantitatively by convergence of observed vortex centers in the
velocity field to the ground truth. We assume this ground truth to be the high
resolution simulation in [45] based on $1,048,576=1024\times 1024$ cells,
resulting in a total of eight vortices, which are also shown in Fig. 14. The
centers of the realized vortices are measured as intersection points of zero
velocity contours for $x$\- and $y$-velocity.

The measured centers by the reference [45] and in our experiment are listed in
Tab. V in Appendix F, where we also specify the Euclidean distance $e$ of the
individual vortex centers to the reference. We designate the L2-norm $e_{L2}$
of these $e$ to be the quantitative performance measure. The characteristics
of the compared meshes are summarized in Tab. IV.

The initial mesh has a uniform background resolution of $64\times 64$ cells
with two thin refinements towards the walls. These refinements are obtained
when one has no prior knowledge on the flow field, apart from that the
temperature and velocity gradient at the wall will be large. This base line
mesh contains $8608$ cells. This resolution turns out to be insufficient as
the vortices 5 and 6 are not reproduced at all. Here it is to note, that this
resolution is very low compared to the $40,000$ cells in [46] and $1,048,576$
cells in [45].

In the following process these cells are redistributed in space and over the
refinement levels to improve the accuracy of the computed flow field. We
alternate in a loop between the following steps:

  1. 1.

Given the current mesh we run the simulation to convergence, which generates a
dataset $(X^{(i)},y^{(i)})$ with mesh-coordinates $X^{(i)}\sim p^{(i)}$ and
absolute velocity $y^{(i)}$.

  2. 2.

Based on $(X^{(i)},y^{(i)})$ we run SABER for $10$ repetitions, randomly
drawing $6000$ training and $2000$ validation samples, giving an averaged
bandwidth function $\sigma_{\text{\tiny{SABER}}}^{(i)}$.

  3. 3.

We propose the new sampling density

$p^{(i+1)}(x)\propto\frac{p^{(i)}(x)^{0.4}}{\sigma_{\text{\tiny{SABER}}}^{(i)}(x)^{1.2}}.$

that is discretized to take values of powers of 4 (each power represents a
refinement) and match the background resolution to not exceed the initial
$8608$ cells.

The hyperparameters of SABER were tuned in the first iteration, giving
$\sigma_{g}=6\cdot 10^{-2}$, $\lambda_{g}=2\cdot 10^{-4}$ and
$\sigma_{j}=10^{\frac{j-4}{8}},1\leq j\leq 9$, and reused subsequently. The
LFC estimate (17) is given by

$\mathfrak{C}(x)=\sigma_{\text{\tiny{SABER}}}(x)^{-1}p^{(i)}(x)^{-\frac{1}{6}}.$

The proposed sampling density according to (6) is then:

$p_{\text{\tiny{Opt}}}(x)\propto\left[\mathfrak{C}(x)^{2}p^{(i)}(x)\right]^{\frac{6}{10}}=\frac{p^{(i)}(x)^{0.4}}{\sigma_{\text{\tiny{SABER}}}^{(i)}(x)^{1.2}}.$

Fig. 15 shows the evolution of the cell distribution and the corresponding
velocity fields.

Figure 15: Mesh evolution for the square cavity with differentially heated side walls. TABLE IV: Vortex centers. Compare Fig. 14 for the reference. Mesh | Resolution | # Cells | # Vortices | $e_{L2}$
---|---|---|---|---
Reference | $1024\times 1024$ | $1,048,576$ | $8$ | —
base line | $64\times 64$ | $8608$ | $6$ | 0.0766
base line | $128\times 128$ | $34,432$ | $7$ | 0.0200
1st iteration | $65\times 65$ | $8512$ | $6$ | 0.0724
2nd iteration | $53\times 53$ | $8602$ | $6$ | 0.0521
3rd iteration | $98\times 98$ | $34,360$ | $8$ | 0.0150

The background resolution in the first iteration is increased to $65\times 65$
and many cells from the fine level on top, bottom and left wall are shifted to
the intermediate level. In the second iteration the resolution in the
background was reduced to $53\times 53$ and both intermediate and fine level
get more cells in the region of the vortices. With each iteration, the
quantitative performance increases as can be seen by a decrease in $e_{L2}$.
Note that the present procedure only reorganizes a fix budget of cells, such
that no perfect solution is expected. While the remaining vortices don’t
emerge, the structure close to them becomes more distinct. When calculating
the third iteration the redistribution has almost converged, only increasing
resolution in the area of the vortices a bit further.

When doubling the background resolution (see Fig. 16 in Appendix F), all $8$
vortices emerge, based on the mesh from the third iteration, whereas vortex 6
remains absent, using the base line. Here, we also observe a qualitative
performance gain. Hence, the present approach for cell redistribution improves
the overall quality of the simulation.

We would like to emphasize that this experiment is a proof-of-concept, since
it was not clear from the start that the machine learning concept of active
learning can be applied successfully for mesh adaptation in numeric
simulation. Future work on this will have to apply a more rigid analysis and a
comparison to classical mesh refinement techniques that are based on local
features of the flow field (e.g. gradients) or error estimators, where one
refines regions of leading numerical error.

##  VI Conclusion

Locally optimal bandwidths can not only be used to raise performance of
inference, but also to gain deep insights into the problem. Their estimation
is a challenging meta inference task due to its indirect formulation of
prediction performance of the underlying task, whose optimization is prone to
overfitting. We proposed a novel local bandwidth model with strong
restrictions on the degrees of freedom, which makes this approach robust.

Based on the locally optimal bandwidth estimates, we deduced an interpretable
measure of local function complexity. With manifold applications of
interpretation in quantum chemistry and active learning in fluid dynamics
simulations, we shed light on local function complexity as a problem intrinsic
property that so far did not receive lots of attention.

Future work will include an extensive study on adaptive mesh selection in
numerical simulations and a deeper analysis of datasets from quantum
chemistry.

## Acknowledgments

The authors would like to thank Stephan Lenz and Martin Geier from the
Institute for computational modeling in civil engineering of the Technische
Universität Braunschweig for the fruitful discussion and active help with the
fluid dynamics simulation.

D. Panknin and T. B. Bui were funded by the Federal Ministry of Education and
Research (BMBF) under the project ALICE II, Autonomous Learning in Complex
Environments (01IB15001B). D. Panknin was also funded by the BMBF project
ALICE III, Autonomous Learning in Complex Environments (01IS18049B) and T. B.
Bui was also funded by the BMBF project AIAx, Machine Learning-driven
Engineering – CAx goes AIAx (01IS18048B). S. Nakajima was funded by the BMBF
project BBDC II, Berlin Big Data Center (01IS18025A). K.-R. Müller was
supported by Institute for Information & Communications Technology Promotion
and funded by the Korea government (MSIT) (No. 2017-0-00451, No.
2017-0-01779). K.-R. Müller also acknowledges support by the Berlin
Mathematics Research Center $\text{MATH}^{+}$ and the BMBF project BZML,
Berlin Center for Machine Learning (01IS18037I).

## References

  * [1] G. R. Lanckriet, N. Cristianini, P. Bartlett, L. E. Ghaoui, and M. I. Jordan, “Learning the kernel matrix with semidefinite programming,” _Journal of Machine learning research_ , vol. 5, no. Jan, pp. 27–72, 2004.
  * [2] R. A. Jacobs, M. I. Jordan, S. J. Nowlan, G. E. Hinton _et al._ , “Adaptive mixtures of local experts.” _Neural computation_ , vol. 3, no. 1, pp. 79–87, 1991.
  * [3] M. I. Jordan and R. A. Jacobs, “Hierarchical mixtures of experts and the em algorithm,” _Neural computation_ , vol. 6, no. 2, pp. 181–214, 1994.
  * [4] K. Pawelzik, J. Kohlmorgen, and K.-R. Müller, “Annealed competition of experts for a segmentation and classification of switching dynamics,” _Neural Computation_ , vol. 8, no. 2, pp. 340–356, 1996.
  * [5] K. Kersting, C. Plagemann, P. Pfaff, and W. Burgard, “Most likely heteroscedastic gaussian process regression,” in _International Conference on Machine Learning (ICML)_ , 2007.
  * [6] G. C. Cawley, N. L. Talbot, and O. Chapelle, “Estimating predictive variances with kernel ridge regression,” in _Machine Learning Challenges Workshop_. Springer, 2005, pp. 56–77.
  * [7] B. Wang and X. Wang, “Bandwidth Selection for Weighted Kernel Density Estimation,” _arXiv e-prints_ , p. arXiv:0709.1616, Sep 2007.
  * [8] M. Mackenzie and A. K. Tieu, “Asymmetric kernel regression,” _IEEE transactions on neural networks_ , vol. 15, no. 2, pp. 276–282, 2004.
  * [9] J. Moody and C. J. Darken, “Fast learning in networks of locally-tuned processing units,” _Neural computation_ , vol. 1, no. 2, pp. 281–294, 1989\.
  * [10] N. Benoudjit, C. Archambeau, A. Lendasse, J. A. Lee, M. Verleysen _et al._ , “Width optimization of the gaussian kernels in radial basis function networks.” in _ESANN_ , vol. 2, 2002, pp. 425–432.
  * [11] D. Zheng, J. Wang, and Y. Zhao, “Non-flat function estimation with a multi-scale support vector regression,” _Neurocomputing_ , vol. 70, no. 1, pp. 420–429, 2006.
  * [12] V. Guigue, A. Rakotomamonjy, and S. Canu, “Kernel basis pursuit,” in _European Conference on Machine Learning_. Springer, 2005, pp. 146–157.
  * [13] S. Ferrari, F. Bellocchio, V. Piuri, and N. A. Borghese, “Multi-scale support vector regression,” in _Neural Networks (IJCNN), The 2010 International Joint Conference on_. IEEE, 2010, pp. 1–7.
  * [14] F. Bellocchio, S. Ferrari, V. Piuri, and N. A. Borghese, “Hierarchical approach for multiscale support vector regression,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 23, no. 9, 2012.
  * [15] R. Caruana, Y. Lou, J. Gehrke, P. Koch, M. Sturm, and N. Elhadad, “Intelligible models for healthcare: Predicting pneumonia risk and hospital 30-day readmission,” in _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_. ACM, 2015, pp. 1721–1730.
  * [16] M. T. Ribeiro, S. Singh, and C. Guestrin, “Why should i trust you?: Explaining the predictions of any classifier,” in _Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining_. ACM, 2016, pp. 1135–1144.
  * [17] S. Lapuschkin, S. Wäldchen, A. Binder, G. Montavon, W. Samek, and K.-R. Müller, “Unmasking clever hans predictors and assessing what machines really learn,” _Nature Communications_ , vol. 10, p. 1096, 2019.
  * [18] “Paper in preparation.”
  * [19] K.-R. Müller, A. Smola, G. Rätsch, B. Schölkopf, J. Kohlmorgen, and V. Vapnik, “Using support vector machines for time series prediction,” _Advances in kernel methods—support vector learning_ , pp. 243–254, 1999\.
  * [20] C. Walder, K. I. Kim, and B. Schölkopf, “Sparse multiscale gaussian process regression,” in _Proceedings of the 25th International Conference on Machine Learning_. ACM, 2008, pp. 1112–1119.
  * [21] J. Yuan, L. Bo, K. Wang, and T. Yu, “Adaptive spherical gaussian kernel in sparse bayesian learning framework for nonlinear regression,” _Expert Systems with Applications_ , vol. 36, no. 2, pp. 3982–3989, 2009.
  * [22] M. Peifer, C. L. FO, S. Paternain, and A. Ribeiro, “Locally adaptive kernel estimation using sparse functional programming,” in _52nd Asilomar Conf. on Signals, Systems, and Computers_. IEEE, 2018, pp. 2022–2026.
  * [23] A. Ye, R. J. Hyndman, and Z. Li, “Local linear multivariate regression with variable bandwidth in the presence of heteroscedasticity,” 2006.
  * [24] M. G. Morgan, M. Henrion, and M. Small, _Uncertainty: a guide to dealing with uncertainty in quantitative risk and policy analysis_. Cambridge university press, 1990.
  * [25] M. Gönen and E. Alpaydin, “Localized multiple kernel regression,” in _Pattern Recognition (ICPR), 2010 20th International Conference on_. IEEE, 2010, pp. 1425–1428.
  * [26] J. Fan and I. Gijbels, “Variable bandwidth and local linear regression smoothers,” _The Annals of Statistics_ , pp. 2008–2036, 1992.
  * [27] C. K. Williams and C. E. Rasmussen, “Gaussian processes for regression,” in _Advances in neural information processing systems_ , 1996, pp. 514–520.
  * [28] S. Sundararajan and S. S. Keerthi, “Predictive approaches for choosing hyperparameters in gaussian processes,” in _Advances in neural information processing systems_ , 2000, pp. 631–637.
  * [29] L. D. Brown, M. Levine _et al._ , “Variance estimation in nonparametric regression via the difference sequence method,” _The Annals of Statistics_ , vol. 35, no. 5, pp. 2219–2232, 2007.
  * [30] B. W. Silverman, _Density estimation for statistics and data analysis_. Routledge, 2018.
  * [31] M. Gönen and E. Alpaydın, “Multiple kernel learning algorithms,” _Journal of machine learning research_ , vol. 12, no. Jul, pp. 2211–2268, 2011.
  * [32] S. E. Yuksel, J. N. Wilson, and P. D. Gader, “Twenty years of mixture of experts,” _IEEE transactions on neural networks and learning systems_ , vol. 23, no. 8, pp. 1177–1193, 2012.
  * [33] V. Tresp, “Mixtures of gaussian processes,” in _Advances in neural information processing systems_ , 2001, pp. 654–660.
  * [34] E. Meeds and S. Osindero, “An alternative infinite mixture of gaussian process experts,” in _Advances in Neural Information Processing Systems_ , 2006, pp. 883–890.
  * [35] C. Yuan and C. Neubauer, “Variational mixture of gaussian process experts,” in _Advances in Neural Information Processing Systems_ , 2009, pp. 1897–1904.
  * [36] Y. Yang and J. Ma, “An efficient em approach to parameter learning of the mixture of gaussian processes,” in _International Symposium on Neural Networks_. Springer, 2011, pp. 165–174.
  * [37] Z. Chen, J. Ma, and Y. Zhou, “A precise hard-cut em algorithm for mixtures of gaussian processes,” in _International Conference on Intelligent Computing_. Springer, 2014, pp. 68–75.
  * [38] S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Schütt, and K.-R. Müller, “Machine learning of accurate energy-conserving molecular force fields,” _Science Advances_ , vol. 3, no. 5, p. e1603015, 2017.
  * [39] S. Chmiela, H. E. Sauceda, K.-R. Müller, and A. Tkatchenko, “Towards exact molecular dynamics simulations with machine-learned force fields,” _Nature Communications_ , vol. 9, no. 1, p. 3887, 2018.
  * [40] W. Pronobis, A. Tkatchenko, and K.-R. Müller, “Many-body descriptors for predicting molecular properties with machine learning: Analysis of pairwise and three-body interactions in molecules,” _Journal of Chemical Theory and Computation_ , vol. 14, no. 6, pp. 2991–3003, 2018.
  * [41] L. Ruddigkeit, R. van Deursen, L. C. Blum, and J.-L. Reymond, “Enumeration of 166 billion organic small molecules in the chemical universe database gdb-17,” _Journal of Chemical Information and Modeling_ , vol. 52, no. 11, pp. 2864–2875, 2012.
  * [42] R. Ramakrishnan, P. O. Dral, M. Rupp, and O. A. Von Lilienfeld, “Quantum chemistry structures and properties of 134 kilo molecules,” _Scientific data_ , vol. 1, p. 140022, 2014.
  * [43] K. T. Schütt, H. E. Sauceda, P.-J. Kindermans, A. Tkatchenko, and K.-R. Müller, “Schnet–a deep learning architecture for molecules and materials,” _The Journal of Chemical Physics_ , vol. 148, no. 24, p. 241722, 2018.
  * [44] G. de Vahl Davis, “Natural convection of air in a square cavity: a bench mark numerical solution,” _International Journal for numerical methods in fluids_ , vol. 3, no. 3, pp. 249–264, 1983.
  * [45] J. Vierendeels, B. Merci, and E. Dick, “Benchmark solutions for the natural convective heat transfer problem in a square cavity with large horizontal temperature differences,” _International Journal of Numerical Methods for Heat & Fluid Flow_, vol. 13, no. 8, pp. 1057–1078, 2003.
  * [46] S. Lenz, M. Krafczyk, M. Geier, S. Chen, and Z. Guo, “Validation of a two-dimensional gas-kinetic scheme for compressible natural convection on structured and unstructured meshes,” _International Journal of Thermal Sciences_ , vol. 136, pp. 299–315, 2019.
  * [47] K. Xu, “A gas-kinetic bgk scheme for the navier–stokes equations and its connection with artificial dissipation and godunov method,” _Journal of Computational Physics_ , vol. 171, no. 1, pp. 289–335, 2001.
  * [48] K. Xu, M. Mao, and L. Tang, “A multidimensional gas-kinetic bgk scheme for hypersonic viscous flow,” _Journal of Computational Physics_ , vol. 203, no. 2, pp. 405–421, 2005.
  * [49] S. Lenz, M. Geier, and M. Krafczyk, “An explicit gas kinetic scheme algorithm on non-uniform cartesian meshes for gpgpu architectures,” _Computers & Fluids_, vol. 186, pp. 58–73, 2019.

##  Appendix A Notation

Let $A,B\in\operatorname*{\mathds{R}}^{m_{1}\times n_{1}}$,
$C\in\operatorname*{\mathds{R}}^{m_{2}\times n_{2}}$ and
$v\in\operatorname*{\mathds{R}}^{n}$.

  * •

$A\otimes C=\begin{bmatrix}C_{11}A&\cdots&C_{1n_{2}}A\\\ \vdots&&\vdots\\\
C_{m_{2}1}A&\cdots&C_{m_{2}n_{2}}A\end{bmatrix}$ is the Kronecker product

  * •

$A\odot B=[A_{ij}B_{ij}]_{ij}$ is the Hadamard (pointwise) product

  * •

In analogy $A\oslash B=[A_{ij}/B_{ij}]_{ij}$ is the pointwise division

  * •

$\operatorname*{\bf{diag}}(v)_{ij}=\begin{cases}v_{i},i=j\\\
0,\text{else}\end{cases}$ is a matrix with $v$ on its diagonal

  * •

$\mathbbmss{1}_{n}=\begin{bmatrix}1&\cdots&1\end{bmatrix}^{\top}\in\operatorname*{\mathds{R}}^{n}$
the length-n vector of ones

  * •

$\mathcal{I}_{n}=\operatorname*{\bf{diag}}(\mathbbmss{1}_{n})$ the identity
matrix in $\operatorname*{\mathds{R}}^{n\times n}$

  * •

vec is the vectorization of a matrix via columnwise concatenation

##  Appendix B Mixture of Experts

In the ME paradigm we assume an observation to be generated by one of several
processes with distinct behavior. The probability that a specific process
drives the observation varies over the input space. If the supports of the
processes, that is, where the respective generating probability is not
negligible, overlap at most moderately, then we can cluster the input space
with respect to the generating processes. In that sense, assigning the most
probable generating process to an input can be formulated as a hidden
classification problem.

Let $(\bm{X},\bm{Y})$ be the training set of size $n$ and let $\bm{Z}$ be
hidden multinomial variables

$\bm{Z}_{ij}=\begin{cases}1,\text{process }j\text{ generated }y_{i}\\\
0,\text{else}\end{cases}.$

indicating the generating one of $L$ potential processes. Finally, let
$\Theta=\\{\theta_{g},\theta_{e_{1}},\ldots,\theta_{e_{L}}\\}$ be the set of
all model parameters for describing both, the expert models
$\operatorname*{\mathds{P}}(\bm{Y}|\bm{X},\theta_{e_{j}})$ when assuming
process $j$ to drive the observations, and the probability
$\operatorname*{\mathds{P}}(\bm{Z}|\bm{X},\theta_{g})$ that specific processes
drive the observations. Then we can write the observation likelihood as

$\operatorname*{\mathds{P}}(\bm{Y}|\bm{X},\Theta)=\sideset{}{{}_{\bm{Z}}}{\sum}\operatorname*{\mathds{P}}(\bm{Y}|\bm{X},\bm{Z},\theta_{e_{1}},\ldots,\theta_{e_{L}})\operatorname*{\mathds{P}}(\bm{Z}|\bm{X},\theta_{g}).$

With respect to maximizing the observation likelihood, the optimal $\Theta$
can be estimated iteratively via an EM-procedure:

In the E-step, we estimate the posterior distribution of the hidden variable
$\bm{Z}^{(t)}\sim\operatorname*{\mathds{P}}(\bm{Z}|\bm{X},\bm{Y},\Theta^{(t)})$
given the current set of parameters $\Theta^{(t)}$. In the M-step, we update
the parameters such that the augmented log-likelihood is maximized

$\Theta^{(t+1)}\leftarrow\operatorname*{\bf{argmax}}_{\Theta}\sideset{}{{}_{\bm{Z}^{(t)}}}{\operatorname*{\mathds{E}}}\log\operatorname*{\mathds{P}}(\bm{Y},\bm{Z}^{(t)}|\bm{X},\Theta).$

If the observations $\bm{Y}$ are mutually independent, conditionally on
$\bm{X}$, then the EM-procedure can be simplified: In the E-step, we can
factorize the hidden variable posterior

$\displaystyle\bm{P}^{(t)}_{ij}$
$\displaystyle=\operatorname*{\mathds{P}}(\bm{Z}_{ij}=1|\bm{X},\bm{Y},\Theta^{(t)})$
$\displaystyle=\frac{\operatorname*{\mathds{P}}(\bm{Z}_{ij}=1|x_{i},\theta^{(t)}_{g})\operatorname*{\mathds{P}}(y_{i}|x_{i},\theta^{(t)}_{e_{j}})}{\sideset{}{{}_{k=1}^{L}}{\sum}\operatorname*{\mathds{P}}(\bm{Z}_{i,k}=1|x_{i},\theta^{(t)}_{g})\operatorname*{\mathds{P}}(y_{i}|x_{i},\theta^{(t)}_{e_{k}})},$

describing the belief that process $j$ drives the $i^{th}$ observation, given
the current set of hyperparameters $\Theta^{(t)}$.

The M-step can be reformulated as

$\displaystyle\Theta^{(t+1)}\leftarrow\operatorname*{\bf{argmax}}_{\Theta}$
$\displaystyle\underbrace{\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{L}\bm{P}^{(t)}_{ij}\log\operatorname*{\mathds{P}}(\bm{Z}_{ij}=1|x_{i},\theta_{g})}_{R_{g}(\theta_{g})}$
$\displaystyle+\sum\limits_{j=1}^{L}\underbrace{\sum\limits_{i=1}^{n}\bm{P}^{(t)}_{ij}\log\operatorname*{\mathds{P}}(y_{i}|x_{i},\theta_{e_{j}})}_{R_{j}(\theta_{e_{j}})}.$

This form reveals, that each expert and the gate can be optimized
independently in the M-step with respect to their own set of parameters. Note
that predictions of the model take the form

$\widehat{y}(x)=\sideset{}{{}_{z}}{\operatorname*{\mathds{E}}}\operatorname*{\mathds{P}}(y|x,\Theta)=\sideset{}{{}_{j=1}^{L}}{\sum}\operatorname*{\mathds{E}}(y|x,\theta_{e_{j}})\operatorname*{\mathds{P}}(z_{j}=1|x,\theta_{g}).$

##  Appendix C Test predictions of GPR and MKLR

The GPR predictions (7) are summarized in Alg. 4.

Algorithm 4
$\widehat{\bm{Y}}^{*}\leftarrow\text{predict}_{\text{GPR}}(\bm{X}^{*}|\bm{X},\bm{Y},v,\theta_{e})$

1:[-0.4ex]0.490.4pt4pt 3pt

2:Input

3:0.490.4pt4pt 3pt

4:Training data
$\bm{X}=[x_{1},\ldots,x_{n}]^{\top}\in\operatorname*{\mathds{R}}^{n\times
d},\bm{Y}=[y_{1},\ldots,y_{n}]^{\top}\in\operatorname*{\mathds{R}}^{n}$

5:Local noise variance $v$ and GPR hyperparameters
$\theta_{e}=(m,\lambda_{f},\sigma_{\varepsilon}^{2},\Sigma)$

6:Test input data
$\bm{X}^{*}=[x_{1}^{*},\ldots,x_{T}^{*}]^{\top}\in\operatorname*{\mathds{R}}^{T\times
d}$

7:[-0.4ex]0.490.4pt4pt 3pt

8:Output

9:0.490.4pt4pt 3pt

10:Predictions $\widehat{\bm{Y}}^{*}\in\operatorname*{\mathds{R}}^{T}$

11:[-0.4ex]0.490.4pt4pt 3pt

12:Procedure

13:0.490.4pt4pt 3pt

14:Set
$\Sigma_{\varepsilon}^{2}=\operatorname*{\bf{diag}}\left(v(\bm{X})\right)$

15:Estimate
$\bm{C}=\sigma_{\varepsilon}^{2}(\exp(\lambda_{f})\bm{K}^{\Sigma}(\bm{X})+\Sigma_{\varepsilon}^{2})$

16:Solve $\bm{C}\alpha=\bm{Y}-m\mathbbmss{1}_{n}$ for $\alpha$

17:Predict $\widehat{y}^{*}\leftarrow
m\mathbbmss{1}_{T}+\sigma_{\varepsilon}^{2}\exp(\lambda_{f})\bm{K}^{\Sigma}(\bm{X}^{*},\bm{X})\alpha$

We would like to summarize the MKLR predictions according to (10) in Alg. 5.

Algorithm 5
$\bm{Q}^{*}\leftarrow\text{predict}_{\text{MKLR}}(\bm{X}^{*}|\bm{X},\bm{\omega},\theta_{g})$

1:[-0.4ex]0.490.4pt4pt 3pt

2:Input

3:0.490.4pt4pt 3pt

4:Training inputs
$\bm{X}=[x_{1},\ldots,x_{n}]^{\top}\in\operatorname*{\mathds{R}}^{n\times d}$

5:MKLR hyperparameters $\theta_{g}=(\Sigma_{g},\lambda_{g})$ and parameters
$\bm{\omega}=[\bm{\alpha}^{\top},\bm{b}^{\top}]^{\top}$

6:Test input data
$\bm{X}^{*}=[x_{1}^{*},\ldots,x_{T}^{*}]^{\top}\in\operatorname*{\mathds{R}}^{T\times
d}$

7:[-0.4ex]0.490.4pt4pt 3pt

8:Output

9:0.490.4pt4pt 3pt

10:Test class assignments $\bm{Q}^{*}\in[0,1]^{T\times L}$

11:[-0.4ex]0.490.4pt4pt 3pt

12:Procedure

13:0.490.4pt4pt 3pt

14:Set
$\bm{q}^{*}=[\exp(\bm{K}^{\Sigma_{g}}(\bm{X}^{*},\bm{X})\bm{\alpha}+\mathbbmss{1}_{T}\bm{b}),\mathbbmss{1}_{T}]$

15:Normalize
$\bm{Q}^{*}\leftarrow\operatorname*{\bf{diag}}(\bm{q}^{*}\mathbbmss{1}_{L})^{-1}\bm{q}^{*}$

16:return $\bm{Q}^{*}$

##  Appendix D Gradient descent of the GPR objective

Recall the objective (9) of the GPR model,

$\ell_{\text{e}}(\bm{X},\bm{Y},\theta_{e},P)=\sum\limits_{i=1}^{n}\frac{P_{i}}{2}\left[\frac{looR_{i}^{2}}{\sigma^{*}_{-i}(x_{i})}+\log(\sigma^{*}_{-i}(x_{i}))\right],$

where $looR_{i}$ and $\sigma^{*}_{-i}(x_{i})$ were the leave-one-out residuals
and predictive variances based on the restricted training set
$(\bm{X}_{-i},\bm{Y}_{-i})$. Define
$\overline{\bm{D}}=\sigma_{\varepsilon}^{2}\bm{C}^{-1}$ and
$W=\operatorname*{\bf{diag}}(\overline{\bm{D}})$. We can rewrite
$[\sigma^{*}_{-i}(x_{i})]_{i}=\sigma_{\varepsilon}^{2}\oslash W$ and
$looR=looY-m\cdot loo\mathbbmss{1}$ with
$looY=\left[\overline{\bm{D}}\bm{Y}\right]\oslash W$ and
$loo\mathbbmss{1}=\left[\overline{\bm{D}}\mathbbmss{1}\right]\oslash W$, such
that:

$\displaystyle\ell_{\text{e}}(\bm{X},\bm{Y},\theta_{e},P)=\sum\limits_{i=1}^{n}\frac{P_{i}}{2}\left[\frac{looR_{i}^{2}}{\sigma_{\varepsilon}^{2}}W_{i}-\log(W_{i})+\log(\sigma_{\varepsilon}^{2})\right]$

Then

$\frac{d\ell_{\text{e}}(\bm{X},\bm{Y},\theta_{e},P)}{d\sigma_{\varepsilon}^{2}}=0\Leftrightarrow\sigma_{\varepsilon}^{2}\leftarrow\frac{\mathbbmss{1}_{n}^{\top}[looR\odot
P\odot W\odot looR]}{\mathbbmss{1}_{n}^{\top}P}$

Therefore the optimal noise variance estimate relates to the weighted mean
squared residual, which is intuitive. Furthermore

$\frac{d\ell_{\text{e}}(\bm{X},\bm{Y},\theta_{e},P)}{dm}=0\Leftrightarrow
m\leftarrow\frac{\mathbbmss{1}_{n}^{\top}[loo\mathbbmss{1}\odot P\odot W\odot
looY]}{\mathbbmss{1}_{n}^{\top}[loo\mathbbmss{1}\odot P\odot W\odot
loo\mathbbmss{1}]}.$

Therefore the optimal mean estimate relates to the weighted mean output. The
remaining parameters are contained in $\overline{\bm{D}}$. Define

$A_{p}=-\frac{d\overline{\bm{D}}}{dp}=\overline{\bm{D}}\frac{d\bm{D}}{dp}\overline{\bm{D}}.$

Then

$\displaystyle\frac{d\ell_{\text{e}}(\bm{X},\bm{Y},\theta_{e},P)}{dp}=-\frac{1}{\sigma_{\varepsilon}^{2}}(looR\odot
P)^{\top}A_{p}(\bm{Y}-m\mathbbmss{1})$
$\displaystyle\quad\quad\,+\frac{1}{2\sigma_{\varepsilon}^{2}}\mathbbmss{1}^{\top}[P\odot\operatorname*{\bf{diag}}(A_{p})\odot[\sigma_{\varepsilon}^{2}\oslash
W+looR\odot looR]].$

It is $A_{\lambda_{f}}=\overline{\bm{D}}(I-\overline{\bm{D}})$. The gradient
with respect to $\Sigma$ is the only one depending on the concrete kernel
choice. We give the explicit gradient for the Gaussian kernel: In the
isotropic case $\Sigma=\exp(\sigma^{x})\mathcal{I}_{d}$ it is
$A_{\sigma^{x}}=\exp(-2\sigma^{x})\left[d(\bm{X})\odot(\bm{D}-I)\right]$,
where $d(\bm{X})_{ij}=(x_{i}-x_{j})^{2}$ is the squared distance matrix. In
the diagonal case
$\Sigma=\operatorname*{\bf{diag}}(\exp(\sigma^{x}_{1}),\ldots,\exp(\sigma^{x}_{d}))$
it is
$A_{\sigma^{x}_{k}}=\exp(-2\sigma^{x}_{k})\left[d^{k}(\bm{X})\odot(\bm{D}-I)\right]$,
where $d^{k}(\bm{X})_{ij}=(x^{k}_{i}-x^{k}_{j})^{2}$ is the coordinate-wise
squared distance matrix.

With this we perform gradient descent for $\lambda_{f}$ and $\Sigma$ until
convergence, applying the closed-form updates for $m$ and
$\sigma_{\varepsilon}^{2}$ in each iteration.

##  Appendix E Minimizing the MKLR objective

We summarized the efficient training of MKLR in Alorithm 6 with explicit
calculation formulas described in this appendix.

Algorithm 6
$\bm{\omega}\leftarrow\text{train}_{\text{MKLR}}(\bm{X},\bm{P},\theta_{g})$

1:[-0.4ex]0.490.4pt4pt 3pt

2:Input

3:0.490.4pt4pt 3pt

4:Training inputs
$\bm{X}=[x_{1},\ldots,x_{n}]^{\top}\in\operatorname*{\mathds{R}}^{n\times d}$

5:Training bandwidths assignment probabilities $\bm{P}\in[0,1]^{n\times L}$

6:MKLR hyperparameters $\theta_{g}=(\Sigma_{g},\lambda_{g})$

7:[-0.4ex]0.490.4pt4pt 3pt

8:Output

9:0.490.4pt4pt 3pt

10:Optimal MKLR parameters $\bm{\omega}$

11:[-0.4ex]0.490.4pt4pt 3pt

12:Procedure

13:0.490.4pt4pt 3pt

14:Initialize $\bm{\alpha}=0$,
$\bm{b}=\log\left(\frac{1}{n}\mathbbmss{1}_{n}^{\top}\bm{P}_{-L}\right)$

15:Estimate the kernel $\bm{K}^{\Sigma_{g}}(\bm{X})$ and
$B=\frac{1}{2}\left(\mathcal{I}_{L-1}-\frac{1}{L}\mathbbmss{1}_{L-1}\mathbbmss{1}_{L-1}^{\top}\right)$

16:Estimate the eigendecompositions
$\bm{K}^{\Sigma_{g}}(\bm{X})=U_{g}\operatorname*{\bf{diag}}(\Lambda_{g})U_{g}^{\top}$
and $B=V\operatorname*{\bf{diag}}(D)V^{\top}$

17:Virtually define $\triangleright$ Don’t calculate explicitly!
$\widetilde{\bm{M}}^{-1}=(V\otimes
U_{g})[\operatorname*{\bf{diag}}(D\otimes\Lambda_{g})+\lambda_{g}\mathcal{I}_{n(L-1)}]^{-1}(V^{\top}\otimes
U_{g}^{\top})$

18:Estimate
$\bm{N}^{-1}=(\mathcal{I}_{L-1}\otimes\mathbbmss{1}_{n}^{\top})\widetilde{\bm{M}}^{-1}(B\otimes\mathbbmss{1}_{n})$
$\triangleright$ see (19)

19:Solve for $\bm{N}$

20:while $\bm{\alpha}$ not converged do

21: Set
$\bm{q}=\exp(\bm{K}^{\Sigma_{g}}(\bm{X})\bm{\alpha}+\mathbbmss{1}_{n}\bm{b})$

22: Normalize
$\bm{Q}=\operatorname*{\bf{diag}}(\mathbbmss{1}_{n}+\bm{q}\mathbbmss{1}_{L-1})^{-1}\bm{q}$

23: Set $\bm{R}=\bm{Q}-\bm{P}_{-L}+\lambda_{g}\bm{\alpha}$

24:
$\textbf{{vec}}(\Delta\bm{b})=\bm{N}\left((\mathcal{I}_{L-1}\otimes\mathbbmss{1}^{\top})\widetilde{\bm{M}}^{-1}\textbf{{vec}}(\bm{R})-\bm{\alpha}^{\top}\mathbbmss{1}\right)$
$\triangleright$ see (20)

25:
$\textbf{{vec}}(\Delta\bm{\alpha})=\widetilde{\bm{M}}^{-1}\textbf{{vec}}(\bm{R}-(\Delta\bm{b}B)\otimes\mathbbmss{1}_{n})$
$\triangleright$ see (21)

26: Update $\bm{\alpha}=\bm{\alpha}-\Delta\bm{\alpha}$,
$\bm{b}=\bm{b}-\Delta\bm{b}$

27:return $\bm{\omega}=[\bm{\alpha}^{\top},\bm{b}^{\top}]^{\top}$

Recall the objective (11) of the MKLR model,

$\displaystyle\ell_{\text{g}}(\bm{\omega},$
$\displaystyle\bm{X},\bm{P},\theta_{g})=\lambda_{g}\sum\limits_{j=1}^{L-1}\alpha_{j}^{\top}\bm{K}^{\Sigma_{g}}(\bm{X})\alpha_{j}$
$\displaystyle+\sum\limits_{i=1}^{n}\left[\log\left(1+\sum\limits_{j=1}^{L-1}\bm{q}_{ij}\right)-\sum\limits_{j=1}^{L-1}\bm{P}_{ij}\log(\bm{q}_{ij})\right].$

Let
$\bm{Q}=[Q_{1},\ldots,Q_{L-1}]=\operatorname*{\bf{diag}}(\mathbbmss{1}_{n}+\bm{q}\mathbbmss{1}_{L-1})^{-1}\bm{q}$
be the normalized form of $\bm{q}$. Then

$\nabla_{\text{g}}(\bm{\omega})=\frac{d\ell_{\text{g}}(\bm{\omega},\bm{X},\bm{P},\theta_{g})}{d\textbf{{vec}}(\bm{\omega})}^{\top}=(\mathcal{I}_{L-1}\otimes\widetilde{\bm{K}})\textbf{{vec}}(\bm{S}),$

where

$\widetilde{\bm{K}}=\begin{bmatrix}\bm{K}&0\\\
\mathbbmss{1}^{\top}&-\lambda_{g}\end{bmatrix}$

and $\bm{S}=\left[S_{1},\ldots,S_{L-1}\right]$ with
$S_{j}=\begin{bmatrix}Q_{j}-P_{j}+\lambda_{g}\alpha_{j}\\\
\mathbbmss{1}^{\top}\alpha_{j}\end{bmatrix}$.

Note that there is no closed form solution for
$\nabla_{\text{g}}(\bm{\omega})=0$. Fortunately, the objective is convex such
that we can approximate the solution via _iteratively reweighted least
squares_ , performing _Newton-Raphson_ updates on the gradient to find its
unique root via

$\textbf{{vec}}(\bm{\omega})\leftarrow\textbf{{vec}}(\bm{\omega})-\bm{H}_{\text{g}}(\bm{\omega})^{-1}\nabla_{\text{g}}(\bm{\omega}),$

where the Hessian is given by

$\displaystyle\bm{H}_{\text{g}}(\bm{\omega})=\frac{d^{2}\ell_{\text{g}}(\bm{\omega},\bm{X},\bm{P},\theta_{g})}{d^{2}\textbf{{vec}}(\bm{\omega})}=\mathcal{I}_{L-1}\otimes\begin{bmatrix}\lambda_{g}\bm{K}&0\\\
0&0\end{bmatrix}$ $\displaystyle+$
$\displaystyle\mathcal{I}_{L-1}\otimes\begin{bmatrix}\bm{K}\\\
\mathbbmss{1}^{\top}\end{bmatrix}\left[\operatorname*{\bf{diag}}(\textbf{{vec}}(\bm{Q}))-\mathcal{Q}\mathcal{Q}^{\top}\right]\mathcal{I}_{L-1}\otimes[\bm{K},\mathbbmss{1}],$

with the block-vector

$\mathcal{Q}=\begin{bmatrix}\operatorname*{\bf{diag}}(Q_{1})\\\ \vdots\\\
\operatorname*{\bf{diag}}(Q_{L-1})\end{bmatrix}.$

Solving $\bm{H}_{\text{g}}(\bm{\omega})^{-1}\nabla_{\text{g}}(\bm{\omega})$ in
each iteration is however a computational bottleneck. In the multiclass case,
solving one such system of $n(L-1)$ linear equations may already be
infeasible. To circumvent both problems, we instead use the bound optimization
technique according to [Böhning, 1992], where the Hessian is replaced by an
upper bound positive definite matrix
$\bm{B}\geq\bm{H}_{\text{g}}(\bm{\omega})$ in the sense of _Loewner ordering_
, that is, $\bm{B}-\bm{H}_{\text{g}}(\bm{\omega})$ is positive semi-definite.
Such an upper bound is given by

$\bm{B}=B\otimes\begin{bmatrix}\bm{K}\\\
\mathbbmss{1}^{\top}\end{bmatrix}[\bm{K},\mathbbmss{1}]+\mathcal{I}_{L-1}\otimes\begin{bmatrix}\lambda_{g}\bm{K}&0\\\
0&0\end{bmatrix},$

where

$B=\frac{1}{2}\left(\mathcal{I}_{L-1}-\frac{1}{L}\mathbbmss{1}\mathbbmss{1}^{\top}\right).$

Then the new update rule

$\bm{\omega}\leftarrow\bm{\omega}-\Delta\bm{\omega}$

is still guaranteed to converge, where
$\textbf{{vec}}(\Delta\bm{\omega})=\bm{B}^{-1}\nabla_{\text{g}}(\bm{\omega})$.
We can rewrite

$\bm{B}=(\mathcal{I}_{L-1}\otimes\widetilde{\bm{K}})\bm{M},$

with

$\bm{M}=B\otimes\begin{bmatrix}\bm{K}&\mathbbmss{1}\\\
0&0\end{bmatrix}+\mathcal{I}_{L-1}\otimes\begin{bmatrix}\lambda_{g}\mathcal{I}_{n}&0\\\
\mathbbmss{1}^{\top}&0\end{bmatrix}.$

Therefore we need to solve
$\bm{M}\textbf{{vec}}(\Delta\bm{\omega})=\textbf{{vec}}(\bm{S})$. For this we
split the system of linear equations into two groups, the rows that correspond
to the update of $\bm{\alpha}$ and $\bm{b}$, respectively: Denote

$\Delta\bm{\omega}=\begin{bmatrix}\Delta\bm{\alpha}\\\
\Delta\bm{b}\end{bmatrix}$

and let $\widetilde{\bm{M}}=B\otimes\bm{K}+\lambda_{g}\mathcal{I}_{n(L-1)}$
and $\bm{R}=\bm{Q}-\bm{P}_{-L}+\lambda_{g}\bm{\alpha}$. Then

$\displaystyle\textbf{{vec}}(\bm{R})=[\widetilde{\bm{M}},B\otimes\mathbbmss{1}]\begin{bmatrix}\textbf{{vec}}(\Delta\bm{\alpha})\\\
\textbf{{vec}}(\Delta\bm{b})\end{bmatrix}$
$\displaystyle\quad\quad\quad=\widetilde{\bm{M}}\textbf{{vec}}(\Delta\bm{\alpha})+(B\otimes\mathbbmss{1})\textbf{{vec}}(\Delta\bm{b})$
$\displaystyle\Leftrightarrow$
$\displaystyle\textbf{{vec}}(\Delta\bm{\alpha})=\widetilde{\bm{M}}^{-1}\textbf{{vec}}(\bm{R}-(\Delta\bm{b}B)\otimes\mathbbmss{1})$

The second set of equations are fulfilled if

$\displaystyle\bm{\alpha}^{\top}\mathbbmss{1}=\mathcal{I}_{L-1}\otimes[\mathbbmss{1}^{\top},0]\textbf{{vec}}(\Delta\bm{\omega})=(\mathcal{I}_{L-1}\otimes\mathbbmss{1}^{\top})\textbf{{vec}}(\Delta\bm{\alpha}).$

Plugging in the first into the second equation we can solve for
$\Delta\bm{b}$:

$\displaystyle\bm{\alpha}^{\top}\mathbbmss{1}=(\mathcal{I}_{L-1}\otimes\mathbbmss{1}^{\top})\widetilde{\bm{M}}^{-1}\textbf{{vec}}(\bm{R}-(\Delta\bm{b}B)\otimes\mathbbmss{1})$
$\displaystyle\Leftrightarrow$
$\displaystyle\textbf{{vec}}(\Delta\bm{b})=\bm{N}\left((\mathcal{I}_{L-1}\otimes\mathbbmss{1}^{\top})\widetilde{\bm{M}}^{-1}\textbf{{vec}}(\bm{R})-\bm{\alpha}^{\top}\mathbbmss{1}\right),$

where we have defined

$\bm{N}^{-1}=(\mathcal{I}_{L-1}\otimes\mathbbmss{1}^{\top})\widetilde{\bm{M}}^{-1}(B\otimes\mathbbmss{1}).$

Note that since $\bm{N}$ is an $(L-1)^{2}$-matrix its inversion is fast. The
advantage over explicit _Newton-Raphson_ updates is as follows:

First of all, since $\widetilde{\bm{M}}$ does not depend on $\bm{\omega}$ we
can calculate its inverse once and reuse it over all iterations. Moreover we
can exploit the structure of $\widetilde{\bm{M}}$, being sum of a Kronecker
product of two symmetric positive semi-definite matrices and a multiple of the
identity. Therefore we can reduce the complexity of its inversion to the
eigendecomposition of the two matrices
$B=V\operatorname*{\bf{diag}}(D)V^{\top}$ and
$\bm{K}^{\Sigma_{g}}=U_{g}\operatorname*{\bf{diag}}(\Lambda_{g})U_{g}^{\top}$:

$\widetilde{\bm{M}}^{-1}=(V\otimes
U_{g})[\operatorname*{\bf{diag}}(D\otimes\Lambda_{g})+\lambda_{g}\mathcal{I}_{n(L-1)}]^{-1}(V^{\top}\otimes
U_{g}^{\top})$

Finally, due to its form we neither need to calculate
$\widetilde{\bm{M}}^{-1}$ explicitly nor we require matrix-vector
multiplication in the $n(L-1)$-dimensional space. This is because for matrices
$A_{1}\in\operatorname*{\mathds{R}}^{n_{1}\times m_{1}}$,
$A_{2}\in\operatorname*{\mathds{R}}^{n_{2}\times m_{2}}$,
$A_{3}\in\operatorname*{\mathds{R}}^{m_{2}\times m_{1}}$ we have

$(A_{1}\otimes
A_{2})\textbf{{vec}}(A_{3})=\textbf{{vec}}(A_{2}A_{3}A_{1}^{\top})$

and for $B_{1}\in\operatorname*{\mathds{R}}^{m_{1}\times
n_{3}},B_{2}\in\operatorname*{\mathds{R}}^{m_{2}\times n_{4}}$ we have

$(A_{1}\otimes A_{2})(B_{1}\otimes B_{2})=(A_{1}B_{1}\otimes A_{2}B_{2}).$

For example we can calculate $\bm{N}^{-1}$ efficiently via

$\displaystyle\bm{N}^{-1}=$
$\displaystyle(V\otimes(\mathbbmss{1}_{n}^{\top}U_{g}))[\operatorname*{\bf{diag}}(D\otimes\Lambda_{g})+\lambda_{g}\mathcal{I}_{n(L-1)}]^{-1}$
$\displaystyle\cdot\,((V^{\top}B)\otimes U_{g}^{\top}\mathbbmss{1}_{n}).$ (19)

Then the updates are given by

$\displaystyle\textbf{{vec}}(\Delta\bm{b})$
$\displaystyle=\bm{N}\Big{(}-\bm{\alpha}^{\top}\mathbbmss{1}_{n}$ (20)
$\displaystyle+V\left[\left(U_{g}^{\top}\bm{R}V\right)\oslash(\Lambda_{g}D^{\top}+\lambda_{g}\mathbbmss{1}_{n}\mathbbmss{1}_{L-1}^{\top})\right]^{\top}U_{g}^{\top}\mathbbmss{1}_{n}\Big{)}$

and

$\displaystyle\Delta\bm{\alpha}=U_{g}\big{[}$
$\displaystyle\left(U_{g}^{\top}(\bm{R}-(\Delta\bm{b}B)\otimes\mathbbmss{1})V\right)$
(21)
$\displaystyle\oslash(\Lambda_{g}D^{\top}+\lambda_{g}\mathbbmss{1}_{n}\mathbbmss{1}_{L-1}^{\top})\big{]}V^{\top}.$

##  Appendix F Supplemental

TABLE V: Vortex centers. Compare Fig. 14 for the reference. | single resolution | double resolution
---|---|---
| Reference | base line | 1st iteration | 2nd iteration | 3rd iteration | base line
Resolution | $1024\times 1024$ | $64\times 64$ | $65\times 65$ | $53\times 53$ | $98\times 98$ | $128\times 128$
# Cells | $1,048,576$ | $8608$ | $8512$ | $8602$ | $34360$ | $34432$
Vortex |  |  |  |  |  |
1 | x | 0.0989 | 0.1064 | 0.1001 | 0.1008 | 0.0992 | 0.1005
y | 0.5016 | 0.4708 | 0.4533 | 0.4868 | 0.4975 | 0.4949
e | — | 0.0317 | 0.0483 | 0.0149 | 0.0041 | 0.0069
2 | x | 0.1639 | 0.1391 | 0.1491 | 0.1553 | 0.1643 | 0.1611
y | 0.8211 | 0.7830 | 0.7987 | 0.8029 | 0.8181 | 0.8175
e | — | 0.0455 | 0.0268 | 0.0201 | 0.0030 | 0.0045
3 | x | 0.2798 | 0.2810 | 0.2751 | 0.2855 | 0.2787 | 0.2790
y | 0.4556 | 0.4336 | 0.4386 | 0.4472 | 0.4524 | 0.4512
e | — | 0.0220 | 0.0176 | 0.0102 | 0.0034 | 0.0044
4 | x | 0.8009 | 0.8298 | 0.8324 | 0.8353 | 0.8109 | 0.8123
y | 0.0419 | 0.0055 | 0.0123 | 0.0142 | 0.0330 | 0.0295
e | — | 0.0465 | 0.0432 | 0.0442 | 0.0134 | 0.0169
5 | x | 0.8040 | — | — | — | 0.8052 | 0.8058
y | 0.4301 | — | — | — | 0.4288 | 0.4273
e | — | — | — | — | 0.0018 | 0.0033
6 | x | 0.8807 | — | — | — | 0.8792 | —
y | 0.1765 | — | — | — | 0.1748 | —
e | — | — | — | — | 0.0023 | —
7 | x | 0.9207 | 0.9301 | 0.9228 | 0.9227 | 0.9203 | 0.9190
y | 0.0793 | 0.0819 | 0.0804 | 0.0835 | 0.0804 | 0.0829
e | — | 0.0098 | 0.0024 | 0.0047 | 0.0012 | 0.0040
8 | x | 0.9256 | 0.9224 | 0.9244 | 0.9248 | 0.9253 | 0.9242
y | 0.3909 | 0.3837 | 0.3895 | 0.3882 | 0.3904 | 0.3905
e | — | 0.0079 | 0.0018 | 0.0028 | 0.0005 | 0.0014
$e_{L2}$ | — | 0.0766 | 0.0724 | 0.0521 | 0.0150 | 0.0200

(a) Base line mesh.

(b) Adapted mesh in third iteration.

Figure 16: Initial and final mesh of the fluid dynamics experiment in double
resolution.