# Variational Hierarchical Mixtures for
Probabilistic Learning of Inverse Dynamics

Hany Abdulsamad, Peter Nickl, Pascal Klink, and Jan Peters H. Abdulsamad is
with the Department of Electrical Engineering and Automation at Aalto
University, Finland. P. Nickl is with the RIKEN Center for Advanced
Intelligence Project, Japan. P. Klink, and J. Peters are with the Department
of Computer Science at the Technical University of Darmstadt, Germany.

###### Abstract

Well-calibrated probabilistic regression models are a crucial learning
component in robotics applications as datasets grow rapidly and tasks become
more complex. Unfortunately, classical regression models are usually either
probabilistic kernel machines with a flexible structure that does not scale
gracefully with data or deterministic and vastly scalable automata, albeit
with a restrictive parametric form and poor regularization. In this paper, we
consider a probabilistic hierarchical modeling paradigm that combines the
benefits of both worlds to deliver computationally efficient representations
with inherent complexity regularization. The presented approaches are
probabilistic interpretations of local regression techniques that approximate
nonlinear functions through a set of local linear or polynomial units.
Importantly, we rely on principles from Bayesian nonparametrics to formulate
flexible models that adapt their complexity to the data and can potentially
encompass an infinite number of components. We derive two efficient
variational inference techniques to learn these representations and highlight
the advantages of hierarchical infinite local regression models, such as
dealing with non-smooth functions, mitigating catastrophic forgetting, and
enabling parameter sharing and fast predictions. Finally, we validate this
approach on large inverse dynamics datasets and test the learned models in
real-world control scenarios.

###### Index Terms:

Hierarchical Local Regression, Generative Models, Dirichlet Process Mixtures,
Inverse Dynamics Control

##  1 Introduction

Principled data-driven, adaptive, and incremental learning is desirable in
domains in which datasets are dynamic and accumulate slowly over time. For
example, robots must build models of their dynamics and the environment as
they interact with the world. Moreover, these models must be computationally
efficient during both learning and evaluation. In the case of general-purpose
robots, these models must incorporate different modalities of continuous and
discrete stochastic random variables and possibly incorporate heteroscedastic
noise [1, 2]. Predominant and successful regression techniques, such as
Gaussian process regression (GPR) [3], artificial neural networks (ANNs) [4],
and local regression (LR) [5], have a mixed set of properties that are useful
in different scenarios.

Gaussian process regression offers a principled Bayesian treatment that
enables continual and incremental learning. Nonetheless, the vanilla
formulation of GPR [3] suffered from many drawbacks that have been gradually
addressed by recent research. Some of these drawbacks are the functional
smoothness assumption [6, 7, 8], limitations when scaling to large datasets
[9, 10, 11, 12, 13, 14] and difficulties modeling heteroscedasticity [15, 16,
17].

Figure 1: Left, gap data learned with infinite local regression (ILR). The top
plot depicts the mean prediction (red) on the training data (dots) and the
true mean function (dashed). The blue area represents the predictive
uncertainty of two standard deviations. The predictive uncertainty in areas
lacking training data is large, and the mean prediction falls back to the
prior. Right, the CMB dataset learned by ILR. The top figure depicts the mean
prediction (red) with a predictive uncertainty of three standard deviations
(shaded blue). ILR captures the heteroscedastic spread of the data with a
handful of local regression models. The bottom plots show the activation of
the models over the input space.

On the other hand, artificial neural networks have proven themselves as very
powerful easy-to-train universal approximators. They are, however, still
susceptible to over-parameterization [18] and catastrophic forgetting [19].
Moreover, despite major progress on the front of Bayesian neural networks
(BNNs) [20, 21, 22, 23, 24, 25, 26], new evidence suggests that issues
regarding the accuracy of uncertainty quantification still need to be tackled
[27, 28].

Finally, local regression methods have had great success in the domain of
robotics and control [29, 30, 31, 32], because of their flexibility, ability
to model hard nonlinearities, and to incorporate new data online. More
generally, local regression is a family of generative mixture of experts (MoE)
techniques that take a basis-function approach to model the input density and
automatically induce local model responsibilities [33, 34, 35, 36], see Figure
1. In contrast, discriminative MoEs rely on an explicitly input-conditioned
gating to choose the local expert [37, 38, 39].

Two categories of LR exist [40], _lazy_ learners that maintain all seen data
points in memory [29, 31], and _memoryless_ learners that compress data by
constructing basis functions in the input space and fitting and storing
locally parameterized regression models [35]. Prominent examples of the latter
include receptive field weighted regression (RFWR) [30] and locally weighted
projection regression (LWPR) [32]. However, these methods are often difficult
to tune as they possess many hyperparameters.

A limited attempt at a Bayesian treatment of LR is made in [41] by
constructing local nonparametric kernels and placing gamma priors on the
kernel widths to alleviate the need to tune the basis functions. This approach
leads to a localized GP formulation that needs to retain the training data in
memory, again leading to the computational issues of vanilla GPR. Local
Gaussian regression (LGR) is a further Bayesian generalization of LR [42]. The
authors treat the local models in a Bayesian framework and couple them via the
loss function that reinforces global coordination. Nonetheless, both
approaches rely on heuristics and thresholds for adding and pruning local
models and fall short of formulating a generative model over input and output.

Following this introduction, we believe that local regression with a proper
generative Bayesian treatment can potentially serve as a powerful general-
purpose function approximator. Moreover, as a probabilistic and efficient
representation, it can drive many low-level applications in control and
robotics that favor fast predictions and do not require deep representations.

We introduce two probabilistic graphical mixture models for local regression
in the upcoming sections. The first, infinite local regression (ILR) [43], is
a generative formulation that relies on the paradigm of Bayesian
nonparametrics (BNP) [44] to automatically grow the mixture size based on
observed data. This technique ultimately results in a general formulation of
related methods that alleviates the need for any heuristic considerations.
However, despite the effectiveness of ILR, and like other local regression
techniques that rely on locally linear or polynomial approximations, it
maintains a one-to-one correspondence between the activations and local
regression units. This effect limits the model’s capacity to share parameters
across the input space and often forces the generation of duplicate
components, needlessly increasing the overall number of parameters. To address
this limitation, we introduce hierarchical infinite local regression (HILR), a
multi-level development of ILR that enables multi-modal activations of the
same regression component, giving the model a structure that allows sharing of
regression parameters across repeating local patterns in the data. This
architecture increases the flexibility of the representation and contributes
towards its compression.

For learning these models, we derive two general variational Bayes (VB)
schemes [45] that efficiently infer the posterior parameters and overcome the
need for computationally heavy sampling methods. We benchmark the models on a
range of tasks highlighting their strengths, such as dealing with
heteroscedasticity, non-continuous functions, and multi-modal activation.
Additionally, we test on large real-world high-dimensional datasets to
benchmark learning the inverse dynamics of robotic manipulators. Most
importantly, we deploy an instance of ILR to perform inverse dynamics control
on a real Barrett-WAM manipulator.

Alternative Bayesian extensions of the generative mixture of experts exist in
the literature. A prominent area of research focused on Gaussian processes
mixtures, in which local components are modeled by separate GP s [46, 47].
Given that a single Gaussian process is an excellent approximator of nonlinear
trends, the motivation for constructing such experts is not to improve the
quality of approximation but rather to reduce the computational complexity and
memory requirements of GP s by associating every component only with a slice
of the overall data, thus limiting the negative effects of cubic scaling.
While this motivation is understandable from a computational standpoint, it is
often unclear how to decompose a nonlinear function into a mixture of
_nonlinear_ experts, especially when each expert has the representational
flexibility of a GP. In contrast, models like ILR and HILR rely on locally
linear experts that offer a natural unit of local approximation.

Finally, the infinite mixtures used in ILR and HILR to regularize model
complexity are rooted in Bayesian nonparametrics. Here we reference
influential work on Markov chain Monte Carlo (MCMC) techniques for Bayesian
nonparametric density estimation [48, 49, 50, 51], which developed the first
seeds of Bayesian inference for Dirichlet processes [52, 53]. These concepts
inspired comparable infinite mixture regression models. However, these
attempts exclusively relied on expensive Gibbs sampling algorithms [54, 55,
56, 57]. Instead, we focus on developing efficient deterministic variational
inference (VI) algorithms that dramatically improve the practical aspects of
training and deploying Bayesian finite and infinite mixture models [58, 59].
For a detailed discussion on the scalability of different paradigms of
Bayesian inference, we refer the reader to [60].

This paper is organized as follows. In Section 2, we review Bayesian linear
regression, finite and infinite mixture models, and variational inference.
Section 3 presents the infinite local regression (ILR)model and a variational
Bayes inference scheme [43]. In Section 4, we introduce the hierarchical
infinite local regression (HILR), an extension of ILR, that incorporates
multi-modal activations and parameter sharing. Finally, in Section 5, we test
ILR and HILR on a set of simulated and real-world benchmarks.

##  2 Preliminaries

This section introduces related concepts, such as Bayesian linear regression,
Bayesian mixture models, Dirichlet processes, and variational inference.

###  2.1 Bayesian Linear Regression

We start by discussing the Bayesian treatment of a single component of a
Bayesian local regression model, namely Bayesian linear regression [61]. The
conditional data likelihood takes a feature vector
$\bm{\mathbf{x}}\in\mathds{R}^{m}$ as a random input variable and returns a
random response variable $\bm{\mathbf{y}}\in\mathds{R}^{d}$ according to a
linear mapping
$\bm{\mathbf{A}}\mathrel{\mathop{\ordinarycolon}}\mathds{R}^{m}\to\mathds{R}^{d}$,
a bias vector $\bm{\mathbf{c}}$, and additive zero-mean noise with a precision
matrix $\bm{\mathbf{V}}$

$\bm{\mathbf{y}}=\bm{\mathbf{A}}\bm{\mathbf{x}}+\bm{\mathbf{c}}+\bm{\mathbf{e}},\quad\bm{\mathbf{e}}\sim\operatorname{N}(\bm{\mathbf{0}},\bm{\mathbf{V}}).$

For a fully Bayesian treatment, we consider all parameters of this model to be
random variables on which we place proper conjugate priors. In this case, we
place matrix-normal (MN) and normal-Wishart (NW) priors on the matrix
$\bm{\mathbf{A}}$, the bias coefficient $\bm{\mathbf{c}}$, and precision
matrix $\bm{\mathbf{V}}$

$p(\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})=\operatorname{MN}(\bm{\mathbf{A}}\,|\,\bm{\mathbf{M}},\bm{\mathbf{V}},\bm{\mathbf{K}})\operatorname{N}(\bm{\mathbf{c}}\,|\,\bm{\mathbf{\theta}},\rho\bm{\mathbf{V}})\operatorname{W}(\bm{\mathbf{V}}\,|\,\bm{\mathbf{\Phi}},\eta),$

where $\bm{\mathbf{M}}$, the mean of $\bm{\mathbf{A}}$, is a $d\times m$
matrix and $\bm{\mathbf{V}}$ and $\bm{\mathbf{K}}$ are $d\times d$ and
$m\times m$ that serve as row and column precision matrices of
$\bm{\mathbf{A}}$, respectively. The mean $\bm{\mathbf{\theta}}$ is an
$m$-dimensional vector, and the scalar $\rho$ modulates the amplitude of the
precision. Finally, the parameters of the Wishart distribution are the
$d\times d$ positive definite scale matrix $\bm{\mathbf{\Phi}}$ and the
degrees of freedom $\eta$. Due to the conjugate nature of the priors, the
joint posteriors
$p(\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}}\,|\,\mathcal{D})$ are
matrix-normal and normal-Wishart distributions, conditioned on the data of $N$
independent and identically distributed data pairs
$\mathcal{D}=\\{(\bm{\mathbf{x}}_{1},\bm{\mathbf{y}}_{1}),\ldots,(\bm{\mathbf{x}}_{N},\bm{\mathbf{y}}_{N})\\}$.

###  2.2 Bayesian Finite Mixture Models

Gaussian mixture models (GMMs) are hierarchical latent variable models with
universal approximation capabilities for arbitrary continuous densities. This
insight is of central interest when connected to density estimation for local
regression models, which are themselves universal nonlinear function
approximators [5]. A finite $K$-component Gaussian mixture of a random
variable $\bm{\mathbf{x}}$ is a weighted linear combination of densities

$\vspace{-0.1cm}p(\bm{\mathbf{x}}\,|\,\bm{\mathbf{\theta}})=\sum_{k=1}^{K}p(\bm{\mathbf{z}}=k\,|\,\bm{\mathbf{\pi}})p(\bm{\mathbf{x}}\,|\,\bm{\mathbf{\theta}}_{k})=\sum_{k=1}^{K}\pi_{k}\operatorname{N}(\bm{\mathbf{x}}\,|\,\bm{\mathbf{\mu}}_{k},\bm{\mathbf{\Lambda}}_{k}),$

with $K$ unique mean vectors $\bm{\mathbf{\mu}}_{k}$ and precision matrices
$\bm{\mathbf{\Lambda}}_{k}$. The latent quantity $\bm{\mathbf{z}}$ is a one-
hot random variable distributed according to a categorical distribution
$p(\bm{\mathbf{z}})=\operatorname{Cat}(\bm{\mathbf{\pi}})$, governed by the
weights $\bm{\mathbf{\pi}}=\\{\pi_{1},\dots,\pi_{K}\\}$ that satisfy
$0\leq\pi_{k}\leq 1$ and $\sum_{k=1}^{K}\pi_{k}=1$.

The Bayesian extension of this model [58] introduces a conjugate normal-
Wishart prior on the means and precision matrices
$(\bm{\mathbf{\mu}}_{k},\bm{\mathbf{\Lambda}}_{k})\sim\operatorname{NW}(\bm{\mathbf{\lambda}})$,
where $\bm{\mathbf{\lambda}}$ contains the hyperparameters. Furthermore, a
conjugate Dirichlet prior, with a concentration parameter $\alpha$, is placed
on the mixing weights
$\bm{\mathbf{\pi}}\sim\operatorname{Dir}(\bm{\mathbf{\alpha}})$. This Bayesian
perspective has proven effective in regularizing the shortcomings of GMMs by
allowing superfluous components to fall back onto their priors instead of
severely over-fitting to small clusters. This effect can be understood as
sparsification bias over $K$ [62, 63].

###  2.3 Dirichlet Process and Stick-Breaking

Scaling the finite mixtures in Section 2.2 to infinite components requires
placing a nonparametric Dirichlet process (DP) prior on the parameters of the
mixture. A DP is a distribution over probability measures $\operatorname{G}$.
We write $\operatorname{G}\sim\operatorname{DP}(\alpha,\operatorname{H})$,
where $\alpha$ is the concentration parameter and $\operatorname{H}$ is the
base measure [64, 65]. Intuitively, a Dirichlet process is a distribution over
distributions, meaning each draw $\operatorname{G}$ is itself a distribution.
The base distribution $\operatorname{H}$ is the mean of the DP, and the
concentration parameter $\alpha$ can be interpreted as an inverse variance.
The larger $\alpha$ is, the smaller the variance, and the process concentrates
more of its mass around the mean distribution $\operatorname{H}$.

We will rely on the stick-breaking construction [66] of a DP as an algorithmic
realization. Stick-breaking delivers an infinite sequence of mixture weights
$\pi_{k}$ of an infinite mixture model from the stochastic process

$\pi_{k}=s_{k}\prod_{l=1}^{k-1}\left(1-s_{l}\right),\quad
s_{k}\sim\operatorname{Beta}(1,\alpha).$

This process is sometimes denoted as
$\bm{\mathbf{\pi}}\sim\operatorname{GEM}(\alpha)$ [64]. The stick-breaking
procedure describes how the random variables $s_{k}$, representing stick
lengths, are drawn from a beta distribution and combined to obtain the mixture
weights $\pi_{k}$. If the concentration parameter $\alpha$ increases, the
magnitude of the mixing weights $\pi_{k}$ decreases on average, and the number
of possible active components increases.

This representation of DPs can be used to replace the priors placed on the
finite Gaussian mixture model [59]. In such a setting, the base
$\operatorname{H}$ is a normal-Wishart distribution, and the sampled measure
$\operatorname{G}\sim\operatorname{DP}(\alpha,\operatorname{NW})$ is a draw of
an unbounded number of parameters
$(\bm{\mathbf{\mu}}_{k},\bm{\mathbf{\Lambda}}_{k})\sim\operatorname{NW}(\bm{\mathbf{\lambda}})$
for an infinite number of clusters, associated with an infinite number of
weights $\pi_{k}$ generated by the stick-breaking process. These draws from a
Dirichlet process are discrete with probability one, which leads to the
clustering effect of the DP. Eventually, the same parameters will be sampled
over and over, forcing the associated data points to cluster.

###  2.4 Variational Inference of Structured Models

MCMC [67], and VI [68] have become the two main approaches for (approximate)
probabilistic inference in graphical models. While MCMC constructs a
stochastic sampling process that converges to the posterior, VI formulates the
inference task as a deterministic optimization problem. Despite its reliance
on a coarser functional posterior approximation, VI is often preferred in
settings with large number of parameters. Moreover, the deterministic nature
of VI circumvents the issue of label switching in Monte-Carlo-based posterior
inference of mixture models.

In a nutshell, in variational inference, a typically intractable posterior is
approximated by a tractable functional distribution $q(\bm{\mathbf{\beta}})$
that minimizes the Kullback-Leibler divergence (KL) to true posterior
$p(\bm{\mathbf{\beta}}\,|\,\mathcal{D})$

$q^{*}(\bm{\mathbf{\beta}})=\operatorname*{arg\,min}\quad\operatorname{KL}(q(\bm{\mathbf{\beta}})\,||\,p(\bm{\mathbf{\beta}}\,|\,\mathcal{D))},$

where $\mathcal{D}$ is observed data and the vector $\bm{\mathbf{\beta}}$
subsumes both the parameters $\bm{\mathbf{\theta}}$ and hidden indicators
$\bm{\mathbf{z}}$ in structured models. Note that the KL in the above
objective is mode-seeking, meaning it will lock on one mode of a possibly
multi-modal posterior. In general, the posterior
$p(\bm{\mathbf{\beta}}\,|\,\mathcal{D})$ is unknown, because the normalizer
$p(\mathcal{D})$ is not tractable. In consequence, the KL cannot be minimized
directly but rather optimized via a related objective that is equal up to the
constant term equivalent to the evidence $p(\mathcal{D})$

$\displaystyle\operatorname{KL}(q\,||\,p)$
$\displaystyle=\mathbb{E}_{q}\left[\log
q(\bm{\mathbf{\beta}})\right]-\mathbb{E}_{q}\left[\log
p(\bm{\mathbf{\beta}}\,|\,\mathcal{D})\right]$
$\displaystyle=\mathbb{E}_{q}\left[\log
q(\bm{\mathbf{\beta}})\right]-\mathbb{E}_{q}\left[\log
p(\mathcal{D},\bm{\mathbf{\beta}})\right]+p(\mathcal{D})$
$\displaystyle=\mathbb{E}_{q}\left[\log
q(\bm{\mathbf{\beta}})\right]-\mathbb{E}_{q}\left[\log
p(\mathcal{D},\bm{\mathbf{\beta}})\right]+\mathrm{const}.$

This modified objective is denoted as the negative evidence lower bound (ELBO)
and can be reformulated to take the traditional maximization objective of VI
algorithms

$\displaystyle\mathrm{ELBO}(q)$ $\displaystyle=\mathbb{E}_{q}\left[\log
p(\mathcal{D},\bm{\mathbf{\beta}})\right]-\mathbb{E}_{q}\left[\log
q(\bm{\mathbf{\beta}})\right]$ $\displaystyle=\mathbb{E}_{q}\left[\log
p(\mathcal{D}\,|\,\bm{\mathbf{\beta}})\right]+\mathbb{E}_{q}\left[\log
p(\bm{\mathbf{\beta}})\right]-\mathbb{E}_{q}\left[\log
q(\bm{\mathbf{\beta}})\right]$ $\displaystyle=\mathbb{E}_{q}\left[\log
p(\mathcal{D}\,|\,\bm{\mathbf{\beta}})\right]-\operatorname{KL}(q(\bm{\mathbf{\beta}})\,||\,p(\bm{\mathbf{\beta}})).$

The choice of the approximate posterior distribution $q(\bm{\mathbf{\beta}})$
is open. In this paper, we focus on variational Bayes methods that rely on the
(structured) mean-field assumption as a general recipe for maximizing the ELBO
[69, 45]. This approximation requires that the posterior factorizes over the
set of the hidden variables
$q(\bm{\mathbf{\beta}})=\prod_{i=1}^{M}q_{i}(\beta_{i})$. It is emphasized
that no other assumptions are made about $q(\bm{\mathbf{\beta}})$. The
resulting posterior will be determined solely by the assumed likelihood and
priors.

Specifically, we follow the scheme of variational Bayes expectation-
maximization (VBEM) as a probabilistic generalization of expectation-
maximization (EM). This approach constitutes a coordinate ascent scheme that
iteratively optimizes the ELBO for individual factors of the approximate
posterior $q(\bm{\mathbf{\beta}})$ while holding the others constant

$\ln q_{j}(\beta_{j})=\mathbb{E}_{q_{i\neq j}}\left[\log
p(\mathcal{D},\bm{\mathbf{\beta}})\right]+\mathrm{const}.$

A more practical version of this optimization can be achieved by using
stochastic variational inference (SVI), a batched stochastic gradient ascent
approach. In the case of the conjugate exponential family, SVI not only
facilitates scalability over large datasets but also resembles a natural
gradient ascent algorithm on the ELBO with favorable convergence properties
[70].

[width=400pt, svgpath=figures/plates/]svg-inkscape/ilr_horiz_svg-tex.pdf_tex

Figure 2: A unified plate notation for infinite mixtures of Bayesian local
regression. Assuming Gaussian and linear Gaussian densities, the basis
parameters $(\bm{\mathbf{\mu}}_{k},\bm{\mathbf{\Lambda}}_{k})$ are sampled
from a normal-Wishart distribution, while the regression parameters
$(\bm{\mathbf{A}}_{k}$, $\bm{\mathbf{c}}_{k})$ and precision matrices
$\bm{\mathbf{V}}_{k}$ are sampled from a matrix-Normal-Wishart, for every
component $k$. The latent variables $\bm{\mathbf{z}}_{n}$ assign every
$\bm{\mathbf{x}}_{n}$ and $\bm{\mathbf{y}}_{n}$ to a component and are drawn
from a categorical distribution parameterized by $\bm{\mathbf{\pi}}$. The
mixture weights $\bm{\mathbf{\pi}}$ are generated by a stick-breaking process
with a concentration parameter $\alpha$.

##  3 The Infinite Mixture of Local Linear Regression Model

Using the previously presented concepts of Bayesian linear regression,
Bayesian mixture models, and Dirichlet processes, we now construct the
Bayesian infinite local regression (ILR)model.

Our approach to tackling the regression problem is a Bayesian joint density
estimation technique over the input $\bm{\mathbf{x}}$ and target
$\bm{\mathbf{y}}$, both conditioned on a latent discrete label
$\bm{\mathbf{z}}$. To avoid the need for a fixed number of components over
$\bm{\mathbf{z}}$, we assume a generative model driven by a stick-breaking
process with a concentration parameter $\alpha$ and a base distribution
$\operatorname{H}$, which is a product of a normal-Wishart, a Matrix-normal,
and Wishart distribution.

The data generation proceeds as depicted in Figure 2. The stick-breaking
process is sampled to generate the categorical weights $\bm{\mathbf{\pi}}$,
mixture activations
$\\{\bm{\mathbf{\mu}}_{k},\bm{\mathbf{\Lambda}}_{k}\\}_{k=1}^{\infty}$, and
regression parameters
$\\{\bm{\mathbf{A}}_{k},\bm{\mathbf{c}}_{k},\bm{\mathbf{V}}_{k}\\}_{k=1}^{\infty}$

$\begin{gathered}\bm{\mathbf{\pi}}(\bm{\mathbf{s}})\sim\operatorname{GEM}(\alpha),\\\
\bm{\mathbf{\Lambda}}_{k}\sim\operatorname{W}(\bm{\mathbf{\Psi}},\nu),\,\,\bm{\mathbf{\mu}}_{k}\sim\operatorname{N}(\bm{\mathbf{m}},\kappa\bm{\mathbf{\Lambda}}_{k}),\\\
\bm{\mathbf{V}}_{k}\sim\operatorname{W}(\bm{\mathbf{\Phi}},\eta),\,\,\bm{\mathbf{A}}_{k}\sim\operatorname{MN}(\bm{\mathbf{M}},\bm{\mathbf{K}},\bm{\mathbf{V}}_{k}),\,\,\bm{\mathbf{c}}_{k}\sim\operatorname{N}(\bm{\mathbf{\theta}},\rho\bm{\mathbf{V}}_{k}),\end{gathered}$

and those parameters are used to generate the one-hot labels
$\bm{\mathbf{z}}_{n}$ and input-output data pairs
$\bm{\mathbf{x}}_{n},\bm{\mathbf{y}}_{n}$

$\displaystyle\bm{\mathbf{z}}_{n}\sim\operatorname{Cat}(\bm{\mathbf{\pi}}(\bm{\mathbf{s}})),\,\,\bm{\mathbf{x}}_{n}\sim\operatorname{N}(\bm{\mathbf{\mu}}_{z_{n}},\bm{\mathbf{\Lambda}}_{z_{n}}),$
$\displaystyle\bm{\mathbf{y}}_{n}\sim\operatorname{N}(\bm{\mathbf{A}}_{z_{n}}\bm{\mathbf{x}}_{n}+\bm{\mathbf{c}}_{z_{n}},\bm{\mathbf{V}}_{z_{n}}).$

Notice that the densities over the input space, parameterized by
$(\bm{\mathbf{\mu}}_{k},\bm{\mathbf{\Lambda}}_{k})$, naturally play the role
of basis functions or so-called receptive fields as in the receptive field
weighted regression (RFWR) [31] and locally weighted projection regression
(LWPR) [32] algorithms.

Given these modeling assumptions, the regression objective can be cast as a
joint density inference problem of the posterior
$p(\bm{\mathbf{Z}},\bm{\mathbf{s}},\bm{\mathbf{\mu}},\allowbreak\bm{\mathbf{\Lambda}},\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}}\,|\,\bm{\mathbf{X}},\bm{\mathbf{Y}})$,
where $\bm{\mathbf{X}}=\\{\bm{\mathbf{x}}_{1},\dots,\bm{\mathbf{x}}_{N}\\}$,
$\bm{\mathbf{Y}}=\\{\bm{\mathbf{y}}_{1},\dots,\bm{\mathbf{y}}_{N}\\}$, and
$\bm{\mathbf{Z}}=\\{\bm{\mathbf{z}}_{1},\dots,\bm{\mathbf{z}}_{N}\\}$. The
result of this inference is then leveraged to perform predictions via
conditioning and marginalization.

###  3.1 Variational Bayes for ILR

For inference, we focus on a variational Bayes EM algorithm that alternates
between a variational expectation step (E-step) and a maximization step
(M-step). Deriving such an algorithm for this model requires pinning down the
following definitions of the likelihood, prior, and posterior.

Complete Data Likelihood. For the general case of multivariate regression with
$m$ inputs and $d$ outputs, we assume the following structured joint
likelihood over the inputs, outputs, and indicator variables

$\displaystyle p(\bm{\mathbf{X}},\bm{\mathbf{Y}},\bm{\mathbf{Z}}\,|\,.)$
$\displaystyle=p(\bm{\mathbf{Z}})\,p(\bm{\mathbf{X}}\,|\,\bm{\mathbf{Z}})\,p(\bm{\mathbf{Y}}\,|\,\bm{\mathbf{X}},\bm{\mathbf{Z}})$
$\displaystyle=\begin{aligned}
&\prod_{n=1}^{N}\operatorname{Cat}(\bm{\mathbf{z}}_{n}\,|\,\bm{\mathbf{\pi}}(\bm{\mathbf{s}}))\\\
&\times\prod_{n=1}^{N}\prod_{k=1}^{\infty}\operatorname{N}(\bm{\mathbf{x}}_{n}\,|\,\bm{\mathbf{\mu}}_{k},\bm{\mathbf{\Lambda}}_{k})^{z_{nk}}\\\
&\times\prod_{n=1}^{N}\prod_{k=1}^{\infty}\operatorname{N}(\bm{\mathbf{y}}_{n}\,|\,\bm{\mathbf{A}}_{k}\bm{\mathbf{x}}_{n}+\bm{\mathbf{c}}_{k},\bm{\mathbf{V}}_{k})^{z_{nk}},\end{aligned}$

where the dimensions of all quantities follow the notation of Bayesian linear
regression.

Infinite Conjugate Prior. We construct a factorized conjugate infinite mixture
prior

$p(\bm{\mathbf{s}},\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}},\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})=\begin{aligned}
&p(\bm{\mathbf{s}})\,p(\bm{\mathbf{\mu}}|\bm{\mathbf{\Lambda}})\,p(\bm{\mathbf{\Lambda}})\\\
&\times
p(\bm{\mathbf{A}}|\bm{\mathbf{V}})\,p(\bm{\mathbf{c}}|\bm{\mathbf{V}})\,p(\bm{\mathbf{V}}).\end{aligned}$

This prior assumes that the cluster means $\bm{\mathbf{\mu}}_{k}$ and
precision matrices $\bm{\mathbf{\Lambda}}_{k}$ are sampled from normal-Wishart
distributions

$p(\bm{\mathbf{\mu}}\,|\,\bm{\mathbf{\Lambda}})\,p(\bm{\mathbf{\Lambda}})=\prod_{k=1}^{\infty}\operatorname{N}(\bm{\mathbf{\mu}}_{k}\,|\,\bm{\mathbf{m}}_{0},\kappa_{0}\bm{\mathbf{\Lambda}}_{k})\operatorname{W}(\bm{\mathbf{\Lambda}}_{k}\,|\,\bm{\mathbf{\Psi}}_{0},\nu_{0}),$

while matrix-normal-Wishart and normal-Wishart priors are placed on the
regression coefficients $(\bm{\mathbf{A}}_{k},\bm{\mathbf{c}}_{k})$ and the
precision matrices $\bm{\mathbf{V}}_{k}$

$\displaystyle
p(\bm{\mathbf{A}}\,|\,\bm{\mathbf{V}})p(\bm{\mathbf{c}}\,|\,\bm{\mathbf{V}})p(\bm{\mathbf{V}})=\prod_{k=1}^{\infty}$
$\displaystyle\operatorname{MN}(\bm{\mathbf{A}}_{k}\,|\,\bm{\mathbf{M}}_{0},\bm{\mathbf{K}}_{0},\bm{\mathbf{V}}_{k})$
$\displaystyle\times\operatorname{N}(\bm{\mathbf{c}}_{k}\,|\,\bm{\mathbf{\theta}}_{0},\rho_{0}\bm{\mathbf{V}}_{k})\operatorname{W}(\bm{\mathbf{V}}_{k}\,|\,\bm{\mathbf{\Phi}}_{0},\eta_{0}).$

The parameters $\pi_{k}$ are generated by a stick-breaking process
$\pi_{k}(\bm{\mathbf{s}})=s_{k}\prod_{l=1}^{k-1}(1-s_{l})$, where the stick
lengths $\bm{\mathbf{s}}$ are independently beta distributed

$p(\bm{\mathbf{s}})=\prod_{k=1}^{\infty}\operatorname{Beta}(s_{k}\,|\,1,\alpha_{0}).$

Truncated Mean-Field Posterior. We rely on a structured mean-field
approximation of the posterior [69] that factorizes between the labels
$q(\bm{\mathbf{Z}})$ and the remaining parameters
$q(\bm{\mathbf{s}},\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}},\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})$,
thus automatically leading to the decomposition $p(.\,|\,\mathcal{D})\approx
q(\bm{\mathbf{Z}})\,q(\bm{\mathbf{s}})\,q(\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})\,q(\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})$.
Further, we follow [59] by allowing a truncation of the posterior while
maintaining an infinite prior, so that $q(s_{\scriptscriptstyle
K}\\!=\\!1)=1$, implying that $\pi_{k}=0$ for $k>K$

$p(.\,|\,\mathcal{D})\approx\begin{aligned}
&\prod_{n=1}^{N}\operatorname{Cat}(\bm{\mathbf{z}}_{n}\,|\,\bm{\mathbf{r}}_{n})\prod_{k=1}^{K-1}\operatorname{Beta}(s_{k}\,|\,\gamma_{k},\alpha_{k})\\\
&\times\prod_{k=1}^{K}\operatorname{N}(\bm{\mathbf{\mu}}_{k}\,|\,\bm{\mathbf{m}}_{k},\kappa_{k}\bm{\mathbf{\Lambda}}_{k})\operatorname{W}(\bm{\mathbf{\Lambda}}_{k}\,|\,\bm{\mathbf{\Psi}}_{k},\nu_{k})\\\
&\times\begin{aligned}
\prod_{k=1}^{K}&\operatorname{MN}(\bm{\mathbf{A}}_{k}\,|\,\bm{\mathbf{M}}_{k},\bm{\mathbf{K}}_{k},\bm{\mathbf{V}}_{k})\\\
&\times\operatorname{N}(\bm{\mathbf{c}}_{k}\,|\,\bm{\mathbf{\theta}}_{k},\rho_{k}\bm{\mathbf{V}}_{k})\operatorname{W}(\bm{\mathbf{V}}_{k}\,|\,\bm{\mathbf{\Phi}}_{k},\eta_{k}),\end{aligned}\end{aligned}$

where $\bm{\mathbf{r}}_{n}$ are the expected responsibilities of the mixture.
During evaluation, the truncation threshold $K$ is chosen to be very high and
is seldom reached.

Variational Expectation Step. In the E-step, the responsibilities are computed
by following the recipe of VBEM in Section 2.4. The responsibilities are the
variational parameters of the posterior categorical

$\displaystyle\log q(\bm{\mathbf{Z}})$ $\displaystyle=\begin{aligned}
&\mathbb{E}_{q(\bm{\mathbf{s}})}\left[\log
p(\bm{\mathbf{Z}}\,|\,\bm{\mathbf{\pi}}(\bm{\mathbf{s}}))\right]+\mathbb{E}_{q(\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})}\left[\log
p(\bm{\mathbf{X}}\,|\,\bm{\mathbf{Z}})\right]\\\
&+\mathbb{E}_{q(\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})}\left[\log
p(\bm{\mathbf{Y}}\,|\,\bm{\mathbf{X}},\bm{\mathbf{Z}})\right]+\operatorname{const}\end{aligned}$
$\displaystyle=\sum_{n=1}^{N}\sum_{k=1}^{K}z_{nk}\log r_{nk}.$

The expectations associated with the data likelihoods of $\bm{\mathbf{x}}_{n}$
and $\bm{\mathbf{y}}_{n}$ can be straightforwardly computed [71]. However, the
expectations associated with infinite-dimensional categorical require more
careful consideration [59]. We provide the necessary details in the appendix.

Variational Maximization Step. The M-step updates the remaining variational
distributions given an approximation of the categorical posterior as follows

$\displaystyle\log q(\bm{\mathbf{s}})$
$\displaystyle=\mathbb{E}_{q(\bm{\mathbf{Z}})}\left[\log
p(\bm{\mathbf{Z}}\,|\,\bm{\mathbf{\pi}}(\bm{\mathbf{s}}))\right]+\log
p(\bm{\mathbf{s}})+\operatorname{const},$ $\displaystyle\log
q(\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})$
$\displaystyle=\mathbb{E}_{q(\bm{\mathbf{Z}})}\left[\log
p(\bm{\mathbf{X}}\,|\,\bm{\mathbf{Z}})\right]+\log
p(\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})+\operatorname{const},$
$\displaystyle\log q(\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})$
$\displaystyle=\begin{aligned} &\mathbb{E}_{q(\bm{\mathbf{Z}})}\left[\log
p(\bm{\mathbf{Y}}\,|\,\bm{\mathbf{X}},\bm{\mathbf{Z}})\right]\\\ &+\log
p(\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})+\operatorname{const}.\end{aligned}$

Each of these updates reflects a conjugate computation of $K$ log-posterior
densities given a log-prior and a log-likelihood weighted by the
responsibilities $r_{nk}=\mathbb{E}[z_{nk}]$. We provide more details and
general computational recipes based on conjugacy rules in the appendix.

###  3.2 Posterior Predictive Distribution

For predicting the function value $\hat{\bm{\mathbf{y}}}$ conditioned on a
test query $\hat{\bm{\mathbf{x}}}$, we marginalize the likelihood over the
posterior parameters to get the joint posterior predictive density. To make
the marginalization tractable, we replace the true posterior with our
approximate variational posterior $q(.\,|\,\mathcal{D})$ inferred under a
training dataset $\mathcal{D}$. The conditional predictive for a single
component $\bm{\mathbf{z}}=k$ is a conditional multivariate Student’s
t-distribution of the form

$\displaystyle
p(\hat{\bm{\mathbf{y}}}\,|\,\hat{\bm{\mathbf{x}}},\hat{\bm{\mathbf{z}}}\\!=\\!k,\mathcal{D})$
$\displaystyle=\mathbb{E}_{q(\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})}\left[p(\hat{\bm{\mathbf{y}}}\,|\,\hat{\bm{\mathbf{x}}},\bm{\mathbf{A}}_{k},\bm{\mathbf{c}}_{k},\bm{\mathbf{V}}_{k})\right]$
$\displaystyle=\operatorname{T}\left(\bm{\mathbf{M}}_{k}\hat{\bm{\mathbf{x}}}+\bm{\mathbf{\theta}}_{k},a_{k}\bm{\mathbf{\Phi}}_{k},\eta_{k}+1\right),$

where we have defined

$a_{k}=1-\bm{\mathbf{u}}^{\top}\left(\bm{\mathbf{L}}_{k}+\bm{\mathbf{u}}\bm{\mathbf{u}}^{\top}\right)^{-1}\bm{\mathbf{u}},$

with $\bm{\mathbf{u}}=[\hat{\bm{\mathbf{x}}},\,1]^{\top}$ and
$\bm{\mathbf{L}}_{k}=\operatorname{Block}\left(\bm{\mathbf{K}}_{k},\rho_{k}\right)$.

Additionally, the joint activation of a component $k$ is a Student’s
t-distribution weighted by the expected categorical probability under the
posterior stick-breaking process

$\displaystyle
p(\hat{\bm{\mathbf{x}}},\hat{\bm{\mathbf{z}}}\\!=\\!k\,|\,\mathcal{D})$
$\displaystyle\propto\mathbb{E}_{q(\bm{\mathbf{s}})}\left[p(\hat{\bm{\mathbf{z}}}\\!=\\!k\,|\,\bm{\mathbf{\pi}}(\bm{\mathbf{s}}))\right]\mathbb{E}_{q(\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})}\left[p(\hat{\bm{\mathbf{x}}}\,|\,\bm{\mathbf{\mu}}_{k},\bm{\mathbf{\Lambda}}_{k})\right]$
$\displaystyle=\begin{aligned}
&\frac{\gamma_{k}}{\gamma_{k}+\alpha_{k}}\prod_{l=1}^{k-1}\left(1-\frac{\gamma_{l}}{\gamma_{l}+\alpha_{l}}\right)\\\
&\times\operatorname{T}\left(\bm{\mathbf{\mu}}_{k},\frac{\kappa_{k}}{1+\kappa_{k}}\bm{\mathbf{\Psi}}_{k},\nu_{k}+1\right).\end{aligned}$

These $K$-activation probabilities enable two prediction techniques. A mode-
prediction, where the most likely active component is selected and used to
perform prediction with the corresponding linear regression model, or a mean-
prediction, that averages the predictions of all components weighted by their
activation probabilities.

###  3.3 Computational Complexity

We calculate the training-time computational cost to be
$\mathcal{O}(NK(d+m)^{3})$, which can be straightforwardly reduced to
$\mathcal{O}(LK(d+m)^{3})$ by applying stochastic updates [70], where $L$ is
the batch size. This result shows linear scalability with the data, which is
considerably more efficient than simple variants of GPR. The test-time
complexity of a mean prediction is $\mathcal{O}(K(d^{3}+dm))$, which combines
the input membership query and the linear matrix transformation for every
model $k$. This computation is, in contrast to GPR, independent of the
training data size, hence the advantage of memoryless locally-parametric
representations during real-time critical applications.

[width=450pt, svgpath=figures/plates/]svg-inkscape/hilr_horiz_svg-tex.pdf_tex

Figure 3: A unified plate notation for hierarchical infinite tied mixtures of
Bayesian local regression models. This model outlines a two-level architecture
that allows sharing of parameters between single components in order to
compress the representation. It can be interpreted as a local regression model
with multi-modal activations. Each unit of the upper level is itself a mixture
of local regression models that share the same slope $\bm{\mathbf{A}}_{m}$ and
output precision $\bm{\mathbf{V}}_{m}$. Each of these $m$ different slopes can
be activated at $k$ unique lower-level input regions centered around
$\bm{\mathbf{\mu}}_{mk}$ and tied via a shared input precision
$\bm{\mathbf{\Lambda}}_{m}$. Both the upper- and lower-level mixtures are
governed by independent Dirichlet process priors.

##  4 The Hierarchical Infinite Mixture of Local Linear Regression Model

The local regression model presented in Section 3 offers a very flexible and
well-regularized alternative to previously developed approaches [30, 32, 42].
However, like similar local regression representations, it cannot account for
shift-invariance across the input space. This drawback can cause it to
generate duplicate regression components to account for similar local function
trends across disconnected regions of the input space. This limitation results
from the hierarchical design that directly couples activations and local
function units via a one-to-one correspondence and enforces a uni-modal
activation per regression component. This coupling can be observed in the
definition of the likelihood in Section 3, where the activation and the local
regression units share the same assignment variable.

It stands to reason that ILR does not offer enough flexibility and hinders
parameter sharing between components. Therefore, we present a modified
formulation of ILR that explicitly accounts for shift-invariance in the input
space and provides the freedom to create regression units with multi-modal,
theoretically infinitely-modal, activations, if needed. The resulting model,
hierarchical infinite local regression (HILR), is an infinite mixture over
infinite mixtures that shares similarities with existing representations
developed for hierarchical clustering [72, 73, 74].

We start by describing the generative process of the hierarchical mixture over
mixtures as depicted in Figure 3. HILR consists of two levels. At the upper
level, a meta Dirichlet process generates the stick-breaking weights
$\bm{\mathbf{\omega}}$, the meta-activations
$\\{\bm{\mathbf{\tau}}_{m},\bm{\mathbf{\Lambda}}_{m}\\}_{m=1}^{\infty}$, and
the shared slope and output precision matrices
$\\{\bm{\mathbf{A}}_{m},\bm{\mathbf{V}}_{m}\\}_{m=1}^{\infty}$

$\begin{gathered}\bm{\mathbf{\omega}}(\bm{\mathbf{t}})\sim\operatorname{GEM}(\beta),\\\
\bm{\mathbf{\Lambda}}_{m}\sim\operatorname{W}(\bm{\mathbf{\Psi}},\nu),\,\,\bm{\mathbf{\tau}}_{m}\sim\operatorname{N}(\bm{\mathbf{m}},\lambda\bm{\mathbf{\Lambda}}_{m}),\\\
\bm{\mathbf{V}}_{m}\sim\operatorname{W}(\bm{\mathbf{\Phi}},\eta),\,\,\bm{\mathbf{A}}_{m}\sim\operatorname{MN}(\bm{\mathbf{M}},\bm{\mathbf{K}},\bm{\mathbf{V}}_{m}),\end{gathered}$

where $\beta$ is the concentration parameter of the upper-level DP and
$\bm{\mathbf{t}}$ are the generated stick lengths.

HILR’s ability to account for shift-invariance is a result of allowing multi-
modal activations. That way, every regression unit can be associated with
multiple local regions in the input space. This structure is realized by
endowing every upper-level component $m$ with its own local DP. These lower-
level DPs, with a concentration parameter $\alpha$, generate the weights
$\bm{\mathbf{\pi}}_{m}$, the activation centers
$\\{\bm{\mathbf{\mu}}_{mk}\\}_{k=1}^{\infty}$ and the shift coefficients
$\\{\bm{\mathbf{c}}_{mk}\\}_{k=1}^{\infty}$

$\displaystyle\bm{\mathbf{\pi}}_{m}(\bm{\mathbf{s}}_{m})\sim\operatorname{GEM}(\alpha),$
$\displaystyle\bm{\mathbf{\mu}}_{mk}\sim\operatorname{N}(\bm{\mathbf{\tau}}_{m},\kappa\bm{\mathbf{\Lambda}}_{m}),\,\,\bm{\mathbf{c}}_{mk}\sim\operatorname{N}(\bm{\mathbf{\theta}},\rho\bm{\mathbf{V}}_{m}).$

We assume here that the lower-level centers $\bm{\mathbf{\mu}}_{mk}$ and
coefficients $\bm{\mathbf{c}}_{mk}$ are tied via their respective upper-level
precision matrices $\bm{\mathbf{\Lambda}}_{m}$ and $\bm{\mathbf{V}}_{m}$.

Finally, given all necessary parameters, the upper- and lower-level labels
$\bm{\mathbf{h}}_{n},\bm{\mathbf{z}}_{n}$ and the data pairs
$\bm{\mathbf{x}}_{n},\bm{\mathbf{y}}_{n}$ are sampled according to the
likelihoods

$\displaystyle\bm{\mathbf{h}}_{n}\sim\operatorname{Cat}(\bm{\mathbf{\omega}}(\bm{\mathbf{t}})),\,\,\bm{\mathbf{z}}_{n}\sim\operatorname{Cat}(\bm{\mathbf{\pi}}(\bm{\mathbf{s}}),\bm{\mathbf{h}}_{n}),$
$\displaystyle\bm{\mathbf{x}}_{n}\sim\operatorname{N}(\bm{\mathbf{\mu}}_{h_{n},z_{n}},\bm{\mathbf{\Lambda}}_{h_{n}}),\,\,\bm{\mathbf{y}}_{n}\sim\operatorname{N}(\bm{\mathbf{A}}_{h_{n}}\bm{\mathbf{x}}_{n}+\bm{\mathbf{c}}_{h_{n},z_{n}},\bm{\mathbf{V}}_{h_{n}}).$

Notice how the upper-level generative process resembles the structure of ILR.
However, in this model, the meta-activations
$(\bm{\mathbf{\tau}}_{m},\bm{\mathbf{\Lambda}}_{m})$ induce another hierarchy
over the lower-level multi-modal activation centers $\bm{\mathbf{\mu}}_{mk}$.
Moreover, the upper-level DP only accounts for the shared slope of the local
regression units $\bm{\mathbf{A}}_{m}$, whereas the shift coefficients
$\bm{\mathbf{c}}_{mk}$ are induced at the lower level, as they are influenced
by the individual activation centers $\bm{\mathbf{\mu}}_{mk}$.

Similar to the procedure in Section 3, we aim to infer a posterior over the
parameters of this model given a set of inputs
$\bm{\mathbf{X}}=\\{\bm{\mathbf{x}}_{1},\dots,\bm{\mathbf{x}}_{N}\\}$ and
outputs $\bm{\mathbf{Y}}=\\{\bm{\mathbf{y}}_{1},\dots,\bm{\mathbf{y}}_{N}\\}$.
To that end, we derive a structured variational Bayes expectation-maximization
algorithm for this model. We start the derivation by defining the complete
data likelihood, the infinite prior, and the mean-field posterior
factorization.

Complete Data Likelihood. The likelihood model is a two-level precision-tied
joint density over the observations $(\bm{\mathbf{X}},\bm{\mathbf{Y}})$ and
the one-hot upper- and lower-labels $(\bm{\mathbf{H}},\bm{\mathbf{Z}})$

$\displaystyle p(.)$
$\displaystyle=p(\bm{\mathbf{H}})\,p(\bm{\mathbf{Z}}\,|\,\bm{\mathbf{H}})\,p(\bm{\mathbf{X}}\,|\,\bm{\mathbf{H}},\bm{\mathbf{Z}})\,p(\bm{\mathbf{Y}}\,|\,\bm{\mathbf{H}},\bm{\mathbf{Z}},\bm{\mathbf{X}})$
$\displaystyle=\begin{aligned}
&\prod_{n=1}^{N}\operatorname{Cat}(\bm{\mathbf{h}}_{n}\,|\,\bm{\mathbf{\omega}}(\bm{\mathbf{t}}))\\\
&\times\prod_{n=1}^{N}\operatorname{Cat}(\bm{\mathbf{z}}_{n}\,|\,\bm{\mathbf{\pi}}(\bm{\mathbf{s}}),\bm{\mathbf{h}}_{n})\\\
&\times\prod_{n=1}^{N}\prod_{m=1}^{\infty}\prod_{k=1}^{\infty}\operatorname{N}(\bm{\mathbf{x}}_{n}\,|\,\bm{\mathbf{\mu}}_{mk},\bm{\mathbf{\Lambda}}_{m})^{z_{nk}\times
h_{nm}}\\\
&\times\prod_{n=1}^{N}\prod_{m=1}^{\infty}\prod_{k=1}^{\infty}\operatorname{N}(\bm{\mathbf{y}}_{n}\,|\,\bm{\mathbf{A}}_{m}\bm{\mathbf{x}}_{n}+\bm{\mathbf{c}}_{mk},\bm{\mathbf{V}}_{m})^{z_{nk}\times
h_{nm}},\end{aligned}$

where the conditioning of the lower-level labels $\bm{\mathbf{z}}_{n}$ on the
upper-level labels $\bm{\mathbf{h}}_{n}$ manifests as a one-hot selector of
the corresponding lower-level categoricals

$\displaystyle\operatorname{Cat}(\bm{\mathbf{z}}_{n}\,|\,\bm{\mathbf{\pi}}(\bm{\mathbf{s}}),\bm{\mathbf{h}}_{n})$
$\displaystyle=\prod_{m=1}^{\infty}\operatorname{Cat}(\bm{\mathbf{z}}_{n}\,|\,\bm{\mathbf{\pi}}_{m}(\bm{\mathbf{s}}_{m}))^{h_{nm}}$
$\displaystyle=\prod_{m=1}^{\infty}\prod_{k=1}^{\infty}\pi_{mk}^{z_{nk}\times
h_{nm}}.$

Infinite Conjugate Prior. We assume a factorized two-level precision-tied
conjugate infinite mixture prior

$p(\bm{\mathbf{t}},\bm{\mathbf{s}},\bm{\mathbf{\mu}},\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}},\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})=\begin{aligned}
&p(\bm{\mathbf{t}})\,p(\bm{\mathbf{s}})\\\ &\times
p(\bm{\mathbf{\mu}}\,|\,\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})\,p(\bm{\mathbf{\tau}}\,|\,\bm{\mathbf{\Lambda}})\,p(\bm{\mathbf{\Lambda}})\\\
&\times
p(\bm{\mathbf{A}}\,|\,\bm{\mathbf{V}})\,p(\bm{\mathbf{c}}|\bm{\mathbf{V}})\,p(\bm{\mathbf{V}}).\end{aligned}$

The meta-activation prior is a normal-Wishart distribution over the meta-
centers $\bm{\mathbf{\tau}}_{m}$ and precision matrices
$\bm{\mathbf{\Lambda}}_{m}$

$p(\bm{\mathbf{\tau}}\,|\,\bm{\mathbf{\Lambda}})p(\bm{\mathbf{\Lambda}})=\prod_{m=1}^{\infty}\operatorname{N}(\bm{\mathbf{\tau}}_{m}\,|\,\bm{\mathbf{m}}_{0},\lambda_{0}\bm{\mathbf{\Lambda}}_{m})\operatorname{W}(\bm{\mathbf{\Lambda}}_{m}\,|\,\bm{\mathbf{\Psi}}_{0},\nu_{0}),$

while the activation centers $\bm{\mathbf{\mu}}_{mk}$ are sampled from a
conditional normal distribution

$p(\bm{\mathbf{\mu}}\,|\,\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})=\prod_{m=1}^{\infty}\prod_{k=1}^{\infty}\operatorname{N}(\bm{\mathbf{\mu}}_{mk}\,|\,\bm{\mathbf{\tau}}_{m},\kappa_{0}\bm{\mathbf{\Lambda}}_{m}).$

The mappings $\bm{\mathbf{A}}_{m}$ and precision matrices
$\bm{\mathbf{V}}_{m}$ are sampled form a matrix-normal-Wishart

$p(\bm{\mathbf{A}}\,|\,\bm{\mathbf{V}})p(\bm{\mathbf{V}})=\prod_{m=1}^{\infty}\operatorname{MN}(\bm{\mathbf{A}}_{m}\,|\,\bm{\mathbf{M}}_{0},\bm{\mathbf{K}}_{0},\bm{\mathbf{V}}_{m})\operatorname{W}(\bm{\mathbf{V}}_{m}\,|\,\bm{\mathbf{\Phi}}_{0},\eta_{0}).$

while the biases $\bm{\mathbf{c}}_{mk}$ are drawn from a $K$-tied conditional
normal distribution

$p(\bm{\mathbf{c}}\,|\,\bm{\mathbf{V}})=\prod_{m=1}^{\infty}\prod_{k=1}^{\infty}\operatorname{MN}(\bm{\mathbf{c}}_{mk}\,|\,\bm{\mathbf{\theta}}_{0},\rho_{0}\bm{\mathbf{V}}_{m}).$

Finally, the stick-breaking priors $p(\bm{\mathbf{t}}$, $\bm{\mathbf{s}})$
follow the definitions from Section 3

$\displaystyle p(\bm{\mathbf{t}})$
$\displaystyle=\prod_{m=1}^{\infty}\operatorname{Beta}(t_{m}\,|\,1,\beta_{0}),$
$\displaystyle p(\bm{\mathbf{s}})$
$\displaystyle=\prod_{m=1}^{\infty}\prod_{k=1}^{\infty}\operatorname{Beta}(s_{mk}\,|\,1,\alpha_{0}).$

Truncated Mean-Field Posterior. We assume a structured decomposition of the
posterior that leads to conjugate computation while maintaining the
dependencies between the discrete labels, the input activations, and the
regression parameters, respectively. Moreover, we apply the truncation scheme
from [59] to establish a tractable posterior approximation while maintaining
an infinite prior

$\displaystyle p(.\,|\,\mathcal{D})$ $\displaystyle\approx
q(\bm{\mathbf{H}})\,q(\bm{\mathbf{Z}}\,|\,\bm{\mathbf{H}})\,q(\bm{\mathbf{t}})\,q(\bm{\mathbf{s}})\,q(\bm{\mathbf{\mu}},\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})\,q(\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})$
$\displaystyle=\begin{aligned}
&\prod_{n=1}^{N}\operatorname{Cat}(\bm{\mathbf{h}}_{n}\,|\,\bm{\mathbf{g}}_{n})\prod_{n=1}^{N}\operatorname{Cat}(\bm{\mathbf{z}}_{n}\,|\,\bm{\mathbf{r}}_{n},\bm{\mathbf{h}}_{n})\\\
&\times\prod_{m=1}^{M-1}\operatorname{Beta}(t_{m}\,|\,\delta_{m},\beta_{m})\\\
&\times\prod_{m=1}^{M-1}\prod_{k=1}^{K-1}\operatorname{Beta}(s_{mk}\,|\,\gamma_{mk},\alpha_{mk})\\\
&\times\prod_{m=1}^{M}\operatorname{N}(\bm{\mathbf{\tau}}_{m}\,|\,\bm{\mathbf{m}}_{m},\lambda_{m}\bm{\mathbf{\Lambda}}_{m})\operatorname{W}(\bm{\mathbf{\Lambda}}_{m}\,|\,\bm{\mathbf{\Psi}}_{m},\nu_{m})\\\
&\times\prod_{m=1}^{M}\operatorname{MN}(\bm{\mathbf{A}}_{m}\,|\,\bm{\mathbf{M}}_{m},\bm{\mathbf{K}}_{m},\bm{\mathbf{V}}_{m})\operatorname{W}(\bm{\mathbf{V}}_{m}\,|\,\bm{\mathbf{\Phi}}_{m},\eta_{m})\\\
&\times\prod_{m=1}^{M}\prod_{k=1}^{K}\operatorname{N}(\bm{\mathbf{\mu}}_{mk}\,|\,\bm{\mathbf{\tau}}_{m},\kappa_{mk}\bm{\mathbf{\Lambda}}_{m})\operatorname{N}(\bm{\mathbf{c}}_{mk}\,|\,\bm{\mathbf{\theta}}_{mk},\rho_{mk}\bm{\mathbf{V}}_{k}),\end{aligned}$

where $\bm{\mathbf{g}}_{n}$ and $\bm{\mathbf{r}}_{n}$ are the upper- and
lower-level posterior responsibilities, respectively.

Variational Expectation Step. The E-step computes the joint posterior
categorical over joint labels $\bm{\mathbf{H}}$ and $\bm{\mathbf{Z}}$

$\displaystyle\log q(\bm{\mathbf{Z}}\,|\,\bm{\mathbf{H}})$
$\displaystyle=\begin{aligned} &\mathbb{E}_{q(\bm{\mathbf{s}})}\left[\log
p(\bm{\mathbf{Z}}\,|\,\bm{\mathbf{H}})\right]\\\
&+\mathbb{E}_{q(\bm{\mathbf{\mu}},\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})}\left[\log
p(\bm{\mathbf{X}}\,|\,\bm{\mathbf{H}},\bm{\mathbf{Z}})\right]\\\
&+\mathbb{E}_{q(\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})}\left[\log
p(\bm{\mathbf{Y}}\,|\,\bm{\mathbf{H}},\bm{\mathbf{Z}},\bm{\mathbf{X}})\right]+\operatorname{const}\end{aligned}$
$\displaystyle=\sum_{n=1}^{N}\sum_{m=1}^{M}\sum_{k=1}^{K}h_{nm}z_{nmk}\log
r_{nmk},$ $\displaystyle\log q(\bm{\mathbf{H}})$
$\displaystyle=\mathbb{E}_{q(\bm{\mathbf{t}})}\left[\log
p(\bm{\mathbf{H}})\right]+\log
q(\bm{\mathbf{Z}}\,|\,\bm{\mathbf{H}})+\operatorname{const}$
$\displaystyle=\sum_{n=1}^{N}\sum_{m=1}^{M}h_{nm}\log g_{nm},$

where these expectations can computed in a similar fashion to Section 3. More
details can be found in the appendix.

Variational Maximization Step. The M-step updates the variational gating,
activations, and regression parameters

$\displaystyle\log q(\bm{\mathbf{t}})$
$\displaystyle=\mathbb{E}_{q(\bm{\mathbf{H}})}\left[\log
p(\bm{\mathbf{H}})\right]+\log p(\bm{\mathbf{t}})+\operatorname{const},$
$\displaystyle\log q(\bm{\mathbf{s}})$
$\displaystyle=\mathbb{E}_{q(\bm{\mathbf{H}},\bm{\mathbf{Z}})}\left[\log
p(\bm{\mathbf{Z}}\,|\,\bm{\mathbf{H}})\right]+\log
p(\bm{\mathbf{s}})+\operatorname{const},$ $\displaystyle\log
q(\bm{\mathbf{\mu}})$ $\displaystyle=\begin{aligned}
&\mathbb{E}_{q(\bm{\mathbf{H}},\bm{\mathbf{Z}},\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})}\left[\log
p(\bm{\mathbf{X}}\,|\,\bm{\mathbf{H}},\bm{\mathbf{Z}})\right]\\\
&+\mathbb{E}_{q(\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})}\left[\log
p(\bm{\mathbf{\mu}}\,|\,\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})\right]+\operatorname{const},\end{aligned}$
$\displaystyle\log q(\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})$
$\displaystyle=\begin{aligned}
&\mathbb{E}_{q(\bm{\mathbf{H}},\bm{\mathbf{Z}},\bm{\mathbf{\mu}})}\left[\log
p(\bm{\mathbf{X}}\,|\,\bm{\mathbf{H}},\bm{\mathbf{Z}})\right]+\log
p(\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})\\\
&+\mathbb{E}_{q(\bm{\mathbf{\mu}})}\left[\log
p(\bm{\mathbf{\mu}}\,|\,\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})\right]+\operatorname{const},\end{aligned}$
$\displaystyle\log q(\bm{\mathbf{c}})$ $\displaystyle=\begin{aligned}
&\mathbb{E}_{q(\bm{\mathbf{H}},\bm{\mathbf{Z}},\bm{\mathbf{A}},\bm{\mathbf{V}})}\left[\log
p(\bm{\mathbf{Y}}\,|\,\bm{\mathbf{H}},\bm{\mathbf{Z}},\bm{\mathbf{X}})\right]\\\
&+\mathbb{E}_{q(\bm{\mathbf{V}})}\left[\log
p(\bm{\mathbf{c}}\,|\,\bm{\mathbf{V}})\right]+\operatorname{const},\end{aligned}$
$\displaystyle\log q(\bm{\mathbf{A}},\bm{\mathbf{V}})$
$\displaystyle=\begin{aligned}
&\mathbb{E}_{q(\bm{\mathbf{H}},\bm{\mathbf{Z}},\bm{\mathbf{c}})}\left[\log
p(\bm{\mathbf{Y}}\,|\,\bm{\mathbf{H}},\bm{\mathbf{Z}},\bm{\mathbf{X}})\right]+\log
p(\bm{\mathbf{A}},\bm{\mathbf{V}})\\\
&+\mathbb{E}_{q(\bm{\mathbf{c}})}\left[\log
p(\bm{\mathbf{c}}\,|\,\bm{\mathbf{V}})\right]+\operatorname{const}.\end{aligned}$

As previously stated, these updates resemble posterior computations weighted
by the responsibilities $g_{nm}=\mathbb{E}[h_{nm}]$ and
$r_{nmk}=\mathbb{E}[z_{nmk}]$. The appendix provides further details.

###  4.1 Posterior Predictive Distribution

Prediction with HILR is akin to that with ILR, as described in Section 3.2. We
briefly state the conditional predictive for a component $\bm{\mathbf{h}}=m$
and an activation $\bm{\mathbf{z}}=k$ as

$\displaystyle
p(\hat{\bm{\mathbf{y}}}\,|\,\hat{\bm{\mathbf{x}}},\hat{\bm{\mathbf{h}}}\\!=\\!m,\hat{\bm{\mathbf{z}}}\\!=\\!k,\mathcal{D})$
$\displaystyle=\mathbb{E}_{q(\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})}\left[p(\hat{\bm{\mathbf{y}}}\,|\,\hat{\bm{\mathbf{x}}},\bm{\mathbf{A}}_{m},\bm{\mathbf{c}}_{mk},\bm{\mathbf{V}}_{m})\right]$
$\displaystyle=\operatorname{T}\left(\bm{\mathbf{M}}_{m}\hat{\bm{\mathbf{x}}}+\bm{\mathbf{\theta}}_{mk},a_{mk}\bm{\mathbf{\Phi}}_{m},\eta_{m}+1\right),$

where the computation of $a_{mk}$ is similar to that in Section 3.2. Further,
the weight of the $k-$th activation of the $m-$th component is computed as
follows

$\displaystyle
p(\hat{\bm{\mathbf{x}}},\hat{\bm{\mathbf{h}}}\\!=\\!m,\hat{\bm{\mathbf{z}}}\\!=\\!k\,|\,\mathcal{D})$
$\displaystyle\propto\begin{aligned}
&\mathbb{E}_{q(\bm{\mathbf{t}})}\left[p(\hat{\bm{\mathbf{h}}}\\!=\\!m\,|\,\bm{\mathbf{\omega}}(\bm{\mathbf{t}}))\right]\\\
&\times\mathbb{E}_{q(\bm{\mathbf{s}}_{m})}\left[p(\hat{\bm{\mathbf{z}}}\\!=\\!k\,|\,\bm{\mathbf{\pi}}_{m}(\bm{\mathbf{s}}_{m}))\right]\\\
&\times\mathbb{E}_{q(\bm{\mathbf{\mu}},\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})}\left[p(\hat{\bm{\mathbf{x}}}\,|\,\bm{\mathbf{\mu}}_{mk},\bm{\mathbf{\Lambda}}_{m})\right]\end{aligned}$
$\displaystyle=\begin{aligned}
&\frac{\delta_{m}}{\delta_{m}+\beta_{m}}\prod_{l=1}^{m-1}\left(1-\frac{\delta_{l}}{\delta_{l}+\beta_{l}}\right)\\\
&\times\frac{\gamma_{mk}}{\gamma_{mk}+\alpha_{mk}}\prod_{l=1}^{k-1}\left(1-\frac{\gamma_{ml}}{\gamma_{ml}+\alpha_{ml}}\right)\\\
&\times\operatorname{T}\left(\bm{\mathbf{\mu}}_{mk},\frac{\kappa_{mk}}{1+\kappa_{mk}}\bm{\mathbf{\Psi}}_{m},\nu_{m}+1\right).\end{aligned}$

##  5 Empirical Evaluation

We evaluate the presented models on a range of tasks. Our goals are (1) to
highlight the advantages of ILR and HILR, such as dealing with out-of-
distribution predictions, recovering an input-dependent noise function,
hierarchical gating, sharing parameters, and the ability to perform Bayesian
sequential updates, (2) to benchmark the models on high dimensional datasets
from real robots, and (3) to deploy the models in a real-world scenario to
further empirically demonstrate its validity. A reference implementation can
be found under https://github.com/hanyas/mimo.

Out-of-distribution Uncertainty. In Figure 1 (left), we apply ILR on a
synthetic dataset with two large gaps. We observe how the predictive
uncertainty strongly reflects the lack of training data in these regions and
how the mean prediction falls back to the prior values. This example
highlights the model’s reasonable uncertainty quantification. The out-of-
distribution behavior of ILR is strongly influenced by the discrete gating and
a query’s activation probability that jointly define the overall membership
weights.

Heteroscedastic Noise. We test on two different problems with input-dependent
noise, the cosmic microwave background (CMB) [75], Figure 1 (right), and a
synthetic dataset from a stochastic sinc function
$y(x)=\operatorname{sinc}(x)+\epsilon$, where the noise $\epsilon$ is
distributed according to zero-mean normal with a standard deviation
$\sigma_{\epsilon}(x)=0.05+0.2(1+\sin(2x))/(1+e^{-0.2x})$. In Figure 4 (left),
we see that ILR can approximate the nonlinear functions well. In particular,
the heteroscedastic noise function is recovered in great detail.

Figure 4: Left, a challenging heteroscedastic example of a sinc function
heavily overlayed with input-dependent noise. The figure shows the mean
prediction (red) on the training data (dots) and the true mean function
(dashed black) corrupted by noise (dashed green). The blue dashed lines
represent the recovered noise process. Right, hierarchical local regression
with HILR using parameter sharing in shift-invariant functions. The top figure
shows the _mode-prediction_ (red) along with two standard deviations of
predictive uncertainty (shaded blue). The bottom plot highlights the multi-
modal activation, leading to shared slope information over non-adjacent
regions.

Figure 5: Left and middle, learning discontinuous functions with ILR. The top
figures show the mode-prediction (red) and two standard deviations confidence
(shaded blue). The left example is a simple step function that can be captured
with linear features, while in the middle, we use a polynomial transformation
of the input for more flexibility. Right, we tackle inverse mapping problems
with ILR. This example includes scattered data that maps the input
$\bm{\mathbf{x}}$ to multiple output values $\bm{\mathbf{y}}$. A
discriminative modeling approach fails in these scenarios, as it tries to
capture the ambiguous mean of the function
$f\mathrel{\mathop{\ordinarycolon}}\bm{\mathbf{x}}\rightarrow\bm{\mathbf{y}}$.
By approximating the joint density over both input and output and using mode-
prediction, ILR can reconstruct these non-unique relations via local linear
approximations. The bottom plots show the activation over the input space.

Hierarchical Parameter Sharing. In Figure 4 (right), we test HILR’s ability to
share slope parameters via multi-modal activations. We consider a dataset
stemming from a periodic triangle signal overlayed with additive noise. HILR
decides to activate two upper- and two lower-level regions to match the data
structure, despite having more degrees of freedom at each level.

Discontinuous (polynomial) Functions. In Figure 5 (left), a combination of
step functions is fitted using mode-prediction, as described in Section 3.2.
By using polynomial input features, more expressive local models can be
realized. Figure 5 (middle) depicts an example of cubic regressors, which are
still linear in the parameters, fitted to data sampled from noisy cubic
polynomials.

Inverse Mapping. One crucial advantage of generative over discriminative
modeling is the ability to deal with non-unique inverse mapping problems. Such
scenarios arise when the same input can be mapped to multiple output values.
Joint modeling of the input-output data allows for flexible conditioning and
alleviates the directional graph constraints. In Figure 5 (right), we show a
simple example of how ILR can learn these mappings, adapted from [64].

Figure 6: Bayesian sequential updates. Mean (red) and a two standard
deviations interval (shaded blue) of the predictive distribution fitted to
sequentially arriving data (three batches) from the chirp dataset (gray dots).
For the second and third plots, the posterior fitted to the previous batches
is used as a prior to perform a Bayesian sequential update. There is no
catastrophic forgetting of previously learned knowledge, and in regions with
no data, the prediction falls back to the prior.

Bayesian Sequential Updates. In Figure 6, we demonstrate a sequential learning
problem. Data from a chirp signal arrives in three batches so that the
posterior of previous batches serves as the prior for the next learning phase.
ILR successfully captures the data trend, and no significant catastrophic
forgetting is observed. The approximation errors due to the mean-field
posterior factorization assumption have little influence because the posterior
updates are localized in the input domain.

Robot Inverse Dynamics. Next, we use ILR and HILR to learn the inverse
dynamics of anthropomorphic manipulators. The dynamics of these manipulators
are governed by the following equation

$\bm{\mathbf{u}}=\bm{\mathbf{M}}(\bm{\mathbf{q}})\ddot{\bm{\mathbf{q}}}+\bm{\mathbf{C}}(\bm{\mathbf{q}},\dot{\bm{\mathbf{q}}})+\bm{\mathbf{G}}(\bm{\mathbf{q}})+\epsilon(\bm{\mathbf{q}},\dot{\bm{\mathbf{q}}},\ddot{\bm{\mathbf{q}}}),$

where $(\bm{\mathbf{q}},\dot{\bm{\mathbf{q}}},\ddot{\bm{\mathbf{q}}})$ are
joint angles, velocities and accelerations, and $\bm{\mathbf{u}}$ are torques.
$\bm{\mathbf{M}}(\bm{\mathbf{q}})$ is the inertia matrix,
$\bm{\mathbf{C}}(\bm{\mathbf{q}},\dot{\bm{\mathbf{q}}})$ are the Coriolis and
centripetal forces, and $\bm{\mathbf{G}}(\bm{\mathbf{q}})$ is the gravity
force.
$\epsilon(\bm{\mathbf{q}},\dot{\bm{\mathbf{q}}},\ddot{\bm{\mathbf{q}}})$ are
general unmodelled nonlinearities such as sticktion/friction and hydraulic and
tendon/cable dynamics, that motivate a data-driven approach to learn the
mapping
$(\bm{\mathbf{q}},\dot{\bm{\mathbf{q}}},\ddot{\bm{\mathbf{q}}})\rightarrow\bm{\mathbf{u}}$.
Later, we use the learned ILR model for online inverse dynamics control.

TABLE I: Accuracy on the SARCOS dataset | MSE | NMSE | Experts
---|---|---|---
ILR | $0.480\pm 0.030$ | $0.003\,40\pm 0.0002$ | $1700$
HILR | $0.530\pm 0.040$ | $0.003\,90\pm 0.0003$ | $1450$
LGR∗ | $8.600\pm 0.000$ | $0.050\,00\pm 0.0000$ | $7000$
LWPR | $2.600\pm 0.030$ | $0.018\,00\pm 0.0002$ | $32\,000$
GPR | N/A | N/A | -
SGPR | $0.850\pm 0.003$ | $0.006\,00\pm 0.000\,008$ | -
rBCM | $0.452\pm 0.005$ | $0.002\,60\pm 0.000\,030$ | $315$
gPoE | $0.460\pm 0.006$ | $0.003\,00\pm 0.000\,075$ | $315$

We use the mean squared error (MSE), normalized mean squared error (NMSE), and
the number of local experts as evaluation criteria. These measures cover the
prediction accuracy as well as the complexity of the learned model. We compare
to popular (probabilistic) methods such as local Gaussian regression
(LGR)[42], locally weighted projection regression (LWPR)[32], Gaussian process
regression (GPR)[3] and sparse Gaussian process regression (SGPR)[10], and two
scalable Gaussian process product of experts: the robust Bayesian committee
machine (rBCM) [12] and the generalized product of experts (gPoE) [11].

We benchmark the prediction accuracy of all regression techniques on a high-
dimensional dataset collected from a 7-DoF (degrees of freedom)
anthropomorphic SARCOS arm [32]. The dataset consists of $44484$ training
points and 4449 test cases. Overall there are 21 input variables
$(\bm{\mathbf{q}},\dot{\bm{\mathbf{q}}},\ddot{\bm{\mathbf{q}}})$ mapping to 7
motor torques $\bm{\mathbf{u}}$. We also benchmark on an inverse dynamics
dataset from a 4-DoF Barrett-WAM manipulator, mapping from a 12-D to 4-D
space. This dataset contains $25000$ training and $5000$ test pairs.

Table I and Table II list the results for both datasets. We report the MSE,
NMSE, and the number of active models over all joints, averaged over five
seeds. We compute the mean and standard deviation for every cell in the table,
except for LGR ∗, because of the unreasonable training times achieved while
using the authors’ code. When evaluating GPR on the SARCOS dataset, we faced
GPU-memory constraints $($32\text{\,}\mathrm{GB}$)$, and we have discarded
this evaluation. For rBCM and gPoE, we assigned an expert to every 1000 data
points, repeated for every output dimension.

The results show that ILR and HILR outperform the related local regression
methods LWPR and LGR, both in terms of prediction accuracy and number of
active experts. Nonetheless, GPR is still the gold standard when the kernel
size is within memory limits. Interestingly, the figures also indicate that
ILR and HILR are competitive with sparse Gaussian process regression (SGPR)
and the two product of experts rBCM and gPoE. Finally, the evidence reveals
that HILR tends to activate roughly $10-15\%$ fewer components than ILR. This
observation hints that HILR may be taking advantage of shift-invariance
patterns in the data and avoiding duplicate regression units. However, this
hypothesis is hard to validate due to the data’s high dimensionality.

TABLE II: Accuracy on the Barrett-WAM dataset | MSE | NMSE | Experts
---|---|---|---
ILR | $0.290\pm 0.050$ | $0.0070\pm 0.0005$ | $1350$
HILR | $0.310\pm 0.065$ | $0.0080\pm 0.0006$ | $1110$
LGR∗ | $0.770\pm 0.000$ | $0.0170\pm 0.0000$ | $3270$
LWPR | $1.000\pm 0.150$ | $0.0370\pm 0.0100$ | $2900$
GPR | $0.100\pm 0.001$ | $0.0023\pm 0.000\,01$ | -
SGPR | $0.180\pm 0.005$ | $0.0063\pm 0.000\,02$ | -
rBCM | $0.380\pm 0.035$ | $0.0190\pm 0.001\,80$ | $100$
gPoE | $0.340\pm 0.013$ | $0.0160\pm 0.000\,60$ | $100$

Real Inverse Dynamics Control. Finally, we demonstrate the validity of the
learned dynamics in a real-robot control scenario. In this evaluation, we
tackle two experiments. In the first, we focus on the Bayesian sequential
learning aspect, whereas the second highlights our model’s ability to deal
with massive amounts of data. Both experiments deploy an ILR model on the
Barrett-WAM to track an 8-shaped trajectory in the $xy$-plane of the end-
effector.

For the first experiment, we collect $30000$ training samples (roughly $1$
minute) from multiple trajectories with different velocity profiles and
Bayesian sequential learning over $15$ batches for multiple seeds. Figure 7
depicts the progression of the learning process. We then select the best model
and perform model-based control to track held-out test trajectories with
unseen velocity profiles. ILR predicts the feed-forward torques while
supported by a low-gain PD-controller for safety considerations. We compare
this controller to one with access to the analytical dynamics and the same
low-gain PD-controller and to a _model-free_ PD-controller. Figure 8 depicts
the tracking performance of the different controllers on two test
trajectories.

In the second experiment, we learn a model covering a larger region of the
state-action space. We generate a larger Barrett-WAM dataset consisting of
$150000$ training examples (roughly 5 minutes). We test the model learned by
ILR on held-out test trajectories and compare it to the same controller
constellation as in the first experiment. As benchmarking criteria, we
evaluate the MSE with respect to the desired trajectory and the mean torque
contribution of the PD-controller to the overall control signal. The rationale
is that a good inverse dynamics model will consistently produce a low MSE
without relying on the PD-controller’s assistance. Table III shows the
benchmarks on three test trajectories. The results indicate that ILR
significantly improves tracking performance with minimal contribution from the
PD-controller. Moreover, we are able to consistently achieve a prediction
frequency of $2000\text{\,}\mathrm{Hz}$ during both tasks, although the
Barrett-WAM robot requires only $500\text{\,}\mathrm{Hz}$.

TABLE III: Tracking error and torque contributed by the PD-controller during the Barrett-WAM real robot task. |  | PD | Analytic+PD | ILR+PD
---|---|---|---|---
T1 | MSE | $2.33\times 10^{-2}$ | $2.16\times 10^{-2}$ | $1.03\times 10^{-3}$
PD-Torque | $8.25\times 10^{0}$ | $7.12\times 10^{0}$ | $1.40\times 10^{0}$
T2 | MSE | $2.60\times 10^{-2}$ | $2.55\times 10^{-2}$ | $9.17\times 10^{-4}$
PD-Torque | $8.71\times 10^{0}$ | $7.41\times 10^{0}$ | $1.33\times 10^{0}$
T3 | MSE | $2.94\times 10^{-2}$ | $3.08\times 10^{-2}$ | $8.96\times 10^{-4}$
PD-Torque | $9.38\times 10^{0}$ | $8.06\times 10^{0}$ | $1.38\times 10^{0}$
Figure 7: 8-Shaped trajectory learning. Bayesian sequential updates on a
dataset collected from a Barrett-WAM. We plot the NMSE for five different
seeds on accumulated data over the number of batches. The NMSE consistently
improves with new data, and no catastrophic forgetting is observed.

Figure 8: 8-shaped trajectory tracking on the Barrett-WAM. We compare three
controllers on two test trajectories (blue), a low-gain PD (black), a low-gain
PD + feed-forward torques from an analytical model (red), and a low-gain PD +
feed-forward torques from ILR (green). The results indicate that ILR delivers
the best tracking performance.

##  6 Conclusion

In this paper, we presented two probabilistic hierarchical local regression
models, ILR and HILR, and derived an efficient variational inference technique
for data-driven learning. These representations are based on infinite mixtures
from Bayesian nonparametrics. We situate our contributions as the next
iteration in a large family of local linear regression techniques such as
RFWR, LWPR, and LGR. We have shown that placing Dirichlet process priors on
Bayesian mixtures of local regression units can regularize model complexity
with minor loss in performance and without relying on heuristics. Moreover, we
have highlighted the advantage of the two-level architecture of HILR in
compressing the learned models and reducing the number of active experts.
Empirical evaluation indicates that the models outperform LWPR and LGR, and
are competitive with sparse GPR and products of experts. Finally, we have
empirically confirmed the practicality of this approach for online inverse
dynamics control on a Barrett-WAM robot.

Nonetheless, these presented concepts still suffer from multiple drawbacks.
The mean-field assumption is a source of significant errors in posterior
inference. Collapsed formulations of Dirichlet process priors promise better
approximations [76]. In addition, Bayesian mixture models are influenced by a
large number of hyperparameters, which cannot be directly optimized via
empirical Bayes [77], leading to lower predictive performance when compared to
optimized GPR. Nonetheless, the ELBO offers an objective for hyperparameter
optimization. Naive gradient-based techniques have proven too brittle due to
their reliance on Euclidean distance metrics. A natural-gradient approach
appears to be a suitable alternative in the future.

Further development of hierarchical local regression may focus on treating ILR
and HILR as layers in a multi-layered representation. This extension would
allow the models to benefit from intermediate nonlinear projections into high
dimensional spaces that have proven powerful in deep neural networks. Another
practical consideration is to incorporate physical inductive biases such as
inverse dynamics [78] to facilitate learning meaningful quantities.

## Acknowledgments

This work has received funding from the European Unions Horizon 2020 research
and innovation program under grant agreement No. 640554 (SKILLS4ROBOTS).

## Appendix A

In the following, we present more detailed derivations of the E-step and
M-step of ILR and HILR.

### E-Step of Infinite Linear Regression

$\displaystyle\log q(\bm{\mathbf{Z}})$ $\displaystyle=\begin{aligned}
&\mathbb{E}_{q(\bm{\mathbf{s}})}\left[\log
p(\bm{\mathbf{Z}}|\bm{\mathbf{\pi}}(\bm{\mathbf{s}}))\right]\\\
&+\mathbb{E}_{q(\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})}\left[\log
p(\bm{\mathbf{X}}|\bm{\mathbf{Z}})\right]\\\
&+\mathbb{E}_{q(\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})}\left[\log
p(\bm{\mathbf{Y}}|\bm{\mathbf{Z}},\bm{\mathbf{X}})\right]+\operatorname{const}\end{aligned}$
$\displaystyle=\begin{aligned}
&\operatorname{const}+\sum_{n=1}^{N}\Big{[}\mathbb{E}_{q(\bm{\mathbf{s}})}\left[\log\operatorname{Cat}(\bm{\mathbf{z}}_{n}|\bm{\mathbf{\pi}}(\bm{\mathbf{s}}))\right]\\\
&\\!+\\!\sum_{k=1}^{K}z_{nk}\mathbb{E}_{q(\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})}\left[\log\operatorname{N}(\bm{\mathbf{x}}_{n}|\bm{\mathbf{\mu}}_{k},\bm{\mathbf{\Lambda}}_{k})\right]\\\
&\\!+\\!\sum_{k=1}^{K}z_{nk}\mathbb{E}_{q(\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})}\left[\log\operatorname{N}(\bm{\mathbf{y}}_{n}|\bm{\mathbf{A}}_{k}\bm{\mathbf{x}}_{n}+\bm{\mathbf{c}}_{k},\bm{\mathbf{V}}_{k})\right]\Big{]}\end{aligned}$
$\displaystyle=\sum_{n=1}^{N}\sum_{k=1}^{K}z_{nk}\log r_{nk},$

where

$\displaystyle\mathbb{E}_{q(\bm{\mathbf{s}})}\left[\log\operatorname{Cat}(\bm{\mathbf{z}}_{n}|\bm{\mathbf{\pi}}(\bm{\mathbf{s}}))\right]$
$\displaystyle\\!=\\!\begin{aligned}
&\sum_{k=1}^{K}z_{nk}\mathbb{E}_{q(s_{k})}\left[\log s_{k}\right]\\\
&+\sum_{k=1}^{K}z_{nk}\mathbb{E}_{q(s_{k})}\left[\sum_{l=1}^{k-1}\log(1-s_{l})\right].\end{aligned}$

The individual expectations of the _log-sticks_ under the beta posteriors are
given by

$\displaystyle\mathbb{E}_{q(s_{k})}\left[\log s_{k}\right]$
$\displaystyle=\Psi\left(\gamma_{k}\right)-\Psi\left(\gamma_{k}+\alpha_{k}\right),$
$\displaystyle\mathbb{E}_{q(s_{k})}\left[\log(1-s_{k})\right]$
$\displaystyle=\Psi\left(\alpha_{k}\right)-\Psi\left(\gamma_{k}+\alpha_{k}\right),$

where $\Psi$ is the Digamma function.

### M-Step of Infinite Linear Regression

$\displaystyle\log q(\bm{\mathbf{s}})$
$\displaystyle=\mathbb{E}_{q(\bm{\mathbf{Z}})}\left[\log
p(\bm{\mathbf{Z}}|\bm{\mathbf{\pi}}(\bm{\mathbf{s}}))\right]+\log
p(\bm{\mathbf{s}})+\operatorname{const}$ $\displaystyle=\begin{aligned}
&\sum_{n=1}^{N}\sum_{k=1}^{K}r_{nk}\log\left[s_{k}\prod_{l=1}^{k-1}(1-s_{l})\right]\\\
&+\sum_{k=1}^{K-1}\log\operatorname{Beta}(s_{k}|1,\alpha_{0})+\operatorname{const},\end{aligned}$
$\displaystyle\log q(\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})$
$\displaystyle=\mathbb{E}_{q(\bm{\mathbf{Z}})}\left[\log
p(\bm{\mathbf{X}}|\bm{\mathbf{Z}})\right]+\log
p(\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})+\operatorname{const}$
$\displaystyle=\begin{aligned}
&\sum_{n=1}^{N}\sum_{k=1}^{K}r_{nk}\log\operatorname{N}(\bm{\mathbf{x}}_{n}|\bm{\mathbf{\mu}}_{k},\bm{\mathbf{\Lambda}}_{k})\\\
&+\sum_{k=1}^{K}\log\operatorname{N}(\bm{\mathbf{\mu}}_{k}|\bm{\mathbf{m}}_{0},\kappa_{0}\bm{\mathbf{\Lambda}}_{k})\\\
&+\sum_{k=1}^{K}\log\operatorname{W}(\bm{\mathbf{\Lambda}}_{k}|\bm{\mathbf{\Psi}}_{0},\nu_{0})+\operatorname{const},\end{aligned}$
$\displaystyle\log q(\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})$
$\displaystyle=\begin{aligned} &\mathbb{E}_{q(\bm{\mathbf{Z}})}\left[\log
p(\bm{\mathbf{Y}}|\bm{\mathbf{Z}},\bm{\mathbf{X}})\right]\\\ &+\log
p(\bm{\mathbf{A}},\bm{\mathbf{c}},\bm{\mathbf{V}})+\operatorname{const}\end{aligned}$
$\displaystyle=\begin{aligned}
&\sum_{n=1}^{N}\sum_{k=1}^{K}r_{nk}\log\operatorname{N}(\bm{\mathbf{y}}_{n}|\bm{\mathbf{A}}_{k}\bm{\mathbf{x}}_{n}+\bm{\mathbf{c}}_{k},\bm{\mathbf{V}}_{k})\\\
&+\sum_{k=1}^{K}\log\operatorname{MN}(\bm{\mathbf{A}}_{k}|\bm{\mathbf{M}}_{0},\bm{\mathbf{K}}_{0},\bm{\mathbf{V}}_{k})\\\
&+\sum_{k=1}^{K}\log\operatorname{N}(\bm{\mathbf{c}}_{k}|\bm{\mathbf{\theta}}_{0},\rho_{0}\bm{\mathbf{V}}_{k})\\\
&+\sum_{k=1}^{K}\log\operatorname{W}(\bm{\mathbf{V}}_{k}|\bm{\mathbf{\Phi}}_{0},\eta_{0})+\operatorname{const},\end{aligned}$

where $r_{nk}=\mathbb{E}_{q(\bm{\mathbf{Z}})}\left[z_{nk}\right]$ are the
expected responsibilities computed in the E-step. Further details on
performing these updates via general exponential family recipes are provided
in the upcoming sections.

### M-Step of Hierarchical Infinite Linear Regression

$\displaystyle\log q(\bm{\mathbf{t}})\begin{aligned}
&=\mathbb{E}_{q(\bm{\mathbf{H}})}\left[\log
p(\bm{\mathbf{H}}|\bm{\mathbf{\omega}}(\bm{\mathbf{t}}))\right]+\log
p(\bm{\mathbf{t}})+\operatorname{const}\\\ &=\begin{aligned}
&\sum_{n=1}^{N}\sum_{m=1}^{M}\hat{g}\log\left[t_{m}\prod_{l=1}^{m-1}(1-t_{l})\right]\\\
&+\sum_{m=1}^{M}\log\operatorname{Beta}(t_{m}|1,\beta_{0})+\operatorname{const},\end{aligned}\end{aligned}$
$\displaystyle\log q(\bm{\mathbf{s}})\begin{aligned}
&=\mathbb{E}_{q(\bm{\mathbf{H}},\bm{\mathbf{Z}})}\left[\log
p(\bm{\mathbf{Z}}|\bm{\mathbf{H}})\right]+\log
p(\bm{\mathbf{s}})+\operatorname{const}\\\ &=\begin{aligned}
&\sum_{n=1}^{N}\sum_{m=1}^{M}\sum_{k=1}^{K}\hat{g}\,\hat{r}\log\left[s_{mk}\prod_{l=1}^{k-1}(1-s_{ml})\right]\\\
&+\sum_{m=1}^{M}\sum_{k=1}^{K}\log\operatorname{Beta}(s_{mk}|1,\alpha_{0})+\operatorname{const},\end{aligned}\end{aligned}$
$\displaystyle\log q(\bm{\mathbf{\mu}})\begin{aligned} &=\begin{aligned}
&\mathbb{E}_{q(\bm{\mathbf{H}},\bm{\mathbf{Z}},\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})}\left[\log
p(\bm{\mathbf{X}}|\bm{\mathbf{H}},\bm{\mathbf{Z}})\right]\\\
&+\mathbb{E}_{q(\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})}\left[\log
p(\bm{\mathbf{\mu}}|\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})\right]+\operatorname{const}\end{aligned}\\\
&=\begin{aligned}
&\sum_{n=1}^{N}\sum_{m=1}^{M}\sum_{k=1}^{K}\hat{g}\,\hat{r}\,\mathbb{E}_{q(\bm{\mathbf{\Lambda}})}\left[\log\operatorname{N}(\bm{\mathbf{x}}_{n}|\bm{\mathbf{\mu}}_{mk},\bm{\mathbf{\Lambda}}_{m})\right]\\\
&+\sum_{m=1}^{M}\sum_{k=1}^{K}\mathbb{E}_{q(\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})}\left[\log\operatorname{N}(\bm{\mathbf{\mu}}_{mk}|\bm{\mathbf{\tau}}_{m},\kappa_{0}\bm{\mathbf{\Lambda}}_{m})\right]\\\
&+\operatorname{const},\end{aligned}\end{aligned}$ $\displaystyle\log
q(\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})\begin{aligned} &=\begin{aligned}
&\mathbb{E}_{q(\bm{\mathbf{H}},\bm{\mathbf{Z}},\bm{\mathbf{\mu}})}\left[\log
p(\bm{\mathbf{X}}|\bm{\mathbf{H}},\bm{\mathbf{Z}})\right]+\log
p(\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})\\\
&+\mathbb{E}_{q(\bm{\mathbf{\mu}})}\left[\log
p(\bm{\mathbf{\mu}}|\bm{\mathbf{\tau}},\bm{\mathbf{\Lambda}})\right]+\operatorname{const}\end{aligned}\\\
&=\begin{aligned}
&\sum_{n=1}^{N}\sum_{m=1}^{M}\sum_{k=1}^{K}\hat{g}\,\hat{r}\,\mathbb{E}_{q(\bm{\mathbf{\mu}})}\left[\log\operatorname{N}(\bm{\mathbf{x}}_{n}|\bm{\mathbf{\mu}}_{mk},\bm{\mathbf{\Lambda}}_{m})\right]\\\
&+\sum_{m=1}^{M}\log\operatorname{N}(\bm{\mathbf{\tau}}_{m}|\bm{\mathbf{m}}_{0},\lambda_{0}\bm{\mathbf{\Lambda}}_{m})\\\
&+\sum_{m=1}^{M}\log\operatorname{W}(\bm{\mathbf{\Lambda}}_{m}|\bm{\mathbf{\Psi}}_{0},\nu_{0})\\\
&+\sum_{m=1}^{M}\sum_{k=1}^{K}\mathbb{E}_{q(\bm{\mathbf{\mu}})}\left[\log\operatorname{N}(\bm{\mathbf{\mu}}_{mk}|\bm{\mathbf{\tau}}_{m},\kappa_{0}\bm{\mathbf{\Lambda}}_{m})\right]\\\
&+\operatorname{const},\end{aligned}\end{aligned}$ $\displaystyle\log
q(\bm{\mathbf{c}})\begin{aligned} &=\begin{aligned}
&\mathbb{E}_{q(\bm{\mathbf{H}},\bm{\mathbf{Z}},\bm{\mathbf{A}},\bm{\mathbf{V}})}\left[\log
p(\bm{\mathbf{Y}}|\bm{\mathbf{H}},\bm{\mathbf{Z}},\bm{\mathbf{X}})\right]\\\
&+\mathbb{E}_{q(\bm{\mathbf{V}})}\left[\log
p(\bm{\mathbf{c}}|\bm{\mathbf{V}})\right]+\operatorname{const}\end{aligned}\\\
&=\begin{aligned}
&\sum_{n=1}^{N}\sum_{m=1}^{M}\sum_{k=1}^{K}\hat{g}\,\hat{r}\,\\\
&\times\mathbb{E}_{q(\bm{\mathbf{A}},\bm{\mathbf{V}})}\left[\log\operatorname{N}(\bm{\mathbf{y}}_{n}|\bm{\mathbf{A}}_{m}\bm{\mathbf{x}}_{n}+\bm{\mathbf{c}}_{mk},\bm{\mathbf{V}}_{m})\right]\\\
&+\sum_{m=1}^{M}\sum_{k=1}^{K}\mathbb{E}_{q(\bm{\mathbf{V}})}\left[\operatorname{MN}(\bm{\mathbf{c}}_{mk}|\bm{\mathbf{\theta}}_{0},\rho_{0}\bm{\mathbf{V}}_{m})\right]\\\
&+\operatorname{const},\end{aligned}\end{aligned}$ $\displaystyle\log
q(\bm{\mathbf{A}},\bm{\mathbf{V}})\begin{aligned} &=\begin{aligned}
&\mathbb{E}_{q(\bm{\mathbf{H}},\bm{\mathbf{Z}},\bm{\mathbf{c}})}\left[\log
p(\bm{\mathbf{Y}}|\bm{\mathbf{H}},\bm{\mathbf{Z}},\bm{\mathbf{X}})\right]+\log
p(\bm{\mathbf{A}},\bm{\mathbf{V}})\\\
&+\mathbb{E}_{q(\bm{\mathbf{c}})}\left[\log
p(\bm{\mathbf{c}}|\bm{\mathbf{V}})\right]+\operatorname{const}\end{aligned}\\\
&=\begin{aligned}
&\sum_{n=1}^{N}\sum_{m=1}^{M}\sum_{k=1}^{K}\hat{g}\,\hat{r}\,\\\
&\times\mathbb{E}_{q(\bm{\mathbf{c}})}\left[\log\operatorname{N}(\bm{\mathbf{y}}_{n}|\bm{\mathbf{A}}_{m}\bm{\mathbf{x}}_{n}+\bm{\mathbf{c}}_{mk},\bm{\mathbf{V}}_{m})\right]\\\
&+\sum_{m=1}^{M}\sum_{k=1}^{K}\mathbb{E}_{q(\bm{\mathbf{c}})}\left[\operatorname{MN}(\bm{\mathbf{c}}_{mk}|\bm{\mathbf{\theta}}_{0},\rho_{0}\bm{\mathbf{V}}_{m})\right]\\\
&+\operatorname{const},\end{aligned}\end{aligned}$

where the quantities
$\hat{g}\,=\,g_{nm}\,=\,\mathbb{E}_{q(\bm{\mathbf{H}})}\left[h_{nm}\right]$
and $\hat{r}=r_{nmk}=\mathbb{E}_{q(\bm{\mathbf{Z}})}\left[z_{nmk}\right]$.
After computing the necessary expectations, these updates correspond to the
posterior update recipes described in later appendix sections.

## Appendix B

Our work mainly considers random variables with probability density functions
belonging to the exponential family. The unified minimal parameterization of
this class of distributions lends itself to convenient and efficient posterior
computation when paired with conjugate priors.

We assume the natural form for a probability density of a random variable
$\bm{\mathbf{x}}$

$f(\bm{\mathbf{x}}|\bm{\mathbf{\eta}})=h(\bm{\mathbf{x}})\exp\left[\bm{\mathbf{\eta}}\cdot\bm{\mathbf{t}}(\bm{\mathbf{x}})-a(\bm{\mathbf{\eta}})\right],$

where $h(\bm{\mathbf{x}})$ is the base measure, $\bm{\mathbf{\eta}}$ are the
natural parameters, $\bm{\mathbf{t}}(\bm{\mathbf{x}})$ are the sufficient
statistics and $a(\bm{\mathbf{\eta}})$ is the log-partition function, or log-
normalizer. Following the same notation, a conjugate prior
$g(\bm{\mathbf{\eta}}|\bm{\mathbf{\lambda}})$ to the likelihood
$f(\bm{\mathbf{x}}|\bm{\mathbf{\eta}})$ has the form

$g(\bm{\mathbf{\eta}}|\bm{\mathbf{\lambda}})=h(\bm{\mathbf{\eta}})\exp\left[\bm{\mathbf{\lambda}}\cdot\bm{\mathbf{t}}(\bm{\mathbf{\eta}})-a(\bm{\mathbf{\lambda}})\right],$

with prior sufficient statistics
$\bm{\mathbf{t}}(\bm{\mathbf{\eta}})=[\bm{\mathbf{\eta}},\,-a(\bm{\mathbf{\eta}})]^{\top}$
and hyperparameters
$\bm{\mathbf{\lambda}}=[\bm{\mathbf{\alpha}},\,\bm{\mathbf{\beta}}]^{\top}$.
By applying Bayes’ rule, we can directly infer the posterior
$q(\bm{\mathbf{\eta}}|\bm{\mathbf{x}})$

$\displaystyle q(\bm{\mathbf{\eta}}|\bm{\mathbf{x}})$ $\displaystyle\propto
f(\bm{\mathbf{x}}|\bm{\mathbf{\eta}})g(\bm{\mathbf{\eta}}|\bm{\mathbf{\lambda}})$
$\displaystyle\propto\exp\left[\bm{\mathbf{\rho}}(\bm{\mathbf{x}},\bm{\mathbf{\lambda}})\cdot\bm{\mathbf{t}}(\bm{\mathbf{\eta}})-a(\bm{\mathbf{\rho}})\right],$

where the posterior natural parameters
$\bm{\mathbf{\rho}}(\bm{\mathbf{x}},\bm{\mathbf{\lambda}})$ are a function of
the likelihood sufficient statistics $\bm{\mathbf{t}}(\bm{\mathbf{x}})$ and
prior hyperparameters $[\bm{\mathbf{\alpha}},\,\bm{\mathbf{\beta}}]$

$\bm{\mathbf{\rho}}(\bm{\mathbf{x}},\bm{\mathbf{\lambda}})=\left[\bm{\mathbf{\alpha}}+\bm{\mathbf{t}}(\bm{\mathbf{x}}),~{}\bm{\mathbf{\beta}}+\bm{\mathbf{1}}\right]^{\top}.$

The structure of the resulting posterior reveals a simple recipe for data-
driven inference. By moving into the natural space, the posterior parameters
are computed by combining the prior hyperparameters with the likelihood’s
sufficient statistics and log-partition function. By definition, every
exponential family distribution has a minimal natural parameterization that
leads to a unique decomposition of these quantities [79].

## Appendix C

We present an outline of all M-step updates. We use an adapted form of the
exponential natural parameterization, as it offers a clear methodology on how
to derive and implement such updates for all relevant distributions.

### Gaussian with Normal-Wishart Prior

A weighted Gaussian likelihood over a random variable
$\bm{\mathbf{x}}\in\mathds{R}^{d}$ has the following precision-based
parameterization

$\displaystyle p(\bm{\mathbf{X}}|\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})$
$\displaystyle=\prod_{n=1}^{N}\operatorname{N}(\bm{\mathbf{x}}_{n}|\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})^{w_{n}}$
$\displaystyle\propto\exp\left\\{\begin{bmatrix}\bm{\mathbf{\Lambda}}\bm{\mathbf{\mu}}\\\\[5.0pt]
\bm{\mathbf{\mu}}^{\top}\bm{\mathbf{\Lambda}}\bm{\mathbf{\mu}}\\\\[5.0pt]
\bm{\mathbf{\Lambda}}\\\\[5.0pt]
\log|\bm{\mathbf{\Lambda}}|\end{bmatrix}\\!\cdot\\!\begin{bmatrix}\sum_{n=1}^{N}w_{n}\bm{\mathbf{x}}_{n}\\\\[5.0pt]
-\frac{1}{2}\sum_{n=1}^{N}w_{n}\\\\[5.0pt]
-\frac{1}{2}\sum_{n=1}^{N}w_{n}\bm{\mathbf{x}}_{n}\bm{\mathbf{x}}_{n}^{\top}\\\\[5.0pt]
\frac{1}{2}\sum_{n=1}^{N}w_{n}\end{bmatrix}\right\\},$

where $w_{n}$ are the importance weights. The conjugate prior
$p(\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})$ is a normal-Wishart distribution
with zero mean

$\displaystyle p(\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})$
$\displaystyle=\operatorname{N}(\bm{\mathbf{\mu}}|\bm{\mathbf{0}},\kappa_{0}\bm{\mathbf{\Lambda}})\operatorname{W}(\bm{\mathbf{\Lambda}}|\bm{\mathbf{\Psi}}_{0},\nu_{0})$
$\displaystyle\propto\exp\left\\{\begin{bmatrix}\bm{\mathbf{0}}\\\\[5.0pt]
-\frac{1}{2}\kappa_{0}\\\\[5.0pt]
-\frac{1}{2}\bm{\mathbf{\Psi}}_{0}^{-1}\\\\[5.0pt]
\frac{1}{2}(\nu_{0}-d)\end{bmatrix}\\!\cdot\\!\begin{bmatrix}\bm{\mathbf{\Lambda}}\bm{\mathbf{\mu}}\\\\[5.0pt]
\bm{\mathbf{\mu}}^{\top}\bm{\mathbf{\Lambda}}\bm{\mathbf{\mu}}\\\\[5.0pt]
\bm{\mathbf{\Lambda}}\\\\[5.0pt]
\log|\bm{\mathbf{\Lambda}}|\end{bmatrix}\right\\}.$

The resulting posterior $q(\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})$ is also a
normal-Wishart

$\displaystyle q(\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})$
$\displaystyle\\!=\\!\operatorname{N}(\bm{\mathbf{\mu}}|\bm{\mathbf{m}},\kappa\bm{\mathbf{\Lambda}})\operatorname{W}(\bm{\mathbf{\Lambda}}|\bm{\mathbf{\Psi}},\nu)$
$\displaystyle\\!\propto\exp\left\\{\\!\begin{bmatrix}\sum_{n=1}^{N}w_{n}\bm{\mathbf{x}}_{n}\\\\[5.0pt]
-\frac{1}{2}(\kappa_{0}+\sum_{n=1}^{N}w_{n})\\\\[5.0pt]
-\frac{1}{2}(\bm{\mathbf{\Psi}}_{0}^{-1}+\sum_{n=1}^{N}w_{n}\bm{\mathbf{x}}_{n}\bm{\mathbf{x}}_{n}^{\top})\\\\[5.0pt]
\frac{1}{2}(\nu_{0}-d+\sum_{n=1}^{N}w_{n})\end{bmatrix}\\!\cdot\\!\begin{bmatrix}\bm{\mathbf{\Lambda}}\bm{\mathbf{\mu}}\\\\[5.0pt]
\bm{\mathbf{\mu}}^{\top}\bm{\mathbf{\Lambda}}\bm{\mathbf{\mu}}\\\\[5.0pt]
\bm{\mathbf{\Lambda}}\\\\[5.0pt]
\log|\bm{\mathbf{\Lambda}}|\end{bmatrix}\\!\right\\}.$

The resulting posterior parameters are

$\displaystyle\kappa=\kappa_{0}+\sum_{n=1}^{N}w_{n},$
$\displaystyle\bm{\mathbf{m}}=1/\kappa\sum_{n=1}^{N}w_{n}\bm{\mathbf{x}}_{n},$
$\displaystyle\nu=\nu_{0}+\sum_{n=1}^{N}w_{n},$
$\displaystyle\bm{\mathbf{\Psi}}=(\bm{\mathbf{\Psi}}_{0}^{-1}+\sum_{n=1}^{N}w_{n}\bm{\mathbf{x}}_{n}\bm{\mathbf{x}}_{n}^{\top}-\kappa\,\bm{\mathbf{m}}\,\bm{\mathbf{m}}^{\top})^{-1}.$

### Linear-Gaussian with Matrix-Normal-Wishart Prior

A weighted linear-Gaussian likelihood takes a random input variable
$\bm{\mathbf{x}}\in\mathds{R}^{d}$ and returns a response random variable
$\bm{\mathbf{y}}\in\mathds{R}^{m}$ according to a linear mapping
$\bm{\mathbf{A}}\mathrel{\mathop{\ordinarycolon}}\mathds{R}^{d}\to\mathds{R}^{m}$

$\displaystyle
p(\bm{\mathbf{Y}}|\bm{\mathbf{X}},\bm{\mathbf{A}},\bm{\mathbf{V}})$
$\displaystyle=\prod_{n=1}^{N}\operatorname{N}(\bm{\mathbf{y}}_{n}|\bm{\mathbf{x}}_{n},\bm{\mathbf{A}},\bm{\mathbf{V}})^{w_{n}}$
$\displaystyle\propto\exp\left\\{\begin{bmatrix}\bm{\mathbf{V}}\bm{\mathbf{A}}\\\\[5.0pt]
\bm{\mathbf{A}}^{\top}\bm{\mathbf{V}}\bm{\mathbf{A}}\\\\[5.0pt]
\bm{\mathbf{V}}\\\\[5.0pt]
\log|\bm{\mathbf{V}}|\end{bmatrix}\cdot\begin{bmatrix}\bm{\mathbf{Y}}\bm{\mathbf{W}}\bm{\mathbf{X}}^{\top}\\\\[5.0pt]
-\frac{1}{2}\bm{\mathbf{X}}\bm{\mathbf{W}}\bm{\mathbf{X}}^{\top}\\\\[5.0pt]
-\frac{1}{2}\bm{\mathbf{Y}}\bm{\mathbf{W}}\bm{\mathbf{Y}}^{\top}\\\\[5.0pt]
\frac{1}{2}\sum_{n=1}^{N}w_{n}\end{bmatrix}\right\\},$

where $w_{n}$ are the weights and $\bm{\mathbf{W}}=\mathrm{diag}(w_{n})$ is
the diagonal weight matrix. The data matrices $\bm{\mathbf{X}}$ and
$\bm{\mathbf{Y}}$ are of size $d\times N$ and $m\times N$, respectively. The
conjugate prior $p(\bm{\mathbf{A}},\bm{\mathbf{V}})$ is a matrix-normal-
Wishart with zero mean

$\displaystyle p(\bm{\mathbf{A}},\bm{\mathbf{V}})$
$\displaystyle=\operatorname{N}(\bm{\mathbf{A}}|\bm{\mathbf{0}},\bm{\mathbf{V}},\bm{\mathbf{K}}_{0})\operatorname{W}(\bm{\mathbf{V}}|\bm{\mathbf{\Psi}}_{0},\nu_{0})$
$\displaystyle\propto\exp\left\\{\begin{bmatrix}\bm{\mathbf{0}}\\\\[5.0pt]
-\frac{1}{2}\bm{\mathbf{K}}_{0}\\\\[5.0pt]
-\frac{1}{2}\bm{\mathbf{\Psi}}_{0}^{-1}\\\\[5.0pt]
\frac{1}{2}(\nu_{0}-m-1+d)\end{bmatrix}\cdot\begin{bmatrix}\bm{\mathbf{V}}\bm{\mathbf{A}}\\\\[5.0pt]
\bm{\mathbf{A}}^{\top}\bm{\mathbf{V}}\bm{\mathbf{A}}\\\\[5.0pt]
\bm{\mathbf{V}}\\\\[5.0pt] \log|\bm{\mathbf{V}}|\end{bmatrix}\right\\}.$

The posterior $q(\bm{\mathbf{\mu}},\bm{\mathbf{\Lambda}})$ is matrix-normal-
Wishart

$\displaystyle q(\bm{\mathbf{A}},\bm{\mathbf{V}})$
$\displaystyle=\operatorname{N}(\bm{\mathbf{A}}|\bm{\mathbf{M}},\bm{\mathbf{V}},\bm{\mathbf{K}})\operatorname{W}(\bm{\mathbf{V}}|\bm{\mathbf{\Psi}},\nu)$
$\displaystyle\\!\propto\\!\exp\left\\{\\!\begin{bmatrix}\bm{\mathbf{Y}}\bm{\mathbf{W}}\bm{\mathbf{X}}^{\top}\\\\[5.0pt]
-\frac{1}{2}(\bm{\mathbf{K}}_{0}+\bm{\mathbf{X}}\bm{\mathbf{W}}\bm{\mathbf{X}}^{\top})\\\\[5.0pt]
-\frac{1}{2}(\bm{\mathbf{\Psi}}_{0}^{-1}+\bm{\mathbf{Y}}\bm{\mathbf{W}}\bm{\mathbf{Y}}^{\top})\\\\[5.0pt]
\frac{1}{2}(\nu_{0}-m-1+d+\sum_{n=1}^{N}w_{n})\end{bmatrix}\\!\cdot\\!\begin{bmatrix}\bm{\mathbf{V}}\bm{\mathbf{A}}\\\\[5.0pt]
\bm{\mathbf{A}}^{\top}\bm{\mathbf{V}}\bm{\mathbf{A}}\\\\[5.0pt]
\bm{\mathbf{V}}\\\\[5.0pt] \log|\bm{\mathbf{V}}|\end{bmatrix}\\!\right\\}.$

The standard-form posterior parameters are

$\displaystyle\bm{\mathbf{K}}=\bm{\mathbf{K}}_{0}+\bm{\mathbf{X}}\bm{\mathbf{W}}\bm{\mathbf{X}}^{\top},$
$\displaystyle\bm{\mathbf{M}}=\bm{\mathbf{Y}}\bm{\mathbf{W}}\bm{\mathbf{X}}^{\top}\bm{\mathbf{K}}^{-1},$
$\displaystyle\nu=\nu_{0}+\sum_{n=1}^{N}w_{n},$
$\displaystyle\bm{\mathbf{\Psi}}=(\bm{\mathbf{\Psi}}_{0}^{-1}+\bm{\mathbf{Y}}\bm{\mathbf{W}}\bm{\mathbf{Y}}^{\top}-\bm{\mathbf{M}}\,\bm{\mathbf{K}}\,\bm{\mathbf{M}}^{\top})^{-1}.$

### Infinite Categorical with Stick-Breaking Prior

We assume one-hot random variable $\bm{\mathbf{z}}$ of infinite size with a
categorical likelihood

$\displaystyle p(\bm{\mathbf{Z}}|\bm{\mathbf{\pi}}(\bm{\mathbf{s}}))$
$\displaystyle=\prod_{n=1}^{N}\operatorname{Cat}(\bm{\mathbf{z}}_{n}|\bm{\mathbf{\pi}}(\bm{\mathbf{s}}))^{w_{nk}}$
$\displaystyle=\prod_{n=1}^{N}\prod_{k=1}^{\infty}\left[\left(s_{k}\prod_{l=1}^{k-1}\left(1-s_{l}\right)\right)^{z_{nk}}\right]^{w_{nk}},$

where $w_{n}$ are the weights. The prior is an infinitely factorized beta
distribution over the stick lengths $s_{k}$

$\displaystyle p(\bm{\mathbf{s}})$
$\displaystyle=\prod_{k=1}^{\infty}\operatorname{Beta}(s_{k}|\gamma_{0},\alpha_{0})$
$\displaystyle\propto\prod_{k=1}^{\infty}s_{k}^{\gamma_{0}-1}\left(1-s_{k}\right)^{\alpha_{0}-1},$

and the posterior is also a factorized beta density truncated up to a level
$K$ [59]

$q(\bm{\mathbf{s}})=\prod_{k=1}^{K}\operatorname{Beta}(s_{k}|\gamma_{k},\alpha_{k}),$

where

$\displaystyle\gamma_{k}$ $\displaystyle=\gamma_{0}+\sum_{n=1}^{N}w_{nk},$
$\displaystyle\alpha_{k}$
$\displaystyle=\alpha_{0}+\sum_{n=1}^{N}\sum_{l=k+1}^{K}w_{nl}.$

## References

  * [1] E. Todorov, “Stochastic optimal control and estimation methods adapted to the noise characteristics of the sensorimotor system,” _Neural Computation_ , 2005.
  * [2] D. Büchler, R. Calandra, B. Schölkopf, and J. Peters, “Control of musculoskeletal systems using learned dynamics models,” _IEEE Robotics and Automation Letters_ , 2018.
  * [3] C. E. Rasmussen and C. K. I. Williams, _Gaussian Processes for Machine Learning_ , 2006.
  * [4] I. Goodfellow, Y. Bengio, and A. Courville, _Deep Learning_ , 2016.
  * [5] L. Wasserman, _All of Nonparametric Statistics_ , 2006.
  * [6] R. Calandra, J. Peters, C. E. Rasmussen, and M. P. Deisenroth, “Manifold Gaussian processes for regression,” in _International Joint Conference on Neural Networks_ , 2016.
  * [7] A. G. Wilson, Z. Hu, R. Salakhutdinov, and E. P. Xing, “Deep kernel learning,” in _International Conference on Artificial Intelligence and Statistics_ , 2016.
  * [8] H. Salimbeni and M. Deisenroth, “Doubly stochastic variational inference for deep Gaussian processes,” in _Advances in Neural Information Processing Systems_ , 2017.
  * [9] R. Herbrich, N. D. Lawrence, and M. Seeger, “Fast sparse Gaussian process methods: The informative vector machine,” in _Advances in Neural Information Processing Systems_ , 2003.
  * [10] M. Titsias, “Variational learning of inducing variables in sparse Gaussian processes,” in _International Conference on Artificial Intelligence and Statistics_ , 2009.
  * [11] Y. Cao and D. J. Fleet, “Generalized product of experts for automatic and principled fusion of Gaussian process predictions,” _arXiv preprint arXiv:1410.7827_ , 2014.
  * [12] M. Deisenroth and J. W. Ng, “Distributed Gaussian processes,” in _International Conference on Machine Learning_ , 2015.
  * [13] M. Bauer, M. van der Wilk, and C. E. Rasmussen, “Understanding probabilistic sparse Gaussian process approximations,” in _Advances in Neural Information Processing Systems_ , 2016.
  * [14] A. G. d. G. Matthews, “Scalable Gaussian process inference using variational methods,” Ph.D. dissertation, University of Cambridge, 2017.
  * [15] Q. V. Le, A. J. Smola, and S. Canu, “Heteroscedastic Gaussian process regression,” in _International Conference on Machine Learning_ , 2005.
  * [16] K. Kersting, C. Plagemann, P. Pfaff, and W. Burgard, “Most likely heteroscedastic Gaussian process regression,” in _International Conference on Machine Learning_ , 2007.
  * [17] H. Liu, Y.-S. Ong, and J. Cai, “Large-scale heteroscedastic regression via Gaussian process,” _IEEE Transactions on Neural Networks and Learning Systems_ , 2020.
  * [18] J. Frankle and M. Carbin, “The lottery ticket hypothesis: Finding sparse, trainable neural networks,” in _International Conference on Learning Representations_ , 2019.
  * [19] M. McCloskey and N. J. Cohen, “Catastrophic interference in connectionist networks: The sequential learning problem,” in _Psychology of Learning and Motivation_ , 1989.
  * [20] R. Neal, “Bayesian learning for neural networks,” Ph.D. dissertation, University of Toronto, 1994.
  * [21] C. Blundell, J. Cornebise, K. Kavukcuoglu, and D. Wierstra, “Weight uncertainty in neural network,” in _International Conference on Machine Learning_ , 2015.
  * [22] B. Lakshminarayanan, A. Pritzel, and C. Blundell, “Simple and scalable predictive uncertainty estimation using deep ensembles,” in _Advances in Neural Information Processing Systems_ , 2017.
  * [23] M. E. Khan, D. Nielsen, V. Tangkaratt, W. Lin, Y. Gal, and A. Srivastava, “Fast and scalable Bayesian deep learning by weight-perturbation in Adam,” _International Conference on Machine Learning_ , 2018.
  * [24] S. Sun, G. Zhang, J. Shi, and R. Grosse, “Functional variational Bayesian neural networks,” _International Conference on Learning Representations_ , 2019.
  * [25] J. Watson, J. A. Lin, P. Klink, J. Pajarinen, and J. Peters, “Latent derivative Bayesian last layer networks,” in _International Conference on Artificial Intelligence and Statistics_ , 2021.
  * [26] E. Daxberger, A. Kristiadi, A. Immer, R. Eschenhagen, M. Bauer, and P. Hennig, “Laplace redux-effortless Bayesian deep learning,” _Advances in Neural Information Processing Systems_ , 2021.
  * [27] F. Wenzel, K. Roth, B. Veeling, J. Swiatkowski, L. Tran, S. Mandt, J. Snoek, T. Salimans, R. Jenatton, and S. Nowozin, “How good is the Bayes posterior in deep neural networks really?” _International Conference on Machine Learning_ , 2020.
  * [28] A. Y. Foong, D. R. Burt, Y. Li, and R. E. Turner, “On the expressiveness of approximate inference in Bayesian neural networks,” _Advances in Neural Information Processing Systems_ , 2020.
  * [29] C. G. Atkeson, A. W. Moore, and S. Schaal, “Locally weighted learning for control,” _Artificial Intelligence Review_ , 1997.
  * [30] S. Schaal and C. G. Atkeson, “Constructive incremental learning from only local information,” _Neural Computation_ , 1998.
  * [31] S. Schaal, C. G. Atkeson, and S. Vijayakumar, “Scalable techniques from nonparametric statistics for real time robot learning,” _Applied Intelligence_ , 2002.
  * [32] S. Vijayakumar, A. D’Souza, and S. Schaal, “Incremental online learning in high dimensions,” _Neural Computation_ , 2005.
  * [33] J. Moody and C. J. Darken, “Fast learning in networks of locally-tuned processing units,” _Neural Computation_ , 1989.
  * [34] L. Xu, M. Jordan, and G. E. Hinton, “An alternative model for mixtures of experts,” _Advances in Neural Information Processing Systems_ , 1994.
  * [35] O. Nelles and R. Isermann, “Basis function networks for interpolation of local linear models,” _Conference on Decision and Control_ , 1996.
  * [36] P. Moerland, “Classification using localized mixtures of experts,” 1999.
  * [37] R. A. Jacobs, M. I. Jordan, S. J. Nowlan, and G. E. Hinton, “Adaptive mixtures of local experts,” _Neural Computation_ , 1991.
  * [38] M. I. Jordan and R. A. Jacobs, “Hierarchical mixtures of experts and the EM algorithm,” _Neural Computation_ , 1994.
  * [39] C. E. Rasmussen and Z. Ghahramani, “Infinite mixtures of Gaussian process experts,” in _Advances in Neural Information Processing Systems_ , 2002.
  * [40] J. A. Ting, F. Meier, S. Vijayakumar, and S. Schaal, “Locally weighted regression for control,” in _Encyclopedia of Machine Learning and Data Mining_ , 2010.
  * [41] J. A. Ting, M. Kalakrishnan, S. Vijayakumar, and S. Schaal, “Bayesian kernel shaping for learning control,” in _Advances in Neural Information Processing Systems_ , 2009.
  * [42] F. Meier, P. Hennig, and S. Schaal, “Incremental local Gaussian regression,” in _Advances in Neural Information Processing Systems_ , 2014\.
  * [43] H. Abdulsamad, P. Nickl, P. Klink, and J. Peters, “A variational infinite mixture for probabilistic inverse dynamics learning,” _IEEE International Conference on Robotics and Automation_ , 2021.
  * [44] N. L. Hjort, C. Holmes, P. Müller, and S. G. Walker, _Bayesian Nonparametrics_ , 2010.
  * [45] M. J. Beal, “Variational algorithms for approximate Bayesian inference,” Ph.D. dissertation, University College London, 2003.
  * [46] E. Meeds and S. Osindero, “An alternative infinite mixture of Gaussian process experts,” in _Advances in Neural Information Processing Systems_ , 2006.
  * [47] C. Yuan and C. Neubauer, “Variational mixture of Gaussian process experts,” in _Advances in Neural Information Processing Systems_ , 2009.
  * [48] M. D. Escobar and M. West, “Bayesian density estimation and inference using mixtures,” _Journal of the American Statistical Association_ , 1995.
  * [49] R. M. Neal, “Markov chain sampling methods for Dirichlet process mixture models,” _Journal of Computational and Graphical Statistics_ , 2000.
  * [50] H. Ishwaran and L. F. James, “Gibbs sampling methods for stick-breaking priors,” _Journal of the American Statistical Association_ , 2001.
  * [51] C. E. Rasmussen, “The infinite Gaussian mixture model,” in _Advances in Neural Information Processing Systems_ , 1999.
  * [52] D. Blackwell, J. B. MacQueen _et al._ , “Ferguson distributions via pólya urn schemes,” _The Annals of Statistics_ , 1973.
  * [53] T. S. Ferguson, “A Bayesian analysis of some nonparametric problems,” _The Annals of Statistics_ , 1973.
  * [54] P. Mueller, A. Erkanli, and M. West, “Bayesian curve fitting using multivariate normal mixtures,” _Biometrika_ , 1996.
  * [55] B. Shahbaba and R. M. Neal, “Nonlinear models using Dirichlet process mixtures,” _The Journal of Machine Learning Research_ , 2009.
  * [56] L. Hannah, D. M. Blei, and W. B. Powell, “Dirichlet process mixtures of generalized linear models,” _The Journal of Machine Learning Research_ , 2011\.
  * [57] C. Gadd, S. Wade, and A. Boukouvalas, “Enriched mixtures of generalised Gaussian process experts,” in _International Conference on Artificial Intelligence and Statistics_ , 2020.
  * [58] H. Attias, “A variational Bayesian framework for graphical models,” in _Advances in Neural Information Processing Systems_ , 2000.
  * [59] D. M. Blei and M. I. Jordan, “Variational inference for Dirichlet process mixtures,” _Bayesian Analysis_ , 2006.
  * [60] E. Angelino, M. J. Johnson, and R. Adams, “Patterns of scalable Bayesian inference,” _Foundations and Trends in Machine Learning_ , 2016.
  * [61] T. Minka, “Bayesian linear regression,” 2000, accessed: 2020-07-28. [Online]. Available: https://tminka.github.io/papers/linear.html
  * [62] M. J. Beal and Z. Ghahramani, “Variational Bayesian learning of directed graphical models with hidden variables,” _Bayesian Analysis_ , 2006.
  * [63] J. Rousseau and K. Mengersen, “Asymptotic behaviour of the posterior distribution in overfitted mixture models,” _Journal of the Royal Statistical Society_ , 2011.
  * [64] K. P. Murphy, _Machine Learning: A Probabilistic Perspective_ , 2012.
  * [65] Y. W. Teh, “Dirichlet process,” in _Encyclopedia of Machine Learning and Data Mining_ , 2010.
  * [66] J. Sethuraman, “A constructive definition of Dirichlet priors,” _Statistica Sinica_ , 1994.
  * [67] S. Brooks, A. Gelman, G. Jones, and X.-L. Meng, _Handbook of Markov chain Monte Carlo_ , 2011.
  * [68] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe, “Variational inference: A review for statisticians,” _Journal of the American Statistical Association_ , 2017.
  * [69] M. Opper and D. Saad, _Advanced Mean Field Methods: Theory and Practice_ , 2001.
  * [70] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley, “Stochastic variational inference,” _Journal of Machine Learning Research_ , 2013.
  * [71] C. M. Bishop, _Pattern Recognition and Machine Learning_ , 2006.
  * [72] H. Z. Yerebakan, B. Rajwa, and M. Dundar, “The infinite mixture of infinite Gaussian mixtures,” in _Advances in Neural Information Processing Systems_ , 2014.
  * [73] T. V. Nguyen, D. Phung, X. Nguyen, S. Venkatesh, and H. Bui, “Bayesian nonparametric multilevel clustering with group-level contexts,” in _International Conference on Machine Learning_ , 2014.
  * [74] V. Huynh, D. Q. Phung, S. Venkatesh, X. Nguyen, M. D. Hoffman, and H. H. Bui, “Scalable nonparametric bayesian multilevel clustering.” in _Conference on Uncertainty in Artificial Intelligence_ , 2016.
  * [75] C. L. Bennett, M. Halpern, G. Hinshaw, N. Jarosik, A. Kogut, M. Limon, S. S. Meyer, L. Page, D. N. Spergel, G. S. Tucker _et al._ , “First year Wilkinson microwave anisotropy probe (WMAP) observations: Preliminary maps and basic results,” _The Astrophysical Journal Supplement Series_ , 2003\.
  * [76] K. Kurihara, M. Welling, and Y. W. Teh, “Collapsed variational Dirichlet process mixture models.” in _International Joint Conferences on Artificial Intelligence_ , 2007.
  * [77] J. S. Maritz and T. Lwin, “Empirical Bayes methods,” _Monographs on Statistics and Applied Probability_ , 1989.
  * [78] D. Nguyen-Tuong and J. Peters, “Using model knowledge for learning inverse dynamics,” in _IEEE international conference on robotics and automation_ , 2010.
  * [79] M. J. Wainwright and M. I. Jordan, “Graphical models, exponential families, and variational inference,” _Foundations and Trends in Machine Learning_ , 2008.