\thanks{The project this report is based on was supported with funds from the German Federal Ministry of Education and Research under project number 01\,IS\,18049\,A.}

The data association problem is concerned with separating data coming from different generating processes, for example when data comes from different data sources, contain significant noise, or exhibit multimodality.

Underpinning our approach is the use of Gaussian process priors to encode the structure of both the data and the data associations.

We present an efficient learning scheme based on doubly stochastic variational inference and discuss how it can be applied to deep Gaussian process priors.

Real-world data often include multiple operational regimes of the considered system, for example a wind turbine or gas turbine~\parencite{hein_benchmark_2017}.

Building a truthful model of such data requires learning two separate models and correctly associating the observed data to each of the dynamical regimes.

Estimating a model in this scenario is often referred to as a \emph{data association problem}~\parencite{Bar-Shalom:1987, Cox93areview}, where we consider the data to have been generated by a mixture of processes and we are interested in factorising the data into these components.

\Cref{fig:choicenet_data} shows an example of faulty sensor data, where sensor readings are disturbed by uncorrelated and asymmetric noise.

Applying standard machine learning approaches to such data can lead to model pollution, where the expressive power of the model is used to explain noise instead of the underlying signal.

A data association problem consisting of two generating processes, one of which is a signal we wish to recover and one is an uncorrelated noise process.

Early approaches to explaining data using multiple generative processes are based on separating the input space and training local expert models explaining easier subtasks~\parencite{jacobs_adaptive_1991,tresp_mixtures_2001, rasmussen_infinite_2002}.

Another approach is presented in~\parencite{bishop_mixture_1994}, where the multimodal regression tasks are interpreted as a density estimation problem.

In contrast, we are interested in a generative process, where data at the same location in the input space could have been generated by a number of global independent processes.

Inherently, the data association problem is ill-posed and requires assumptions on both the underlying functions and the association of the observations.

In~\parencite{lazaro-gredilla_overlapping_2012} the authors place Gaussian process (GP) priors on the different generative processes which are assumed to be relevant globally.

A drawback is that the model cannot give a posterior estimate about the relevance of the different generating processes at different locations in the input space.

Another approach in~\parencite{bodin_latent_2017} expands this model by allowing interdependencies between the different generative processes and formulating the association problem as an inference problem on a latent space and a corresponding covariance function.

However, in this approach the number of components is a free parameter and is prone to overfitting, as the model has no means of turning off components.

In this paper, we formulate a Bayesian model for the data association problem.

Underpinning our approach is the use of GP priors which encode structure both on the functions and the associations themselves, allowing us to incorporate the available prior knowledge about the proper factorization into the learning problem.

The use of GP priors allows us to achieve principled regularization without reducing the solution space leading to a well-regularized learning problem.

Importantly, we simultaneously solve the association problem for the training data taking both inputs and outputs into account while also obtaining posterior belief about the relevance of the different generating processes in the input space.

Our model can describe non-stationary processes in the sense that a different number of processes can be activated in different locations in the input space.

In \cref{sec:model}, we propose the data association with Gaussian processes model (DAGP).

In \cref{sec:variational_approximation}, we present an efficient learning scheme via a variational approximation which allows us to simultaneously train all parts of our model via stochastic optimization and show how the same learning scheme can be applied to deep GP priors.

We demonstrate our model on a noise separation problem, an artificial multimodal data set, and a multi-regime regression problem based on the cart-pole benchmark in \cref{sec:experiments}.

The assignment $\mat{a_n}$ to a latent function is driven by input-dependent weights $\mat{\alpha_n^{\pix{k}}}$ which encode the relevance of the different functions at $\mat{x_n}$.

The different parts of the model are determined by the hyperparameters $\mat{\theta}, \mat{\sigma}$ (yellow) and variational parameters $\mat{u}$ (blue).

The data association with Gaussian processes (DAGP) model assumes that there exist $K$ independent functions $\Set*{f^{\pix{k}}}_{k=1}^K$, which generate pairs of observations $\D=\Set*{(\mat{x_n}, \mat{y_n})}_{n=1}^N$.

The assignment of the $\nth{n}$ data point to one of the functions is specified by the indicator vector $\mat{a_n}\in\Set*{0, 1}^K$, which has exactly one non-zero entry.

For notational conciseness, we follow the GP related notation in \parencite{hensman_scalable_2015} and collect all $N$ inputs as $\mat{X}=\left(\mat{x_1}, \ldots, \mat{x_N}\right)$ and all outputs as $\mat{Y}=\left(\mat{y_1}, \ldots, \mat{y_N}\right)$.

We further denote the $\nth{k}$ latent function value associated with the $\nth{n}$ data point as $\rv{f_n^{\pix{k}}}=\Fun{f^{\pix{k}}}{\mat{x_n}}$ and collect them as $\mat{F^{\pix{k}}}=\left(\rv{f_1^{\pix{k}}}, \ldots, \rv{f_N^{\pix{k}}}\right)$ and $\mat{F}=\left(\mat{F^{\pix{1}}}, \ldots, \mat{F^{\pix{K}}}\right)$.

We refer to the $\nth{k}$ entry in $\mat{a_n}$ as $a_n^{\pix{k}}$ and denote $\mat{A}=\left(\mat{a_1}, \ldots, \mat{a_N}\right)$.

Given this notation, the marginal likelihood of DAGP can be separated into the likelihood, the latent function processes, and the assignment process and is given by,

$\Prob*{\mat{F}\given\mat{X}}=\prod_{k=1}^K \Gaussian*{\mat{F^{\pix{k}}}\given\Fun*{\mu^{\pix{k}}}{\mat{X}}, \Fun*{\K^{\pix{k}}}{\mat{X}, \mat{X}}}$ with mean function $\mu^{\pix{k}}$ and kernel $\K^{\pix{k}}$.

First, we assume that the $\mat{a_n}$ are drawn independently from multinomial distributions with logit parameters $\mat{\alpha_n}=\left(\alpha_n^{\pix{1}}, \ldots, \alpha_n^{\pix{K}}\right)$.

One approach to specify $\mat{\alpha_n}$ is to assume them to be known a priori and to be equal for all data points~\parencite{lazaro-gredilla_overlapping_2012}.

Specifically, we assume that there is a relationship between the location in the input space $\mathbf{x}$ and the associations.

By placing independent GP priors on $\mat{\alpha^{\pix{k}}}$, we can encode our prior knowledge of the associations by the choice of covariance function

The prior on the assignments $\mat{A}$ is given by marginalizing the $\mat{\alpha^{\pix{k}}}$, which, when normalized, parametrize a batch of multinomial distributions,

Modelling the relationship between the input and the associations allows us to efficiently model data, which, for example, is unimodal in some parts of the input space and bimodal in others.

Since the GPs of the $\mat{\alpha^{\pix{k}}}$ use a zero mean function, our prior assumption is a uniform distribution of the different generative processes everywhere in the input space.

If inference on the $\mat{a_n}$ reveals that, say, all data points at similar positions in the input space can be explained by the same $\nth{k}$ process, the belief about $\mat{\alpha}$ can be adjusted to make a non-uniform distribution favorable at this position, thereby increasing the likelihood via $\Prob*{\mat{A}\given\mat{X}}$.

This mechanism introduces an incentive for the model to use as few functions as possible to explain the data and importantly allows us to predict a relative importance of these functions when calculating the posterior of the new observations $\mat{x_\ast}$.

\Cref{fig:dynamic_graphical_model} shows the resulting graphical model, which divides the generative process for every data point in the application of the latent functions on the left side and the assignment process on the right side.

The interdependencies between the data points are introduced through the GP priors on $\rv{f_n^{\pix{k}}}$ and $\rv{\alpha_n^{\pix{k}}}$ and depend on the hyperparameters $\mat{\theta}=\Set*{\mat{\theta^{\pix{k}}}, \mat{\theta_\alpha^{\pix{k}}}, \sigma^{\pix{k}}}_{k=1}^K$.

Since many real word systems are inherently hierarchical, prior knowledge can often be formulated more easily using composite functions~\parencite{kaiser_bayesian_2018}.

This bound can be sampled efficiently and allows us to optimize both the models for the different processes $\Set*{f^{\pix{k}}}_{k=1}^K$ and our belief about the data assignments $\Set*{\mat{a_n}}_{n=1}^N$ simultaneously using stochastic optimization.

As first introduced by~\textcite{titsias_variational_2009}, we augment all GP in our model using sets of $M$ inducing points $\mat{Z^{\pix{k}}}=\left(\mat{z_1^{\pix{k}}}, \ldots, \mat{z_M^{\pix{k}}}\right)$ and their corresponding function values $\mat{u^{\pix{k}}}=\Fun*{f^{\pix{k}}}{\mat{Z^{\pix{k}}}}$, the inducing variables.

We collect them as $\mat{Z}=\Set*{\mat{Z^{\pix{k}}}, \mat{Z_\alpha^{\pix{k}}}}_{k=1}^K$ and $\mat{U}=\Set*{\mat{u^{\pix{k}}}, \mat{u_\alpha^{\pix{k}}}}_{k=1}^K$.

Taking the function $f^{\pix{k}}$ and its corresponding GP as an example, the inducing variables $\mat{u^{\pix{k}}}$ are jointly Gaussian with the latent function values $\mat{F^{\pix{k}}}$ of the observed data by the definition of GPs.

We follow~\parencite{hensman_gaussian_2013} and choose the variational approximation $\Variat*{\mat{F^{\pix{k}}}, \mat{u^{\pix{k}}}}=\Prob*{\mat{F^{\pix{k}}}\given\mat{u^{\pix{k}}}, \mat{X}, \mat{Z^{\pix{k}}}}\Variat*{\mat{u^{\pix{k}}}}$ with $\Variat*{\mat{u^{\pix{k}}}}=\Gaussian*{\mat{u^{\pix{k}}}\given\mat{m^{\pix{k}}}, \mat{S^{\pix{k}}}}$.

This formulation introduces the set $\Set*{\mat{Z^{\pix{k}}}, \mat{m^{\pix{k}}}, \mat{S^{\pix{k}}}}$ of variational parameters indicated in~\cref{fig:dynamic_graphical_model}.

A central assumption of this approximation is that given enough well-placed inducing variables $\mat{u^{\pix{k}}}$, they are a sufficient statistic for the latent function values $\mat{F^{\pix{k}}}$.

This implies conditional independence of the $\mat{f_n^{\pix{k}}}$ given $\mat{u^{\pix{k}}}$ and $\mat{X}$.

Given a set of assignments $\mat{A}$, this factorization along the data points is preserved in our model due to the assumed independence of the different functions in~\cref{eq:true_marginal_likelihood}.

Following the ideas of doubly stochastic variational inference (DSVI) presented by~\textcite{salimbeni_doubly_2017} in the context of deep GPs, we maintain these correlations between different parts of the model while assuming factorization of the variational distribution.

An important property of the variational bound for DSVI~\parencite{salimbeni_doubly_2017} is that taking samples for single data points is straightforward and can be implemented efficiently.

Specifically, for some $k$ and $n$, samples $\mat{\hat{f}_n^{\pix{k}}}$ from $\Variat*{\mat{f_n^{\pix{k}}}}$ are independent of all other parts of the model and can be drawn using samples from univariate unit Gaussians using reparametrizations~\parencite{kingma_variational_2015,rezende_stochastic_2014}.

Note that it would not be necessary to sample from the different processes, since $\Variat*{\mat{F^{\pix{k}}}}$ can be computed analytically~\parencite{hensman_gaussian_2013}.

However, we apply the sampling scheme to the optimization of both the assignment processes $\mat{\alpha}$ and the assignments $\mat{A}$ as for $\mat{\alpha}$, the analytical propagation of uncertainties through the $\softmax$ renormalization and multinomial likelihoods is intractable but can easily be evaluated using sampling.

We optimize $\Ell_{\text{DAGP}}$ to simultaneously recover maximum likelihood estimates of the hyperparameters $\mat{\theta}$, the variational parameters $\Set*{\mat{Z}, \mat{m}, \mat{S}}$, and assignments $\mat{A}$.

For every $n$, we represent the belief about $\mat{a_n}$ as a $K$-dimensional discrete distribution $\Variat*{\mat{a_n}}$.

This distribution models the result of drawing a sample from $\Multinomial*{\mat{a_n}\given\Fun{\softmax}{\mat{\alpha_n}}}$ during the generation of the data point $(\mat{x_n}, \mat{y_n})$.

Since we want to optimize $\Ell_{\text{DAGP}}$ using (stochastic) gradient descent, we need to employ a continuous relaxation to gain informative gradients of the bound with respect to the binary (and discrete) vectors $\mat{a_n}$.

One straightforward way to relax the problem is to use the current belief about $\Variat*{\mat{a_n}}$ as parameters for a convex combination of the $\mat{f_n^{\pix{k}}}$, that is, to approximate $\mat{f_n}\approx\sum_{k=1}^K \Variat*{\mat{a_n^{\pix{k}}}}\mat{\hat{f}_n^{\pix{k}}}$.

Explaining data points as mixtures of the different generating processes violates the modelling assumption that every data point was generated using exactly one function but can substantially simplify the learning problem.

Instead of directly using $\Variat*{\mat{a_n}}$ to combine the $\mat{f_n^{\pix{k}}}$, we first draw a sample $\mat{\hat{a}_n}$ from a concrete random variable as suggested by~\textcite{maddison_concrete_2016}, parameterized by $\Variat*{\mat{a_n}}$.

Based on a temperature parameter $\lambda$, a concrete random variable enforces sparsity but is also continuous and yields informative gradients using automatic differentiation.

Samples from a concrete random variable are unit vectors and for $\lambda\to0$ their distribution approaches a discrete distribution.

Our approximate evaluation of the bound in \cref{eq:variational_bound} during optimization has multiple sources of stochasticity, all of which are unbiased.

First, we approximate the expectations using Monte Carlo samples $\mat{\hat{f}_n^{\pix{k}}}$, $\mat{\hat{\alpha}_n^{\pix{k}}}$, and $\mat{\hat{a}_n}$.

And second, the factorization of the bound along the data allows us to use mini-batches for optimization~\parencite{salimbeni_doubly_2017, hensman_gaussian_2013}.

The predictive posteriors of the $K$ functions $\Variat*{\mat{f_\ast^{\pix{k}}}\given\mat{x_\ast}}$ are given by $K$ independent shallow GPs and can be calculated analytically~\parencite{hensman_gaussian_2013}.

Samples from the predictive density over $\Variat*{\mat{a_\ast}\given\mat{x_\ast}}$ can be obtained by sampling from the GP posteriors $\Variat*{\mat{\alpha_\ast^{\pix{k}}}\given\mat{x_\ast}}$ and renormalizing the resulting vector $\mat{\alpha_\ast}$ using the $\softmax$-function.

The distribution $\Variat*{\mat{a_\ast}\given\mat{x_\ast}}$ reflects the model's belief about how many and which of the $K$ generative processes are relevant at the test location $\mat{x_\ast}$ and their relative probability.

Analogously to~\parencite{salimbeni_doubly_2017}, our new prior assumption about the $\nth{k}$ latent function values $\Prob*{\mat{F^{\prime\pix{k}}}\given\mat{X}}$ is given by,

Similar to the single-layer case, we introduce sets of inducing points $\mat{Z_l^{\prime\pix{k}}}$ and a variational distribution over their corresponding function values $\Variat*{\mat{u_l^{\prime\pix{k}}}}=\Gaussian*{\mat{u_l^{\prime\pix{k}}}\given\mat{m_l^{\prime\pix{k}}}, \mat{S_l^{\prime\pix{k}}}}$.

We collect the latent multi-layer function values as $\mat{F^\prime}=\Set{\mat{F_l^{\prime\pix{k}}}}_{k=1,l=1}^{K,L}$ and corresponding $\mat{U^\prime}$ and assume an extended variational distribution,

As the $\nth{n}$ marginal of the $\nth{L}$ layer depends only on the $\nth{n}$ marginal of all layers above sampling from them remains straightforward~\parencite{salimbeni_doubly_2017}.

To calculate the first term, samples have to be propagated through the deep GP structures.

This extended bound thus has complexity $\Fun*{\Oh}{NM^2LK}$ to evaluate in the general case and complexity $\Fun*{\Oh}{NM^2\cdot\Fun{\max}{L, K}}$ if the assignments $\mat{a_n}$ take binary values.

We use an implementation of DAGP in TensorFlow~\parencite{tensorflow2015-whitepaper} based on GPflow~\parencite{matthews_gpflow_2017} and the implementation of DSVI~\parencite{salimbeni_doubly_2017}.

All models can solve standard regression problems and yield unimodal predictive distributions or, in case of multi-layer perceptrons (MLP), a single point estimate.

Mixture density networks (MDN)~\parencite{bishop_mixture_1994} and the infinite mixtures of Gaussian processes (RGPR)~\parencite{rasmussen_infinite_2002} model yield multi-modal posteriors through mixtures with many components but do not solve an association problem.

Similarly, Bayesian neural networks with added latent variables (BNN+LV)~\parencite{depeweg_learning_2016} represent such a mixture through a continuous latent variable.

Both the overlapping mixtures of Gaussian processes (OMGP)~\parencite{lazaro-gredilla_overlapping_2012} model and DAGP explicitly model the data association problem and yield independent models for the different generating processes.

However, OMGP assumes global relevance of the different modes.

In contrast, DAGP infers a spacial posterior of this relevance.

Additional to a joint posterior explaining the data, DAGP also gives insight into the relative importance of the different processes in different parts of the input space.

DAGP is able to explicitly recover the changing number of modes in a data set.

\emph{Separate models for independent generating processes avoid model pollution.}

We apply DAGP to a one-dimensional regression problem with uniformly distributed asymmetric outliers in the training data.

We use a task proposed by~\textcite{choi_choicenet_2018} where we sample $x \in[-3, 3]$ uniformly and apply the function $\Fun{f}{x}=(1-\delta)(\Fun{\cos}{\sfrac{\pi}{2}\cdot x}\Fun{\exp}{-(\sfrac{x}{2})^2}+\gamma)+\delta\cdot\epsilon$, where $\delta\sim\Fun{\Ber}{\lambda}$, $\epsilon\sim\Fun{\Uni}{-1, 3}$ and $\gamma\sim\Gaussian{0, 0.15^2}$.

To avoid pathological solutions for high outlier ratios, we add a prior to the likelihood variance of the first process, which encodes our assumption that there actually is a signal in the training data.

The model proposed in~\parencite{choi_choicenet_2018}, called ChoiceNet (CN), is a specific neural network structure and inference algorithm to deal with corrupted data.

Since we can encode the same prior knowledge about the signal and noise processes in both OMGP and DAGP, the results of the two models are comparable:

For low outlier rates, they correctly identify the outliers and ignore them, resulting in a predictive posterior of the signal equivalent to standard GP regression without outliers.

In the special case of 0\,\% outliers, the models correctly identify that the process modelling the noise is not necessary, thereby simplifying to standard GP regression.

For high outlier rates, stronger prior knowledge about the signal is required to still identify it perfectly.

\Cref{fig:choicenet} shows the DAGP posterior for an outlier rate of 60\,\%.

While the function has still been identified well, some of the noise is also explained using this process, thereby introducing slight errors in the predictions.

The DAGP posterior on an artificial data set with bimodal and trimodal parts.

The joint predictions (top) are mixtures of four Gaussians weighed by the assignment probabilities $\mat{\alpha}$ (bottom).

The weights are represented via the opacity of the modes.

The model has learned that the mode $k =2$ is irrelevant, that the mode $k =1$ is only relevant around the interval $[0, 5]$.

Outside this interval, the mode $k =3$ is twice as likely as the mode $k =4$.

The concrete assignments $\mat{a}$ (middle) of the training data show that the mode $k =1$ is only used to explain observations where the training data is trimodal.

We uniformly sample 350 data points in the interval $x \in[-2\pi, 2\pi]$ and obtain $y_1=\Fun{\sin}{x}+\epsilon$, $y_2=\Fun{\sin}{x}-2\Fun{\exp}{-\sfrac{1}{2}\cdot(x-2)^2}+\epsilon$ and $y_3=-1-\sfrac{3}{8\pi}\cdot x +\sfrac{3}{10}\cdot\Fun*{\sin}{2x}+\epsilon$ with additive independent noise $\epsilon\sim\Gaussian*{0, 0.005^2}$.

The resulting data set $\D=\Set{\left( x, y_1\right), \left( x, y_2\right), \left( x, y_3\right)}$ is trimodal in the interval $[0, 5]$ and is otherwise bimodal with one mode containing double the amount of data than the other.

We use squared exponential kernels as priors for both the $f^{\pix{k}}$ and $\alpha^{\pix{k}}$ and $25$ inducing points in every GP.

The figure shows the posterior belief about the assignments $\mat{A}$ and illustrates that DAGP recovered that it needs only three of the four available modes to explain the data.

Importantly, DAGP does not only cluster the data with respect to the generating processes but also infers a factorization of the input space with respect to the relative importance of the different processes.

The model has disabled the mode $k =2$ in the complete input space and has learned that the mode $k =1$ is only relevant in the interval $[0, 5]$ where the three enabled modes each explain about a third of the data.

Outside this interval, the model has learned that one of the modes has about twice the assignment probability than the other one, thus correctly reconstructing the true generative process.

The DAGP is implicitly incentivized to explain the data using as few modes as possible through the likelihood term of the inferred $\mat{a_n}$ in \cref{eq:variational_bound}.

At $x =-10$ the inferred modes and assignment processes start reverting to their respective priors away from the data.

We evaluate the two inferred sub-models for the default system (center) and short-pole system (right).

We provide gray baseline comparisons with BNN+LV and GPR models which cannot solve the data assignment problem.

BNN+LV yields joint predictions which cannot be separated into sub-models.

Specialized GPR models trained the individual training sets give a measure of the possible performance if the data assignment problem would be solved perfectly.