\documentstyle[jair,twoside,11pt,theapa]{article}
\input epsf %% at top of file
\input psfig

\newenvironment{defs}%
{\begin{list}{}{\setlength{\itemindent}{-\leftmargin}\setlength{\itemsep}{0pt}\setlength{\parsep}{0pt}}}%
{\end{list}}

\jairheading{4}{1996}{129-145}{11/95}{3/96}
\ShortHeadings{Active Learning with Statistical Models}
{Cohn, Ghahramani \& Jordan}
\firstpageno{129}

\begin{document}

\title{Active Learning with Statistical Models}
\author{\name David A.~Cohn \email cohn@harlequin.com \\
	\name Zoubin Ghahramani \email zoubin@cs.toronto.edu \\
	\name Michael I.~Jordan \email jordan@psyche.mit.edu \\
       \addr Center for Biological and Computational Learning \\
	Dept.~of Brain and Cognitive Sciences \\
	Massachusetts Institute of Technology \\
	Cambridge, MA 02139 USA}

\maketitle

\begin{abstract}
For many types of machine learning algorithms, one can compute the
statistically ``optimal'' way to select training data.  In this paper,
we review how optimal data selection techniques have been used with
feedforward neural networks.  We then show how the same principles may
be used to select data for two alternative, statistically-based
learning architectures: mixtures of Gaussians and locally weighted
regression.  While the techniques for neural networks are
computationally expensive and approximate, the techniques for mixtures
of Gaussians and locally weighted regression are both efficient and
accurate.  Empirically, we observe that the optimality criterion
sharply decreases the number of training examples the learner needs in
order to achieve good performance.
\end{abstract}

\section{Introduction}
The goal of machine learning is to create systems that can improve
their performance at some task as they acquire experience or data. In
many natural learning tasks, this experience or data is gained
interactively, by taking actions, making queries, or doing
experiments. Most machine learning research, however, treats the
learner as a passive recipient of data to be processed. This
``passive'' approach ignores the fact that, in many situations, the
learner's most powerful tool is its ability to act, to gather data,
and to influence the world it is trying to understand.  Active
learning is the study of how to use this ability effectively.

Formally, active learning studies the closed-loop phenomenon of a
learner selecting actions or making queries that influence what data
are added to its training set.  Examples include selecting joint
angles or torques to learn the kinematics or dynamics of a robot arm,
selecting locations for sensor measurements to identify and locate
buried hazardous wastes, or querying a human expert to classify an
unknown word in a natural language understanding problem.

When actions/queries are selected properly, the data requirements for
some problems decrease drastically, and some NP-complete learning
problems become polynomial in computation time
\cite{angluin88,baum-lang91a}.  In practice, active learning offers
its greatest rewards in situations where data are expensive or
difficult to obtain, or when the environment is complex or dangerous.
In industrial settings each training point may take days to gather and
cost thousands of dollars; a method for optimally selecting these
points could offer enormous savings in time and money.

There are a number of different goals which one may wish to achieve
using active learning. One is optimization, where the learner performs
experiments to find a set of inputs that maximize some response
variable. An example of the optimization problem would be finding the
operating parameters that maximize the output of a steel mill or candy
factory. There is an extensive literature on optimization, examining
both cases where the learner has some prior knowledge of the
parameterized functional form and cases where the learner has no such
knowledge; the latter case is generally of greater interest to machine
learning practitioners. The favored technique for this kind of
optimization is usually a form of response surface methodology
\cite{box-draper87}, which performs experiments that guide
hill-climbing through the input space.

A related problem exists in the field of adaptive control, where one
must learn a control policy by taking actions. In control problems,
one faces the complication that the value of a specific action may not
be known until many time steps after it is taken.  Also, in control
(as in optimization), one is usually concerned with the performing
well {\em during} the learning task and must trade of exploitation of
the current policy for exploration which may improve it. The subfield
of dual control \cite{feldbaum65} is specifically concerned with
finding an optimal balance of exploration and control while learning.

In this paper, we will restrict ourselves to examining the problem of
supervised learning: based on a set of potentially noisy training
examples ${\cal D} = \{(x_i, y_i)\}_{i=1}^m$, where $x_i \in X$ and
$y_i \in Y$, we wish to learn a general mapping $X \rightarrow Y$.  In
robot control, the mapping may be $state \times action \rightarrow
new~state$; in hazard location it may be $sensor~reading \rightarrow
target~position$.  In contrast to the goals of optimization and
control, the goal of supervised learning is to be able to efficiently
and accurately predict $y$ for a given $x$.

In active learning situations, the learner itself is responsible for
acquiring the training set. Here, we assume it can iteratively select
a new input $\tilde{x}$ (possibly from a constrained set), observe the
resulting output $\tilde{y}$, and incorporate the new example
$(\tilde{x},\tilde{y})$ into its training set. This contrasts with
related work by Plutowski and White \citeyear{pluto-white93}, which is
concerned with filtering an existing data set.  In our case,
$\tilde{x}$ may be thought of as a query, experiment, or action,
depending on the research field and problem domain. The question we
will be concerned with is how to choose which $\tilde{x}$ to try next.

There are many heuristics for choosing $\tilde{x}$, including choosing
places where we don't have data \cite{whitehead91}, where we perform
poorly \cite{linden-webber93}, where we have low confidence
\cite{thrun-moeller92}, where we expect it to change our model
\cite{cal90,cal94}, and where we previously found data that resulted
in learning \cite{schmid-storck93}.  In this paper we will consider
how one may select $\tilde{x}$ in a statistically ``optimal'' manner
for some classes of machine learning algorithms. We first briefly
review how the statistical approach can be applied to neural networks,
as described in earlier work \cite{mackay92,cohn94a}. Then, in
Sections~\ref{mixture_section} and \ref{lwr_section} we consider two
alternative, statistically-based learning architectures: mixtures of
Gaussians and locally weighted regression.
Section~\ref{experiment_section} presents the empirical results of
applying statistically-based active learning to these architectures.
While optimal data selection for a neural network is computationally
expensive and approximate, we find that optimal data selection for the
two statistical models is efficient and accurate.

\section{Active Learning -- A Statistical Approach}
\label{stats_section}
We begin by defining $P(x,y)$ to be the unknown joint distribution
over $x$ and $y$, and $P(x)$ to be the known marginal distribution of
$x$ (commonly called the {\em input distribution}).  We denote the
learner's output on input $x$, given training set $\cal D$ as
$\hat{y}(x;{\cal D})$.\footnote{We present our equations in the
univariate setting. All results in the paper apply equally to the
multivariate case.} We can then write the expected error of the
learner as follows:
\begin{eqnarray}
\int_x E_T\left[\left(\hat{y}(x;{\cal D}) - y(x)\right)^2 | x \right]
	P(x) dx, \label{full-error-eqn}
\end{eqnarray}
where $E_T[\cdot]$ denotes expectation over $P(y|x)$ and over training
sets ${\cal D}$.  The expectation inside the integral may be 
decomposed as follows \cite{gbd92}:
\begin{eqnarray}
E_T\left[\left(\hat{y}(x;{\cal D}) - y(x)\right)^2 | x \right] & = &
E\left[\left(y(x) - E[y|x]\right)^2 \right]
	\label{decomp-eqn}\\
&& + \left(E_{\cal D}\left[\hat{y}(x;{\cal D})\right] - E[y|x]\right)^2
	\nonumber\\
&& + E_{\cal D}\left[\left(\hat{y}(x;{\cal D})
		- E_{\cal D}[\hat{y}(x;{\cal D})]\right)^2 \right] 
	\nonumber
\end{eqnarray}
where $E_{\cal D}[\cdot]$ denotes the expectation over training sets
${\cal D}$ and the remaining expectations on the right-hand side are
expectations with respect to the conditional density $P(y|x)$.  It is
important to remember here that in the case of active learning, the
distribution of $\cal D$ may differ substantially from the joint
distribution $P(x,y)$.

The first term in Equation~\ref{decomp-eqn} is the variance of $y$
given $x$ --- it is the {\em noise} in the distribution, and does not
depend on the learner or on the training data.  The second term is the
learner's {\em squared bias}, and the third is its {\em variance};
these last two terms comprise the mean squared error of the learner
with respect to the regression function $E[y|x]$.  When the second
term of Equation~\ref{decomp-eqn} is zero, we say that the learner is
{\em unbiased}. We shall assume that the learners considered in this
paper are approximately unbiased; that is, that their squared bias is
negligible when compared with their overall mean squared error.  Thus
we focus on algorithms that minimize the learner's error by minimizing
its variance:
\begin{eqnarray}
\sigma^2_{\hat{y}} \equiv \sigma^2_{\hat{y}}(x) = E_{\cal D}\left[\left(\hat{y}(x;{\cal D})
		- E_{\cal D}[\hat{y}(x;{\cal D})]\right)^2 \right].
	\label{passive-var-eqn}
\end{eqnarray}
(For readability, we will drop the explicit dependence on $x$ and
${\cal D}$ --- unless denoted otherwise, $\hat{y}$ and
$\sigma^2_{\hat{y}}$ are functions of $x$ and $\cal D$.)  In an active
learning setting, we will have chosen the $x$-component of our
training set $\cal D$; we indicate this by rewriting
Equation~\ref{passive-var-eqn} as
\begin{eqnarray*}
\sigma^2_{\hat{y}} = \left<\left(\hat{y} - \left<\hat{y}\right>\right)^2 
	\right>,
\end{eqnarray*}
where $\left<\cdot\right>$ denotes $E_{\cal D}[\cdot]$ given a fixed
$x$-component of $\cal D$. When a new input $\tilde{x}$ is selected
and queried, and the resulting $(\tilde{x}, \tilde{y})$ added to the
training set, $\sigma^2_{\hat{y}}$ should change. We will denote the
expectation (over values of $\tilde{y}$) of the learner's new variance
as
\begin{eqnarray}
\left<\tilde{\sigma}^2_{\hat{y}}\right> =
	 E_{{\cal D} \cup (\tilde{x},
		 \tilde{y})}\left[\sigma^2_{\hat{y}} | \tilde{x}\right].
	\label{new-var-eqn}
\end{eqnarray}

\subsection{Selecting Data to Minimize Learner Variance}
In this paper we consider algorithms for active learning which select
data in an attempt to minimize the value of
Equation~\ref{new-var-eqn}, integrated over $X$. Intuitively, the
minimization proceeds as follows: we assume that we have an estimate
of $\sigma^2_{\hat{y}}$, the variance of the learner at $x$. If, for
some new input $\tilde{x}$, we knew the conditional distribution
$P(\tilde{y} | \tilde{x})$, we could compute an estimate of the
learner's new variance at $x$ given an additional example at
$\tilde{x}$. While the true distribution $P(\tilde{y} | \tilde{x})$ is
unknown, many learning architectures let us approximate it by giving
us estimates of its mean and variance. Using the estimated
distribution of $\tilde{y}$, we can estimate
$\left<\tilde{\sigma}^2_{\hat{y}}\right>$, the expected variance of
the learner after querying at $\tilde{x}$.

Given the estimate of $\left<\tilde{\sigma}^2_{\hat{y}}\right>$, which
applies to a given $x$ and a given query $\tilde{x}$, we must
integrate $x$ over the input distribution to compute the integrated
average variance of the learner.  In practice, we will compute a Monte
Carlo approximation of this integral, evaluating
$\left<\tilde{\sigma}^2_{\hat{y}}\right>$ at a number of {\em
reference points} drawn according to $P(x)$. By querying an
$\tilde{x}$ that minimizes the average expected variance over the
reference points, we have a solid statistical basis for choosing new
examples.

\subsection{Example: Active Learning with a Neural Network}
In this section we review the use of techniques from Optimal
Experiment Design (OED) to minimize the estimated variance of a neural
network \cite{fedorov72,mackay92,cohn94a}.  We will assume we have
been given a learner $\hat{y} = f_{\hat{w}}()$, a training set ${\cal
D} = \{(x_i, y_i)\}_{i=1}^m$ and a parameter vector estimate $\hat{w}$
that maximizes some likelihood measure given ${\cal D}$. If, for
example, one assumes that the data were produced by a process whose
structure matches that of the network, and that noise in the process
outputs is normal and independently identically distributed, then the
negative log likelihood of $\hat{w}$ given ${\cal D}$ is proportional
to
\[
S^2 = \frac{1}{m}\sum_{i=1}^m \left(y_i - \hat{y}(x_i)\right)^2.
\]
The maximum likelihood estimate for $\hat{w}$ is that which minimizes
$S^2$.

The estimated output variance of the network is
\[
\sigma^2_{\hat{y}} \approx S^2
\left( \frac{\partial \hat{y}(x)}{\partial w} \right)^T
\left( \frac{\partial^2 S^2}{\partial w^2} \right)^{-1}
\left( \frac{\partial \hat{y}(x)}{\partial w} \right),
	\mbox{  \cite{mackay92}}
\]
where the true variance is approximated by a second-order Taylor
series expansion around $S^2$. This estimate makes the assumption that
$\partial \hat{y}/\partial w$ is locally linear.  Combined with the
assumption that $P(y|x)$ is Gaussian with constant variance for all
$x$, one can derive a closed form expression for
$\left<\tilde{\sigma}^2_{\hat{y}}\right>$. See Cohn \citeyear{cohn94a}
for details.

In practice, $\partial \hat{y}/\partial w$ may be highly nonlinear,
and $P(y|x)$ may be far from Gaussian; in spite of this, empirical
results show that it works well on some problems \cite{cohn94a}.  It
has the advantage of being grounded in statistics, and is optimal
given the assumptions.  Furthermore, the expectation is differentiable
with respect to $\tilde{x}$. As such, it is applicable in continuous
domains with continuous action spaces, and allows hillclimbing to find
the $\tilde{x}$ that minimizes the expected model variance.

For neural networks, however, this approach has many disadvantages.
In addition to relying on simplifications and assumptions which hold
only approximately, the process is computationally expensive.
Computing the variance estimate requires inversion of a $|w| \times
|w|$ matrix for each new example, and incorporating new examples into
the network requires expensive retraining. Paass and Kindermann
\citeyear{paass-kindermann95} discuss a Markov-chain based sampling
approach which addresses some of these problems. In the rest of this
paper, we consider two ``non-neural'' machine learning architectures
that are much more amenable to optimal data selection.

\section{Mixtures of Gaussians}
\label{mixture_section}
The mixture of Gaussians model is a powerful estimation and prediction
technique with roots in the statistics literature
\cite{titterington85}; it has, over the last few years, been adopted
by researchers in machine learning
\shortcite{csktfs88,nowlan91,specht91,ghahramani-jordan94}.  The model
assumes that the data are produced by a mixture of $N$ multivariate
Gaussians $g_i$, for $i = 1,...,N$ (see Figure~\ref{mixture_fig}).

In the context of learning from random examples, one begins by
producing a joint density estimate over the input/output space $X
\times Y$ based on the training set $\cal D$.  The EM algorithm
\cite{dlr77} can be used to efficiently find a locally optimal fit of
the Gaussians to the data. It is then straightforward to compute
$\hat{y}$ given $x$ by conditioning the joint distribution on $x$ and
taking the expected value.
\begin{figure}[hbt]
\centerline{\psfig{figure=mixture_fig.ps,height=2in}}
\vspace{-0.2in}
\caption{\label{mixture_fig}
Using a mixture of Gaussians to compute $\hat{y}$. The Gaussians model
the data density. Predictions are made by mixing the conditional
expectations of each Gaussian given the input $x$.}
\end{figure}

One benefit of learning with a mixture of Gaussians is that there is
no fixed distinction between inputs and outputs --- one may specify
any subset of the input-output dimensions, and compute expectations on
the remaining dimensions. If one has learned a forward model of the
dynamics of a robot arm, for example, conditioning on the outputs
automatically gives a model of the arm's inverse dynamics.  With the
mixture model, it is also straightforward to compute the mode of the
output, rather than its mean, which obviates many of the problems of
learning direct inverse models \cite{ghahramani-jordan94}.

For each Gaussian $g_i$ we will denote the input/output means as
$\mu_{x,i}$ and $\mu_{y,i}$ and variances and covariances as
$\sigma^2_{x,i}$, $\sigma^2_{y,i}$ and $\sigma_{xy,i}$ respectively.
We can then express the probability of point $(x,y)$, given $g_i$ as
\begin{eqnarray}
P(x,y | i) = \frac{1}{2\pi\sqrt{|\Sigma_i|}}
	\exp\left[-\frac{1}{2}({\bf x} - \mu_i)^T\Sigma_i^{-1}({\bf x} - \mu_i)
	\right]
	\label{joint-eqn}
\end{eqnarray}
where we have defined
\begin{eqnarray*}
\begin{array}{lcr}
{\bf x} = \left[\begin{array}{c} x \\ y \end{array}\right] &
{\bf \mu}_i = \left[\begin{array}{c} \mu_{x,i} \\ \mu_{y,i} \end{array}\right]
&
{\bf \Sigma}_i = \left[\begin{array}{ll} \sigma^2_{x,i} & \sigma_{xy,i} \\
			 \sigma_{xy,i} & \sigma^2_{y,i} \end{array}\right].
\end{array}
\end{eqnarray*}
In practice, the true means and variances will be unknown, but can be
estimated from data via the EM algorithm.  The (estimated) conditional
variance of $y$ given $x$ is then
\[\sigma^2_{y|x,i} = \sigma^2_{y,i} -
\frac{\sigma^2_{xy,i}}{\sigma^2_{x,i}}.
\]
and the conditional expectation $\hat{y}_i$ and variance
$\sigma^2_{\hat{y},i}$ given $x$ are:
\begin{eqnarray}
\hat{y}_i = \mu_{y,i} + \frac{\sigma_{xy,i}}{\sigma^2_{x,i}}(x - \mu_{x,i}),& &
\sigma^2_{\hat{y},i} = \frac{\sigma^2_{y|x,i}}{n_i}
	\left(1 + \frac{(x - \mu_{x,i})^2}{\sigma^2_{x,i}}\right).
	\label{mixture_est_eqn}
\end{eqnarray}
Here, $n_i$ is the amount of ``support'' for the Gaussian $g_i$ in the 
training data. It can be computed as
\[
n_i = \sum^m_{j=1} \frac{P(x_j, y_j | i)}{\sum_{k=1}^N P(x_j, y_j | k)}.
\]

The expectations and variances in Equation~\ref{mixture_est_eqn} are
mixed according to the probability that $g_i$ has of being responsible
for $x$, prior to observing $y$:
\[
h_i \equiv h_i(x) = \frac{P(x | i)}{\sum_{j=1}^N P(x | j)},
\]
where
\begin{eqnarray}
P(x | i) = \frac{1}{\sqrt{2 \pi \sigma^2_{x,i}}}
	\exp\left[-\frac{(x - \mu_{x,i})^2}{2\sigma^2_{x,i}}\right].
	\label{marginal-eqn}
\end{eqnarray}

For input $x$ then, the conditional expectation
$\hat{y}$ of the resulting mixture and its variance may be written:
\begin{eqnarray*}
\hat{y} = \sum_{i=1}^N h_i \: \hat{y}_i, & &
\sigma^2_{\hat{y}} = \sum_{i=1}^N \frac{h^2_i\sigma^2_{y|x,i}}{n_i}
	\left(1 + \frac{(x - \mu_{x,i})^2}{\sigma^2_{x,i}}\right),
\end{eqnarray*}
where we have assumed that the $\hat{y}_i$ are independent in
calculating $\sigma^2_{\hat{y}}$.  Both of these terms can be computed
efficiently in closed form. It is also worth noting that
$\sigma^2_{\hat{y}}$ is only one of many variance measures we might be
interested in. If, for example, our mapping is stochastically
multivalued (that is, if the Gaussians overlapped significantly in the
$x$ dimension), we may wish our prediction $\hat{y}$ to reflect the
most likely $y$ value.  In this case, $\hat{y}$ would be the mode, and
a preferable measure of uncertainty would be the (unmixed) variance of
the individual Gaussians.

\subsection{Active Learning with a Mixture of Gaussians}
In the context of active learning, we are assuming that the input
distribution $P(x)$ is known. With a mixture of Gaussians, one
interpretation of this assumption is that we know $\mu_{x,i}$ and
$\sigma^2_{x,i}$ for each Gaussian. In that case, our application of
EM will estimate only $\mu_{y,i}$, $\sigma^2_{y,i}$, and
$\sigma_{xy,i}$.

Generally however, knowing the input distribution will not correspond
to knowing the actual $\mu_{x,i}$ and $\sigma^2_{x,i}$ for each
Gaussian. We may simply know, for example, that $P(x)$ is uniform, or
can be approximated by some set of sampled inputs. In such cases, we
must use EM to estimate $\mu_{x,i}$ and $\sigma^2_{x,i}$ in addition
to the parameters involving $y$. If we simply estimate these values
from the training data, though, we will be estimating the joint
distribution of $P(\tilde{x},y|i)$ instead of $P(x,y|i)$. To obtain a
proper estimate, we must correct Equation~\ref{joint-eqn} as
follows:
\begin{eqnarray}
P(x,y|i) = P(\tilde{x},y|i)\frac{P(x|i)}{P(\tilde{x}|i)}. \label{correction-eqn}
\end{eqnarray}
Here, $P(\tilde{x}|i)$ is computed by applying
Equation~\ref{marginal-eqn} given the mean and $x$ variance of the
training data, and $P(x|i)$ is computed by applying the same equation
using the mean and $x$ variance of a set of reference data drawn
according to $P(x)$.

\medskip
If our goal in active learning is to minimize variance, we should
select training examples $\tilde{x}$ to minimize
$\left<\tilde{\sigma}^2_{\hat{y}}\right>$. With a mixture of
Gaussians, we can compute $\left<\tilde{\sigma}^2_{\hat{y}}\right>$
efficiently. The model's estimated distribution of $\tilde{y}$ given
$\tilde{x}$ is explicit:

\[
P(\tilde{y}|\tilde{x}) = 
	\sum_{i=1}^N \tilde{h}_iP(\tilde{y}|\tilde{x}, i)
	= \sum_{i=1}^N \tilde{h}_i 
	{\cal N}(\hat{y}_i(\tilde{x}), \sigma^2_{y|x,i}(\tilde{x})),
\]
where $\tilde{h}_i \equiv h_i(\tilde{x})$, and ${\cal N}(\mu,
\sigma^2)$ denotes the normal distribution with mean $\mu$ and variance
$\sigma^2$.  Given this, we can model the change in each $g_i$
separately, calculating its expected variance given a new point
sampled from $P(\tilde{y}|\tilde{x}, i)$ and weight this change by
$\tilde{h}_i$.  The new expectations combine to form the learner's new
expected variance
\begin{eqnarray}
\left<\tilde{\sigma}^2_{\hat{y}}\right> & = & \sum_{i=1}^N 
	 \frac{h^2_i\left<\tilde{\sigma}^2_{y|x,i}\right>}{n_i + \tilde{h}_i}
   \left(1 + \frac{(x - \mu_{x,i})^2}{\sigma^2_{x,i}}\right)
	\label{mix-exp}
\end{eqnarray}
where the expectation can be computed exactly in closed form:
\begin{eqnarray*}
\left<\tilde{\sigma}^2_{y,i}\right> =
	\frac{n_i\sigma^2_{y,i}}{n_i+\tilde{h}_i}
	+ \frac{n_i\tilde{h}_i  \left(\sigma^2_{y|\tilde{x},i} 
		+ (\hat{y}_i(\tilde{x}) - \mu_{y,i})^2\right)}
	     {(n_i+\tilde{h}_i)^2},
&&
\left<\tilde{\sigma}^2_{y|x,i}\right> =
   \left<\tilde{\sigma}^2_{y,i}\right> - 
	\frac{\left<\tilde{\sigma}^2_{xy,i}\right>}
	     {\sigma^2_{x,i}},
\\
\left<\tilde{\sigma}^{\mbox{ }}_{xy,i}\right> = 
	\frac{n_i\sigma_{xy,i}}
	     {n_i+\tilde{h}_i}
	+  \frac{n_i\tilde{h}_i  (\tilde{x} - \mu_{x,i})
			     (\hat{y}_i(\tilde{x}) - \mu_{y,i})}
			       {(n_i+\tilde{h}_i)^2},
&&
\left<\tilde{\sigma}^2_{xy,i}\right>  = 
   \left<\tilde{\sigma}_{xy,i}\right>^2
    + \frac{n_i^2\tilde{h}_i^2\sigma^2_{y|\tilde{x},i}
	(\tilde{x} - \mu_{x,i})^2}
	     {(n_i+\tilde{h}_i)^4}.
\end{eqnarray*}
If, as discussed earlier, we are also estimating $\mu_{x,i}$ and
$\sigma^2_{x,i}$, we must take into account the effect of the new
example on those estimates, and must replace $\mu_{x,i}$ and
$\sigma^2_{x,i}$ in the above equations with
\begin{eqnarray*}
\tilde{\mu}_{x,i} = \frac{n_i\mu_{x,i} + \tilde{h}_i\tilde{x}}
                {n_i + \tilde{h}_i}, & &
\tilde{\sigma}^2_{x,i} = \frac{n_i\sigma^2_{x,i}}
                            {n_i+\tilde{h}_i}
                        + \frac{n_i\tilde{h}_i  (\tilde{x} - \mu_{x,i})^2}
                               {(n_i+\tilde{h}_i)^2}.
\end{eqnarray*}

We can use Equation~\ref{mix-exp} to guide active learning. By
evaluating the expected new variance over a reference set given
candidate $\tilde{x}$, we can select the $\tilde{x}$ giving the lowest
expected model variance. Note that in high-dimensional spaces, it may
be necessary to evaluate an excessive number of candidate points to
get good coverage of the potential query space. In these cases, it is
more efficient to differentiate Equation~\ref{mix-exp} and hillclimb
on $\partial\left<\tilde{\sigma}^2_{\hat{y}}\right>/\partial
\tilde{x}$ to find a locally maximal $\tilde{x}$. See, for example,
\cite{cohn94a}.

\section{Locally Weighted Regression}
\label{lwr_section}
Model-based methods, such as neural networks and the mixture of
Gaussians, use the data to build a parameterized model. After
training, the model is used for predictions and the data are generally
discarded. In contrast, ``memory-based'' methods are non-parametric
approaches that explicitly retain the training data, and use it each
time a prediction needs to be made. Locally weighted regression (LWR)
is a memory-based method that performs a regression around a point of
interest using only training data that are ``local'' to that point.
One recent study demonstrated that LWR was suitable for real-time
control by constructing an LWR-based system that learned a difficult
juggling task \cite{schaal-atkeson94}.

\begin{figure}[bht]
\centerline{\psfig{figure=lwr_fig.ps,height=2in}}
\vspace{-0.2in}
\caption{\label{lwr_fig}
In locally weighted regression, points are weighted by proximity to
the current $x$ in question using a kernel. A regression is then
computed using the weighted points.}
\end{figure}

We consider here a form of locally weighted regression that is a
variant of the LOESS model \cite{cdg88}. The LOESS model performs a
linear regression on points in the data set, weighted by a kernel
centered at $x$ (see Figure~\ref{lwr_fig}). The kernel shape is a
design parameter for which there are many possible choices: the
original LOESS model uses a ``tricubic'' kernel; in our experiments we
have used a Gaussian
\[
h_i(x) \equiv h(x - x_i) = \exp(-k(x - x_i)^2),
\]
where $k$ is a smoothing parameter. In Section~\ref{setting-k} we will
describe several methods for automatically setting $k$.

For brevity, we will drop the argument $x$ for $h_i(x)$, and define $n
= \sum_i h_i$. We can then write the estimated means and covariances
as:
\begin{eqnarray*}
& &\mu_x = \frac{\sum_i h_i x_i}{n}, \mbox{ } 
\sigma^2_{x} = \frac{\sum_i h_i  (x_i - \mu_x)^2}{n}, \mbox{ } 
\sigma_{xy} = \frac{\sum_i h_i  (x_i - \mu_x)(y_i - \mu_y)}{n}  \\
& &\mu_y = \frac{\sum_i h_i y_i}{n},  \mbox{ }
\sigma^2_{y} = \frac{\sum_i h_i  (y_i - \mu_y)^2}{n},  \mbox{ } 
\sigma^2_{y|x} = \sigma^2_{y} - \frac{\sigma^2_{xy}}{\sigma^2_{x}}.
\end{eqnarray*}
We use the data covariances to express the conditional expectations
and their estimated variances:
\begin{eqnarray}
	\hat{y} =  \mu_y + \frac{\sigma_{xy}}{\sigma^2_x}(x - \mu_x),  & &
\sigma^2_{\hat{y}} = \frac{\sigma^2_{y|x}}{n^2}
	\left(\sum_i h_i^2 +
	\frac{(x - \mu_x)^2}{\sigma_x^2}
	\sum_i h_i^2\frac{(x_i - \mu_x)^2}{\sigma_x^2}\right)	
\end{eqnarray}

\label{smooth_sec}
\begin{figure}[t]
\centerline{\psfig{figure=smooth_fig.ps,height=2in}}
\vspace{-0.2in}
\caption{\label{smooth_fig}
The estimator variance is minimized when the kernel includes as many
training points as can be accommodated by the model. Here the linear
LOESS model is shown. Too large a kernel includes points that degrade
the fit; too small a kernel neglects points that increase confidence
in the fit.}
\end{figure}
\subsection{Setting the Smoothing Parameter $k$}
\label{setting-k}
There are a number of ways one can set $k$, the smoothing parameter.
The method used by Cleveland et al. \citeyear{cdg88} is to set
$k$ such that the reference point being predicted has a predetermined
amount of support, that is, $k$ is set so that $n$ is close to some
target value. This has the disadvantage of requiring assumptions about
the noise and smoothness of the function being learned. Another
technique, used by Schaal and Atkeson \citeyear{schaal-atkeson94},
sets $k$ to minimize the crossvalidated error on the training set. A
disadvantage of this technique is that it assumes the distribution of
the training set is representative of $P(x)$, which it may not be in
an active learning situation. A third method, also described by Schaal
and Atkeson \citeyear{schaal-atkeson94}, is to set $k$ so as to
minimize the estimate of $\sigma^2_{\hat{y}}$ at the reference
points. As $k$ decreases, the regression becomes more global.  The
total weight $n$ will increase (which decreases $\sigma^2_{\hat{y}}$),
but so will the conditional variance $\sigma^2_{y|x}$ (which increases
$\sigma^2_{\hat{y}}$).  At some value of $k$, these two quantities
will balance to produce a minimum estimated variance (see
Figure~\ref{smooth_fig}). This estimate can be computed for arbitrary
reference points in the domain, and the user has the option of using
either a different $k$ for each reference point or a single global $k$
that minimizes the average $\sigma^2_{\hat{y}}$ over all reference
points.  Empirically, we found that the variance-based method gave the
best performance.

\subsection{Active Learning with Locally Weighted Regression}
As with the mixture of Gaussians, we want to select $\tilde{x}$ to
minimize $\left<\tilde{\sigma}^2_{\hat{y}}\right>$. To do this, we
must estimate the mean and variance of $P(\tilde{y}|\tilde{x})$. With
locally weighted regression, these are explicit: the mean is
$\hat{y}(\tilde{x})$ and the variance is $\sigma^2_{y|\tilde{x}}$.
The estimate of $\left<\tilde{\sigma}^2_{\hat{y}}\right>$ is also
explicit. Defining $\tilde{h}$ as the weight assigned to $\tilde{x}$
by the kernel we can compute these expectations exactly in closed form.
For the LOESS model, the learner's expected new variance is
\begin{eqnarray}
 \left<\tilde{\sigma}^2_{\hat{y}}\right> =  
	\frac{\left<\tilde{\sigma}^2_{y|x}\right>}{(n + \tilde{h})^2}
	\left[\sum_i h_i^2 + \tilde{h}^2 + 
	\frac{(x - \tilde{\mu}_x)^2}{\tilde{\sigma}_x^2}
	\left(\sum_i h_i^2\frac{(x_i - \tilde{\mu}_x)^2}
		{\tilde{\sigma}_x^2}
	 + \tilde{h}^2\frac{(\tilde{x} - \tilde{\mu}_x)^2}
		{\tilde{\sigma}_x^2}\right)\right].\label{lwr-exp} & &
\end{eqnarray}
Note that, since $\sum_i h_i^2 (x_i - \mu_x)^2 = \sum_i h_i^2 x^2_i +
\mu^2_x\sum_i h_i^2 - 2 \mu_x\sum_i h_i^2 x_i$, the new expectation of
Equation~\ref{lwr-exp} may be efficiently computed by caching the
values of $\sum_i h_i^2 x^2_i$ and $\sum_i h_i^2 x_i$. This obviates
the need to recompute the entire sum for each new candidate point.
The component expectations in Equation~\ref{lwr-exp} are computed as
follows:
\begin{eqnarray*}
\left<\tilde{\sigma}^2_{y|x}\right> =
   \left<\tilde{\sigma}^2_{y}\right> - 
	\frac{\left<\tilde{\sigma}^2_{xy}\right>}
	     {\tilde{\sigma}^2_{x}},
& &
\left<\tilde{\sigma}^2_{y}\right> =
	\frac{n\sigma^2_{y}}
	     {n+\tilde{h}}
	+ \frac{n\tilde{h}  \left(\sigma^2_{y|\tilde{x}} + (\hat{y}(\tilde{x}) - \mu_{y})^2\right)}
	     {(n+\tilde{h})^2},
\\
\tilde{\mu}_{x} = \frac{n\mu_{x} + \tilde{h}\tilde{x}}
		{n + \tilde{h}},
& &
\left<\tilde{\sigma}_{xy}\right> = 
	\frac{n\sigma_{xy}}
	     {n+\tilde{h}}
	+  \frac{n\tilde{h}  (\tilde{x} - \mu_{x})
			     (\hat{y}(\tilde{x}) - \mu_{y})}
			       {(n+\tilde{h})^2},
\\
\tilde{\sigma}^2_{x} = \frac{n\sigma^2_{x}}
			    {n+\tilde{h}}
			+ \frac{n\tilde{h}  (\tilde{x} - \mu_{x})^2}
			       {(n+\tilde{h})^2},
& &
\left<\tilde{\sigma}^2_{xy}\right>  = 
   \left<\tilde{\sigma}_{xy}\right>^2
    + \frac{n^2\tilde{h}^2\sigma^2_{y|\tilde{x}}(\tilde{x} - \mu_{x})^2}
	     {(n+\tilde{h})^4}.
\end{eqnarray*}
Just as with the mixture of Gaussians, we can use the expectation in
Equation~\ref{lwr-exp} to guide active learning.

\section{Experimental Results}
\label{experiment_section}
For an experimental testbed, we used the ``Arm2D'' problem described
by Cohn \citeyear{cohn94a}. The task is to learn the kinematics of
a toy 2-degree-of-freedom robot arm (see Figure~\ref{arm_fig}). The
inputs are joint angles $(\Theta_1, \Theta_2)$, and the outputs are
the Cartesian coordinates of the tip $(X_1, X_2)$. One of the implicit
assumptions of both models described here is that the noise is
Gaussian in the output dimensions. To test the robustness of the
algorithm to this assumption, we ran experiments using no noise, using
additive Gaussian noise in the outputs, and using additive Gaussian
noise in the inputs. The results of each were comparable; we report
here the results using additive Gaussian noise in the inputs. Gaussian
input noise corresponds to the case where the arm effectors or joint
angle sensors are noisy, and results in non-Gaussian errors in the
learner's outputs. The input distribution $P(x)$ is assumed to be
uniform.
\begin{figure}[hbt]
\centerline{\psfig{figure=arm_fig.ps,height=2in}}
\vspace{-0.2in}
\caption{\label{arm_fig} The arm kinematics problem. The learner
attempts to predict tip position given a set of joint angles
$(\theta_1, \theta_2)$.}
\end{figure}

We compared the performance of the variance-minimizing criterion by
comparing the {\em learning curves} of a learner using the criterion
with that of one learning from random samples.  The learning curves
plot the mean squared error and variance of the learner as its
training set size increases. The curves are created by starting with
an initial sample, measuring the learner's mean squared error or
estimated variance on a set of ``reference'' points (independent of the
training set), selecting and adding a new example to the training set,
retraining the learner on the augmented set, and repeating.

On each step, the variance-minimizing learner chose a set of 64
unlabeled reference points drawn from input distribution $P(x)$.  It
then selected a query $\tilde{x} = (\theta_1, \theta_2)$ that it
estimated would minimize $\left<\tilde{\sigma}^2_{y|x}\right>$ over
the reference set.  In the experiments reported here, the best
$\tilde{x}$ was selected from another set of 64 ``candidate'' points
drawn at random on each iteration.\footnote{As described earlier, we
could also have selected queries by hillclimbing on
$\partial\left<\tilde{\sigma}^2_{y|x}\right>/\partial\tilde{x}$; in
this low dimensional problem it was more computationally efficient to
consider a random candidate set.}


\subsection{Experiments with Mixtures of Gaussians}
With the mixtures of Gaussians model, there are three design
parameters that must be considered --- the number of Gaussians, their
initial placement, and the number of iterations of the EM algorithm.
We set these parameters by optimizing them on the learner using random
examples, then used the same settings on the learner using the
variance-minimization criterion. Parameters were set as follows:
Models with fewer Gaussians have the obvious advantage of requiring
less storage space and computation.  Intuitively, a small model should
also have the advantage of avoiding overfitting, which is thought to
occur in systems with extraneous parameters. Empirically, as we
increased the number of Gaussians, generalization improved
monotonically with diminishing returns (for a fixed training set size
and number of EM iterations).  The test error of the larger models
generally matched that of the smaller models on small training sets
(where overfitting would be a concern), and continued to decrease on
large training sets where the smaller networks ``bottomed out.'' We
therefore preferred the larger mixtures, and report here our results
with mixtures of 60 Gaussians.  We selected initial placement of the
Gaussians randomly, chosen uniformly from the smallest hypercube
containing all current training examples.  We arbitrarily chose the
identity matrix as an initial covariance matrix.  The learner was
surprisingly sensitive to the number of EM iterations. We examined a
range of 5 to 40 iterations of the EM algorithm per step. Small
numbers of iterations (5-10) appear insufficent to allow convergence
with large training sets, while large numbers of iterations (30-40)
degraded performance on small training sets. An ideal training regime
would employ some form of regularization, or would examine the degree
of change between iterations to detect convergence; in our
experiments, however, we settled on a fixed regime of 20 iterations
per step.
\begin{figure}[hbt]
\centerline{
	\parbox{3.0in}{\psfig{figure=em_noise_var.ps,width=3in}}
	\parbox{3.0in}{\psfig{figure=em_noise.ps,width=3in}}}
\vspace{-0.25in}
\caption{\label{em_figs} Variance and MSE learning curves for mixture
of 60 Gaussians trained on the Arm2D domain. Dotted lines denote
standard error for average of 10 runs, each started with one initial
random example.}
\end{figure}

Figure~\ref{em_figs} plots the variance and MSE learning curves for a
mixture of 60 Gaussians trained on the Arm2D domain with 1\% input
noise added. The estimated model variance using the
variance-minimizing criterion is significantly better than that of the
learner selecting data at random. The mean squared error, however,
exhibits even greater improvement, with an error that is consistently
$1/3$ that of the randomly sampling learner.

\subsection{Experiments with LOESS Regression}
With LOESS, the design parameters are the the size and shape of the
kernel. As described earlier, we arbitrarily chose to work with a
Gaussian kernel; we used the variance-based method for automatically
selecting the kernel size.
\begin{figure}[hbt]
\centerline{
	\parbox{3in}{\psfig{figure=new_lwr_noise_var.ps,width=3in}}
	\parbox{3in}{\psfig{figure=new_lwr_noise.ps,width=3in}}}
\vspace{-0.25in}
\caption{\label{lwr_figs}
Variance and MSE learning curves for LOESS model trained on the Arm2D
domain. Dotted lines denote standard error for average of 60 runs,
each started with a single initial random example.}
\end{figure}

In the case of LOESS, both the variance and the MSE of the learner
using the variance-minimizing criterion are significantly lower than
those of the learner selecting data randomly.  It is worth noting that
on the Arm2D domain, this form of locally weighted regression also
significantly outperforms both the mixture of Gaussians and the neural
networks discussed by Cohn \citeyear{cohn94a}.

\subsection{Computation Time}
One obvious concern about the criterion described here is its
computational cost. In situations where obtaining new examples may
take days and cost thousands of dollars, it is clearly wise to expend
computation to ensure that those examples are as useful as possible.
In other situations, however, new data may be relatively inexpensive,
so the computational cost of finding optimal examples must be
considered.

Table~\ref{time_table} summarizes the computation times for the two
learning algorithms discussed in this paper.\footnote{The times
reported are ``per reference point'' and ``per candidate per reference
point''; overall time must be computed from the number of candidates
and reference points examined. In the case of the LOESS model, for
example, with 100 training points, 64 reference points and 64
candidate points, the time required to select an action would be $(58
+ 0.16 \times 100)\times 4096 \mu$seconds, or about 0.3 seconds.} Note
that, with the mixture of Gaussians, training time depends linearly on
the number of examples, but prediction time is
independent. Conversely, with locally weighted regression, there is no
``training time'' per se, but the cost of additional examples accrues
when predictions are made using the training set.

While the training time incurred by the mixture of Gaussians may make
it infeasible for selecting optimal action learning actions in
realtime control, it is certainly fast enough to be used in many
applications.  Optimized, parallel implementations will also enhance
its utility.\footnote{It is worth mentioning that approximately half
of the training time for the mixture of Gaussians is spent computing
the correction factor in Equation~\ref{correction-eqn}. Without the
correction, the learner still computes $P(y|x)$, but does so by
modeling the training set distribution rather than the reference
distribution. We have found however, that for the problems examined,
the performance of such ``uncorrected'' learners does not differ
appreciably from that of the ``corrected'' learners.}  Locally
weighted regression is certainly fast enough for many control
applications, and may be made faster still by optimized, parallel
implementations. It is worth noting that, since the prediction speed
of these learners depends on their training set size, optimal data
selection is doubly important, as it creates a parsimonious training
set that allows faster predictions on future points.

\begin{table}
\begin{center}
\begin{tabular}{|l|l|l|l|} \hline
	& Training & Evaluating Reference & Evaluating Candidates \\ \hline
Mixture	& $3.9 + 0.05m$ sec & 15000 $\mu$sec      & 1300 $\mu$sec \\ 
LOESS	& -                 & $92+9.7m$ $\mu$sec   &  $58+0.16m$ $\mu$sec \\
\hline
\end{tabular}
\end{center}
\caption{\label{time_table}
Computation times on a Sparc 10 as a function of training set size
$m$.  Mixture model had 60 Gaussians trained for 20 iterations.
Reference times are per reference point; candidate times are per
candidate point per reference point.}
\end{table}

\section{Discussion}
Mixtures of Gaussians and locally weighted regression are two
statistical models that offer elegant representations and efficient
learning algorithms. In this paper we have shown that they also offer
the opportunity to perform active learning in an efficient and
statistically correct manner. The criteria derived here can be
computed cheaply and, for problems tested, demonstrate good predictive
power. In industrial settings, where gathering a single data point may
take days and cost thousands of dollars, the techniques described here
have the potential for enormous savings.

In this paper, we have only considered function approximation
problems. Problems requiring classification could be handled
analogously with the appropriate models. For learning classification
with a mixture model, one would select examples so as to maximize
discriminability between Gaussians; for locally weighted regression,
one would use a logistic regression instead of the linear one
considered here \cite{weisberg85}.

Our future work will proceed in several directions. The most important
is active bias minimization. As noted in Section~\ref{stats_section},
the learner's error is composed of both bias and variance. The
variance-minimizing strategy examined here ignores the bias component,
which can lead to significant errors when the learner's bias is
non-negligible.  Work in progress examines effective ways of measuring
and optimally eliminating bias \cite{cohn95}; future work will examine
how to jointly minimize both bias and variance to produce a criterion
that truly minimizes the learner's expected error.

Another direction for future research is the derivation of variance-
(and bias-) minimizing techniques for other statistical learning
models. Of particular interest is the class of models known as
``belief networks'' or ``Bayesian networks'' \cite{pearl88,hgc94}.
These models have the advantage of allowing inclusion of domain
knowledge and prior constraints while still adhering to a
statistically sound framework. Current research in belief networks
focuses on algorithms for efficient inference and learning; it would
be an important step to derive the proper criteria for learning {\em
actively} with these models.

%\newpage
\appendix
\section{Notation}
\begin{tabular}{ll}
& General \\ \hline
$X$		& input space \\
$Y$		& output space \\
$x$		& an arbitrary point in the input space \\
$y$		& true output value corresponding to input $x$ \\
$\hat{y}$	& predicted output value corresponding to input $x$ \\
$x_i$		& ``input'' part of example $i$ \\
$y_i$		& ``output'' part of example $i$ \\
$m$		& the number of examples in the training set \\
$\tilde{x}$	& specified input of a query \\
$\tilde{y}$	& the (possibly not yet known) output of query $\tilde{x}$ \\
$\sigma^2_{\hat{y}}$ & estimated variance of $\hat{y}$ \\ 
$\tilde{\sigma}^2_{\hat{y}}$ & new variance of $\hat{y}$, after
example $(\tilde{x}, \tilde{y})$ has been added \\ 
$\left<\tilde{\sigma}^2_{\hat{y}}\right>$ & the expected value of
	 $\tilde{\sigma}^2_{\hat{y}}$ \\
$P(x)$		& the (known) natural distribution over $x$ \\ 
& \\
& Neural Network \\ \hline
$w$		& a weight vector for a neural network \\
$\hat{w}$	& estimated ``best'' $w$ given a training set \\
$f_{\hat{w}}()$	& function computed by neural network given $\hat{w}$ \\
$S^2$		& average estimated noise in data, used as an estimate
		for $\sigma^2_y$  \\
& \\
& Mixture of Gaussians \\ \hline
$N$		& total number of Gaussians \\
$g_i$		& Gaussian number $i$ \\
$n_i$		& total point weighting attributed to Gaussian $i$ \\
$\mu_{x,i}$	& estimated $x$ mean of Gaussian $i$ \\
$\mu_{y,i}$	& estimated $y$ mean of Gaussian $i$ \\
$\sigma^2_{x,i}$& estimated $x$ variance of Gaussian $i$ \\
$\sigma^2_{y,i}$& estimated $y$ variance of Gaussian $i$ \\
$\sigma_{xy,i}$ & estimated $xy$ covariance of Gaussian $i$ \\
$\sigma^2_{y|x,i}$ & estimated $y$ variance of Gaussian $i$, given $x$ \\
$P(x, y|i)$	& joint distribution of input-output pair given Gaussian $i$ \\
$P(x|i)$	& distribution $x$ being given Gaussian $i$ \\
$h_i$		& weight of a given point that is attributed to Gaussian $i$ \\
$\tilde{h}_i$	& weight of new point $(\tilde{x}, \tilde{y})$ that is
		attributed to Gaussian $i$ \\
& \\
& Locally Weighted Regression \\ \hline
$k$		& kernel smoothing parameter \\
$h_i$		& weight given to example $i$ by kernel centered at $x$ \\
$n$		& sum of weights given to all points by kernel \\
$\mu_x$		& mean of inputs, weighted by kernel centered at $x$ \\
$\mu_y$		& mean of outputs, weighted by kernel centered at $x$ \\
$\tilde{h}$	& weight of new point $(\tilde{x}, \tilde{y})$ given 
		kernel centered at $x$ \\
\end{tabular}


\acks{David Cohn's current address is: Harlequin, Inc., One Cambridge
Center, Cambridge, MA 02142 USA. Zoubin Ghahramani's current address is:
Department of Computer Science, University of Toronto, Toronto,
Ontario M5S 1A4 CANADA. This work was funded by NSF grant CDA-9309300,
the McDonnell-Pew Foundation, ATR Human Information Processing
Laboratories and Siemens Corporate Research. We are deeply indebted to
Michael Titterington and Jim Kay, whose careful attention and
continued kind help allowed us to make several corrections to an
earlier version of this paper.}

\vskip 0.2in
\bibliography{statmodels}
\bibliographystyle{theapa}

\end{document}


