Probabilistic Regressor Chains with Monte Carlo Methods

TL;DR

This paper introduces a probabilistic approach to regressor chains using Monte Carlo methods, improving multi-output regression by addressing previous limitations and enhancing flexibility and effectiveness.

Contribution

It develops a sequential Monte Carlo scheme for probabilistic regressor chains, overcoming prior greedy inference limitations and expanding applicability in multi-output regression tasks.

Findings

Monte Carlo methods improve regression chain performance

Probabilistic regressor chains are flexible across data types

Enhanced understanding of classifier and regressor chains

Abstract

A large number and diversity of techniques have been offered in the literature in recent years for solving multi-label classification tasks, including classifier chains where predictions are cascaded to other models as additional features. The idea of extending this chaining methodology to multi-output regression has already been suggested and trialed: regressor chains. However, this has so-far been limited to greedy inference and has provided relatively poor results compared to individual models, and of limited applicability. In this paper we identify and discuss the main limitations, including an analysis of different base models, loss functions, explainability, and other desiderata of real-world applications. To overcome the identified limitations we study and develop methods for regressor chains. In particular we present a sequential Monte Carlo scheme in the framework of a probabilistic regressor chain, and we show it can be effective, flexible and useful in several types of data. We place regressor chains in context in general terms of multi-output learning with continuous outputs, and in doing this shed additional light on classifier chains.

Tables

Table 1. Table 1 : Results of 10-fold cross validation.

Dataset	$L$	IR.B	IR.K	RC.B	RC.K	MC.B	MC.D	PF.R/B	PF.N/B
Synth	2	1.03	1.02	1.04	1.05	1.44	1.02	1.51	1.42
Andro	6	0.60	0.24	0.50	0.21	0.77	0.91	1.04	1.02
EDM	2	0.61	0.43	0.60	0.42	0.97	0.80	0.98	1.17
ENB	2	0.10	0.01	0.10	0.01	0.15	0.58	0.16	0.20
Jura	3	0.38	0.35	0.38	0.35	0.57	0.77	0.58	0.50
Slump	3	0.49	0.40	0.49	0.40	0.80	0.85	0.64	0.54
Avg Rank		3.67	2.67	2.33	1.67	6.00	6.33	6.67	6.67

Equations

$$\hat{\mathbf{y}}=[{\hat{y}}_1,\dots ,{\hat{y}}_L]=[h_1(\mathbf{x}),\dots ,h_L(\mathbf{x})]$$

$$\hat{\mathbf{y}}=[h_1(\mathbf{x}),h_2(\mathbf{x},{\hat{y}}_1),\dots ,h_L(\mathbf{x},{\hat{y}}_1,\dots ,{\hat{y}}_{L-1})]$$

$$\hat{\mathbf{y}}=\mathbf{h}(\mathbf{x})$$

$$=\underset{\mathbf{y}\in \mathcal{Y}}{\mathrm{argmax}}\prod _{j=1}^{L}P(y_j|\mathbf{x},y_1,\dots ,y_{j-1}),$$

$$y_j^{(m)}$$

$$w_j^{(m)}$$

$$\hat{\mathbf{y}}$$

$${\hat{y}}_2$$

$$={\mathbf{w}}_2^{\top }[x,{\hat{y}}_j]$$

$$={\mathbf{w}}_2^{\top }[x,w_{1}x]$$

$$=w_{2,1}x+w_{2,2}{w}_{1}x$$

$$=x(w_{2,1}+w_{2,1}{w}_1)=w^{\prime }x$$

$$\hat{\mathbf{\mu }}$$

$${\hat{\mu }}_j$$

$$=\int y_{j}p(y_1,\dots ,y_L|\mathbf{x})d{y}_1\cdots d{y}_L=\int y_j{p}_j(y_j|\mathbf{x})d{y}_j$$

$$p(\mathbf{y}|\mathbf{x})=p(y_1,y_2,\dots ,y_L|\mathbf{x})=p(y_1|\mathbf{x})\prod _{j=2}^{L}p(y_j|\mathbf{x},y_1,\dots ,y_{j-1})$$

$$y_j^{(m)}\sim p(y_j|\mathbf{x},y_1^{(m)},\dots ,y_{j-1}^{(m)}),$$

$${\mathbf{y}}^{(m^{\ast })}=\mathrm{argmax}\left\{p({\mathbf{y}}^{(1)}|\mathbf{x}),\cdots ,p({\mathbf{y}}^{(M)}|\mathbf{x})\right\},$$

$$\hat{\mathbf{y}}=g\left(\{{\mathbf{y}}^{(m)},{\mathbf{w}}^{(m)}{\}}_{m=1}^M\right)$$

$${\mathbb{E}}_{\mathbf{Y}\sim \mathcal{D}}[c(\mathbf{Y},\mathbf{h}(\mathbf{x}))|\mathbf{x}]$$

$$p_j=\mathcal{N}(y_j|{\mathbf{\mu }}_j,{\mathbf{\Sigma }}_j)$$

$${\theta }^{\ast }=\underset{\theta }{\mathrm{argmin}}\text{KL}(q(y_j|\theta )\Vert p_j)$$

$$p(y_j,\mathbf{z})=\frac{1}{N}\sum _{i=1}^N{K}_h(\mathbf{z},{\mathbf{z}}_i)K_h(y_j,y_i)\text{and}p(\mathbf{z})=\frac{1}{N}\sum _{i=1}^N{K}_h(\mathbf{z},{\mathbf{z}}_i)$$

$$\mathcal{M}:\left\{\begin{array}[]{l}f(y_{j}|y_{j-1})\\ \ell(\mathbf{x}_{j}|y_{j})\end{array}\right.$$

$$y_j^{(m)}$$

$$w_j^{(m)}$$

$${\hat{y}}_j=\mathbb{E}[Y_j]$$

$$f(y_j|y_1,\dots ,y_{j-1})$$

$$\ell (\mathbf{x})\propto p(y_j|\mathbf{x},y_1^{(m)},\dots ,y_{j-1}^{(m)})$$

$$\widehat{ESS}({\bar{w}}_{1:M})=\frac{1}{{\sum }_{i=1}^M{\bar{w}}_i^2}\text{or}\widehat{ESS}({\bar{w}}_{1:M})=\frac{1}{\max {\bar{w}}_i}$$

$$w_j^{(m)}=w_{j-1}^{(m)}\frac{\ell (y_j^{(m)}|\mathbf{x},{\tilde{y}}_1^{(m)},\dots ,{\tilde{y}}_{j-1}^{(m)})}{f_j(y_j^{(m)}|{\tilde{y}}_1^{(m)},\dots ,{\tilde{y}}_{j-1}^{(m)})}$$

$$\{{\tilde{y}}_j^{(1)},\dots ,{\tilde{y}}_j^{(M)}\}\leftarrow \{{\tilde{y}}_{j,K}^{(1)},\dots ,{\tilde{y}}_{j,K}^{(M)}\}$$

$$\text{MSE}=\frac{1}{N}\sum _{i=1}^N\sum _{j=1}^L(y_{ij}-{\hat{y}}_{ij}{)}^2\text{and }\text{MAE}=\frac{1}{N}\sum _{i=1}^N\sum _{j=1}^L|y_{ij}-{\hat{y}}_{ij}|$$

$$\text{0/1}\approx \frac{1}{N}\sum _{i=1}^N\{\begin{array}{ll}0& \text{if ∥y^i−yi∥2<c}\\ 1& \text{otherwise}\end{array}$$

$$y_1$$

Full text

1 Introduction

Multi-dimensional data is ever-more present in industrial and scientific contexts. For example, multi-label classification has made a significant impact in the machine learning literature over recent years, where data points are naturally associated with multiple outputs111Distinguished from multi-class classification in that a single output may take multiple values – but only one such value is assigned per output. Rather than a naive method of building one model per output; advanced methods can model the outputs together, resulting in better predictive performance. Although the potential of individual models is periodically revived under particular scenarios, the vast majority of literature proposes joint modeling, and as a result show improvement in both predictive performance and efficiency. A recent review of this area containing many useful references, is given in [41].

An established method in multi-label classification is that of classifier chains, where a model is trained for each label, but estimates of the other models are used as additional features, in a cascaded chain along the target labels. This is exemplified (and contrasted to the naive approach) in Figure 1.

This ‘chaining’ mechanism, although simple in its basic form, has proved successful in multi-label classification and provided dozens of modifications, extensions, and improvements, in the literature (see, e.g., [36, 10, 34, 11, 41, 35, 32] and references therein).

It is flexible, in the sense that any base model case be used to predict each label, and powerful, in the sense that even relatively simple base models lead to greater predictive performance than if they had been used separately.

Given its impact in multi-label classification, the motivation of this paper is to more-thoroughly investigate the use of the chaining approach in the regression context with continuous output variables. This investigation is needed since only greedy inference can be applied directly, but even in this case the application of chaining to the regression context is not straightforward and the effectiveness of such application is not yet widely explored.

After summarizing background work (Section 2), we provide a rigorous discussion and development of probabilistic regressor chains, including a survey of possible approaches (Sections 3 and 4). Following this, we develop a sequential Monte Carlo scheme (Section 5). This scheme allows sampling of candidate paths through the target space, which is useful for many applications. We implement and test several of these on synthetic and real-world datasets involving multiple continuous outputs (Section 6), the results of which we reflect upon in detailed discussion (Section 7). After looking at connections to related and potential future work (Section 8) we draw conclusions and make recommendations (Section 9).

2 Chain Methods

Given a dataset $\mathcal{D}=\{(\mathbf{x}_{i},\mathbf{y}_{i})\}_{i=1}^{N}$, where $\mathbf{x}_{i}\in\mathcal{X}\subseteq\mathbb{R}^{D}$ and $\mathbf{y}_{i}\in\mathcal{Y}\subseteq\mathbb{R}^{L}$, we are interested building a model that can provide predictions corresponding to $L$ target variables222In cases where it is not clear from context, we will denote $y_{ij}$ as the value of the $j$-th attribute of the $i$-th instance, $\mathbf{\hat{y}}=[\hat{y}_{1},\ldots,\hat{y}_{L}]$, for any given input observation $\mathbf{x}=[x_{1},\ldots,x_{D}]$.

In the multi-output classification scenario [36, 10, 34, 12], the $j$-th output takes some value $y_{j}\in\{1,\ldots,K_{j}\}$ (we may highlight the popular case of binary outputs where each $y_{i}\in\{0,1\}$, known often as multi-label classification) or – in the case of continuous outputs, i.e., multi-output regression – then each $y_{j}\in\mathbb{R}$.

A naive approach is to simply build a separate model for each target independently, such that

[TABLE]

where each model $h_{j}$ has been built from training set $\{(\mathbf{x}_{i},y_{ij})\}_{i=1}^{N}$, as a traditional single-output classifier or regressor according to the domain of $Y_{j}$ (clearly, a regressor, if $Y_{j}\in\mathbb{R}$). The exact class of model (type of regressor) is largely dependent on preference and/or driven by domain assumptions and constraints.

Particularly with regard to multi-output classification, a large volume of literature proposes a plethora of models that out-compete this baseline by modelling the dependencies among outputs; consider [36, 10, 40, 34, 11, 12, 41] and references therein. The multi-output regression literature is smaller in volume (consider, e.g., [3, 38]), being both more recently considered for specialized methods, and being more difficult to obtain better results with such methods – precisely a motivating factor behind this work.

In prior work ([36, 10, 34, 12], and many others), a chaining mechanism was studied for multi-label classification, known as classifier chains. Due to the relatively large amount of work in the literature on this topic, including continued recent interest (e.g., [20, 39] among many others), we argue it is worth studying this mechanism in particular with application to multi-output regression, which we will denote regressor chains.

It is useful to study regressor chains in light of previous developments in classifier chains. In its simplest form, the estimates of other labels are used as feature inputs for the following classifier, thus augmenting the input space along the chain. Given some test instance $\mathbf{x}$, we may obtain an estimate of $\mathbf{\hat{y}}=[\hat{y}_{1},\ldots,\hat{y}_{L}]$ as

[TABLE]

where a recursion takes placed based on $\hat{y}_{j}=h_{j}(\ldots)$.

In the classification context, this is known as a classifier chain and, specifically, one with greedy inference. Note that models $h_{j}$ can be estimated individually and in parallel at training time, similarly to those in Eq. (1), and – likewise – those models may be any class appropriate for each individual target (a binary classifier for binary output, etc).

In this case of greedy inference, a single path/combination is considered. However, note that this $\mathbf{\hat{y}}$ in particular is just one possible prediction, and not necessarily the best. The addition of probabilistic inference was an important development [10, 11], where a search is carried out, typically to obtain a MAP estimate333Although configurations for other losses are also possible; see [10], thus maximizing $0/1$ loss; hence

[TABLE]

where $\mathcal{Y}$ is the space of possible paths over the $L$ outputs and $P(y_{j}|\mathbf{x},\cdot)$ is a pmf associated with the $j$-th model (henceforth we use the abbreviation $P_{j}$ where suitable). This is known as a probabilistic classifier chain. The need for a probabilistic interpretation (i.e., pmfs $P_{j}$) can be addressed by building such models; logistic regression is a common choice.

One may observe that the search space in Eq. (3) is exponential with the number of labels $L$. Indeed, complete inference (all branches of the probability tree are explored) is usually prohibitive. Therefore, one may consider the search space as a probability tree, and conduct a tree search, where each $P_{j}$ provides a weight on each branch. Figure 3 offers some visual intuition.

Many search methods have been applied for this purpose (a survey is given in [26]). In particular we highlight the approach of Monte-Carlo sampling [34] as – unlike many of the search approaches – is directly applicable to regressor chains, and as such is most relevant to the material developed in Section 5 of this work. In this approach, $M$ samples are taken, $m=1,\ldots,M$, and weighted as

[TABLE]

where $w^{(m)}_{0}=1$. An example is shown in Figure 3 where $M=2$. A final prediction is obtained as

[TABLE]

where $m^{*}=\operatornamewithlimits{argmax}_{m\in\{1,\ldots,M\}}w_{L}^{(m)}$ (index of the maximum weight), and where $\mathbf{y}^{(m)}=[y^{(m)}_{1},\ldots,y^{(m)}_{L}]$ is a complete sequence across label indices. Note that complexity is determined by $M\ll 2^{L}$.

It has been noticed that a chaining procedure can be applied in an off-the-shelf manner to the multi-output regression context involving continuous targets, namely with greedy inference; see, e.g., [3]. However, as we highlight and discuss in the following section, even though the application is identical, there are some major differences in the way regressor chains behave and the relative results they obtain (with respect to non-chained methods).

3 The Poor Behaviour of Regressor Chains

Applying greedy inference in chains in the case of regression is – exactly as in classification – a case of each output simply being “plugged in” to the following model as an additional feature. Recall that this simply means that predictions $\hat{y}_{1},\ldots,\hat{y}_{j-1}$ are treated as fixed observations (i.e., and not random variables) when inferring $y_{j}$. Several papers have trialed this approach, e.g., [3]. Unfortunately, predictive performance is not as impressive wrt independent regression models. We explain why.

Recall Eq. (2) from above. A first choice of model in the regression context may be least squares where $\hat{y}_{j}=h_{j}(\mathbf{x})=\mathbf{w}_{j}^{\top}[x_{1},\ldots,x_{D},\hat{y}_{1},\ldots,\hat{y}_{j-1}]$. Plugging this in to the equation, and supposing for simplicity $D=1$ input and $L=2$ outputs, we observe that444Let $\mathbf{w}_{j}=[w_{j,1},w_{j,2}]$

[TABLE]

and thus clearly we see that label predictions $\hat{y}_{1},\ldots,\hat{y}_{j-1}$ are superfluous in predicting $\hat{y}_{j}$; the regressor chain in this context simply performs a series of linear transformations which can naturally be represented as a single transformation.

It may be tempting to argue for application to the case where $x_{j}$ is observed only at time $x_{j}$. In this scenario we presume that inputs are received over time, thus estimating $\hat{y}_{j}=h(x_{1},\ldots,x_{j},\hat{y}_{1},\ldots,\hat{y}_{j-1})$. Unfortunately, we may show that $\hat{y}_{1},\ldots,\hat{y}_{j-1}$ are still not needed wrt $\hat{y}_{j}$ since all available information already comes from $x_{1},\ldots,x_{j}$ directly via earlier labels. This is seen clearly in illustration; Figure 3.

Another fundamental issue is the selection of loss metric. An obvious and popular choice for regression is based on the mean squared error (MSE) loss criterion (as indeed considered in the earlier experimentation of regressor chains in [3, 38] among others).

Under random variables $\bm{Y}=Y_{1},\ldots,Y_{L}$, we see that the minimum MSE (MMSE) estimator under observation $\mathbf{x}$ is given as follows,

[TABLE]

wrt marginal densities $p_{j}$ (homologous to the pmfs of Eq. (4)). Let us emphasise that $\hat{\bm{\mu}}=[\hat{\mu}_{1},\ldots,\hat{\mu}_{L}]$ is a vector of $L$ integrals where

[TABLE]

If we assume a Gaussian distribution $\epsilon_{j}\in\mathcal{N}(0,\sigma_{j}^{2})$ (a natural view of ordinary least squares [17]) then again we may connect to the case discussed at the beginning of this section, where $\hat{\mu}_{1}=\mathbf{w}^{\top}_{1}\mathbf{x}$ and so on for $\hat{\mu}_{1},\ldots,\hat{\mu}_{L}$.

More generally, we are not only interested in obtaining point-wise estimators $\mu_{j}$, but we are also interested in extracting all the statistical information encoded within this posterior $p(\mathbf{y}|\mathbf{x})$, such as uncertainty measures, credible intervals, and quantiles as well. So, our goal is to approximate complex integrals involving this posterior

[TABLE]

Having an estimate of the shape of $p$ then allows us to estimate other values aside from the expected value (i.e., Eq. (8)) such as the median or the mode. We remark that classifier chains typically predicts a mode in the form of a MAP estimate.

Notice, however, that – in spite of using the chain rule to factorize the joint in this manner this estimate (and hence related integrals, Eq. (8)) are generally intractable. This is unlike in probabilistic classifier chains where each pmf $P_{j}$ models a variable taking a finite number of discrete values – inherently more tractable in most cases.

In all cases an exhaustive search in the joint is intractable for a large number of outputs $L$, hence efficient tree search methods such as Monte Carlo sampling (again, refer to Figure 3). However, a tree cannot be formed on continuous output space, and hence no such search is initially possible in the regression context. Discussion of the output space also brings us to identify a further pitfall in regression chains: the possibility of extensive “error propagation”. In classifier chains, this is well known, and easily detectable (e.g., in ). However, the $P_{j}\in[0,1]$ space is inherently limited. On the other hand, under regressior chains, the estimate of greedy inference may completely degenerate and become lost in $p_{j}\in\mathbb{R}$ space, we progress down the chain towards $j=L$.

Therefore we have isolated and discussed some main drawbacks to the application of greedy regressor chains: the type of loss functions often considered in regression (e.g., MSE) totally change the behaviour of the chaining methodology, and with a linear regression scheme, such as least squares, there is no benefit to the chaining mechanism. Even worse, error propagation may be catastrophic.

Having identified the lacking in regressor chains, as compared to their classification counterpart, in the following section we study the idea of probabilistic regressor chains which we later build on to suggest new methods.

4 Probabilistic Regressor Chains with Monte Carlo Search

In probabilistic classifier chains, inference may be cast as a search in the probability tree, as discussed above. Many tree-search methods are applicable. We gave an example of Monte Carlo search which takes $M$ weighted samples across the tree (recall Eq. (4) and Eq. (5)). In the regression context there is no inherent tree to search. We could however create a tree by using particular values $y_{j}^{(m)}\in\mathbb{R}$ as nodes. This requires modeling of the marginals in such a way that we may draw samples

[TABLE]

and in particular these samples should be at an ‘interesting’ part of the space, such as near the modes.

Figure 4 shows an example illustration of samples $\{y_{j}^{(m)}\}$ forming paths and thus a tree with each complete path from root to leave as a possible prediction $\mathbf{y}\in\mathbb{R}^{L}$.

This is an interesting idea, as with a cloud of samples, we can not only obtain a tree and apply existing methods used in the classification context,, but also easily obtain an approximation a MAP estimator as well, as in Eq. (6). Explicitly we can write this as:

[TABLE]

and more generally (estimators other than MAP), a prediction may be given as path

[TABLE]

chosen by function $g$ over a tree with branches defined across nodes $\{\mathbf{y}^{(m)}\}$ and weighted by $\{w^{(m)}_{j}\}$, i.e., holding the marginal branch cost from depth $j-1$ to $j$ by the $m$-th sample. Recall Figure 4. It means that Eq. (11) makes use of the fact that $p(\mathbf{y}^{(m)}|\mathbf{x})=\prod_{j=1}^{L}w^{(m)}_{j}$. In the remainder of this work we occasionally omit the subscript and use $w^{(m)}=p(\mathbf{y}^{(m)}|\mathbf{x})$ as the path cost from root to a leaf, i.e., $w^{(m)}=p(\mathbf{y}^{(m)}|\mathbf{x})$ where $p$ as in Eq. (4) (for notational simplicity).

Suppose we define a cost function $c$ based on path costs $w^{(m)}_{j}$. For example, cost is given whenever crossing a low-density region according to the estimated $p_{j}$s. To minimize expected loss,

[TABLE]

we use our estimation of $p$ (and decomposition as ) and Monte Carlo samples to replace the integral, thus giving something like Eq. (11) (depending on the choice of loss function).

A fundamental consideration is how to estimate each component $p_{j}$ of $p$, and how to draw samples from it. The major issue at stake is that that a powerful non-linear predictor is required to justify chains (as discussed above) yet such a predict may not be ideal for sampling from. We address this issue with a sequential Monte Carlo method, presented in the following section. But let us first look briefly at some available alternatives, all of which could be considered a variation of probabilistic regressor chains.

4.1 Discretization and classification

A first approach to the inference problem is to simply discretize the label space and proceed from a classification perspective. This is not the same as taking samples (which is done at prediction time). As explained above, the classification perspective of chains is easier to justify, and a wealth of methods available. Of course not all problems are suited to discretization of the label space, but often it is suitable and even common practice, for example in the domain of reinforcement learning (see, e.g., [42]).

In this case we are learning $p(y_{j}\in\mathcal{S}^{(k)}_{j}|\mathbf{z})\Leftrightarrow P(y_{j}=k)$ (the RHS relates to Eq. (3)) where bin $\mathcal{S}^{(k)}_{j}\subset\mathbb{R}$; and thus removing the need to estimate any integrals. Clearly it is fundamental to choose a suitable set of bins (which correspond to a finite set of class labels).

4.2 Bayesian regression

From the perspective of a greedy chain (earlier labels are fixed input), we may elabourate $p_{j}$ as Bayesian linear regression:

[TABLE]

with sufficient statistics $\bm{\mu}_{j}$ and $\bm{\Sigma}_{j}$. As a Gaussian, sampling from this distribution (e.g., as in Eq. (10)) is straightforward. However, it based on a linear combination of $\bm{\mu}_{j}=\bm{\theta}^{\top}\mathbf{x}$ and provides a unimodal Gaussian-shaped estimate, not suitable for many types of data.

4.3 Variational inference

We may approximate each555As a reminder, $p_{j}=p(y_{j}|\mathbf{x},y_{1},\ldots,y_{j})$ $p_{j}$ with some other distribution $q(y_{j}|\theta)\approx p_{j}$ as in variational Bayesian methods, which turns inference into an optimization problem

[TABLE]

(minimizing Kullback-Leibler divergence); see, e.g., [1]. We could then sample $y_{j}\sim q(y_{j}|\theta)$ instead. Unlike MCMC methods, this approach does not provide an exact model of $p_{j}$ in the limit (with sufficient samples).

4.4 Density estimation

Let us denote $p_{j}=p(y_{j}|\mathbf{x},y_{1},\ldots,y_{j-1})=p(y_{j}|\mathbf{z})$. We may also model the target density $p(y_{j}|\mathbf{z})=p(y_{j},\mathbf{z})/p(\mathbf{z})$ using a non-parametric method such as a Parzen window (i.e., a kernel density estimate, KDE), where

[TABLE]

for some kernel $K_{h}$ (of bandwidth parameter $h$).

Sampling can be carried out for certain kernels such the Gaussian kernel. Although of course, the general disadvantages of kernel methods apply (quadratic complexity wrt number of examples, and difficulty in incremental modeling).

5 Sequential Monte Carlo Regressor Chains

The methods just described above are adequate to build a probabilistic regressor chain if the resulting approximation of $p_{j}$ is adequate for both sampling (obtaining $y^{(m)}_{j}$) and accurate evaluation (obtaining $w^{(m)}_{j}$). This restricts applicability considerably, and our selection is reduced especially with high-dimensional input observations. For example, many methods may provide an evaluation of a pdf, but not useful sampling.

In this section we build a state-space model for probabilistic regressor chains, separating the sampling and evaluation functions. It is essentially a particle filter (PF, see, e.g., [13]). We make some particular considerations and adaptations for its application in a probabilistic regressor chain.

5.1 The particle filter

To briefly review: a particle filter in a state space model consists of a model

[TABLE]

running over time-steps $j=1,\ldots,L$, encompassing the transition function and observation function, $f$ and $\ell$, respectively666Regarding notation, it is fundamental for those familiar with the literature on particle filters and continuous state-space models, to observe that in this work $y_{j}$ is the state, and $\mathbf{x}_{j}$ is the observation; as in line with the machine learning literature. See [13] for an in-depth introduction and survey.

Under this notation, the vanilla particle filtering method for obtaining a marginal estimation (i.e., the prediction777In machine learning, prediction often refers to an estimation for the current instance, rather than for some future time instance, as often in dynamical systems terminology) for $y_{j}$ can be written as:

[TABLE]

We highlight the strong connection to Eq. (4)–Eq. (6) in Monte Carlo classifier chains. However, in the classification context, the posterior pmf $P$ can be sampled from easily; whereas here, since that might not be the case, we use the auxiliary function $f$ to propose samples at each step in the chain.

There are important differences from typical applications of particle filters, namely 1) in our case the model (Eq. (14)) is learned from the data (i.e., no domain knowledge assumptions), 2) a single observation $\mathbf{x}$ is relevant to an entire sequence of states $y_{1},\ldots,y_{L}$ (that is to say, the isotemporal case), and 3) the full cascade is considered rather than the standard single-order Markovian model.

5.2 Training

There are two functions which we need to learn; as according to Eq. (14): the transition function $f$ and observation function $\ell$. Unlike the vanilla particle filter model described above, we wish to take into account the chain cascade,

[TABLE]

Any suitable method for modeling this density can be considered, as long as we are able to sample from it (the methods given in subsections 4.1–4.4 are applicable here). Except, in this case we may remove $\mathbf{x}$ which simplifies and speeds-up sampling.

For function $\ell$ we do only need to evaluate the function up to some normalizing constant, and thus we have more flexibility in choosing a function to learn. Since we involve observation $\mathbf{x}_{i}$, predictive power is particularly important at this step. Any method giving accurate

[TABLE]

may be considered.

There is no need to impose a particular choice at this stage – we will address options in the experimental investigations in Section 6.

5.3 Inference

Algorithm 1 elaborates our Sequential Monte Carlo (SMC) scheme. It can be specifically designed to measure posterior $p(\mathbf{y}|\mathbf{x})$; namely to approximate the complex integrals involving of type Eq. (3) which represents a MMSE estimator, or other loss functions that we are more interested in such as the mode.

Note that the Effective Sample Size (ESS) approximations is used to decide when the resampling step (and the MCMC/AIS schemes if required) as thus prevent sample degeneration (i.e., error propagation), either

[TABLE]

is typically appropriate [23].

Step 2c in the algorithm is not strictly necessary, nevertheless this step may be useful in cases where the sequential scheme is struggling, and we include a detailed description in A, of the application of $N$ independent MCMC schemes, specifically, $N$ Metropolis-Hastings (MH) methods (at the $j$-th iteration of the SMC algorithm above).

Of course, if a model efficiently meets both constraints for $f$ and $\ell$ (as just mentioned above) we can use the same model for both, thus simplifying Eq. (18) and recovering the approach described in Section 5.

6 Experiments

In this section we compare some of the approaches we discussed, identified, and developed above. Namely, we compare independent regression models (IR), regressor chains with greedy inference (RC) with our Monte Carlo methods, discussed in Section 4 as well as further developed as a Particle Filter in Section 5 (PFC). We compare different base estimators. Recall that PFC takes both a model for sampling and a model for evaluation – not necessarily from the same model class. These models are summarized and denoted as follows:

[TABLE]

such that (in the tables of results, Table 1), IR.B denotes independent Bayesian regression, PF.R/B denotes particle filter chains with a discretized random forest base classifier for sampling, and Bayesian regression for evaluation; and so on. Note that MC and PF are always set to maximize the $0/1$ approximation (by selecting a mode), whereas IR and RC (by default) maximize MSE, i.e., predict the mean. In both cases we consider $M=100$ samples/particles per test example; $\eta=0.1$.

IR and RC are implemented in the well-known Scikit-Learn framework888https://scikit-learn.org. We implemented our novel contributions using the Scikit-Multiflow framework999https://scikit-multiflow.github.io/ [27]; which is based on Scikit-Learn. If not explicitly stated, then default parameters are used.

We carry out an empirical evaluation on synthetic and real-world datasets. The real-world sets are described and referenced in [38]101010Available online: http://mulan.sourceforge.net/datasets-mtr.html; covering a number of real-world applications involving predicting the multi-components of sea-water, residential buildings, concrete pouring, and natural resources; so as to over sufficient variety. The synthetic dataset is described and shown in Figure 6.

We consider the following loss metrics:

[TABLE]

(mean squared error and absolute error, respectively) and denoting

[TABLE]

as an approximation of mean $0/1$ loss (inverse accuracy) which approaches the $0/1$-loss estimate as constant $c$ goes to [math]. This loss is designed to reward models which find a joint mode, which should correspond to a path which is likely to occur in practice (e.g., in the training data).

Averaged results over 10-fold cross validation are provided in Table 1 for the different metrics.

7 Discussion

In this section we discuss empirical results. We look at these results with the goal of drawing conclusions about the behaviour of regressor chains in general, and secondly as justifying both acceptable performance and usability of our proposed methodologies.

The empirical results confirm that greedy regressor chains (RC) shows little to no advantage against independent estimators when a linear base model is used (only a small exception under the Andro dataset, if we are to compare RC.B and IR.B). These findings are completely in line with the analysis so far: classification models involve an inherent non-linearity (such as for example the sigmoid function in logistic regression) which adds predictive power via the chain structure, but this not inherently the case in regression.

We do not need to discuss the power of non-linear modeling for regression, as this is an elementary concept, but it is particularly interesting to observe how well regressor chains with a non-linear base learner (e.g., RC.K) can perform better on average (i.e., in terms of average rank) against its independent counterpart (IR.K). It is in this sense that we begin to find an argument to use regressor chains.

Although regressor chains may be effective, we point out the risk of degeneration of estimates across the chain with poorly regularized base models. Indeed, we found this using standard stochastic gradient descent linear regression, estimates diverged so far (obtaining greater than 1000 MSE) that it there was no point to include them in the table. This effect is well-mentioned in many analyses of classifier chains in the literature under the term of error propagation, but in a multi-output binary classification setting, the posterior $P(y_{j}|\mathbf{x},y_{j-1},\ldots,y_{1})$ such propagation is always constrained between [math] and $1$. However, under RCs the trajectory may become increasingly lost and isolated (from the true path) in $\mathbb{R}^{L}$ space as prediction progresses along the chain. Our choice of Bayesian regression adds some regularization which counters this. Nevertheless it is also one of the issues we took into account with our developments of probabilistic regressor chains.

Results show that acceptable predictive performance can be usually obtained by Monte-Carlo approaches (MC, PF), although on average it does not achieve top results overall. At first glance it seems difficult to justify either of these two methods in either of their configurations, but, we can take a closer look at particular evaluation metrics. We find that when it is important to find modes of the posterior ($0/1$), MC/B and PF./B outperform RC.B about half of the time. We can speculate that this would be similar in a hypothetical comparison of MC./K but the implementation we used (from ScikitLearn) does not offer sampling from this method. Furthermore, such a kernel-based approach has its own disadvantages (quadratic complexity, difficult application to incremental learning, and so on) which are limiting in modern settings of large and dynamic datasets. In addition, we can emphasise the important result on Synth: 1.00 (RC.K) vs 0.51 (PF.R/B) is one of the largest differences in performances obtained (actually 0.50 would be the Bayes optimal result for this data). Finally, and perhaps most importantly: Greedy RC provides almost no form of interpretation.

Results are on the Synth dataset are worth exploring in more detail (shown in Figure 6; results shown in the first row in Table 1; detailed step-by-step illustration of performance in Figure 5). In particular, notice that any path crossing a low density region is unlikely to exist in practice, even though it provides the best estimate under MSE. In any kind of data with a multi-modal density in the output space, the our developed PF approach is highly suited, as it is able to track the path-evolution across different modes along its pass of the chain. Furthermore, it is able to offer interpretation or explainability of results by showing actual paths taken, and the density estimate via the cloud of particles; see, e.g., Figure 5.

It is a pity that the benchmark datasets we used as per related studies, do not seem complex enough to reveal this behaviour, although it is easy to provide examples of real-world possibilities. For example, in the prediction of trajectories for vehicles, it is useful only to estimate paths which do not cross out of transportation axes (for example, a truck does not pass through a river but over one of the bridges that cross it; for example as covered in [25] in the discrete-waypoint scenario). Furthermore, it can be interesting to produce a set of possible trajectories rather than a single estimation. Likewise, in time series forecasting and anomaly detection, we may want to consider different hypotheses, rather than outputting a single estimator of maximum likelihood.

The need for interpretable models is certainly of increasing interest, as machine learning methods are used in more sectors, for example medical and security, where it is often essential to be able to provide detailed explanation of results.

For simplicity, in the Synth data, the observation $x$ does not affect the prediction for $y_{2}$. It means that if that variable was added (in a real-world case) following building and deployment of the model predicting $y_{1}$, it would be more efficient to update the multi-output with a chain than an independent classifier. Of course on this dataset any discussion of efficiency is not warranted, but as the number of labels and instances grows, we can see chain methods operating as a kind of transfer learning, adapting partially pre-trained models by adding links from the outputs of those models. A further discussion is beyond the scope of this work.

8 Related Methods and Additional Considerations

To the best of our knowledge this is the first work treating regression chains in depth, and particularly from a probabilistic point of view with a sequential Monte Carlo approach. Nevertheless, since it is tackling a well-established problem (i.e., multi-output regression) it is equally important to contribute a discussion of related methods in other areas, and additional considerations treated in the classifier-chains literature already, as we do in this section.

8.1 Neural networks

A residual neural network (ResNet, [18]) includes skip layers similar to the way a chain model does, except layer-wise rather than node-wise. It has obtained noteworthy performance on deep learning tasks. The original paper (2015) obtained good performance in image classification, though has also indeed been studied in the context of regression (such as in [21]). Aside from its prevalence in classification, another difference is that in ResNets only the last layer is used for prediction, with other labels simply forming internal nodes. Although, similar multi-label architectures (including, e.g., [5, 31]) have been used which consider each (or at least several nodes/layers) as the targets, i.e., multiple targets. Although this work deals primarily with classification problems, this deep-learning approach tackles well one of the issues faced by regressor chains: needing strong non-linearities to serve as useful representations, as highlighted in our discussion and results. Other neural models have been considered specific to the multi-output regression case, including ensemble settings, e.g., [16].

8.2 Probabilistic models

A closely related approach to the inference in probabilistic regressor chains can be seen as performing a standard regression with noisy inputs [9, 19]. Indeed, suppose that $L=2$, and

[TABLE]

Under the assumption of additive noise, we can write

[TABLE]

where $h_{j}$ are the regression models. For performing a proper inference, all the statistical features of $y_{1}$ should be taken into account (i.e., the uncertainty), hence a regression problem with noisy inputs. Gaussian Processes (GPs) provide a method relevant to this context, e.g., [29, 6], and in particular the idea of warped GPs [37, 22, 7] (a kind of hierarchical GP, thereby giving a kind of ‘depth’). Most of these cases are elaborated where $L=2$. We remark (also with regard to discussion above) that there is of course no clear separating definition between (deep) neural network and probabilistic models.

8.3 State space models

The previous considerations have direct application to the inference and prediction in so-called state-space models. Such models are formed by a transition dynamic equation and an observation equation111111Once more, we have opted for notation where $x$ is the observation,

[TABLE]

where $y_{t}$ is the state at time $t$, and $\epsilon$ is perturbation noise. Note that even if the mappings $h$ and noises $\epsilon$ are known, the inference and prediction in this stochastic system can be interpreted, at each $t$, as a noisy input regression problem since

[TABLE]

showing strong parallels with the noisy regression model above. Furthermore, when the elements of the dynamic system are not deterministically known the problem becomes even more complex. In [8, 2], authors suggest to model $h_{d}$, $h_{o}$ as two GPs.

8.4 Chain order

A natural question which arises is – does the order of the chain affect results? This question has been answered in the classification context (in fact, it does affect results) in various other works, mainly by the use of random ensemble methods (e.g., [36, 40]; each model with a different/random order), hill-climbing approaches (e.g., [12, 34]; the best of a number of trialed orders is taken) and methods based on label dependence (e.g., [26, 30] and references therein). Most of these methods are also directly applicable to the regression case; indeed, the chain-order space and complexity is exactly the same, and thus ensemble and hill-climbing methods need no specific adaptation; label-dependence based methods only need to consider dependence among continuous variables rather than discrete. Thus, there is no reason to suggest why the regression context needs special treatment, and for this reason we do not specifically readdress this consideration in this work.

8.5 Multiple passes

Another aspect to consider is why not pass over the set of variables (i.e., through the chain) multiple times per test instance. This is an interesting proposal and would be an easy extension to our methods: simply iterate over the chain a second time, and plug in training labels $y_{L+1}=y_{1},\ldots,y_{L+L}=y_{L}$ and treat them thenceforth as the target labels of interest. We notice that one can view this as a special case of regressor stacking as described in [3] (in the general case, we simply feed all predictions as inputs to a second model – not necessarily a chain model). On the other hand, time complexity increases significantly with each pass along the chain. One can also see a connection to recurrent neural networks (RNN), as noticed in recent work by [28], unrolled across time.

From the probabilistic perspective, an approach of many passes can be seen as a kind of Gibbs sampling (where the graph is undirected). In fact, an undirected and fully connected network (rather than a ‘chain’) removes the question of label order entirely. This has been developed in the context of conditional dependency networks as presented by [15] (for multi-label classification). We notice that this framework is also applicable to the regression case, whenever sampling from the conditional is possible; namely, (as in [15]) using Gibbs sampling to estimate the marginal mode:

[TABLE]

(having adjusted the equations for continuous target variables) where $M_{0}$ indicates the initial samples during the burn-in time which are not considered. This is relatively straightforward to incorporate in the regression setting. It is, however, worth emphasising that (as found empirically in [34], for example), that the number of samples must be relatively much greater since a burn-in time of $M_{0}>L$ samples is required, and $M\gg L$.

8.6 Other areas of application

Regressor chains can be applied to any problem involving multiple continuous output variables. Time series forecasting is a natural application, e.g., [29]. Previously classifier chains was applied to a discretized version of route prediction in urban traffic modeling [33]; and an application of regressor chains in continuous space would be perhaps even more natural.

One interesting similarity is to the application of Monte Carlo tree search for continuous action spaces in reinforcement learning (a good example of this is given in [42]). Forming a tree to search on top of a continuous space is part of the problem we tackle. We notice that the authors of this cited work use a kernel regression approach which is what we have empirically found to work well in this work on regressor chains. Unlike our study, there is no use of a chain in the sense of a cascade, which is the main focus of this work. Clearly investigation of the use of regressor chains in such related areas could be promising.

9 Conclusions

In this work we have looked at extending the approach of chaining models to the regression context. This was motivated by the fact that, although the idea had attracted initial interest and found potential applicability to multi-output regression tasks, the behaviour of regressor chains was found not to emulate that of classifier chains, thus warranting this in-depth study, to identify, unravel, and overcome the weaknesses of this application.

We identified, discussed, and dealt with several important issues in this regard:

• •

Regressor chains are only useful when non-linearities are used in the base models, and this is equally applicable to the standard multi-output case, and the case where attribute values arrive in sequence

• •

It is difficult to justify regressor chains for MSE estimation, unless paying particular attention to the previous point.

• •

Error propagation (progressive divergence from the true path estimate) can be much more severe in the regression case due to the unbounded output space

• •

Greedy inference in regressor chains, on top of the above points, is unable to provide useful representation or interpretation about the output space.

To our knowledge, this is the first work which looks in-depth at regressor chains in the probabilistic context. We surveyed a number of applicable methods, and developed our own based on sequential Monte Carlo methodology, particularly crafted to tackle these identified issues. We use non-linear base models as particular to the ‘chains’ methodology where base models themselves can be considered as a flexible hyperparameter, guided by aspects of the problem domain, rather than a hardwired design – this, we argue, should not be seen as a ‘nuisance’ parameter but as an attractive feature for adaptation across different areas. We provide a mode-seeking mechanism, rather than only (but as well as, potentially) an estimator for minimum squared error. Error propagation is controlled by a resampling scheme. And interpretation is provided by a tree generated by the point cloud of sample paths, which offer also a description of the underlying density.

Probabilistic regressor chains have peculiarities that distinguish it from other models, such as the full cascade involving all outputs and hyper-parametrization of the base models. However, to properly place regressor chains in context with the wider literature, we also identified and discussed connections to a range of related work which had not been noticed.

We may argue that our analysis not only facilitates understanding the performance of regressor chains that we develop, but also earlier throws more light into the performance considerations of classifier chains, modifications of which are still under active development and recent publication.

Our Monte Carlo methods suggest to be promising especially for tasks where path explainability (i.e., different hypotheses regrading the path taken through output space) is of more importance than outputting the result of a minimum mean squared error estimator. Such tasks include medical research, anomaly detection, de-noising, and missing-value imputation; providing plenty of application lines along which to develop this work further and built additional connections with related areas of the literature. In future work we also intend to explore areas (such as dynamic chains and recurrent models) that are being developed in parallel in the classification context to approach themes relating to chain order and structure.

Appendix A Parallel Metropolis-Hastings (MH) chains

For completeness, we elaborate the algorithm based on [4, 14], starting from the equal weighted set of samples obtained at step 2a of Algorithm 1, $\{\tilde{y}_{j,0}^{(1)},\ldots,\tilde{y}_{j,0}^{(M)}\}$. After $K$ iterations, we obtain the final set of samples, illustrated as follows for the $m$-th particle:

• •

For $k=1,\ldots,K:$

Draw $z\sim q(y|y_{j,k-1}^{(m)})$. 2. 2.

Setting for simplicity $p_{j}(z)=p_{j}(z|\mathbf{x},{\widetilde{y}}_{1:j-1}^{(m)})$, accept the movement $y_{j,k}^{(m)}=z$ with probability

[TABLE]

Otherwise, with probability $1-\alpha$, set $y_{j,k}^{(m)}=y_{j,k-1}^{(m)}$

Therefore final set of samples for an iteration is $\{{\widetilde{y}}_{j}^{(1)},\ldots,{\widetilde{y}}_{j}^{(M)}\}=\{{\widetilde{y}}_{j,K}^{(1)},\ldots,{\widetilde{y}}_{j,K}^{(m)}\}$.