Human-in-the-loop Active Covariance Learning for Improving Prediction in Small Data Sets

TL;DR

This paper introduces a human-in-the-loop method for efficiently eliciting expert knowledge about feature similarities to improve covariance estimation and prediction accuracy in small, high-dimensional datasets.

Contribution

It proposes a novel sequential decision-making approach for eliciting pairwise feature similarities, enhancing covariance learning and predictive performance in small data settings.

Findings

Improved predictive accuracy on simulated data

Enhanced performance on real high-dimensional datasets

Efficient elicitation process reduces expert effort

Abstract

Learning predictive models from small high-dimensional data sets is a key problem in high-dimensional statistics. Expert knowledge elicitation can help, and a strong line of work focuses on directly eliciting informative prior distributions for parameters. This either requires considerable statistical expertise or is laborious, as the emphasis has been on accuracy and not on efficiency of the process. Another line of work queries about importance of features one at a time, assuming them to be independent and hence missing covariance information. In contrast, we propose eliciting expert knowledge about pairwise feature similarities, to borrow statistical strength in the predictions, and using sequential decision making techniques to minimize the effort of the expert. Empirical results demonstrate improvement in predictive performance on both simulated and real data, in high-dimensional linear regression tasks, where we learn the covariance structure with a Gaussian process, based on sequential elicitation.

Tables

Table 1. Table 1: Comparison of the number of required samples and feedbacks to reach a particular MSE level. Values are averages over 50 independent runs. Initial MSE values for Yelp and Amazon are 1.88 1.88 1.88 and 2.26 2.26 2.26 , respectively.

		More Samples		More Feedbacks
	MSE	RND	ACT	SeqElc
Yelp	1.85	5	5	13
	1.80	15	9	60
Amazon	2.20	20	14	55
	2.18	27	17	90

Equations

$$k_{\mathbf{\gamma }}({\mathbf{u}}_i,{\mathbf{u}}_j)=\exp \left(-\frac{{\Vert {\mathbf{u}}_i-{\mathbf{u}}_j\Vert }_{{\mathbf{A}}_{\mathbf{\gamma }}}^2}{2}\right),$$

$$f_{ij}\sim \text{mbox}Bernoulli(q),$$

$$\text{centering}{\mathbb{E}}_{\pi ({\tilde{f}}_{ij}\mid \mathcal{D},\mathcal{F})}\left[\sum _k\text{mbox}KL\left[\pi \left(\tilde{y}\mid {\mathbf{x}}_k,\mathcal{D},\mathcal{F},{\tilde{f}}_{ij}\right)\parallel \pi \left(\tilde{y}\mid {\mathbf{x}}_k,\mathcal{D},\mathcal{F}\right)\right]\right]\text{@add@centering}$$

$$(ij{)}^{\ast }=\arg \underset{(ij)\in R_p}{\max }(\text{mbox}ExpectedInformationGain),$$

$$\mathbf{y}=\mathbf{X\beta }+\mathbf{ϵ},$$

$$\mathbf{\beta }\mid {\sigma }^2,{\tau }^2,\mathbf{\gamma }\sim \mathcal{N}\left(0,{\sigma }^2{\tau }^2\mathbf{K}\right),$$

$${\mathbf{\mu }}^{\ast }$$

$${\mathbf{\Sigma }}^{\ast }$$

$$a^{\ast }$$

$$b^{\ast }$$

$$\pi (\tilde{y}\mid \mathbf{y})$$

$$=\text{mbox}{t}_{2a^{\ast }}\left({\tilde{\mathbf{x}}}^T{\mathbf{\mu }}^{\ast },\frac{b^{\ast }}{a^{\ast }}(1+{\tilde{\mathbf{x}}}^T{\mathbf{\Sigma }}^{\ast }\tilde{\mathbf{x}})\right),$$

$$p(\tilde{y}\mid y)=t_{2a^{\ast }}\left({\tilde{\mathbf{x}}}^T{\mathbf{\mu }}^{\ast },\frac{b^{\ast }}{a^{\ast }}\left(1+{\tilde{\mathbf{x}}}^T{\mathbf{\Sigma }}^{\ast }\tilde{\mathbf{x}}\right)\right)$$

$$p(\tilde{y}\mid y,f)=t_{2a^{\ast }}\left({\tilde{\mathbf{x}}}^T{\mathbf{\mu }}_{new}^{\ast },\frac{b_{new}^{\ast }}{a^{\ast }}\left(1+{\tilde{\mathbf{x}}}^T{\mathbf{\Sigma }}_{new}^{\ast }\tilde{\mathbf{x}}\right)\right).$$

$$\begin{array}{r}p(\tilde{y}\mid y)=\mathcal{N}\left({\tilde{\mathbf{x}}}^T{\mathbf{\mu }}^{\ast },\frac{b^{\ast }}{a^{\ast }-1}\left(1+{\tilde{\mathbf{x}}}^T{\mathbf{\Sigma }}^{\ast }\tilde{\mathbf{x}}\right)\right)\end{array}$$

$$\begin{array}{r}p(\tilde{y}\mid y,f)=\mathcal{N}\left({\tilde{\mathbf{x}}}^T{\mathbf{\mu }}_{new}^{\ast },\frac{b_{new}^{\ast }}{a^{\ast }-1}\left(1+{\tilde{\mathbf{x}}}^T{\mathbf{\Sigma }}_{new}^{\ast }\tilde{\mathbf{x}}\right)\right).\end{array}$$

$$\begin{gathered} \text{mbox}KL\left(p(\tilde{y}\mid y,f)\parallel p(\tilde{y}\mid y)\right)=\log \left(\sqrt{\frac{b^{\ast }\left(1+{\tilde{\mathbf{x}}}^T{\mathbf{\Sigma }}^{\ast }\tilde{\mathbf{x}}\right)}{b_{new}^{\ast }\left(1+{\tilde{\mathbf{x}}}^T{\mathbf{\Sigma }}_{new}^{\ast }\tilde{\mathbf{x}}\right)}}\right)+\frac{(a^{\ast }-1){\left({\tilde{\mathbf{x}}}^T\left({\mathbf{\mu }}_{new}^{\ast }-{\mathbf{\mu }}^{\ast }\right)\right)}^2}{2b^{\ast }\left(1+{\tilde{\mathbf{x}}}^T{\mathbf{\Sigma }}^{\ast }\tilde{\mathbf{x}}\right)} \\ +\frac{b_{new}^{\ast }\left(1+{\tilde{\mathbf{x}}}^T{\mathbf{\Sigma }}_{new}^{\ast }\tilde{\mathbf{x}}\right)}{2b^{\ast }\left(1+{\tilde{\mathbf{x}}}^T{\mathbf{\Sigma }}^{\ast }\tilde{\mathbf{x}}\right)}-\frac{1}{2} \end{gathered}$$

$$\begin{array}{r}p(\mathbf{y})=\int p\left(\mathbf{y}\mid \mathbf{\beta },{\sigma }^2\right)p\left(\mathbf{\beta }\mid {\sigma }^2\right)p({\sigma }^2)d\mathbf{\beta }d{\sigma }^2=\int \mathcal{N}\left(\mathbf{y}\mid \mathbf{X\beta },{\sigma }^2\mathbf{I}\right)\mathcal{N}\left(\mathbf{\beta }\mid 0,{\sigma }^2{\tau }^2\mathbf{K}\right)\text{mbox}IG\left({\sigma }^2\mid a_{\sigma },b_{\sigma }\right)d\mathbf{\beta }d{\sigma }^2\end{array}$$

$$\begin{array}{r}=\int {\left(2\pi {\sigma }^2\right)}^{-\frac{N+M}{2}}\mid {\tau }^2\mathbf{K}{\mid }^{-1/2}\text{mbox}exp\left[-\frac{1}{2{\sigma }^2}\{{\left(\mathbf{y}-\mathbf{X\beta }\right)}^T\left(\mathbf{y}-\mathbf{X\beta }\right)+{\mathbf{\beta }}^T{\left({\tau }^2\mathbf{K}\right)}^{-1}\mathbf{\beta }\}\right]d\mathbf{\beta }\text{mbox}IG\left({\sigma }^2\mid a_{\sigma },b_{\sigma }\right)d{\sigma }^2\end{array}$$

$$\begin{array}{r}=\int {\left(2\pi {\sigma }^2\right)}^{-\frac{N}{2}}\mid \mathbf{I}+\mathbf{X}{\tau }^2{\mathbf{KX}}^T{\mid }^{-1/2}\text{mbox}exp\left[-\frac{1}{2{\sigma }^2}\{{\mathbf{y}}^T{\left(\mathbf{I}+\mathbf{X}{\tau }^2\mathbf{K}{\mathbf{X}}^T\right)}^{-1}\mathbf{y}\}\right]\text{mbox}IG\left({\sigma }^2\mid a_{\sigma },b_{\sigma }\right)d{\sigma }^2\end{array}$$

$$\begin{array}{r}=\int \mathcal{N}\left(\mathbf{y}\mid 0,{\sigma }^2\left(\mathbf{I}+{\tau }^2\mathbf{XK}{\mathbf{X}}^{\mathbf{T}}\right)\right)\text{mbox}IG\left({\sigma }^2\mid a_{\sigma },b_{\sigma }\right)d{\sigma }^2\end{array}$$

$$\begin{array}{r}\int (2\pi {)}^{-\frac{N}{2}}({\sigma }^2{)}^{-\frac{N}{2}}\mid \mathbf{\Sigma }{\mid }^{-1/2}\text{mbox}exp\left(-\frac{1}{2{\sigma }^2}{\mathbf{y}}^T{\mathbf{\Sigma }}^{-1}\mathbf{y}\right)\frac{b_{\sigma }^{a_{\sigma }}}{\Gamma (a_{\sigma })}({\sigma }^2{)}^{-a_{\sigma }-1}\text{mbox}exp(-\frac{b_{\sigma }}{{\sigma }^2})d{\sigma }^2\end{array}$$

$$\begin{array}{r}=(2\pi {)}^{-\frac{N}{2}}\frac{b_{\sigma }^{a_{\sigma }}}{\Gamma (a_{\sigma })}\mid \mathbf{\Sigma }{\mid }^{-1/2}\int ({\sigma }_2{)}^{-\frac{N}{2}-a_{\sigma }-1}\text{mbox}exp\left(-\frac{1}{{\sigma }^2}\left[b_{\sigma }+\frac{1}{2}{\mathbf{y}}^T{\mathbf{\Sigma }}^{-1}\mathbf{y}\right]\right)d{\sigma }^2\end{array}$$

$$\begin{array}{r}=(2\pi {)}^{-\frac{N}{2}}\frac{\Gamma \left(a_{\sigma }+\frac{N}{2}\right)}{\Gamma (a_{\sigma })}{b}_{\sigma }^{a_{\sigma }}\mid \mathbf{\Sigma }{\mid }^{-1/2}{\left(b_{\sigma }+\frac{1}{2}{\mathbf{y}}^T{\mathbf{\Sigma }}^{-1}\mathbf{y}\right)}^{-\left(a_{\sigma }+\frac{N}{2}\right)}\end{array}$$

Full text

1 Introduction

Data sets with a small number of observations $n$ and a large number of variables $p$, a.k.a “small $n$, large $p$”, are common in many applications and pose big challenges for statistical methods Fan and Li (2006). Potential approaches to alleviate these challenges are either to provide more samples, which could be very expensive or even impossible in some applications, or to regularize models using additional sources of information. An example of such an additional source is human expert who may have substantial subject-specific knowledge that can be used to improve data analysis.

In Bayesian statistics, expert-provided information can be formulated as a prior distribution over the parameters to regularize the model. Extracting and incorporating experts’ knowledge into statistical models is known as knowledge elicitation. When designing an elicitation method, two important choices will affect the success. First is what is assumed of the knowledge of the expert: only domain expertise or additional statistical expertise. In Bayesian statistics, most of the existing knowledge elicitation approaches require the domain expert to be expert, or at least knowledgeable, in statistics as well. The second choice is the assumptions made on the content of knowledge that is elicited, which also translates into the structure of the model used for the elicitation. Most existing methods for eliciting importance of variables for prediction tasks assume independence among the parameters to make the elicitation easier, at the cost of ignoring expert knowledge on dependencies.

In this paper, we propose a method for elicitation that avoids both problems: only domain expertise is required, and the underlying assumptions are generalized from independence of the parameters to pairwise dependencies. Technically, we formulate the problem as a sequential decision-making task of choosing which queries to ask to maximally improve predictions with minimal expert effort. The queries will be about pairwise relationships of features, to expand beyond the earlier assumption of their independence.

With experiments on synthetic and real data, we demonstrate that the approach is able to improve the predictive performance by constructing a more informative covariance matrix using expert feedback.

Contribution. In summary, the main contributions of the paper are: (1) A method for learning a prior with a full covariance matrix (without an independence assumption among features) using information provided by a domain expert, (2) an algorithm for actively eliciting expert’s knowledge on pairwise similarities of the features, to maximally improve predictive performance while minimizing the expert’s workload.

2 Related Work

There exists a rich literature on knowledge elicitation for improving the performance of statistical models. A number of studies have proposed to elicit scores or point estimates of an unknown parameter or quantity of interest directly from the expert Garthwaite et al. (2008); Haakma et al. (2011); O’Hagan (1998). These approaches typically either assume that the expert has a high level of statistical training, to be able to assess and formulate prior distributions, or is assisted by a human facilitator in an elicitation session. The main goal is accuracy of the knowledge instead of time-efficiency, which makes the elicitation burdensome. In an attempt to make knowledge elicitation more automatic and resource-efficient, Micallef et al. Micallef et al. (2017) proposed an interactive approach to elicit the expert’s knowledge on the relevance of individual features to improve the predictive accuracy of an underlying model. They used a sequential decision making method, multi-armed bandits, to guide the interaction toward features that would likely benefit from the expert’s input. Daee et al. Daee et al. (2017) formulated knowledge elicitation as a probabilistic inference problem where expert knowledge about individual features is sequentially extracted to improve predictions in a sparse linear regression task. The work was extended by Sundin et al. Sundin et al. (2018) by considering information about the direction of a relevant feature. In all these works, independence among features is assumed in the prior distribution. These approaches could be used to complement our method.

In a similar scenario, Xiao et al. Xiao et al. (2018) proposed a knowledge elicitation approach for extracting expert knowledge on influential relationships between pairs of variables in a directed acyclic graph, i.e., whether a variable $a$ is likely to be up- or down-stream of variable $b$. Their goal, however, is not prediction but to learn the structure of the graph. Afrabandpey et al. Afrabandpey et al. (2017) proposed an interactive method for eliciting knowledge on pairwise feature similarities based on an interactive display and used it to improve predictive performance of linear regression. This work is similar to ours, but differs in two important aspects: First, their method is passive - the queries do not change according to what has been elicited so far - which makes the feedback giving process exhaustive for the expert. Second, the expert feedbacks are only post-hoc connected to the prediction problem through a metric learning process, while our approach adopts expert feedback directly by constructing an informative prior distribution for the underlying model.

Our method builds upon two further works. Krupka and Tishby Krupka and Tishby (2007) studied cases where there are meta-features that capture prior knowledge on features. Using these meta-features, they proposed a method to improve generalization of a classification model by constructing a prior covariance (or a Gaussian process) on the parameters without independence assumptions among the variables. Their approach, however, is not a human-in-the-loop system and does not elicit expert knowledge. Yang et al. Yang et al. (2007) proposed a Bayesian framework for actively learning a similarity metric among samples using pairwise constraints labeled by a human agent. The method aims to improve pairwise distances by querying the unlabeled pairs of samples with the greatest uncertainty in relative distance. In contrast, our approach is designed to optimize the predictive performance of the model using feedback on pairs of features and we query pairs based on their expected improvement on predictive performance.

3 Proposed Approach

We consider the task of learning a probabilistic predictive function $y=h_{\bm{\theta}}(\bm{x})$, parametrized by $\bm{\theta}$, given a data set $\mathcal{D}=\{(y_{i},\bm{x}_{i});i=1,\ldots,n\}$ of $n$ outcomes $y_{i}$ and feature vectors $\bm{x}_{i}=[x_{i1},\ldots,x_{ip}]^{\mathrm{T}}$ pairs. We assume that $n$ is small and cannot be increased, so that it is worthwhile to spend effort on collecting and encoding prior knowledge or other external information. Here, we particularly focus on using knowledge elicitation from domain expert to learn the dependencies between the elements of $\bm{\theta}=[\theta_{1},\ldots,\theta_{L}]^{\mathrm{T}}$, where $L$ is the total number of parameters of the predictive function. This approach is complementary to many of the previous prior and expert knowledge elicitation approaches.

In addition to the data $\mathcal{D}$, we assume the availability of meta-features $\bm{u}_{l}\in\mathbb{R}^{d}$ for each parameter $\theta_{l}$. Our approach to improve the predictive model $h_{\bm{\theta}}$ is then formulated in three components as shown in Figure 1: (1) An auxiliary mapping $g$ from the meta-features $\bm{u}_{l}$ to the parameters $\theta_{l}$, $\theta_{l}=g_{\bm{\gamma}}(\bm{u}_{l})$, with parameters $\bm{\gamma}$. (2) An observation model for expert feedback that encodes the expert knowledge to information about $g$ and consequently to the dependencies between the elements of $\bm{\theta}$. (3) A query algorithm that chooses questions to the expert to optimize the effort of the expert for improving the predictions. Each component is described in more detail in the following sections.

To concretize the approach, consider as an example a sentiment analysis task where the goal is to predict a rating assigned by a reviewer to a product from the text of the review. The training data is the set of all reviews given by different reviewers together with their corresponding ratings. The feature vector $\bm{x}$ describes occurrence of keywords in the review text (bag-of-words). Assume linear regression for prediction, $h_{\bm{\beta}}(\bm{x})=\bm{\beta}^{\mathrm{T}}\bm{x}$, with $\bm{\theta}=\bm{\beta}$ being the regression coefficients. A simple approach to constructing a meta-feature vector $\bm{u}_{l}$ for each keyword $l$ is to take the vector of occurrences of the keyword in the training set reviews. This implies that features that do not occur in the same documents are dissimilar and vice versa. Using a Gaussian process as the auxiliary mapping $g$ then corresponds to learning a full covariance matrix for the prior of the regression coefficients $\bm{\beta}$ based on the meta-features. Knowledge about the pairwise similarity of the role of the keywords in the rating prediction is sequentially queried from the expert. The rationale is that the expert has knowledge about the correlation among regression coefficients of at least some keywords, i.e., either that two keywords positively affect the output (both cause an increase in the rating) or they negatively affect it (decrease the rating). The extracted knowledge will then be incorporated into the model to improve the prior covariance matrix and thus predictions. The case of prediction task with linear regression is described in detail in Section 3.4.

3.1 Modelling Parameter Dependencies

The auxiliary mapping $g$ is used to model the dependencies between the parameters $\bm{\theta}$. We model it with a Gaussian process with a zero mean and covariance function $k_{\bm{\gamma}}$, where $\bm{\gamma}$ is a vector of parameters Krupka and Tishby (2007); Rasmussen and Williams (2006). In particular, this implies that $\bm{\theta}\sim\mathcal{N}(\bm{0},\bm{K})$, where $\bm{K}$ is a covariance matrix (kernel). The element $K_{ij}=k_{\bm{\gamma}}(\bm{u}_{i},\bm{u}_{j})$ describes the prior covariance between the parameters $\theta_{i}$ and $\theta_{j}$ based on the similarity of their meta-features $\bm{u}_{i}$ and $\bm{u}_{j}$ (according to the covariance function).

To illustrate the effect of the auxiliary mapping on learning $\bm{\theta}$, consider a two-dimensional special case, with $\theta_{1}$ and $\theta_{2}$. If the prior correlation implied by $\bm{K}$ is equal to one, then $\theta_{1}=\theta_{2}$; that is, the $g$ maps the two meta-feature vectors $\bm{u}_{1}$ and $\bm{u}_{2}$ to the same point and, effectively, there is only one parameter to be learned from the data set $\mathcal{D}$. If the prior correlation is zero, learning information about one parameter has no effect on the other, and we need to estimate both of the parameters from $\mathcal{D}$. Of course, often we would be somewhere between these two extremes, such that some statistical strength can be shared between learning of the two parameters.

The covariance function $k_{\bm{\gamma}}$ defines many of the properties (smoothness, flexibility, magnitude, etc.) of the auxiliary mapping based on the meta-features $\bm{u}_{i}$ and $\bm{u}_{j}$. Many different covariance functions are possible (see, e.g., Rasmussen and Williams (2006)). We use the Gaussian kernel

[TABLE]

where $\bm{A}_{\bm{\gamma}}\in\mathbb{R}^{d\times d}$ is a diagonal matrix defining the distance metric with elements $\gamma_{l}$, such that $\left\lVert\bm{u}_{i}-\bm{u}_{j}\right\rVert_{\bm{A}_{\bm{\gamma}}}^{2}=\sum_{l}\gamma_{l}(u_{il}-u_{jl})^{2}$. Here, the $\bm{\gamma}$ are inverse length scales: they determine how influential each dimension of the meta-feature is in the distance metric; $\gamma_{i}=0$ means that the $i^{th}$ dimension of the meta-features does not influence the distance (all features are equal w.r.t. this descriptor). Increasing the value from [math] means that it is increasingly important in defining the distance. Inference about $\bm{\gamma}$ is done by defining a half-normal prior distribution $\bm{\gamma}\sim\mathcal{N}^{+}\left(\gamma_{0}\bm{1}_{d\times 1},\delta\bm{I}_{d\times d}\right)$. The $\gamma_{0}$ and $\delta$ are constant values.

The auxiliary mapping $g$ defines a prior on the parameters $\bm{\theta}$, based on the meta features. It would be attractive to learn $g$ from the data $\mathcal{D}$, but that would require more data than available in the ”small $n$, large $p$” case we are targeting. That is why we turn next to learning it based on knowledge elicitation instead.

3.2 Feedback Model

A feedback model is used to incorporate expert’s knowledge into the model. In this work, we consider feedback on the pairwise similarities of the roles of the parameters $\bm{\theta}$ in the prediction. For example, if the $\bm{\theta}$ are linear regression weights, the expert tells the model which pairs of features affect the prediction output similarly (positively and/or negatively). Expert’s feedback will be either “Similar”, “Dissimilar”, or “I don’t know”. The latter will be ignored by the model and has no effect on prior covariance and consequently on the predictions. Similar to Yang et al. (2007), we define the following likelihood for the feedbacks:

[TABLE]

where $f_{ij}$ is the feedback given on the similarity of the features $i$ and $j$ ($f_{ij}=1$ for similar and $f_{ij}=-1$ dissimilar feedback), $q=\frac{1}{1+\mbox{exp}\left(\bm{\gamma}^{T}\bm{w}_{ij}-\xi\right)}$ is the parameter of the Bernoulli distribution, $\bm{w}_{ij}=\left[\bm{u}_{i}-\bm{u}_{j}\right]^{2}$, and $\xi$ is the threshold determining the similarity/dissimilarity of the two features. Since feedbacks are binary, the feedback model is basically defined as a logistic regression of the distances among features $\bm{w}_{ij}$. Two features are more likely to be similar, $f_{ij}=1$, only when their distance is less than the threshold $\xi$ and vice versa. To avoid setting the value of the threshold manually, we define a half-normal prior distribution over it, $\xi\sim\mathcal{N}^{+}\left(m_{\xi},\sigma^{2}_{\xi}\right)$.

3.3 Query Algorithm

Given the set of pairs of features, we would like to query first the similarity of the pair that will give us maximal improvement on predictive performance. Our query algorithm chooses the pair of features that maximizes the expected utility, which is a standard criterion in decision theory. The utility is task dependent - in our prediction task it is natural to choose the utility to be information gain of the predictions. The same rationale was earlier used successfully by Daee et al. (2017), for querying individual variables. The larger the information gain, the greater impact the feedback will have on the posterior predictive distribution. This is defined as the expected Kullback–Leibler divergence (KL) between the posterior predictive distribution before seeing the feedback and after giving the feedback:

[TABLE]

where $\mathcal{F}$ is the set of feedbacks given up to the current step, and the summation goes over the training data. The expectation is over the unknown feedbacks, given everything observed so far. The posterior predictive distribution of the feedback, $\pi\!\left(\tilde{f_{ij}}\mid\mathcal{D},\mathcal{F}\right)$, is obtained using the data and all the feedback given by the expert up to the current iteration.

3.3.1 Computational Complexity

In each iteration, we need to compute the utilities for all pairs of features, except those to which the expert already gave feedback in the previous iterations. The complexity is $O(p^{2})$ which is expensive if $p$ is large. To reduce this complexity, we used the idea presented in Xiong et al. (2015): we construct a smaller pool of the pairs by randomly sampling from the original set of feature pairs, excluding the pairs to which the expert gave feedback previously. We then select the most informative pair from among the pairs in the smaller pool based on their utilities. Denoting the randomly selected pool by $R_{p}$, the objective function for the query algorithm will change to

[TABLE]

where the expected information gain is computed using Equation 3. Although the optimal pair of $R_{p}$ might not necessarily be optimal for the full set, the amount of degradation in the performance is negligible. This is due to the fact that in real data sets with large number of features, there will be significant redundancies among the pairs. The complexity reduces from $O(p^{2})$ to $O(\mid\!R_{p}\!\mid)$, where $\mid\!R_{p}\!\mid$ denotes the cardinality of $R_{p}$.

The above strategy is effective only if the selected pair from $R_{p}$ can be guaranteed to be near optimal with high probability. A pair is near-optimal if it is among the $\epsilon$ top-ranked pairs in the original set, where $\epsilon$ is a small scalar (e.g. 0.001). The near-optimality of the selected pair from the smaller pool can be proven by the following proposition:

Proposition 1. The best pair in $R_{p}$ according to Equation 4 is among the $\epsilon$ top-ranked pairs of the original set with probability $1-(1-\epsilon)^{\mid\!R_{p}\!\mid}$.

Proof. Based on the definition of near-optimal pairs, the probability that a pair does not belong to the near-optimal pairs is $1-\epsilon$. Since the size of the smaller pool is $\mid\!R_{p}\!\mid$, the probability of obtaining a near optimal query is $1-(1-\epsilon)^{\mid\!R_{p}\!\mid}$.

By setting $\mid\!R_{p}\!\mid$ to a reasonable multiple of $p$, e.g. $10p$, we reduce the complexity from $O(p^{2})$ to $O(p)$, while with high probability (probability of $\sim 0.99$ for a moderate size $p=400$), the selected pair from $R_{p}$ will be among the top $0.1$ percent pairs in the original set. To increase computational speed, we implemented computation of the utilities in parallel. This is straightforward since utilities of different pairs are independent.

3.4 Application to Linear Regression

The described approach is general and could, for example, be applied directly to Gaussian process regression where $\bm{\theta}$ would be the predictive function evaluations at the training data points. In this paper, we apply the method to the important case of linear regression with “small $n$, large $p$” data sets, which occurs in a wide array of practical applications. In particular, we take $\bm{\theta}=\bm{\beta}$ to be the regression weights. We assume availability of an expert who has prior knowledge on the pairwise similarities of the role of the features in the prediction task.

The linear regression model is defined as

[TABLE]

where $\bm{X}\in\mathbb{R}^{n\times p}$ is the matrix of predictors, $\bm{y}\in\mathbb{R}^{n\times 1}$ is the vector of all dependent variables, $\epsilon\sim\mathcal{N}(\bm{0},\sigma^{2}\bm{I})$ is the vector of residuals, and $\bm{\beta}\in\mathbb{R}^{p\times 1}$ is the regression coefficient vector to be inferred. With these assumptions, $\bm{y\mid X,\beta},\bm{\sigma}^{2}\sim\mathcal{N}\left(\bm{X\beta},\sigma^{2}\bm{I}\right)$. The goal is to learn the posterior distribution of $\bm{\beta}$ given the training data and the expert feedback. The proposed approach assumes the following prior distribution for $\bm{\beta}$:

[TABLE]

where $\bm{K}$ is the covariance matrix defined by the Gaussian covariance function $k_{\bm{\gamma}}$ and we have also introduced a scalar magnitude parameter $\tau^{2}$. The expert knowledge affects the predictions through its effect on this covariance matrix.

To complete the model formulation, we define inverse gamma distributions for the hyper-parameters $\sigma^{-2}\sim\mbox{Gamma}\left(a_{\sigma},b_{\sigma}\right)$ and $\tau^{-2}\!\sim\!\mbox{Gamma}\left(a_{\tau},b_{\tau}\right)$, where $a_{\sigma}$, $b_{\sigma}$, $a_{\tau}$ and $b_{\tau}$ are constant values. Setting the values of the hyper-parameters is discussed in Section 4. Plate diagram of the model is provided in appendix C.

3.5 Posterior Inference and Prediction

Parameters of the proposed model are $\{\bm{\beta},\sigma^{2},\tau^{2},\bm{\gamma},\xi\}$. Due to the complexity of the model, there is no closed-form posterior distribution. Instead, we obtain maximum a posteriori (MAP) estimates of the parameters by implementing the model in the probabilistic programming language Stan Carpenter et al. (2016) (version 2.17.3, codes in the supplementary material). To reduce computational burden of computing the MAP estimate of $\bm{\beta}$, which requires the inversion of the covariance matrix $\bm{K}$, we marginalized the likelihood over $\bm{\beta}$ and $\sigma^{2}$. This results in a Student’s t-distribution with $2a_{\sigma}$ degrees of freedom for the marginal likelihood, $\bm{y}\sim\mbox{MVSt}_{2a_{\sigma}}\left(0,\frac{b_{\sigma}}{a_{\sigma}}\bm{\Sigma}\right)$, with $\bm{\Sigma}=\bm{I}+\tau^{2}\bm{X}\bm{K}\bm{X}^{T}$. The matrix $\bm{\Sigma}$ is of dimension $n\times n$ which is much smaller than $p\times p$. More details on the marginal likelihood derivation are provided in appendix B. MAP estimate of $\bm{\gamma}$ is used to compute $\bm{K}$. The joint posterior distribution of $\bm{\beta}$ and $\sigma^{2}$ is then obtained as a Normal-Inverse-Gamma distribution, $\mbox{NIG}(\bm{\mu}^{*},\bm{\Sigma}^{*},a^{*},b^{*})$, with parameters

[TABLE]

Finally, predictions for the future input $\tilde{\bm{x}}$ can be done using the posterior predictive distribution:

[TABLE]

where $t_{2a^{*}}$ denotes the univariate Student’s t-distribution with $2a^{*}$ degrees of freedom. Details on computing the KL divergence of Equation 3 for the linear regression model are provided in appendix A.

4 Experimental Results

In this section, the performance of the proposed method is evaluated in several “small $n$, large $p$” regression problems on both simulated and real data with a simulated user. As far as we know, no other method has been proposed for precisely the same task, i.e., improving generalization by constructing a prior distribution using feedback on pairs of features. Therefore, we compare with two natural alternatives:

random query,

• -

non-sequential version of our algorithm, which computes the utilities once, before observing user feedback, and never updates the utilities.

4.1 Synthetic Data

Data generation: We simulate data from a linear regression model, $\bm{y}=\bm{X}\bm{\beta}+\bm{\epsilon}$, with $\bm{\epsilon}\sim\mathcal{N}(\bm{0},\sigma^{2}\bm{I})$ and $\sigma^{2}=5$. The data are generated similarly to Li et al. (2010), with the difference that instead of manually constructing the optimal regression coefficient, $\bm{\beta}_{opt}$, we construct the optimal $\bm{\gamma}$, denoted by $\bm{\gamma}_{opt}$, and the meta-features. They are used to construct the optimal covariance structure, $\bm{K}_{opt}$, using which we sample $\bm{\beta}_{opt}$ from $\mathcal{N}\left(\bm{0},\sigma^{2}\tau^{2}\bm{K}_{opt}\right)$, assuming $\tau^{2}=1$. The idea is to check whether or not the model can learn $\bm{\gamma}_{opt}$ and consequently $\bm{\beta}_{opt}$ using the feedback given by an oracle. We set the number of features to $25$ which are divided into four groups of sizes $10$, $5$, $5$, $5$. Features within each group are highly correlated while features across groups are independent. Meta-features are vectors of size $9$ and dimensions $1$, $2$, $4$, $6$ are assumed to be important for learning the optimum covariance structure, i.e. $\bm{\gamma}_{opt}=\left[1,1,0,1,0,1,0,0,0\right]^{T}$. The optimum covariance matrix constructed using $\bm{\gamma}_{opt}$ and meta-features is shown in Figure 2.a..

Training and test outputs are generated by randomly sampling $\bm{X}_{tr}\in\mathbb{R}^{20\times 25}$ and $\bm{X}_{te}\in\mathbb{R}^{1000\times 25}$ from a standard normal distribution, multiplying them by $\bm{\beta}_{opt}$ and adding normally distributed noise with variance $\sigma^{2}$.

Learning of the model and results: The hyperparameters of the model are $a_{\sigma}=2$ and $b_{\sigma}=7$ for $\sigma^{2}$ and $a_{\tau}=2$ and $b_{\tau}=4$ for $\tau^{2}$, to reflect relatively vague information on the residual variance and scaling parameter, respectively. The prior for $\bm{\gamma}$ is set to $\mathcal{N}^{+}\left(1,0.5\right)$ meaning that a priori we expect almost all dimensions of the meta-features to be equally important, i.e., $\bm{\gamma}_{0}\cong\left[1,1,1,1,1,1,1,1,1\right]^{T}$. In other words, any difference in any dimensions of the meta-feature will affect the (dis)similarity of the features. The hyperparameters of the threshold variable are $\mu_{\xi}=20$ and $\sigma^{2}_{\xi}=10$ to cover wide ranges of values.

The model was first trained to obtain an estimate of the regression coefficients without any feedback ($\bm{\beta}_{0}$). The posterior of the regression coefficients is then used to estimate the utility of each pair (to be used in the next round) and to compute the output of the test data. In each round, the simulated user gives one feedback for a pair queried using three different algorithms, i.e., random, non-sequential (using the utilities computed in round [math]) and sequential (updates utilities of each pair, except those to which the user already gave feedback, using estimated parameters in the previous round). Figures 2.b-d show how $\bm{K}$ evolves using feedback on pairs queried by the sequential model. The total number of possible feedbacks is $\frac{25\times 24}{2}=300$; using only one third of the feedbacks, the sequential model learned almost the optimal covariance structure.

We ran the model $50$ times with different randomly generated data with the setting explained above. The average predictive performances of the model with the different query algorithms are shown in Figure 3. The sequential model obtains the same amount of improvement faster than the random and non-sequential models. Both the random and non-sequential models also perform well, which is due to the simple structure of the data.

4.2 Real Data with Simulated Expert

We tested the proposed method in the task of review ratings prediction from textual reviews, using subsets of Amazon review data and Yelp data. Both data sets contain textual reviews represented as bag-of-words and their corresponding ratings, integer values in the range $1-5$. Features correspond to keywords (unigrams and bigrams). For Amazon, we use the kitchen appliances subset which contains $5149$ reviews Blitzer et al. (2007). For Yelp, we consider $4086$ reviews from the year 2004. In our analysis of both data sets, we removed stop words Salton (1971) and kept features that appeared in at least $100$ reviews. This resulted in data sets containing $340$ features for Amazon and $241$ features for Yelp. From each data set, $1000$ randomly selected reviews were used for training and tesing of the model and the rest as an “expert data set” for constructing the simulated expert. Among the $1000$ reviews for the model, $10\%$ were randomly selected for training and the remaining $90\%$ for testing. This resulted in a training set of size $100\times 340$ for Amazon and of size $100\times 241$ for Yelp, both of which small $n$, large $p$. Training data were normalized to have zero mean and unit variance. As meta-features, we used the transpose of the unnormalized training data. This is the simplest possible meta-feature according to which features (keywords) that are not in the same documents are dissimilar and vice versa. Better meta-features can be defined using word2vec Mikolov et al. (2013) or more advanced representation of the keywords; here we show that the model works fairly well even with such naive meta-features.

4.2.1 Simulated Expert Feedback

To construct the simulated expert, we trained a sparse regression model using LASSO Friedman et al. (2010) on the expert data set. Features with positive regression coefficients are considered pairwise similar. Those with negative coefficients are also similar, while features with zero weight are assumed to be difficult for the simulated expert to evaluate and the feedback will be “I don’t know” for any pair containing at least one of those features. Since our goal is to evaluate the predictive performance of the models and feedback on these features has no effect on models’ performance, we remove all pairs containing any of these features to make the runs faster. For both data sets, there are almost $20$ features in each of the two groups, resulting in $(40\times 39)/2=780$ pairs in total. Using proposition 1 for the sequential model, in each iteration we randomly selected $400$ pairs among which the most informative one will be selected to query to the expert.

The hyperparameters of the model were set to the same values as in the synthetic data, with the only difference that for both data sets we set $\tau^{2}=0.01$, obtained by cross-validation.

4.2.2 Results

Figure 4 shows the predictive performance of query algorithms on each data set. Results are averaged over $50$ runs with randomly selected training, test, and expert data. Query algorithms ask one question at a time and the model was ran for $100$ rounds. For all algorithms and on both data sets, additional knowledge from the simulated expert reduces the prediction error. Yet the amount of reduction depends on the query algorithm. The proposed method in sequential mode obtains fastest improvement. The random model performs better than the non-sequential model. This is expected since the non-sequential model is prone to querying redundant pairs; for instance when there is a cluster of similar features that are highly informative, asking a few of them is sufficient to inform of their grouping. However, the non-sequential model will query all the possible pairs since it does not update the utilities. To compare the effect of learned covariances using expert feedback on regularizing the model, we also compare to the baselines that make a sparsity assumption to regularize the model: LASSO, Ridge and Elastic net. It can be seen that a relatively small number of expert feedback could regularize the model well and achieve good prediction performance compared to the commonly used regularization methods.

4.2.3 Knowledge Elicitation vs. Collecting More Samples

Table 1 compares the number of pairwise feedback required to obtain the same MSE value achieved by adding additional samples (reviews from the expert data set) to the training set. For the latter, we use two strategies: random selection of samples (shown by RND), and active learning strategy (shown by ACT), which selects samples based on maximizing expected information gain (similar to Seeger (2008)). For feedback collection, we adopt our sequential knowledge elicitation approach (shown by SeqElc). All strategies have the same “small $n$, large $p$” setting as starting points with $\bm{X}_{tr}\in\mathbb{R}^{100\times 241}$ for Yelp and $\bm{X}_{tr}\in\mathbb{R}^{100\times 340}$ for Amazon.

The comparison shows the potential of expert knowledge on feature similarities in the prediction when obtaining more samples is impossible or very expensive. According to the table, a particular performance is obtained by a comparable number of expert feedback and additional data. It should be noted that values obtained for the sequential elicitation model are for its best case scenario since we already removed “I don’t know”s which results in a much smaller search pool.

5 Discussion and Conclusion

We proposed a knowledge elicitation approach that incorporates pairwise constraints provided by a human expert for constructing an informative prior distribution. The knowledge elicitation problem is formulated as sequential probabilistic inference which combines expert knowledge with training data. To avoid overwhelming the expert, we proposed a query selection approach that facilitates efficient interaction. To further increase the efficiency of the interaction, we formulated a theoretically motivated subsampling approach. Results for “small $n$, large $p$” problems in simulated and real data with simulated users showed improved prediction accuracy already with a small number of feedbacks.

From our experience, the proposed algorithm works well for moderately high-dimensional problems (e.g., up to 500 predictors). For very high-dimensional problems, with more than thousand predictors, the computations become infeasible due to the large number of pairs. However, the presented approach is general and can be used to complement other existing knowledge elicitation techniques where the user is asked about relevance of features. In this case we can query only pairs of relevant features which significantly improves the feasibility of the model.

Acknowledgments

This work was financially supported by the Academy of Finland (grants 294238, 319264, and 313195), by the Vilho, Yrjö and Kalle Väisälä Foundation of the Finnish Academy of Science and Letters and by Foundation for Aalto University Science and Technology. The authors acknowledged the computational resources provided by the Aalto Science-IT Project.

Appendix A: KL divergence of predictive distributions

In the following formulas, distributions are shown by $p(.)$ in contrast to the main text where we used $\pi(.)$. Let us assume that the predictive distribution before receiving a new feedback $f$ is shown by

[TABLE]

where $\tilde{\bm{x}}$ is the future input with its corresponding output $\tilde{y}$ and $a^{*}$, $b^{*}$, $\bm{\mu}^{*}$ and $\bm{\Sigma}^{*}$ are defined in Section 3.5. The predictive distribution after observing a new feedback $f$ is shown by

[TABLE]

Note that, each feedback updates the point estimates of $\tau^{2}$ and $\bm{K}$ and only those variables depend on these two parameters will be updated, i.e. $b^{*}$, $\bm{\mu}^{*}$, and $\bm{\Sigma}^{*}$. The updated values are shown by $b_{new}^{*}$, $\bm{\mu}_{new}^{*}$, and $\bm{\Sigma}_{new}^{*}$, respectively (Section 3.5).

For the sake of computational simplicity, we relax the above student’s t-distributions to univariate Gaussian distributions. This is possible particularly due the large degree of freedoms of the two distributions (note that $a^{*}=a_{\sigma}+\frac{N}{2}$ where $N$ is the number of samples). Approximations are as follows:

[TABLE]

With these approximations, the KL-divergence of the two Gaussian distributions are obtained as

[TABLE]

Note that since we are using the training data to compute the utilities of the queries, in the implementation, each training data is used as the unobserved input $\tilde{\bm{x}}$ and the KL-divergence is computed as the summation over all the training data.

Appendix B: Marginal likelihood derivation

[TABLE]

To make the following equations shorter, we define $\bm{\Sigma}=\bm{I}+\tau^{2}\bm{XK}\bm{X}^{T}$, then we can write:

[TABLE]

where $\Gamma(.)$ and $\mid.\mid$ represent the Gamma function and the determinant, respectively. The last term is the form of a multivariate student’s t-distribution with $2a_{\sigma}$ degrees of freedom and parameters $\bm{0}$ and $\frac{b_{\sigma}}{a_{\sigma}}\bm{\Sigma}$.

Appendix C: Plate Diagram of the Linear Regression Model

Plate diagram of the prediction model (left) and feedback model (right). Feedbacks $f_{l}$ are sequentially queried from the expert where the total number of feedbacks is equal to the number of feature pairs.

Stan Codes