Optimal Learning for Sequential Decision Making for Expensive Cost Functions with Stochastic Binary Feedbacks

Yingfei Wang Department of Computer Science, Princeton University, Princeton, NJ 08540, yingfei@cs.princeton.edu Chu Wang The Program in Applied and Computational Mathematics, Princeton University, Princeton, NJ 08544 Warren Powell Department of Operations Research and Financial Engineering, Princeton University, Princeton, NJ 08544

Abstract

We consider the problem of sequentially making decisions that are rewarded by “successes" and “failures" which can be predicted through an unknown relationship that depends on a partially controllable vector of attributes for each instance. The learner takes an active role in selecting samples from the instance pool. The goal is to maximize the probability of success in either offline (training) or online (testing) phases. Our problem is motivated by real-world applications where observations are time consuming and/or expensive. We develop a knowledge gradient policy using an online Bayesian linear classifier to guide the experiment by maximizing the expected value of information of labeling each alternative. We provide a finite-time analysis of the estimated error and show that the maximum likelihood estimator based produced by the KG policy is consistent and asymptotically normal. We also show that the knowledge gradient policy is asymptotically optimal in an offline setting. This work further extends the knowledge gradient to the setting of contextual bandits. We report the results of a series of experiments that demonstrate its efficiency.

1 Introduction

There are many real-world optimization tasks where observations are time consuming and/or expensive. This often occurs in experimental sciences where testing different values of controllable parameters of a physical system is expensive. Another example arises in health services, where physicians have to make medical decisions (e.g. a course of drugs, surgery, and expensive tests) and we can characterize an outcome as a success (patient does not need to return for more treatment) or a failure (patient does need followup care such as repeated operations). This puts us in the setting of optimal learning where the number of objective function samples are limited, requiring that we learn from our decisions as quickly as possible. This represents a distinctly different learning environment than what has traditionally been considered using popular policies such as upper confidence bounding (UCB) which have proven effective in fast-paced settings such as learning ad-clicks.

In this paper, we are interested in a sequential decision-making setting where at each time step we choose one of finitely many decisions and observe a stochastic binary reward where each instance is described by various controllable and uncontrollable attributes. The goal is to maximize the probability of success in either offline (training) or online (testing) phases. There are a number of applications that can easily fit into our success/failure model:

•

Producing single-walled nanotubes. Scientists have physical procedures to produce nanotubes. It can produce either single-walled or double walled nanotubes through an unknown relationship with the controllable parameters, e.g. laser poser, ethylene, Hydrogen and pressure. Yet only the single-walled nanotubes are acceptable. The problem is to quickly learn the best parameter values with the highest probability of success (producing single-walled nanotubes).
•

Personalized health care. We consider the problem of how to choose clinical pathways (including surgery, medication and tests) for different upcoming patients to maximize the success of the treatment.
•

Minimizing the default rate for loan applications. When facing borrowers with different background information and credit history, a lending company needs to decide whether to grant a loan, and with what terms (interest rate, payment schedule).
•

Enhancing the acceptance of the admitted students. A university needs to decide which students to admit and how much aid to offer so that the students will then accept the offer of admission and matriculate.

Our work will initially focus on offline settings such as laboratory experiments or medical trials where we are not punished for errors incurred during training and instead are only concerned with the final recommendation after the offline training phases. We then extend our discussion to online learning settings where the goal is to minimize cumulative regret. We also consider problems with partially controllable attributes, which is known as contextual bandits. For example, in the health care problems, we do not have control over the patients (which is represented by a feature vector including demographic characteristics, diagnoses, medical history) and can only choose the medical decision. A university cannot control which students are applying to the university. When deciding whether to grant a loan, the lending company cannot choose the personal information of the borrowers.

Scientists can draw on an extensive body of literature on the classic design of experiments [12, 59, 42] where the goal is to decide which observations to make when fitting a function. Yet in our setting, the decisions are guided by a well-defined utility function (that is, maximize the probability of success). The problem is related to the literature on active learning [49, 55, 18, 50], where our setting is most similar to membership query synthesis where the learner may request labels for any unlabeled instance in the input space to learn a classifier that accurately predicts the labels of new examples. By contrast, our goal is to maximize a utility function such as the success of a treatment. Moreover, the expense of labeling each alternative sharpens the conflicts of learning the prediction and finding the best alternative.

Another similar sequential decision making setting is multi-armed bandit problems [2, 6, 14, 39, 52, 7, 56]. Different belief models have been studied under the name of contextual bandits, including linear models [9] and Gaussian process regression [33]. The focus of bandit work is minimizing cumulative regret in an online setting, while we consider the performance of the final recommendation after an offline training phase. There are recent works to address the problem we describe here by minimizing the simple regret. But first, the UCB type policies [1] are not best suited for expensive experiments. Second, the work on simple regret minimization [28, 27] mainly focuses on real-valued functions and do not consider the problem with stochastic binary feedbacks.

There is a literature on Bayesian optimization [26, 8, 44]. Efficient global optimization (EGO), and related methods such as sequential kriging optimization (SKO) [32, 30] assume a Gaussian process belief model which does not scale to the higher dimensional settings that we consider. Others assume lookup table, or low-dimensional parametric methods, e.g. response surface/surrogate models [24, 31, 46]. The existing literature mainly focuses on real-valued functions and none of these methods are directly suitable for our problem of maximizing the probability of success with binary outcomes. A particularly relevant body of work in the Bayesian optimization literature is the expected improvement (EI) for binary outputs [53]. Yet when EI decides which alternative to measure, it is based on the expected improvement over current predictive posterior distribution while ignoring the potential change of the posterior distribution resulting from the next stochastic measurement (see Section 5.6 of [44] and [30] for detailed explanations).

We investigate a knowledge gradient policy that maximizes the value of information, since this approach is particularly well suited to problems where observations are expensive. After its first appearance for ranking and selection problems [16], KG has been extended to various other belief models (e.g. [41, 43, 57]). Yet there is no KG variant designed for binary classification with parametric belief models. In this paper, we extend the KG policy to the setting of classification problems under a Bayesian classification belief model which introduces the computational challenge of working with nonlinear belief models.

This paper makes the following contributions. 1. Due to the sequential nature of the problem, we develop a fast online Bayesian linear classification procedure based on Laplace approximation for general link functions to recursively predict the response of each alternative. 2. We design a knowledge-gradient type policy for stochastic binary responses to guide the experiment. It can work with any choice of link function $\sigma(\cdot)$ and approximation procedure. 3. Since the knowledge gradient policy adaptively chooses the next sampled points, we provide a finite-time analysis on the estimated error that does not rely on the i.i.d. assumption. 4. We show that the maximum likelihood estimator based on the adaptively sampled points by the KG policy is consistent and asymptotically normal. We show furthermore that the knowledge gradient policy is asymptotic optimal. 5. We extend the KG policy to contextual bandit settings with stochastic binary outcomes.

The paper is organized as follows. In Section 2, we establish the mathematical model for the problem of sequentially maximizing the response under binary outcomes. We give a brief overview of (Bayesian) linear classification in Section 3. We develop an online Bayesian linear classification procedure based on Laplace approximation in Section 4. In Section 5, we propose a knowledge-gradient type policy for stochastic binary responses with or with context information. We theoretically study its finite-time and asymptotic behavior. Extensive demonstrations and comparisons of methods for offline objective are demonstrated in Section 6.

2 Problem formulation

We assume that we have a finite set of alternatives $\bm{x}\in\mathcal{X}=\{\bm{x}_{1},\dots,\bm{x}_{M}\}$ . The observation of measuring each $\bm{x}$ is a binary outcome $y\in\{-1,+1\}$ /{failure, success} with some unknown probability $p(y=+1|\bm{x})$ . The learner sequentially chooses a series of points $(\bm{x}^{0},\dots,\bm{x}^{N-1})$ to run the experiments. Under a limited measurement budget $N$ , the goal of the learner is to recommend an implementation decision $\bm{x}^{N}$ that maximizes $p(y=+1|\bm{x}^{N})$ .

We adopt probabilistic models for classification. Under general assumptions, the probability of success can be written as a link function acting on a linear function of the feature vector

p(y=+1|\bm{x})=\sigma(\bm{w}^{T}\bm{x}).

In this paper, we illustrate the ideas using the logistic link function $\sigma(a)=\frac{1}{1+\text{exp}(-a)}$ and probit link function $\sigma(a)=\Phi(a)=\int_{-\infty}^{a}\mathcal{N}(s|0,1^{2})\text{d}s$ given its analytic simplicity and popularity, but any monotonically increasing function $\sigma:\mathbb{R}\mapsto[0,1]$ can be used. The main difference between the two sigmoid functions is that the logistic function has slightly heavier tails than the normal CDF. Classification using the logistic function is called logistic regression and that using the normal CDF is called probit regression.

We start with a multivariate prior distribution for the unknown parameter vector $\bm{w}$ . At iteration $n$ , we choose an alternative $\bm{x}^{n}$ to measure and observe a binary outcome $y^{n+1}$ assuming labels are generated independently given $\bm{w}$ . Each alternative can be evaluated more than once with potentially different outcomes. Let $\mathcal{D}^{n}=\{(\bm{x}^{i},y^{i+1})\}_{i=0}^{n}$ denote the previous measured data set for any $n=0,\dots,N$ . Define the filtration $(\mathcal{F}^{n})_{n=0}^{N}$ by letting $\mathcal{F}^{n}$ be the sigma-algebra generated by $\bm{x}^{0},y^{1},\dots,\bm{x}^{n-1},y^{n}$ . We use $\mathcal{F}^{n}$ and $\mathcal{D}^{n}$ interchangeably. Note that the notation here is slightly different from the (passive) PAC learning model where the data are i.i.d. and are denoted as $\{(\bm{x}_{i},y_{i})\}$ . Yet in our (adaptive) sequential decision setting, measurement and implementation decisions $\bm{x}^{n}$ are restricted to be $\mathcal{F}^{n}$ -measurable so that decisions may only depend on measurements made in the past. This notation with superscript indexing time stamp is standard, for example, in control theory, stochastic optimization and optimal learning. We use Bayes’ theorem to form a sequence of posterior predictive distributions $p(\bm{w}|\mathcal{D}^{n})$ .

The next lemma states the equivalence of using true probabilities and sample estimates when evaluating a policy. The proof is left in the supplementary material.

Lemma 2.1.

Let $\Pi$ be the set of policies, $\pi\in\Pi$ , and $\bm{x}^{\pi}=\arg\max_{\bm{x}}p(y=+1|\bm{x},\mathcal{D}^{N})$ be the implementation decision after the budget $N$ is exhausted. Then

\mathbb{E}_{\bm{w}}[p(y=+1|\bm{x}^{\pi},\bm{w})]=\mathbb{E}_{\bm{w}}[\max_{\bm{x}}p(y=+1|\bm{x},\mathcal{D}^{N})],

where the expectation is taken over the prior distribution of $\bm{w}$ .

By denoting $\mathcal{X}^{I}$ as an implementation policy for selecting an alternative after the measurement budget is exhausted, then $\mathcal{X}^{I}$ is a mapping from the history $\mathcal{D}^{N}$ to an alternative $\mathcal{X}^{I}(\mathcal{D}^{N})$ . Then as a corollary of Lemma 2.1, we have [44]

\max_{\mathcal{X}^{I}}\mathbb{E}\big{[}p\big{(}y=+1|\mathcal{X}(\mathcal{D}^{N})\big{)}\big{]}=\max_{\bm{x}}p(y=+1|\bm{x},\mathcal{D}^{N}).

In other words, the optimal decision at time $N$ is to go with our final set of beliefs. By the equivalence of using true probabilities and sample estimates when evaluating a policy as stated in Lemma 2.1, while we want to learn the unknown true value $\max_{\bm{x}}p(y=+1|\bm{x})$ , we may write our objective function as

\max_{\pi\in\Pi}\mathbb{E}^{\pi}[\max_{\bm{x}}p(y=+1|\bm{x},\mathcal{D}^{N})].

(1)

3 Background: Linear classification

Linear classification, especially logistic regression, is widely used in machine learning for binary classification [29]. Assume that the probability of success $p(y=+1|\bm{x})$ is a parameterized function $\sigma(\bm{w}^{T}\bm{x})$ and further assume that observations are independently of each other. Given a training set $\mathcal{D}=\{(\bm{x}_{i},y_{i})\}_{i=1}^{n}$ with $\bm{x}_{i}$ a $d$ -dimensional vector and $y_{i}\in\{-1,+1\}$ , the likelihood function $p(\mathcal{D}|\bm{w})$ is $p(\mathcal{D}|\bm{w})=\prod_{i=1}^{n}\sigma(y_{i}\cdot\bm{w}^{T}\bm{x}_{i}).$ The weight vector $\bm{w}$ is found by maximizing the likelihood $p(\mathcal{D}|\bm{w})$ or equivalently, minimizing the negative log likelihood:

\min_{\bm{w}}\sum_{i=1}^{n}-\log(\sigma(y_{i}\cdot\bm{w}^{T}\bm{x}_{i})).

In order to avoid over-fitting, especially when there are a large number of parameters to be learned, $l_{2}$ regularization is often used. The estimate of the weight vector $\bm{w}$ is then given by:

\min_{\bm{w}}\frac{\lambda}{2}\|\bm{w}\|^{2}-\sum_{i=1}^{n}\log(\sigma(y_{i}\cdot\bm{w}^{T}\bm{x}_{i})).

(2)

It can be shown that this log-likelihood function is globally concave in $\bm{w}$ for both logistic regression or probit regression. As a result, numerous optimization techniques are available for solving it, such as steepest ascent, Newton’s method and conjugate gradient ascent.

This logic is suitable for batch learning where we only need to conduct the minimization once to find the estimation of weight vector $\bm{w}$ based on a given batch of training examples $\mathcal{D}$ . Yet due to the sequential nature of our problem setting, observations come one by one as in online learning. After each new observation, if we retrain the linear classifier using all the previous data, we need to re-do the minimization, which is computationally inefficient. In this paper, we instead extend Bayesian linear classification to perform recursive updates with each observation.

A Bayesian approach to linear classification models requires a prior distribution for the weight parameters $\bm{w}$ , and the ability to compute the conditional posterior $p(\bm{w}|\mathcal{D})$ given the observation. Specifically, suppose we begin with an arbitrary prior $p(\bm{w})$ and apply Bayes’ theorem to calculate the posterior: $p(\bm{w}|\mathcal{D})=\frac{1}{Z}p(\mathcal{D}|\bm{w})p(\bm{w}),$ where the normalization constant $Z$ is the unknown evidence. An $l_{2}$ -regularized logistic regression can be interpreted as a Bayesian model with a Gaussian prior on the weights with standard deviation $1/\sqrt{\lambda}$ .

Unfortunately, exact Bayesian inference for linear classifiers is intractable since the evaluation of the posterior distribution comprises a product of sigmoid functions; in addition, the integral in the normalization constant is intractable as well, for both the logistic function or probit function. We can either use analytic approximations to the posterior, or solutions based on Monte Carlo sampling, foregoing a closed-form expression for the posterior. In this paper, we consider different analytic approximations to the posterior to make the computation tractable.

3.1 Online Bayesian probit regression based on assumed Gaussian density filtering

Assumed density filtering (ADF) is a general online learning schema for computing approximate posteriors in statistical models [5, 35, 40, 48]. In ADF, observations are processed one by one, updating the posterior which is then approximated and is used as the prior distribution for the next observation.

For a given Gaussian prior distribution on some latent parameter $\bm{\theta}$ , $p(\bm{\theta})=\mathcal{N}(\bm{\theta}|\bm{\mu},\bm{\Sigma})$ and a likelihood $t(\bm{\theta}):=p(\mathcal{D}|\bm{\theta})$ , the posterior $p(\bm{\theta}|\mathcal{D})$ is generally non-Gaussian,

p(\bm{\theta}|\mathcal{D})=\frac{t(\bm{\theta})p(\bm{\theta})}{\int t(\tilde{\bm{\theta}})p(\tilde{\bm{\theta}})\text{d}\tilde{\bm{\theta}}}.

We find the best approximation by minimizing the Kullback-Leibler (KL) divergence between the true posterior $p(\bm{\theta}|\mathcal{D})$ and the Gaussian approximation. It is well known that when $q(x)$ is Gaussian, the distribution $q(x)$ that minimizes KL $(p(x)||q(x))$ is the one whose first and second moments match that of $p(x)$ . It can be shown that the Gaussian approximation $\hat{q}(\bm{\theta})=\mathcal{N}(\bm{\theta}|\hat{\bm{\mu}},\hat{\bm{\Sigma}})$ found by moment matching is given as:

\hat{\bm{\mu}}=\bm{\mu}+\bm{\Sigma}\bm{g},~~~\hat{\bm{\Sigma}}=\bm{\Sigma}-\bm{\Sigma}(\bm{g}\bm{g}^{T}-2\bm{G})\bm{\Sigma},

(3)

where the vector $\bm{g}$ and the matrix $\bm{G}$ are given by

\bm{g}=\frac{\partial\log\bigg{(}Z(\tilde{\bm{\mu}},\tilde{\bm{\Sigma}})\bigg{)}}{\partial\tilde{\bm{\mu}}}|_{\tilde{\bm{\mu}}=\bm{\mu},\tilde{\bm{\Sigma}}=\bm{\Sigma}},~~~\bm{G}=\frac{\partial\log\bigg{(}Z(\tilde{\bm{\mu}},\tilde{\bm{\Sigma}})\bigg{)}}{\partial\tilde{\bm{\Sigma}}}|_{\tilde{\bm{\mu}}=\bm{\mu},\tilde{\bm{\Sigma}}=\bm{\Sigma}},

and the normalization function $Z$ is defined by

Z(\bm{\mu},\bm{\Sigma}):=\int t(\tilde{\bm{\theta}})\mathcal{N}(\tilde{\bm{\theta}}|\bm{\mu},\bm{\Sigma})\text{d}\tilde{\bm{\theta}}.

For the sake of analytic convenience, we only consider probit regression under assumed Gaussian density filtering. Specifically, the distribution of $\bm{w}$ after $n$ observations is modeled as $p(\bm{w})=\mathcal{N}(\bm{w}|\bm{\mu}^{n},\bm{\Sigma}^{n})$ . The likelihood function for the next available data $(\bm{x},y)$ is $t(\bm{w}):=\Phi(y\bm{w}^{t}\bm{x})$ . Thus we have,

p(\bm{w}|\bm{x},y)\propto\Phi(y\bm{w}^{t}\bm{x})\mathcal{N}(\bm{w}|\bm{\mu}^{n},\bm{\Sigma}^{n}).

Since the convolution of the normal CDF and a Gaussian distribution is another normal CDF, moment matching (3) results in an analytical solution to the Gaussian approximation:

	$\displaystyle\bm{\mu}^{n+1}$	$\displaystyle=$	$\displaystyle\bm{\mu}^{n}+\frac{y\bm{\Sigma}^{n}\bm{x}}{\sqrt{1+\bm{x}^{T}\bm{\Sigma}^{n}\bm{x}}}v\bigg{(}\frac{y\bm{x}^{T}\bm{\mu}^{n}}{\sqrt{1+\bm{x}^{T}\bm{\Sigma}^{n}\bm{x}}}\bigg{)},$		(4)
	$\displaystyle\bm{\Sigma}^{n+1}$	$\displaystyle=$	$\displaystyle\bm{\Sigma}^{n}-\frac{(\bm{\Sigma}^{n}\bm{x})(\bm{\Sigma}^{n}\bm{x})^{T}}{1+\bm{x}^{T}\bm{\Sigma}^{n}\bm{x}}w\bigg{(}\frac{y\bm{x}^{T}\bm{\mu}^{n}}{\sqrt{1+\bm{x}^{T}\bm{\Sigma}^{n}\bm{x}}}\bigg{)},$		(5)

where

v(z):=\frac{\mathcal{N}(z|0,1)}{\Phi(z)}\text{ and }w(z):=v(z)\Big{(}v(z)+z\Big{)}.

In this work, we focus on diagonal covariance matrices $\bm{\Sigma}^{n}$ with $(\sigma_{i}^{n})^{2}$ as the diagonal element due to computational simplicity and its equivalence with $l_{2}$ regularization, resulting in the following update for the posterior parameters:

	$\displaystyle\mu^{n+1}_{i}$	$\displaystyle=$	$\displaystyle\mu^{n}_{i}+\frac{yx_{i}(\sigma_{i}^{n})^{2}}{\tilde{\sigma}}v(\frac{y\bm{x}^{T}\bm{\mu}^{n}}{\tilde{\sigma}}),$		(6)
	$\displaystyle(\sigma_{i}^{n+1})^{2}$	$\displaystyle=$	$\displaystyle(\sigma_{i}^{n})^{2}-\frac{x_{i}^{2}(\sigma^{n}_{i})^{4}}{\tilde{\sigma}^{2}}w(\frac{y\bm{x}^{T}\bm{\mu}^{n}}{\tilde{\sigma}}),$		(7)

where $\tilde{\sigma}^{2}:=1+\sum_{j=1}^{d}(\sigma_{j}^{n})^{2}x_{j}^{2}$ . See, for example, [22] and [10] for successful applications of this online probit regression model in prediction of click-through rates and stream-based active learning.

Due to the popularity of logistic regression and the computational limitations of ADF (on general link functions other than probit function), we develop an online Bayesian linear classification procedure for general link functions to recursively predict the response of each alternative in the next section.

4 Online Bayesian Linear Classification based on Laplace approximation

In this section, we consider the Laplace approximation to the posterior and develop an online Bayesian linear classification schema for general link functions.

4.1 Laplace approximation

Laplace’s method aims to find a gaussian approximation to a probability density defined over a set of continuous variables. It can be obtained by finding the mode of the posterior distribution and then fitting a Gaussian distribution centered at that mode [4]. Specifically, define the logarithm of the unnormalized posterior distribution as

\displaystyle\Psi(\bm{w})=\log p(\mathcal{D}|\bm{w})+\log p(\bm{w}).

Since the logarithm of a Gaussian distribution is a quadratic function, we consider a second-order Taylor expansion to $\Psi$ around its MAP (maximum a posteriori) solution $\hat{\bm{w}}=\arg\max_{\bm{w}}\Psi(\bm{w})$ :

\Psi(\bm{w})\approx\Psi(\hat{\bm{w}})-\frac{1}{2}(\bm{w}-\hat{\bm{w}})^{T}\bm{H}(\bm{w}-\hat{\bm{w}}),

(8)

where $\bm{H}$ is the Hessian of the negative log posterior evaluated at $\hat{\bm{w}}$ :

\bm{H}=-\nabla^{2}\Psi(\bm{w})|_{\bm{w}=\hat{\bm{w}}}.

By exponentiating both sides of Eq. (8), we can see that the Laplace approximation results in a normal approximation to the posterior

p(\bm{w}|\mathcal{D})\approx\mathcal{N}(\bm{w}|\hat{\bm{w}},\bm{H}^{-1}).

(9)

For multivariate Gaussian priors $p(\bm{w})=\mathcal{N}(\bm{w}|\bm{m},\bm{\Sigma})$ ,

\Psi(\bm{w}|\bm{m},\bm{\Sigma})=-\frac{1}{2}(\bm{w}-\bm{m})^{T}\bm{\Sigma}^{-1}(\bm{w}-\bm{m})+\sum_{i=1}^{n}\log(\sigma(y_{i}\cdot\bm{w}^{T}\bm{x}_{i})),

(10)

and the Hessian $\bm{H}$ evaluated at $\hat{\bm{w}}$ is given for both logistic and normal CDF link functions as:

\bm{H}=\bm{\Sigma}^{-1}-\sum_{i=1}^{n}\hat{t}_{i}\bm{x}_{i}\bm{x}_{i}^{T},

(11)

where $\hat{t}_{i}:=\frac{\partial^{2}\log p(y_{i}|\bm{x}_{i},\bm{w})}{\partial f_{i}^{2}}|_{f_{i}=\hat{\bm{w}}^{T}\bm{x}_{i}}$ and $f_{i}=\bm{w}^{T}\bm{x}_{i}$ .

4.2 Online Bayesian linear classification based on Laplace approximation

Starting from a Gaussian prior $\mathcal{N}(\bm{w}|\bm{m}^{0},\bm{\Sigma}^{0})$ , after the first $n$ observations, the Laplace approximated posterior distribution is $p(\bm{w}|\mathcal{D}^{n})\approx\mathcal{N}(\bm{w}|\bm{m}^{n},\bm{\Sigma}^{n})$ according to (9). We formally define the state space $\mathcal{S}$ to be the cross-product of $\mathbb{R}^{d}$ and the space of positive semidefinite matrices. At each time $n$ , our state of knowledge is thus $S^{n}=(\bm{m}^{n},\bm{\Sigma}^{n})$ . Observations come one by one due to the sequential nature of our problem setting. After each new observation, if we retrain the Bayesian classifier using all the previous data, we need to calculate the MAP solution of (10) with $\mathcal{D}=\mathcal{D}^{n}$ to update from $S^{n}$ to $S^{n+1}$ . It is computationally inefficient even with a diagonal covariance matrix. It is better to extend the Bayesian linear classifier to handle recursive updates with each observation.

Here, we propose a fast and stable online algorithm for model updates with independent normal priors (with $\bm{\Sigma}=\lambda^{-1}\bm{I}$ , where $\bm{I}$ is the identity matrix), which is equivalent to $l_{2}$ regularization and which also offers greater computational efficiency [58]. At each time step $n$ , the Laplace approximated posterior $\mathcal{N}(\bm{w}|\bm{m}^{n},\bm{\Sigma}^{n})$ serves as a prior to update the model when the next observation is made. In this recursive way of model updating, previously measured data need not be stored or used for retraining the model. By setting the batch size $n=1$ in Eq. (10) and (11), we have the sequential Bayesian linear model for classification as in Algorithm 1, where $\hat{t}:=\frac{\partial^{2}\log(\sigma(yf))}{\partial f^{2}}|_{f=\hat{\bm{w}}^{T}\bm{x}}$ .

input : Regularization parameter

\lambda>0

m_{j}=0

q_{j}=\lambda

. (Each weight

w_{j}

has an independent prior

\mathcal{N}(w_{j}|m_{j},q_{j}^{-1})

)

for $t=1$ to $T$ do

Get a new point

(\bm{x},y)

Find

\hat{\bm{w}}

as the maximizer of (10):

-\frac{1}{2}\sum_{j=1}^{d}q_{i}(w_{i}-m_{i})^{2}+\log(\sigma(y_{i}\bm{w}^{T}\bm{x}_{i})).

m_{j}=\hat{w}_{j}

Update

q_{i}

according to (11):

q_{j}\leftarrow q_{j}-\hat{t}x^{2}_{j}

end for

Algorithm 1 Online Bayesian linear classification

It is generally assumed that $\log\sigma(\cdot)$ is concave to ensure a unique solution of Eq. (10). It is satisfied by commonly used sigmoid functions for classification problems, including logistic function, probit function, complementary log-log function $\sigma(a)=1-\exp(-\exp(a))$ and log-log function $\exp(-\exp(-a))$ .

We can tap a wide range of convex optimization algorithms including gradient search, conjugate gradient, and BFGS method [60]. But if we set $n=1$ and $\Sigma=\lambda^{-1}\bm{I}$ in Eq. (10), a stable and efficient algorithm for solving

\arg\max_{\bm{w}}-\frac{1}{2}\sum_{j=1}^{d}q_{i}(w_{i}-m_{i})^{2}+\log(\sigma(y\bm{w}^{T}\bm{x}))

(12)

can be obtained as follows. First, taking derivatives with respect to $w_{i}$ and setting $\frac{\partial F}{\partial w_{i}}$ to zero, we have

q_{i}(w_{i}-m_{i})=\frac{yx_{i}\sigma^{\prime}(y\bm{w}^{T}x)}{\sigma(y\bm{w}^{T}x)},~~~~i=1,2,\dots,d.

Defining $p$ as

p:=\frac{\sigma^{\prime}(y\bm{w}^{T}x)}{\sigma(y\bm{w}^{T}x)},

we then have $w_{i}=m_{i}+yp\frac{x_{i}}{q_{i}}.$ Plugging this back into the definition of $p$ to eliminate $w_{i}$ ’s, we get the equation for $p$ :

p=\frac{\sigma^{\prime}(p\sum_{i=1}^{d}x_{i}^{2}/q_{i}+y\bm{m}^{T}x)}{\sigma(p\sum_{i=1}^{d}x_{i}^{2}/q_{i}+y\bm{m}^{T}x)}.

Since $\log(\sigma(\cdot))$ is concave, by its derivative we know the function $\sigma^{\prime}/\sigma$ is monotonically decreasing, and thus the right hand side of the equation decreases as $p$ goes from $0$ to $\infty$ . We notice that the right hand side is positive when $p=0$ and the left hand side is larger than the right hand side when $p=\sigma^{\prime}(y\bm{m}^{T}x)/\sigma(y\bm{m}^{T}x)$ . Hence the equation has a unique solution in interval $[0,\sigma^{\prime}(y\bm{m}^{T}x)/\sigma(y\bm{m}^{T}x)]$ . A simple one dimensional bisection method is sufficient to efficiently find the root $p^{*}$ and thus the solution to the $d$ -dimensional optimization problem (12).

We illustrate and validate this schema using logistic and probit functions. For logistic function $\sigma(\bm{w}^{T}\bm{x})=-\log(1+\exp(-y\bm{w}^{T}\bm{x}))$ , by setting $\partial F/\partial w_{i}=0$ for all $i$ and then by denoting $(1+\exp(y\mathbf{w}^{T}\mathbf{x}))^{-1}$ as $p$ , we have

q_{i}(w_{i}-m_{i})=ypx_{i},~~~~i=1,2,\dots,d,

resulting in the following equation for $p$ :

\frac{1}{p}=1+\exp\Big{(}y\sum_{i=1}^{d}(m_{i}+yp\frac{x_{i}}{q_{i}})x_{i}\Big{)}=1+\exp(y\mathbf{m}^{T}\mathbf{x})\exp\Big{(}y^{2}p\sum_{i=1}^{d}\frac{x_{i}^{2}}{q_{i}}\Big{)}.

It is easy to see that the left hand side decreases from infinity to 1 and the right hand side increases from 1 when $p$ goes from 0 to 1, therefore the solution exists and is unique in $[0,1]$ .

For normal CDF link function $\sigma(\cdot)=\Phi(\cdot)$ , the computation is a little lengthier. First use $\phi(\cdot)$ to denote the standard normal distribution $\mathcal{N}(\cdot|0,1^{2})$ , we have

\frac{\partial F}{\partial w_{i}}=-q_{i}(w_{i}-m_{i})+\frac{yx_{i}\phi(\bm{w}^{T}\bm{x})}{\Phi(y\bm{w}^{T}\bm{x})}.

Let $p=\frac{\phi(\bm{w}^{T}\bm{x})}{\Phi(y\bm{w}^{T}\bm{x})}$ and $\partial F/\partial w_{i}=0$ for all $i$ , and thus we have $w_{i}=m_{i}+yp\frac{x_{i}}{q_{i}}$ . Plugging these into the definition of $p$ we have

p=\frac{\phi\big{(}\sum_{i=1}^{d}(m_{i}+yp\frac{x_{i}}{q_{i}})x_{i}\big{)}}{\Phi\big{(}y\sum_{i=1}^{d}(m_{i}+yp\frac{x_{i}}{q_{i}})x_{i}\big{)}}.

(13)

Define the right-hand-size of the above equation as $g(p)$ . The next lemma shows that $\frac{\text{d}}{\text{d}{p}}g(p)\leq 0$ and thus the right-hand-side of Eq. (13) is non-increasing. Together with the fact that the left hand side increases from 0 to 1, the bisection method can also be used based on Eq. (13). The proof can be found in Appendix 9.

Lemma 4.1.

Define

g(p):=\frac{\phi\big{(}\sum_{i=1}^{d}(m_{i}+yp\frac{x_{i}}{q_{i}})x_{i}\big{)}}{\Phi\big{(}y\sum_{i=1}^{d}(m_{i}+yp\frac{x_{i}}{q_{i}})x_{i}\big{)}}.

We have $\frac{\text{d}}{\text{d}{p}}g(p)\leq 0$ for every $p\in\mathbb{R}$ .

5 Knowledge Gradient Policy for Bayesian Linear Classification Belief Model

We begin with a brief review of the knowledge gradient (KG) for ranking and selection problems, where each of the alternative can be measured sequentially to estimates its unknown underlying expected performance $\mu_{x}$ . The goal is to adaptively allocate alternatives to measure so as to find an implementation decision that has the largest mean after the budget is exhausted. In a Bayesian setting, the performance of each alternative is represented by a (non-parametric) lookup table model of Gaussian distribution. Specifically, by imposing a Gaussian prior $\mathcal{N}(\bm{\mu}|\bm{\theta}^{0},\bm{\Sigma}^{0})$ , the posterior after the first $n$ observations is denoted by $\mathcal{N}(\bm{\mu}|\bm{\theta}^{n},\bm{\Sigma}^{n})$ . At the $n$ th iteration, the knowledge gradient policy chooses its ( $n+1$ )th measurement to maximize the single-period expected increase in value [16]:

\nu_{x}^{\text{KG},n}=\mathbb{E}[\max_{x^{\prime}}\theta_{x^{\prime}}^{n+1}-\max_{x^{\prime}}\theta_{x^{\prime}}^{n}|x^{n}=x,S^{n}].

It enjoys nice properties, including myopic and asymptotic optimality. KG has been extended to various belief models (e.g. [41, 43, 47, 57]). The knowledge gradient can be extended to online problems where we need to maximize cumulative rewards [47],

\nu_{x}^{\text{OLKG},n}=\theta_{x}^{n}+\tau\nu^{\text{KG},n}_{x},

where $\tau$ reflects a planning horizon.

Yet there is no KG variant designed for binary classification with parametric models, primarily because of the computational intractability of dealing with nonlinear belief models. In what follows, we first formulate our learning problem as a Markov decision process and then extend the KG policy for stochastic binary outcomes where, for example, each choice (say, a medical decision) influences the success or failure of a medical outcome.

5.1 Markov decision process formulation

Our learning problem is a dynamic program that can be formulated as a Markov decision process. Define the state space $\mathcal{S}$ as the space of all possible predictive distributions for $\bm{w}$ . By Bayes’ Theorem, the transition function $T$ : $\mathcal{S}\times\mathcal{X}\times\{-1,1\}$ is:

T\bigg{(}q(\bm{w}),\bm{x},y\bigg{)}\propto q(\bm{w})\sigma(y\bm{w}^{T}\bm{x}).

(14)

If we start from a Gaussian prior $\mathcal{N}(\bm{w}|\bm{\mu}^{0},\bm{\Sigma}^{0})$ , after the first $n$ observed data, the approximated posterior distribution is $p(\bm{w}|\mathcal{D}^{n})\approx\mathcal{N}(\bm{w}|\bm{\mu}^{n},\bm{\Sigma}^{n})$ . The state space $\mathcal{S}$ is the cross-product of $\mathbb{R}^{d}$ and the space of positive semidefinite matrices. The transition function for updating the belief state depends on the belief model $\sigma(\cdot)$ and the approximation strategy. For example, for different update equations in Algorithm 1 and (6)(7) under different approximation methods, the transition function can be defined as follows with degenerate state space $\mathcal{S}:=\mathbb{R}^{d}\times[0,\infty)^{d}$ :

Definition 5.1.

The transition function based on online Bayesian classification with Laplace approximation $T$ : $\mathcal{S}\times\mathcal{X}\times\{-1,1\}$ is defined as

T^{\text{L}}\Big{(}(\bm{\mu},\bm{\sigma}^{2}),\bm{x},y\Big{)}=\Big{(}\arg\min_{\bm{w}}\Psi\big{(}\bm{w}|\bm{\mu},\bm{\sigma}^{-2}\big{)},\bm{\sigma}^{-2}+\hat{t}(y)\cdot\textbf{diag}(\bm{x}\bm{x}^{T})\Big{)},

where $\hat{t}(y):=\frac{\partial^{2}\log p(y|\bm{x},\bm{w})}{\partial f^{2}}|_{f=\hat{\bm{w}}^{T}\bm{x}}$ for either logistic or probit functions, $\textbf{diag}(\bm{x}\bm{x}^{T})$ is a column vector containing the diagonal elements of $\bm{x}\bm{x}^{T}$ and $\bm{\sigma}^{-2}$ is understood as a column vector containing $\sigma_{i}^{-2}$ , so that $S^{n+1}=T^{\text{L}}(S^{n},\bm{x},Y^{n+1})$ . $Y^{n+1}$ denotes the unobserved binary random variable at time $n$ .

Definition 5.2.

The transition function based on assumed density filtering $T$ : $\mathcal{S}\times\mathcal{X}\times\{-1,1\}$ is defined as

T^{\text{ADF}}\Big{(}(\bm{\mu},\bm{\sigma}^{2}),\bm{x},y\Big{)}=\Big{(}\bm{\mu}+\frac{y\bm{x}^{T}\bm{\sigma}^{2}}{\tilde{\sigma}}v(\frac{y\bm{x}^{T}\bm{\mu}}{\tilde{\sigma}}),\bm{\sigma}^{2}-\frac{(\bm{x}^{2})^{T}\bm{\sigma}^{4}}{\tilde{\sigma}^{2}}w(\frac{y\bm{x}^{T}\bm{\mu}}{\tilde{\sigma}})\Big{)},

where $\tilde{\sigma}=\sqrt{1+(\bm{x}^{2})^{T}\bm{\sigma}^{2}}$ , $v(z):=\frac{\mathcal{N}(z|0,1)}{\Phi(z)}$ , $w(z):=v(z)\Big{(}v(z)+z\Big{)}$ and $\bm{x}^{2}$ is understood as the column vector containing $x_{i}^{2}$ , so that $S^{n+1}=T^{\text{ADF}}(S^{n},\bm{x},Y^{n+1})$ .

In a dynamic program, the value function is defined as the value of the optimal policy given a particular state $S^{n}$ at time $n$ , and may also be determined recursively through Bellman’s equation. In the case of stochastic binary feedback, the terminal value function $V^{N}:\mathcal{S}\mapsto\mathbb{R}$ is given by (1) as

V^{N}(s)=\max_{\bm{x}}p(y=+1|\bm{x},s),\forall s\in\mathcal{S}.

The dynamic programming principle tells us that the value function at any other time $n=1,\dots,N$ , $V^{n}$ , is given recursively by

V^{n}(s)=\max_{\bm{x}}\mathbb{E}[V^{n+1}(T(s,\bm{x},Y^{n+1}))|\bm{x},s],\forall s\in\mathcal{S}.

Since the curse of dimensionality on the state space $\mathcal{S}$ makes direct computation of the value function intractable, computationally efficient approximate policies need to be considered. A computationally attractive policy for ranking and selection problems is known as the knowledge gradient (KG) [16], which will be extended to handle Bayesian classification models in the next section.

5.2 Knowledge Gradient for Binary Responses

The knowledge gradient of measuring an alternative $\bm{x}$ can be defined as follows:

Definition 5.3.

The knowledge gradient of measuring an alternative $\bm{x}$ while in state $s$ is

\nu_{\bm{x}}^{\text{KG}}(s):=\mathbb{E}\Big{[}V^{N}\Big{(}T(s,\bm{x},Y)\Big{)}-V^{N}(s)|\bm{x},s\Big{]}.

(15)

$V^{N}(s)$ is deterministic given $s$ and is independent of alternatives $\bm{x}$ . Since the label for alternative $\bm{x}$ is not known at the time of selection, the expectation is computed conditional on the current belief state $s=(\bm{\mu},\bm{\Sigma})$ . Specifically, given a state $s=(\bm{\mu},\bm{\Sigma})$ , the outcome $y$ of an alternative $\bm{x}$ is a random variable that follows from a Bernoulli distribution with a predictive distribution

\displaystyle p(y=+1|\bm{x},s)=\int p(y=+1|\bm{x},\bm{w})p(\bm{w}|s)\text{d}\bm{w}=\int\sigma(\bm{w}^{T}\bm{x})p(\bm{w}|s)\text{d}\bm{w}.

(16)

We can calculate the expected value in the next state as

			$\displaystyle\mathbb{E}[V^{N}(T(s,\bm{x},y))]$
		$\displaystyle=$	$\displaystyle p(y=+1\|\bm{x},s)V^{N}\Big{(}T(s,\bm{x},+1)\Big{)}+p(y=-1\|\bm{x},s)V^{N}\Big{(}T(s,\bm{x},-1)\Big{)}$
		$\displaystyle=$	$\displaystyle p(y=+1\|\bm{x},s)\cdot\max_{\bm{x}^{\prime}}p\Big{(}y=+1\|\bm{x}^{\prime},T(s,\bm{x},+1)\Big{)}$
			$\displaystyle+p(y=-1\|\bm{x},s)\cdot\max_{\bm{x}^{\prime}}p\Big{(}y=+1\|\bm{x}^{\prime},T(s,\bm{x},-1)\Big{)}.$

The knowledge gradient policy suggests at each time $n$ selecting the alternative that maximizes $\nu_{\bm{x}}^{\text{KG},n}(s^{n})$ where ties are broken randomly. Because of the errors incurred by approximation and numerical calculation, the tie should be understood as within $\epsilon$ -accuracy. The knowledge gradient policy can work with any choice of link function $\sigma(\cdot)$ and approximation procedures by adjusting the transition function $T(s,x,\cdot)$ accordingly. That is, $T^{\text{L}}$ or $T^{\text{ADF}}$ .

The predictive distribution $\int\sigma(\bm{w}^{T}\bm{x})p(\bm{w}|s)\text{d}\bm{w}$ is obtained by marginalizing with respect to the distribution specified by current belief state $p(\bm{w}|s)=\mathcal{N}(\bm{w}|\bm{\mu},\bm{\Sigma})$ . Denoting $a=\bm{w}^{T}\bm{x}$ and $\delta(\cdot)$ as the Dirac delta function, we have $\sigma(\bm{w}^{T}\bm{x})=\int\delta(a-\bm{w}^{T}\bm{x})\sigma(a)\text{d}a.$ Hence

\int\sigma(\bm{w}^{T}\bm{x})p(\bm{w}|s)\text{d}\bm{w}=\int\sigma(a)p(a)\text{d}a,

where $p(a)=\int\delta(a-\bm{w}^{T}\bm{x})p(\bm{w}|s)\text{d}\bm{w}.$ Since the delta function imposes a linear constraint on $\bm{w}$ and $p(\bm{w}|s)=\mathcal{N}(\bm{w}|\bm{\mu},\bm{\sigma}^{2})$ is Gaussian, the marginal distribution $p(a)$ is also Gaussian. We can evaluate $p(a)$ by calculating the mean and variance of this distribution [4]. We have

	$\displaystyle\mu_{a}$	$\displaystyle=$	$\displaystyle\mathbb{E}[a]=\int p(a)a\text{ d}a=\int p(\bm{w}\|s)\bm{w}^{T}\bm{x}\text{ d}\bm{w}=\bm{\mu}^{T}\bm{x},$
	$\displaystyle\sigma_{a}^{2}$	$\displaystyle=$	$\displaystyle\text{Var}[a]=\int p(\bm{w}\|s)\big{(}(\bm{w}^{T}\bm{x})^{2}-(\bm{\mu}^{T}\bm{x})^{2}\big{)}\text{ d}\bm{w}=\sum_{j=1}^{d}\sigma_{j}^{2}x_{j}^{2}.$

Thus $\int\sigma(\bm{w}^{T}\bm{x})p(\bm{w}|s)\text{d}\bm{w}=\int\sigma(a)p(a)\text{d}a=\int\sigma(a)\mathcal{N}(a|\mu_{a},\sigma^{2}_{a})\text{d}a.$

For probit function $\sigma(a)=\Phi(a)$ , the convolution of a Gaussian and a normal CDF can be evaluated analytically. Thus for probit regression, the predictive distribution can be solved exactly as $p(y=+1|\bm{x},s)=\Phi(\frac{\mu_{a}}{\sqrt{1+\sigma^{2}_{a}}}).$ Yet, the convolution of a Gaussian with a logistic sigmoid function cannot be evaluated analytically. We apply the approximation $\sigma(a)\approx\Phi(\alpha a)$ with $\alpha=\pi/8$ (see [3, 51]), leading to the following approximation for the convolution of a logistic sigmoid with a Gaussian

p(y=+1|\bm{x},s)=\int\sigma(\bm{w}^{T}\bm{x})p(\bm{w}|s)\text{d}\bm{w}\approx\sigma(\kappa(\sigma^{2}_{a})\mu_{a}),

where $\kappa(\sigma^{2})=(1+\pi\sigma^{2}/8)^{-1/2}$ .

Because of the one-step look ahead, the KG calculation can also benefit from the online recursive update of the belief either from ADF or from online Bayesian linear classification based on Laplace approximation. We summarize the decision rules of the knowledge gradient policy at each iteration under different sigmoid functions and different approximation methods in Algorithm 2, 3, and 4, respectively.

input :

m_{j}

q_{j}

(Each weight

w_{j}

has an independent prior

\mathcal{N}(w_{j}|m_{j},q_{j}^{-1})

)

for $\bm{x}$ in $\mathcal{X}$ do

Let

\Psi(\bm{w},y)=-\frac{1}{2}\sum_{j=1}^{d}q_{j}(w_{j}-m_{j})^{2}-\log(1+\exp(-y\bm{w}^{T}\bm{x}))

Use bisection method to find

\hat{\bm{w}}_{+}=\arg\max_{\bm{w}}\Psi(\bm{w},+1),~\hat{\bm{w}}_{-}=\arg\max_{\bm{w}}\Psi(\bm{w},-1)

\mu=\bm{m}^{T}\bm{x}

\sigma^{2}=\sum_{j=1}^{d}q_{j}^{-1}x_{j}^{2}

Let

\sigma(a)=\big{(}1+\exp(-a)\big{)}^{-1}

\mu_{+}(\bm{x}^{\prime}):=\hat{\bm{w}}_{+}^{T}\bm{x}^{\prime}

\mu_{-}(\bm{x}^{\prime}):=\hat{\bm{w}}_{-}^{T}\bm{x}^{\prime}

\sigma^{2}_{\pm}(\bm{x}^{\prime}):=\sum_{j=1}^{d}\Big{(}q_{j}+\sigma(\hat{\bm{w}}^{T}_{\pm}\bm{x})(1-\sigma(\hat{\bm{w}}_{\pm}^{T}\bm{x}))x^{2}_{j}\Big{)}^{-1}(x^{\prime}_{j})^{2}

\tilde{\nu}_{\bm{x}}=\sigma(\kappa(\sigma^{2})\mu)\cdot\max_{\bm{x^{\prime}}}\sigma\big{(}\kappa(\sigma_{+}^{2}(\bm{x}^{\prime}))\mu_{+}(\bm{x^{\prime}})\big{)}+\sigma(-\kappa(\sigma^{2})\mu)\cdot\max_{\bm{x^{\prime}}}\sigma\big{(}\kappa(\sigma_{-}^{2}(\bm{x}^{\prime}))\mu_{-}(\bm{x^{\prime}})\big{)}

end for

\bm{x}^{\text{KG}}\in\arg\max_{\bm{x}}\tilde{\nu}_{x}

Algorithm 2 Knowledge Gradient Policy for Logistic Model based on Laplace approximation

input :

m_{j}

q_{j}

(Each weight

w_{j}

has an independent prior

\mathcal{N}(w_{j}|m_{j},q_{j}^{-1})

)

for $\bm{x}$ in $\mathcal{X}$ do

Let

\Psi(\bm{w},y)=-\frac{1}{2}\sum_{j=1}^{d}q_{i}(w_{i}-m_{i})^{2}+\log(\Phi(y\bm{w}^{T}\bm{x}))

Use bisection method to find

\hat{\bm{w}}_{+}=\arg\max_{\bm{w}}\Psi(\bm{w},+1),\hat{\bm{w}}_{-}=\arg\max_{\bm{w}}\Psi(\bm{w},-1)

\mu=\bm{m}^{T}\bm{x}

\sigma^{2}=\sum_{j=1}^{d}q_{j}^{-1}x_{j}^{2}

\mu_{+}(\bm{x}^{\prime}):=\hat{\bm{w}}_{+}^{T}\bm{x}^{\prime}

\mu_{-}(\bm{x}):=\hat{\bm{w}}_{-}^{T}\bm{x}^{\prime}

\sigma^{2}_{\pm}(\bm{x}^{\prime}):=\sum_{j=1}^{d}\Big{(}q_{j}+(\frac{\mathcal{N}(\hat{\bm{w}}^{T}_{\pm}\bm{x}|0,1)^{2}}{\Phi(\hat{\bm{w}}^{T}_{\pm}\bm{x})^{2}}+\frac{\hat{\bm{w}}^{T}_{\pm}\bm{x}\mathcal{N}(\hat{\bm{w}}^{T}_{\pm}\bm{x}|0,1)}{\Phi(\hat{\bm{w}}^{T}_{\pm}\bm{x})})x^{2}_{j}\Big{)}^{-1}(x^{\prime}_{j})^{2}

\tilde{\nu}_{\bm{x}}=\Phi(\frac{\mu}{\sqrt{1+\sigma^{2}}})\cdot\max_{\bm{x^{\prime}}}\Phi(\frac{\mu_{+}(\bm{x}^{\prime})}{\sqrt{1+\sigma_{+}^{2}(\bm{x}^{\prime})}})+\Phi(-\frac{\mu}{\sqrt{1+\sigma^{2}}})\cdot\max_{\bm{x^{\prime}}}\Phi(\frac{\mu_{-}(\bm{x}^{\prime})}{\sqrt{1+\sigma_{-}^{2}(\bm{x}^{\prime})}})

end for

\bm{x}^{\text{KG}}\in\arg\max_{\bm{x}}\tilde{\nu}_{x}

Algorithm 3 Knowledge Gradient Policy for Probit Model based on Laplace approximation

input :

m_{j}

q_{j}

(Each weight

w_{j}

has an independent prior

\mathcal{N}(w_{j}|m_{j},q_{j}^{-1})

)

\sigma_{j}^{2}=1/q_{j}

for $\bm{x}$ in $\mathcal{X}$ do

Define

v(z):=\frac{\mathcal{N}(z|0,1)}{\Phi(z)}\text{ and }w(z):=v(z)\Big{(}v(z)+z\Big{)}

\mu=\bm{m}^{T}\bm{x}

\sigma^{2}=\sum_{j=1}^{d}\sigma_{j}^{2}x_{j}^{2}

m_{+j}=m_{j}+\frac{x_{j}\sigma_{j}^{2}}{\sqrt{1+\sigma^{2}}}v(\frac{\bm{m}^{T}\bm{x}}{\sqrt{1+\sigma^{2}}})

m_{-j}=m_{j}-\frac{x_{j}\sigma_{j}^{2}}{\sqrt{1+\sigma^{2}}}v(-\frac{\bm{m}^{T}\bm{x}}{\sqrt{1+\sigma^{2}}})

\sigma_{+j}^{2}=\sigma_{j}^{2}-\frac{x_{j}^{2}\sigma_{j}^{4}}{\sqrt{1+\sigma^{2}}^{2}}w(\frac{\bm{m}^{T}\bm{x}}{\sqrt{1+\sigma^{2}}})

\sigma_{-j}^{2}=\sigma_{j}^{2}-\frac{x_{j}^{2}\sigma_{j}^{4}}{\sqrt{1+\sigma^{2}}^{2}}w(-\frac{\bm{m}^{T}\bm{x}}{\sqrt{1+\sigma^{2}}})

\mu_{+}(\bm{x}^{\prime}):=\bm{m}_{+}^{T}\bm{x}^{\prime}

\mu_{-}(\bm{x}):=\bm{m}_{-}^{T}\bm{x}^{\prime}

\sigma_{+}^{2}(\bm{x}^{\prime}):=\sum_{j=1}^{d}\sigma_{+j}^{2}(x^{\prime}_{j})^{2}

\sigma_{-}^{2}(\bm{x}^{\prime}):=\sum_{j=1}^{d}\sigma_{-j}^{2}(x^{\prime}_{j})^{2}

\tilde{\nu}_{\bm{x}}=\Phi(\frac{\mu}{\sqrt{1+\sigma^{2}}})\cdot\max_{\bm{x^{\prime}}}\Phi(\frac{\mu_{+}(\bm{x}^{\prime})}{\sqrt{1+\sigma_{+}^{2}(\bm{x}^{\prime})}})+\Phi(-\frac{\mu}{\sqrt{1+\sigma^{2}}})\cdot\max_{\bm{x^{\prime}}}\Phi(\frac{\mu_{-}(\bm{x}^{\prime})}{\sqrt{1+\sigma_{-}^{2}(\bm{x}^{\prime})}})

end for

\bm{x}^{\text{KG}}\in\arg\max_{\bm{x}}\tilde{\nu}_{x}

Algorithm 4 Knowledge Gradient Policy for Probit Model based on assumed density filtering

We next present the following finite-time bound on the the mean squared error (MSE) of the estimated weight for Bayesian logistic regression. The proof is in the supplement. Since the learner plays an active role in selecting the measurements, the bound does not make the i.i.d. assumption of the examples which differs from the PAC learning bound. Without loss of generality, we assume $\|\bm{x}\|_{2}\leq 1$ , $\forall{\bm{x}\in\mathcal{X}}$ . In our finite-time analysis, we would like to highlight the difference between the ideal n-step estimate $\bm{w}^{n}=\arg\max_{\bm{w}}\Psi(\bm{w}|\bm{m}^{0},\Sigma^{0})$ with the prior distribution $p(\bm{w}^{*})=\mathcal{N}(\bm{w}^{*}|\bm{m}^{0},\bm{\Sigma}^{0})$ , and its approximation $\hat{\bm{w}}^{n}$ . The approximation error may come from the Laplace approximation we adopt when an explicit form is not available, or from accumulation of numerical error.

Theorem 5.4.

Let $\mathcal{D}^{n}$ be the $n$ measurements produced by the KG policy. Then with probability $P_{d}(M)$ , the expected error of $\bm{w}^{n}$ is bounded as

\mathbb{E}_{\bm{y}\sim\mathcal{B}(\mathcal{D}^{n},\bm{w}^{*})}||\bm{w}^{n}-\bm{w}^{*}||_{2}\leq\frac{C_{min}+\lambda_{min}\big{(}\bm{\Sigma}^{-1}\big{)}}{2},

where the distribution $\mathcal{B}(\mathcal{D}^{n},\bm{w}^{*})$ is the vector Bernoulli distribution $Pr(y^{i}=+1)=\sigma((\bm{w}^{*})^{T}\bm{x}^{i})$ , $P_{d}(M)$ is the probability of a d-dimensional standard normal random variable appears in the ball with radius $M=\frac{1}{16}\frac{\lambda_{min}^{2}}{\sqrt{\lambda_{max}}}$ and $C_{min}=\lambda_{min}\Big{(}\frac{1}{n}\sum_{i=1}\sigma((\bm{w}^{*})^{T}\bm{x}^{i})\big{(}1-\sigma((\bm{w}^{*})^{T}\bm{x}^{i})\big{)}\bm{x}^{i}(\bm{x}^{i})^{T}\Big{)}.$ The same finite time bound holds for $\hat{\bm{w}}^{n}$ as long as the approximation satisfies

\Psi(\hat{\bm{w}}^{n}|\bm{m}^{0},\Sigma^{0})\geq\Psi(\bm{w}^{n}|\bm{m}^{0},\Sigma^{0})-\frac{1}{32}\lambda^{3}_{min}\big{(}\bm{\Sigma}^{-1}\big{)}.

(17)

In the special case where $\bm{\Sigma}^{0}=\lambda^{-1}\bm{I}$ , we have $\lambda_{max}=\lambda_{min}=\lambda$ and $M=\frac{\lambda^{3/2}}{16}$ . The bound holds with higher probability $P_{d}(M)$ with larger $\lambda$ . This is natural since a larger $\lambda$ represents a normal distribution with narrower bandwidth, resulting in a more concentrated $\bm{w}^{*}$ around $\bm{m}^{0}$ .

Condition (17) quantifies how well the approximation should be carried out in order for the bound to hold. It places a condition on the optimality of the log-likelihood value instead of the distance to the global maximizer. This is a friendly condition since the dependence of the maximizer over the objective function is generally not continuous. The approximation error could come from either the Laplace approximation or the numerical calculation in practice, which calls for experiments to further verity the performance of the algorithm. Indeed, we conduct extensive experiments in Section 6 to demonstrate the behaviors of the proposed algorithm together with bench mark methods.

If the goal of the learner is to maximize the cumulative successes as in bandit settings, rather than finding the final recommendation after the offline training phases, the online knowledge gradient policy can be modified by deleting $X^{\text{KG},n}$ at each time step $n$ as:

X^{\text{KG},n}(S^{n})=\arg\max_{\bm{x}}p(y=+1|x,S^{n})+\tau\nu_{\bm{x}}^{\text{KG}}(S^{n}),

(18)

where $\tau$ reflects a planning horizon.

5.3 Behavior and Asymptotic Optimality

In this section, we study theoretically the behavior of the KG policy, especially in the limit as the number of measurements $N$ grows large. For the purposes of the theoretical analysis, we do not approximate the predictive posterior distribution. We use logistic function as the sigmoid link function throughout this section. Yet the theoretical results can be generalized to other link functions. We begin by showing the positive value of information (benefits of measurement).

Proposition 5.5 (Benefits of measurement).

The knowledge gradient of measuring any alternative $\bm{x}$ while in any state $s\in\mathcal{S}$ is nonnegative, $\mu_{\bm{x}}^{\text{KG}}(s)\geq 0.$ The state space $\mathcal{S}$ is the space of all possible predictive distributions for $\bm{w}$ .

The next lemma shows that the asymptotic probability of success of each alternative is well defined.

Lemma 5.6.

For any alternative $\bm{x}$ , $p_{\bm{x}}^{n}$ converges almost surely to a random variable $p_{\bm{x}}^{\infty}=\mathbb{E}\big{[}\sigma(\bm{w}^{T}\bm{x})|\mathcal{F}^{\infty}\big{]}$ , where $p_{\bm{x}}^{n}$ is short hand notation for $p(y=+1|\bm{x},\mathcal{F}^{n})=\mathbb{E}\big{[}\sigma(\bm{w}^{T}\bm{x})|\mathcal{F}^{n}\big{]}.$

Proof.

ProofSince $|\sigma(\bm{w}^{T}\bm{x})|\leq 1$ , the definition of $p_{\bm{x}}^{n}$ implies that $p_{\bm{x}}^{n}$ is a bounded martingale and hence converges. ∎

The rest of this section shows that this limiting probability of success of each alternative is one in which the posterior is consistent and thus the KG is asymptotically optimal. We also show that as the number of measurements grows large, the maximum likelihood estimator $\bm{w}_{\text{MLE}}$ based on the alternatives measured by the KG policy is consistent and asymptotically normal. The next proposition states that if we have measured an alternative infinitely many times, there is no benefit to measuring it one more time. This is a key step for establishing the consistency of the KG policy and the MLE. The proof is similar to that by [15] with additional mathematical steps for Bernoulli distributed random variables. See Appendix 11.2 for details.

Proposition 5.7.

If the policy $\pi$ measures alternative $\bm{x}$ infinitely often almost surely, then the value of information of that alternative $\nu_{\bm{x}}(\mathcal{F}^{\infty})=0$ almost surely under policy $\pi$ .

Without loss of generality, we assume $\|\bm{x}\|_{2}\leq 1$ and that the $d\times d$ matrix $P$ formed by $(\bm{x}_{1},\bm{x}_{2},\dots,\bm{x}_{d})$ is invertible. The next theorem states the strong consistency and asymptotic normality of the maximum likelihood estimator $\bm{w}^{n}_{\text{MLE}}$ (e.g. with $\lambda=0$ ) based on KG’s sequential selection of alternatives by verification of the following regularity conditions:

(C₁)

The exogenous variables are uniformly bounded.
(C₂)

Let $\lambda_{1n}$ and $\lambda_{dn}$ be respectively the smallest and the largest eigenvalue of the Fisher information of the first $n$ observations $\bm{F}_{n}(\bm{w}^{*})$ . There exists $C$ such that $\lambda_{kN}/\lambda_{1N}<C$ , $\forall N$ .
(C₃)

The smallest eigenvalue of the Fisher information is divergent, $\lambda_{1N}\rightarrow+\infty.$

Theorem 5.8.

The sequence of $\bm{w}^{n}_{\text{MLE}}$ based on KG’s sequential selection of alternatives converges almost surely to the true value $\bm{w}^{*}$ and is asymptotically normal:

\bm{F}_{n}^{\frac{1}{2}}(\bm{w}^{n}_{\text{MLE}}-\bm{w}^{*})\xrightarrow{d}\mathcal{N}(0,\bm{I}).

Proof.

Proof We first prove that for any alternative $\bm{x}$ , it will be measured infinitely many times almost surely. We will prove it by contradiction. If this is not the case, then there exists a time $T$ such that for any $n>T$ ,

\mu^{\text{KG},n}_{\bm{x}}<\max_{\bm{x}\in\mathcal{X}}\mu^{\text{KG},n}_{\bm{x}}-\epsilon.

(19)

This is because otherwise the difference between the KG value of $\bm{x}$ and the maximum KG value will be smaller than $\epsilon$ for infinitely many times, and thus the probability of not measuring $\bm{x}$ after $T$ will be 0.

In addition, since the KG value is always non-negative, we have $\max_{\bm{x}\in\mathcal{X}}\mu^{\text{KG},n}_{\bm{x}}>\epsilon$ for each $n>T$ . Notice that $\mathcal{X}$ is a finite set, then it follows that there exists an alternative $\bm{x}^{\prime}$ such that the following happens infinitely many times:

\mu^{\text{KG},n}_{\bm{x}^{\prime}}=\max_{\bm{x}\in\mathcal{X}}\mu^{\text{KG},n}_{\bm{x}}.

(20)

Therefore, $\bm{x}^{\prime}$ is measured infinitely many times almost surely. However, we have proved in Proposition 5.7 that $\mu^{\text{KG},n}_{\bm{x}^{\prime}}$ goes to 0 zero as long as we measure $\bm{x}^{\prime}$ infinitely many times, which contradicts (20). The contradiction shows that our original assumption that $\bm{x}$ only being measured finite times is incorrect. As a consequence, we have prove that, almost surely, $\bm{x}$ will be measured infinitely many times.

Since our proof is for arbitrary $\bm{x}$ , we actually proved that every alternative will be measured infinitely many times, which immediately leads to $\max_{\bm{x}\in\mathcal{X}}\mu^{\text{KG},n}_{\bm{x}}\rightarrow 0$ . Therefore the algorithm will eventually pick the alternative uniformly at random. Hence it satisfies the conditions (C₂) and (C₃), leading to the strong consistency and asymptotic normality [21, 13, 25, 11].

In particular, we can prove that the smallest eigenvalue of the Fisher’s information goes to infinity in a simple way. Without loss of generality, assume that the alternatives $\bm{x}_{1},\bm{x}_{2},\dots,\bm{x}_{d}$ are in general linear position, which means that the $d\times d$ matrix $P:=(\bm{x}_{1},\bm{x}_{2},\dots,\bm{x}_{d})$ is invertible. We use $\mathcal{X}^{\prime}$ to denote the set of these $d$ alternatives and denote $P^{-1}$ by $Q$ . Then $Q\bm{x}_{t}\bm{x}_{t}^{T}Q^{T}$ is a matrix whose every element equals 0 except that its $k$ -th diagonal element equals one.

We use $\bm{F}_{n}$ to denote the Fisher’s information at time $n$ . Then from the fact that the matrix of type $\bm{x}\bm{x}^{T}$ is positive semi-definite, it is straightforward to see that

\bm{F}_{n}(\bm{w})\geq\min_{\bm{x}\in\mathcal{X}}(1-\sigma(\bm{x}^{T}\bm{w}))\sigma(\bm{x}^{T}\bm{w})\sum_{1\leq t\leq n,\bm{x}_{t}\in\mathcal{X}^{\prime}}\bm{x}_{t}\bm{x}_{t}^{T},

where the constant $\min_{\bm{x}\in\mathcal{X}}(1-\sigma(\bm{x}^{T}\bm{w}))\sigma(\bm{x}^{T}\bm{w})>0$ since $\mathcal{X}$ is a finite set. Now define matrix $R_{n}$ as

R_{n}:=\sum_{1\leq t\leq n,\bm{x}_{t}\in\mathcal{X}^{\prime}}Q\bm{x}_{t}\bm{x}_{t}^{T}Q^{T},

which is a diagonal matrix whose $i$ -th element equals the times that $\bm{x}_{i}$ is estimated. Since $\bm{x}_{k}$ is measured infinitely many times for $1\leq k\leq d$ , then each diagonal element of $R_{n}$ goes to infinity. Now notice that $\sum_{1\leq t\leq n,\bm{x}_{t}\in\mathcal{X}^{\prime}}\bm{x}_{t}\bm{x}_{t}^{T}=Q^{-1}R_{n}(Q^{-1})^{T}$ is congruent to $R_{n}$ and $Q^{-1}$ is a constant matrix, then it follows that the smallest eigenvalue of $\sum_{1\leq t\leq n,\bm{x}_{t}\in\mathcal{X}^{\prime}}\bm{x}_{t}\bm{x}_{t}^{T}$ , and hence the smallest eigenvalue of $\bm{F}_{n}$ , goes to infinity.

Similarly, by using the congruent argument, we only need to show that almost surely for any $k$ , there exists a constant $C$ such that $n_{k}(t)/t\geq C$ , where $n_{k}(t)$ is the time that $n_{k}(t)$ is measured in the first $t$ measurements. Since the algorithm eventually selects alternatives uniformly at random, without loss of generality, we assume $n_{k}(t)$ is the num of $t$ i.i.d Bernoulli random variable with probability $1/m$ , then it is clear to check that $\mathbb{E}[n_{k}(t)]=t/m$ , and $Var[n_{k}(t)]=t(m-1)/m^{2}$ . Let $L$ be the event that $\underline{\lim}_{t\rightarrow\infty}n_{k}(t)/t<1/2m$ . Notice that $L$ implies that there exists a sequence $t_{1}<t_{2}<t_{3}<\dots$ such that $n_{k}(t_{s})/t_{s}<1/2m$ . However, by Chebyshev’s inequality, we have

\mathbb{P}\left[\frac{n_{k}(t_{s})}{t_{s}}<\frac{1}{2m}\right]\leq\mathbb{P}\left[\left|n_{k}(t_{s})-\frac{t_{s}}{m}\right|\geq\frac{t_{s}}{2m}\right]\leq\frac{4(m-1)}{t_{s}}.

Notice that the above inequality holds for any $t_{s}$ , which means $\mathbb{P}[L]=0$ and thus $\underline{\lim}_{t\rightarrow\infty}n_{k}(t)/t\geq 1/2m$ almost surely, which completes the proof. ∎

After establishing the consistency and asymptotic normality for $\bm{w}_{\text{MLE}}^{n}$ , for any $\lambda>0$ as the inverse of the variance of the prior Gaussian distribution, the asymptotic bias of the estimator $\bm{w}^{n}_{\lambda}$ is as follows [36]:

\bm{E}[\bm{w}^{n}_{\lambda}-\bm{w}^{*}]=-2\lambda\{\bm{F}_{n}(\bm{w}^{*})+2\lambda\bm{I}\}^{-1}\bm{w}^{*},

and the asymptotic variance of $\bm{w}^{n}_{\lambda}$ is $\{\bm{F}_{n}(\bm{w}^{*})+2\lambda\bm{I}\}^{-1}\bm{F}_{n}(\bm{w}^{*})\{\bm{F}_{n}(\bm{w}^{*})+2\lambda\bm{I}\}^{-1}.$

Finally, we show in the next theorem that given the opportunity to measure infinitely often, for any given neighborhood of $\bm{w}^{*}$ , the probability that the posterior distribution lies in this neighborhood goes to 1 and KG will discover which one is the true best alternative. The detailed proof can be found in Appendix 11.3.

Theorem 5.9 (Asymptotic optimality).

For any true parameter value $\bm{w}^{*}\in\mathbb{R}^{d}$ and any normal prior distribution with positive definite covariance matrix, under the knowledge gradient policy, the posterior is consistent and the KG policy is asymptotically optimal: $\arg\max_{\bm{x}}\mathbb{E}^{\text{KG}}[\sigma(\bm{w}^{T}\bm{x})|\mathcal{F}^{\infty}]=\arg\max_{\bm{x}}\sigma((\bm{w^{*}})^{T}\bm{x})$ .

5.4 Knowledge Gradient for contextual bandits

In many of the applications, the “success” or “failure” can be predicted though an unknown relationship that depends on a partially controllable vector of attributes for each instance. For example, the patient attributes are given which we do not have control over. We can only choose a medical decision, and the outcome is based on both the patient attributes and the selected medical decision.

This problem is known as a contextual bandit [34, 37, 33]. At each round $n$ , the learner is presented with a context vector $\bm{c}_{n}$ (patient attributes) and a set of actions $\bm{a}\in\mathcal{A}$ (medical decisions). After choosing an alternative $x$ , we observe an outcome $y^{n+1}$ . The goal is to find a policy that selects actions such that the cumulative reward is as large as possible over time. In contrast, the previous section considers a context-free scenario with the goal of maximizing the probability of success after the offline training phase so that the error incurred during the training is not punished.

Since there are very few patients with the same characteristics, it is unlikely that we can work with one type of patients that have the same characteristics to produce statistically reliable performance measurements on medical outcomes. The contextual bandit formulation provides the potential to use parametric belief models that allow us to learn relationships across a wide range of patients and medical decisions and personalize health care based on the patient attributes.

Each of the past observations are made of triplets $(\bm{c}^{n},\bm{a}^{n},y^{n+1})$ . Under Bayesian linear classification belief models, the binomial outcome $p(y=+1|\bm{c},\bm{a})$ is predicted though $\sigma\big{(}F(\bm{c},\bm{a})\big{)}$ ,

F(\bm{c},\bm{a})=\theta+\bm{\alpha}^{T}\bm{\varphi_{\bm{c}}}+\bm{\beta}^{T}\bm{\psi_{\bm{a}}},

(21)

where suppose each context $\bm{c}$ and action $\bm{x}$ is represented by feature vectors $\bm{\varphi_{\bm{c}}}$ and $\bm{\psi_{\bm{a}}}$ , respectively. At each round $n$ , the model updates can be slightly modified based on the observation triplet $(\bm{c}^{n},\bm{a}^{n},y^{n+1})$ by treating $1||\bm{\varphi_{\bm{c}}}||\bm{\psi_{\bm{a}}}$ as the alternative $\bm{x}$ , where $\bm{u}||\bm{v}$ denotes the concatenation of the two vectors $\bm{u}$ and $\bm{v}$ .

The knowledge gradient $\nu_{\bm{a}|\bm{c}}^{\text{KG}}(s)$ of measuring an action $\bm{a}$ given context $\bm{c}$ can be defined as follows:

Definition 5.10.

The knowledge gradient of measuring an action $\bm{a}$ given a context $\bm{c}$ while in state $s$ is

\nu_{\bm{a}|\bm{c}}^{\text{KG}}(s):=\mathbb{E}\Big{[}V^{N}\Big{(}T\big{(}s,1||\bm{\varphi_{\bm{c}}}||\bm{\psi_{\bm{a}}},y\big{)}\Big{)}-V^{N}(s)|\bm{c},\bm{a},s\Big{]}.

(22)

The calculation of $\nu_{\bm{a}|\bm{c}}^{\text{KG}}(s)$ can be modified based on Section 5.2 by replacing $\bm{x}$ as $1||\bm{\varphi_{\bm{c}}}||\bm{\psi_{\bm{a}}}$ throughout.

Since the objective in the contextual bandit problems is to maximize the cumulative number of successes, the knowledge gradient policy developed in Section 5.2 for stochastic binary feedback can be easily extended to online learning following the “stop-learning” (SL) policy adopted by [47]. The action $X^{\text{KG},n}_{\bm{c}}(s^{n})$ that is chosen by KG at time $n$ given a context $\bm{c}$ and a knowledge state $s^{n}$ can be obtained as:

X^{\text{KG},n}_{\bm{c}}(S^{n})=\arg\max_{\bm{a}}p(y=+1|\bm{c},\bm{a},S^{n})+\tau\nu_{\bm{a}|\bm{c}}^{\text{KG},n}(S^{n}),

where $\tau$ reflects a planning horizon.

It should be noted that if the context are constant over time, $X^{\text{KG},n}_{\bm{c}}(S^{n})$ is degenerated to the case of context-free multi-armed bandit problems as discussed in Eq. (18).

6 Experimental results

In this section, we evaluate the proposed method in offline learning settings where we are not punished for errors incurred during training and only concern with the final recommendation after the offline training phases.

We experiment with both synthetic datasets and the UCI machine learning repository [38] which includes classification problems drawn from settings including sonar, glass identification, blood transfusion, survival, breast cancer (wpbc), planning relax and climate model failure. We first analyze the behavior of the KG policy and then compare it to state-of-the-art learning algorithms. On synthetic datasets, we randomly generate a set of $M$ $d$ -dimensional alternatives $\bm{x}$ from $[-3,3]$ . At each run, the stochastic binary labels are simulated using a $d+1$ -dimensional weight vector $\bm{w}^{*}$ which is sampled from the prior distribution $w^{*}_{i}\sim\mathcal{N}(0,\lambda)$ . The $+1$ label for each alternative $\bm{x}$ is generated with probability $\sigma(w^{*}_{0}+\sum_{j=1}^{d}w^{*}_{d}x_{d})$ . For each UCI dataset, we use all the data points as the set of alternatives with their original attributes. We then simulate their labels using a weight vector $\bm{w}^{*}$ . This weight vector could have been chosen arbitrarily, but it was in fact a perturbed version of the weight vector trained through logistic regression on the original dataset.

6.1 Behavior of the KG policy

Refer to caption — Figure 1: The scatter plots illustrate the KG values at 1-4 iterations from left to right with both the color and the size reflecting the magnitude. The star, the red square and pink circle indicate the true best alternative, the alternative to be selected and the implementation decision, respectively.

To better understand the behavior of the KG policy, Fig. 1 shows the snapshot of the KG policy at each iteration on a $2$ -dimensional synthetic dataset and a $3$ -dimensional synthetic dataset in one run. Fig. 1 shows the snapshot on a 2-dimensional dataset with 200 alternatives. The scatter plots show the KG values with both the color and the size of the point reflecting the KG value of the corresponding alternative. The star denotes the true alternative with the largest response. The red square is the alternative with the largest KG value. The pink circle is the implementation decision that maximizes the response under current estimation of $\bm{w}^{*}$ if the budget is exhausted after that iteration.

It can be seen from the figure that the KG policy finds the true best alternative after only three measurements, reaching out to different alternatives to improve its estimates. We can infer from Fig. 1 that the KG policy tends to choose alternatives near the boundary of the region. This criterion is natural since in order to find the true maximum, we need to get enough information about $\bm{w}^{*}$ and estimate well the probability of points near the true maximum which appears near the boundary. On the other hand, in a logistic model with labeling noise, a data $\bm{x}$ with small $\bm{x}^{T}\bm{x}$ inherently brings little information as pointed out by [61]. For an extreme example, when $\bm{x}=\bm{0}$ the label is always completely random for any $\bm{w}$ since $p(y=+1|\bm{w},\bm{0})\equiv 0.5$ . This is an issue when perfect classification is not achievable. So it is essential to label a data with larger $\bm{x}^{T}\bm{x}$ that has the most potential to enhance its confidence non-randomly.

Fig. 2 illustrates the snapshots of the KG policy on a 3-dimensional synthetic dataset with 300 alternatives. It can be seen that the KG policy finds the true best alternative after only 10 measurements. This set of plots also verifies our statement that the KG policy tends to choose data points near the boundary of the region.

Also depicted in Fig. 3 is the absolute class distribution error of each alternative, which is the absolute difference between the predictive probability of class $+1$ under current estimate and the true probability after $6$ iterations for the 2-dimensional dataset and 10 iterations for the 3-dimensional dataset. We see that the probability at the true maximum is well approximated, while moderate error in the estimate is located away from this region of interest.

6.2 Comparison with other policies

Recall that our goal is to maximize the expected response of the implementation decision. We define the Opportunity Cost (OC) metric as the expected response of the implementation decision $\bm{x}^{N+1}:=\arg\max_{\bm{x}}p(y=+1|\bm{x},\bm{w}^{N})$ compared to the true maximal response under weight $\bm{w}^{*}$ :

\text{OC}:=\max_{\bm{x}\in\mathcal{X}}p(y=+1|\bm{x},\bm{w}^{*})-p(y=+1|\bm{x}^{N+1},\bm{w}^{*}).

Note that the opportunity cost is always non-negative and the smaller the better. To make a fair comparison, on each run, all the time- $N$ labels of all the alternatives are randomly pre-generated according to the weight vector $\bm{w}^{*}$ and shared across all the competing policies. We allow each algorithm to sequentially measure N = 30 alternatives.

We compare with the following state-of-the-art active learning and Bayesian optimization policies that are compatible with logistic regression: Random sampling (Random), a myopic method that selects the most uncertain instance each step (MostUncertain), discriminative batch-mode active learning (Disc) [23] with batch size equal to 1, expected improvement (EI) [53] with an initial fit of 5 examples and Thompson sampling (TS) [7]. Besides, as upper confidence bounds (UCB) methods are often considered in bandit and optimization problems, we compare against UCB on the latent function $\bm{w}^{T}\bm{x}$ (UCB) [37] with $\alpha$ tuned to be 1. All the state transitions are based on the online Bayesian logistic regression framework developed in Section 4, while different policies provides different rules for measurement decisions at each iteration. The experimental results are shown in figure 4. In all the figures, the x-axis denotes the number of measured alternatives and the y-axis represents the averaged opportunity cost averaged over 100 runs.

It is demonstrated in Fig. 4 that KG outperforms the other policies in most cases, especially in early iterations, without requiring a tuning parameter. As an unbiased selection procedure, random sampling is at least a consistent algorithm. Yet it is not suitable for expensive experiments where one needs to learn the most within a small experimental budget. MostUncertain and Disc perform quite well on some datasets while badly on others. A possible explanation is that the goal of active leaning is to learn a classifier which accurately predicts the labels of new examples so their criteria are not directly related to maximize the probability of success aside from the intent to learn the prediction. After enough iterations when active learning methods presumably have the ability to achieve a good estimator of $\bm{w}^{*}$ , their performance will be enhanced. Thompson sampling works in general quite well as reported in other literature [7]. Yet, KG has a better performance especially during the early iterations. In the case when an experiment is expensive and only a small budget is allowed, the KG policy, which is designed specifically to maximize the response, is preferred.

We also note that KG works better than EI in most cases, especially in Fig. 4(b), 4(c) and 4(e). Although both KG and EI work with the expected value of information, when EI decides which alternative to measure, it is based on the expected improvement over current predictive posterior distribution while ignoring the potential change of the posterior distribution resulting from the next (stochastic) outcome $y$ . In comparison, KG considers an additional level of expectation over the random (since at the time of decision, we have not yet observed outcome) binary output $y$ .

Finally, KG, EI and Thompson sampling outperform the naive use of UCB policies on the latent function $\bm{w}^{T}\bm{x}$ due to the errors in the variance introduced by the nonlinear transformation. At each time step, the posterior of $\log{\frac{p}{1-p}}$ is approximated as a Gaussian distribution. An upper confidence bound on $\log{\frac{p}{1-p}}$ does not translate to one on $p$ with binary outcomes. In the meantime, KG, EI and Thompson sampling make decisions in the underlying binary outcome probability space and find the right balance of exploration and exploitation.

7 Conclusion

Motivated by real world applications, we consider binary classification problems where we have to sequentially run expensive experiments, forcing us to learn the most from each experiment. With a small budget of measurements, the goal is to learn the classification model as quickly as possible to identify the alternative with the highest probability of success. Due to the sequential nature of this problem, we develop a fast online Bayesian linear classifier for general response functions to achieve recursive updates. We propose a knowledge gradient policy using Bayesian linear classification belief models, for which we use different analytical approximation methods to overcome computational challenges. We further extend the knowledge gradient to the contextual bandit settings. We provide a finite-time analysis on the estimated error, and show that the maximum likelihood estimator based on the adaptively sampled points by the KG policy is consistent and asymptotically normal. We show furthermore that the knowledge gradient policy is asymptotic optimal. We demonstrate its efficiency through a series of experiments.

References

[1] Jean-Yves Audibert and Sébastien Bubeck. Best arm identification in multi-armed bandits. In COLT-23th Conference on Learning Theory-2010, pages 13–p, 2010.
[2] P. Auer, N. Cesa-Bianchi, and P. Fischer. Finite-time analysis of the multiarmed bandit problem. Mach. Learn., 47(2-3):235–256, 2002.
[3] David Barber and Christopher M Bishop. Ensemble learning for multi-layer networks. Advances in neural information processing systems, pages 395–401, 1998.
[4] Christopher M Bishop et al. Pattern recognition and machine learning, volume 4. springer New York, 2006.
[5] Xavier Boyen and Daphne Koller. Tractable inference for complex stochastic processes. In Proceedings of the Fourteenth conference on Uncertainty in artificial intelligence, pages 33–42. Morgan Kaufmann Publishers Inc., 1998.
[6] Sébastien Bubeck and Nicolò Cesa-Bianchi. Regret analysis of stochastic and nonstochastic multi-armed bandit problems. Foundations and Trends in Machine Learning, 2012.
[7] Olivier Chapelle and Lihong Li. An empirical evaluation of thompson sampling. In Advances in neural information processing systems, pages 2249–2257, 2011.
[8] Stephen E Chick. New two-stage and sequential procedures for selecting the best simulated system. Operations Research, 49(5):732–743, 2001.
[9] Wei Chu, Lihong Li, Lev Reyzin, and Robert E Schapire. Contextual bandits with linear payoff functions. In International Conference on Artificial Intelligence and Statistics, pages 208–214, 2011.
[10] Wei Chu, Martin Zinkevich, Lihong Li, Achint Thomas, and Belle Tseng. Unbiased online active learning in data streams. In Proceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 195–203. ACM, 2011.
[11] David Roxbee Cox and David Victor Hinkley. Theoretical statistics. CRC Press, 1979.
[12] M. H. DeGroot. Optimal Statistical Decisions. McGraw-Hill, 1970.
[13] Ludwig Fahrmeir and Heinz Kaufmann. Consistency and asymptotic normality of the maximum likelihood estimator in generalized linear models. The Annals of Statistics, pages 342–368, 1985.
[14] Sarah Filippi, Olivier Cappe, Aurélien Garivier, and Csaba Szepesvári. Parametric bandits: The generalized linear case. In Advances in Neural Information Processing Systems, pages 586–594, 2010.
[15] Peter Frazier, Warren Powell, and Savas Dayanik. The knowledge-gradient policy for correlated normal beliefs. INFORMS journal on Computing, 21(4):599–613, 2009.
[16] Peter I Frazier, Warren B Powell, and Savas Dayanik. A knowledge-gradient policy for sequential information collection. SIAM Journal on Control and Optimization, 47(5):2410–2439, 2008.
[17] David A Freedman. On the asymptotic behavior of bayes’ estimates in the discrete case. The Annals of Mathematical Statistics, pages 1386–1403, 1963.
[18] Yoav Freund, H Sebastian Seung, Eli Shamir, and Naftali Tishby. Selective sampling using the query by committee algorithm. Machine learning, 28(2-3):133–168, 1997.
[19] Subhashis Ghosal, Jayanta K Ghosh, RV Ramamoorthi, et al. Posterior consistency of dirichlet mixtures in density estimation. Ann. Statist, 27(1):143–158, 1999.
[20] Subhashis Ghosal and Anindya Roy. Posterior consistency of gaussian process prior for nonparametric binary regression. The Annals of Statistics, pages 2413–2429, 2006.
[21] Christian Gourieroux and Alain Monfort. Asymptotic properties of the maximum likelihood estimator in dichotomous logit models. Journal of Econometrics, 17(1):83–97, 1981.
[22] Thore Graepel, Joaquin Q Candela, Thomas Borchert, and Ralf Herbrich. Web-scale bayesian click-through rate prediction for sponsored search advertising in microsoft’s bing search engine. In Proceedings of the 27th International Conference on Machine Learning (ICML-10), pages 13–20, 2010.
[23] Yuhong Guo and Dale Schuurmans. Discriminative batch mode active learning. In Advances in neural information processing systems, pages 593–600, 2008.
[24] H-M Gutmann. A radial basis function method for global optimization. Journal of Global Optimization, 19(3):201–227, 2001.
[25] Shelby J Haberman and Shelby J Haberman. The analysis of frequency data, volume 194. University of Chicago Press Chicago, 1974.
[26] Donghai He, Stephen E Chick, and Chun-Hung Chen. Opportunity cost and OCBA selection procedures in ordinal optimization for a fixed number of alternative systems. Systems, Man, and Cybernetics, Part C: Applications and Reviews, IEEE Transactions on, 37(5):951–961, 2007.
[27] Philipp Hennig and Christian J Schuler. Entropy search for information-efficient global optimization. The Journal of Machine Learning Research, 13(1):1809–1837, 2012.
[28] Matthew D Hoffman, Bobak Shahriari, and Nando de Freitas. On correlation and budget constraints in model-based bandit optimization with application to automatic machine learning. In AISTATS, pages 365–374, 2014.
[29] David W Hosmer Jr and Stanley Lemeshow. Applied logistic regression. John Wiley & Sons, 2004.
[30] Deng Huang, Theodore T Allen, William I Notz, and N Zeng. Global optimization of stochastic black-box systems via sequential kriging meta-models. Journal of global optimization, 34(3):441–466, 2006.
[31] Donald R Jones. A taxonomy of global optimization methods based on response surfaces. Journal of global optimization, 21(4):345–383, 2001.
[32] Donald R Jones, Matthias Schonlau, and William J Welch. Efficient global optimization of expensive black-box functions. Journal of Global optimization, 13(4):455–492, 1998.
[33] Andreas Krause and Cheng S Ong. Contextual gaussian process bandit optimization. In Advances in Neural Information Processing Systems, pages 2447–2455, 2011.
[34] John Langford and Tong Zhang. The epoch-greedy algorithm for contextual multi-armed bandits. In Advances in Neural Information Processing Systems, 2008.
[35] Steffen L Lauritzen. Propagation of probabilities, means, and variances in mixed graphical association models. Journal of the American Statistical Association, 87(420):1098–1108, 1992.
[36] Saskia Le Cessie and Johannes C Van Houwelingen. Ridge estimators in logistic regression. Applied statistics, pages 191–201, 1992.
[37] Lihong Li, Wei Chu, John Langford, and Robert E Schapire. A contextual-bandit approach to personalized news article recommendation. In Proceedings of the 19th international conference on World wide web, pages 661–670. ACM, 2010.
[38] M. Lichman. UCI machine learning repository, 2013.
[39] Dhruv Kumar Mahajan, Rajeev Rastogi, Charu Tiwari, and Adway Mitra. Logucb: an explore-exploit algorithm for comments recommendation. In Proceedings of the 21st ACM international conference on Information and knowledge management, pages 6–15. ACM, 2012.
[40] Peter S Maybeck. Stochastic models, estimation, and control, volume 3. Academic press, 1982.
[41] Martijn RK Mes, Warren B Powell, and Peter I Frazier. Hierarchical knowledge gradient for sequential sampling. The Journal of Machine Learning Research, 12:2931–2974, 2011.
[42] D. C. Montgomery. Design and Analysis of Experiments. John Wiley and Sons, 2008.
[43] Diana M Negoescu, Peter I Frazier, and Warren B Powell. The knowledge-gradient algorithm for sequencing experiments in drug discovery. INFORMS Journal on Computing, 23(3):346–363, 2011.
[44] Warren B Powell and Ilya O Ryzhov. Optimal learning. John Wiley & Sons, 2012.
[45] Pradeep Ravikumar, Martin J Wainwright, John D Lafferty, et al. High-dimensional ising model selection using ℓ1-regularized logistic regression. The Annals of Statistics, 38(3):1287–1319, 2010.
[46] Rommel G Regis and Christine A Shoemaker. Constrained global optimization of expensive black box functions using radial basis functions. Journal of Global Optimization, 31(1):153–171, 2005.
[47] Ilya O Ryzhov, Warren B Powell, and Peter I Frazier. The knowledge gradient algorithm for a general class of online learning problems. Operations Research, 60(1):180–195, 2012.
[48] Mehran Sahami, Susan Dumais, David Heckerman, and Eric Horvitz. A bayesian approach to filtering junk e-mail. In Learning for Text Categorization: Papers from the 1998 workshop, volume 62, pages 98–105, 1998.
[49] Andrew I Schein and Lyle H Ungar. Active learning for logistic regression: an evaluation. Machine Learning, 68(3):235–265, 2007.
[50] Burr Settles. Active learning literature survey. University of Wisconsin, Madison, 52(55-66):11, 2010.
[51] David J Spiegelhalter and Steffen L Lauritzen. Sequential updating of conditional probabilities on directed graphical structures. Networks, 20(5):579–605, 1990.
[52] Niranjan Srinivas, Andreas Krause, Sham M Kakade, and Matthias Seeger. Gaussian process optimization in the bandit setting: No regret and experimental design. arXiv preprint arXiv:0912.3995, 2009.
[53] Matthew Tesch, Jeff Schneider, and Howie Choset. Expensive function optimization with stochastic binary outcomes. In Proceedings of The 30th International Conference on Machine Learning, pages 1283–1291, 2013.
[54] Surya T Tokdar and Jayanta K Ghosh. Posterior consistency of logistic gaussian process priors in density estimation. Journal of Statistical Planning and Inference, 137(1):34–42, 2007.
[55] Simon Tong and Daphne Koller. Support vector machine active learning with applications to text classification. The Journal of Machine Learning Research, 2:45–66, 2002.
[56] Yingfei Wang, Hua Ouyang, Chu Wang, Jianhui Chen, Tsvetan Asamov, and Yi Chang. Efficient ordered combinatorial semi-bandits for whole-page recommendation. In AAAI, pages 2746–2753, 2017.
[57] Yingfei Wang, Kristofer G Reyes, Keith A Brown, Chad A Mirkin, and Warren B Powell. Nested-batch-mode learning and stochastic optimization with an application to sequential multistage testing in materials science. SIAM Journal on Scientific Computing, 37(3):B361–B381, 2015.
[58] Yingfei Wang, Chu Wang, and Warren Powell. The knowledge gradient for sequential decision making with stochastic binary feedbacks. In Proceedings of the 33rd international conference on Machine learning, 2016.
[59] G. B. Wetherill and K. D. Glazebrook. Sequential Methods in Statistics. Chapman and Hall, 1986.
[60] Stephen J Wright and Jorge Nocedal. Numerical optimization, volume 2. Springer New York, 1999.
[61] Tong Zhang and F Oles. The value of unlabeled data for classification problems. In Proceedings of the Seventeenth International Conference on Machine Learning,(Langley, P., ed.), pages 1191–1198. Citeseer, 2000.

8 Proof of Lemma 2.1

First notice that for any fixed point $\bm{x}$ ,

$\displaystyle\mathbb{E}_{N}[\text{Pr}(y=+1\|\bm{x},\bm{w})]$	$\displaystyle=$	$\displaystyle\mathbb{E}_{N}[\sigma(\bm{w}^{T}\bm{x})]$
	$\displaystyle=$	$\displaystyle\int\sigma(\bm{w}^{T}\bm{x})\text{Pr}(\bm{w}\|\mathcal{D}^{N})\text{d}\bm{w}$
	$\displaystyle=$	$\displaystyle\text{Pr}(y=+1\|\bm{x},\mathcal{D}^{N}).$

By the tower property of conditional expectations, and since $\bm{x}^{\pi}$ is $\mathcal{F}^{N}$ measurable,

$\displaystyle\mathbb{E}[\text{Pr}(y=+1\|\bm{x}^{\pi},\bm{w})]$	$\displaystyle=$	$\displaystyle\mathbb{E}[\sigma(\bm{w}^{T}\bm{x}^{\pi})]$
	$\displaystyle=$	$\displaystyle\mathbb{E}\mathbb{E}_{N}[\sigma(\bm{w}^{T}\bm{x}^{\pi})]$
	$\displaystyle=$	$\displaystyle\mathbb{E}[\text{Pr}(y=+1\|\bm{x}^{\pi},\mathcal{D}^{N})].$

Then by the definition of $\bm{x}^{\pi}$ , we have $\text{Pr}(y=+1|\bm{x}^{\pi},\mathcal{D}^{N})=\max_{\bm{x}}\text{Pr}[y=+1|\bm{x},\mathcal{D}^{N}].$

9 Proof of Lemma 4.1.

Let $f(p)=\sum_{i=1}^{d}(m_{i}+yp\frac{x_{i}}{q_{i}})x_{i}$ . Then

	$\displaystyle\frac{\text{d}}{\text{d}{p}}g(p)$	$\displaystyle=$	$\displaystyle-\bigg{(}\frac{\phi(f)^{2}}{y\Phi(yf)^{2}}+\frac{f\phi(f)}{\Phi(yf)}\bigg{)}\cdot y\sum_{i=1}^{d}\frac{x_{i}^{2}}{q_{i}}$
		$\displaystyle=$	$\displaystyle-\frac{\phi(f)}{\Phi(yf)^{2}}\sum_{i=1}^{d}\frac{x_{i}^{2}}{q_{i}}\bigg{(}\phi(f)+yf\Phi(yf)\bigg{)}.$

It is obvious that $\frac{\phi(f)}{\Phi(yf)^{2}}\sum_{i=1}^{d}\frac{x_{i}^{2}}{q_{i}}\geq 0$ . We only need to show that $\phi(f)+yf\Phi(yf)\geq 0$ . Since the normal density function $\mathcal{N}(\cdot|0,1)$ is symmetric with respect to the origin and $y\in\{+1,-1\}$ , we have $\phi(f)+yf\Phi(yf)=\phi(yf)+yf\Phi(yf)$ . Define $u(t)=\phi(t)+t\Phi(t)$ . Since $\lim_{t\to-\infty}u(t)=0$ , we show $u(t)\geq 0$ for $t\in[-\infty,\infty]$ by proving that $u^{\prime}(t)\geq 0$ .

First note that $\frac{\text{d}}{\text{d}t}e^{-t^{2}/2}=-te^{-t^{2}/2}$ . So we have $\phi^{\prime}(t)=-t\phi(t)$ . Thus $u^{\prime}(t)=-t\phi(t)+\Phi(t)+t\phi(t)=\Phi(t)\geq 0$ .

10 Proof of Theorem 5.4

For notation simplicity, we use $\bm{\Sigma}$ and $\bm{m}$ instead of $\bm{\Sigma}^{0}$ and $\bm{m}^{0}$ , and use $\bm{x}_{i}$ , $y_{i}$ to denote $\bm{x}^{i}$ and $y^{i}$ in the proof. Consider the objective function $f(\bm{w})$ defined as

f(\bm{w})=\frac{1}{2}(\bm{w}-\bm{m})^{T}\bm{\Sigma}^{-1}(\bm{w}-\bm{m})-\mathbb{E}_{\bm{y}\sim\mathcal{B}(\mathcal{D}^{n},\bm{w}^{*})}\Big{[}\sum_{i=1}^{n}\log\sigma(y_{i}\bm{w}^{T}\bm{x}_{i})\Big{]}.

We use $g(\bm{w})$ and $h(\bm{w})$ to denote the quadratic term and the rest, respectively. By applying the mean value theorem to $h(\bm{w})$ , we have

h(\bm{w})-h(\bm{w}^{*})=(\bm{w}-\bm{w}^{*})^{T}\nabla h(\bm{w}^{*})+\frac{1}{2}(\bm{w}-\bm{w}^{*})^{T}H\big{(}\bm{w}^{*}+\eta(\bm{w}-\bm{w}^{*})\big{)}(\bm{w}-\bm{w}^{*}),

where $H$ is the Hessian of $h(\bm{w})$ . To analyze the first and second order terms in the expansion, we use a similar technique adopted by [45]. For the gradient in the first order term, we have

			$\displaystyle\nabla h(\bm{w}^{})=\mathbb{E}_{\bm{y}\sim\mathcal{B}(\mathcal{D}^{n},\bm{w}^{})}\Big{[}\sum_{i=1}^{n}\big{(}1-\sigma(y_{i}\bm{w}^{T}\bm{x}_{i})\big{)}y_{i}\bm{x}_{i}\Big{]}\bigg{\|}_{\bm{w}=\bm{w}^{*}}$
		$\displaystyle=$	$\displaystyle\sum_{i=1}^{n}\Big{(}\big{(}1-\sigma(\bm{w}^{T}\bm{x}_{i})\big{)}\sigma((\bm{w}^{})^{T}\bm{x}_{i})-\big{(}1-\sigma(-\bm{w}^{T}\bm{x}_{i})\big{)}\big{(}1-\sigma((\bm{w}^{})^{T}\bm{x}_{i})\big{)}\Big{)}\bm{x}_{i}\bigg{\|}_{\bm{w}=\bm{w}^{*}}=0.$

For the second order term, we have

			$\displaystyle H\big{(}\bm{w}^{}+\eta(\bm{w}-\bm{w}^{}))$
		$\displaystyle=$	$\displaystyle\mathbb{E}_{\bm{y}\sim\mathcal{B}(\mathcal{D}^{n},\bm{w}^{})}\Bigg{(}\sum_{i=1}^{n}\sigma(y_{i}(\bm{w}^{}+\eta(\bm{w}-\bm{w}^{}))^{T}\bm{x}_{i})\Big{(}1-\sigma(y_{i}(\bm{w}^{}+\eta(\bm{w}-\bm{w}^{*}))^{T}\bm{x}_{i})\Big{)}\bm{x}_{i}\bm{x}_{i}^{T}\Bigg{)}$
		$\displaystyle=$	$\displaystyle\sum_{i=1}^{n}\sigma((\bm{w}^{}+\eta(\bm{w}-\bm{w}^{}))^{T}\bm{x}_{i})\Big{(}1-\sigma((\bm{w}^{}+\eta(\bm{w}-\bm{w}^{}))^{T}\bm{x}_{i})\Big{)}\bm{x}_{i}\bm{x}_{i}^{T}$
		$\displaystyle=$	$\displaystyle\sum_{i=1}^{n}J_{i}(\bm{w},\eta)\bm{x}_{i}\bm{x}_{i}^{T},$

where $J_{i}(\bm{w},\eta)=\sigma((\bm{w}^{*}+\eta(\bm{w}-\bm{w}^{*}))^{T}\bm{x}_{i})\Big{(}1-\sigma((\bm{w}^{*}+\eta(\bm{w}-\bm{w}^{*}))^{T}\bm{x}_{i})\Big{)}$ .

We expand $J_{i}(\bm{w},\eta)$ to its first order and use mean value theorem again,

$\displaystyle\|J_{i}(\bm{w},\eta)-J_{i}(\bm{w}^{*},\eta)\|$	$\displaystyle=$	$\displaystyle\|\eta(\bm{w}-\bm{w}^{*})^{T}\bm{x}_{i}J_{i}^{\prime}\|$
	$\displaystyle\leq$	$\displaystyle\|\sigma(1-\sigma)(1-2\sigma)\|\|(\bm{w}-\bm{w}^{*})^{T}\bm{x}_{i}\|$
	$\displaystyle\leq$	$\displaystyle\\|\bm{w}-\bm{w}^{*}\\|_{2},$

where we omit the dependence of $\sigma$ on $(\bm{w}^{*}+\eta(\bm{w}-\bm{w}^{*}))^{T}\bm{x}_{i}$ for simplicity and use the fact $\sigma\in(0,1)$ .

Combining the first order and second order analysis and denoting $\|\bm{w}-\bm{w}^{*}\|_{2}$ as $R$ , we have

$\displaystyle h(\bm{w})-h(\bm{w}^{*})$	$\displaystyle\geq$	$\displaystyle\frac{1}{2}\\|\bm{w}-\bm{w}^{}\\|_{2}^{2}\lambda_{min}\Big{(}H\big{(}\bm{w}^{}+\eta(\bm{w}-\bm{w}^{*})\big{)}\Big{)}$
	$\displaystyle\geq$	$\displaystyle\frac{1}{2}R^{2}\lambda_{min}\Big{(}\sum_{i=1}^{n}\left[\sigma(\bm{x}_{i}^{T}\bm{w}^{})\big{(}1-\sigma(\bm{x}_{i}^{T}\bm{w}^{})\big{)}-R\right]\bm{x}i\bm{x}_{i}^{T}\Big{)}$
	$\displaystyle\geq$	$\displaystyle\frac{1}{2}R^{2}\lambda_{min}\Big{(}\sum_{i=1}^{n}\sigma(\bm{x}_{i}^{T}\bm{w}^{})\big{(}1-\sigma(\bm{x}_{i}^{T}\bm{w}^{})\big{)}\bm{x}i\bm{x}_{i}^{T}\Big{)}-\frac{1}{2}R^{3}$

where we use the fact that $\bm{x}_{i}\bm{x}_{i}^{T}$ is positive semi-definite and $\|\bm{x}_{i}\|_{2}\leq 1$ . By using $C_{min}$ to denote the minimal eigenvalue in the last inequality, we have

h(\bm{w})-h(\bm{w}^{*})\geq\frac{1}{2}C_{min}R^{2}-\frac{1}{2}R^{3}.

On the other hand, for the quadratic term $g(\bm{w})$ , we have

	$\displaystyle g(\bm{w})-g(\bm{w}^{*})$	$\displaystyle=$	$\displaystyle(\bm{w}-\bm{w}^{})^{T}\bm{\Sigma}^{-1}(\bm{w}^{}-\bm{m})+\frac{1}{2}(\bm{w}-\bm{w}^{})^{T}\bm{\Sigma}^{-1}(\bm{w}-\bm{w}^{})$
		$\displaystyle\geq$	$\displaystyle-DR\sqrt{\lambda_{max}\big{(}\bm{\Sigma}^{-1}\big{)}}+\frac{1}{2}\lambda_{min}\big{(}\bm{\Sigma}^{-1}\big{)}R^{2},$

where $D=\|\bm{\Sigma}^{-\frac{1}{2}}(\bm{w}^{*}-\bm{m})\|_{2}$ .

Now define function $F(\Delta)$ as $F(\Delta)=f(\bm{w}^{*}+\Delta)-f(\bm{w}^{*})$ , then we have

	$\displaystyle F(\Delta)$	$\displaystyle=$	$\displaystyle g(\bm{w})-g(\bm{w}^{})+h(\bm{w}^{n})-h(\bm{w}^{})$		(23)
		$\displaystyle\geq$	$\displaystyle-\sqrt{\lambda_{max}\big{(}\bm{\Sigma}^{-1}\big{)}}DR+\frac{1}{2}\lambda_{min}\big{(}\bm{\Sigma}^{-1}\big{)}R^{2}+\frac{1}{2}C_{min}R^{2}-\frac{1}{2}R^{3}.$		(23)

From now on we will use the simplified symbols $\lambda_{max}$ and $\lambda_{min}$ instead of $\lambda_{min}\big{(}\bm{\Sigma}^{-1}\big{)}$ and $\lambda_{max}\big{(}\bm{\Sigma}^{-1}\big{)}$ . We consider the sign of $F(\Delta)$ in the case when $\|\Delta\|_{2}=\frac{1}{2}(C_{min}+\lambda_{min})$ and

D\leq\frac{1}{16}\frac{\lambda_{min}^{2}}{\sqrt{\lambda_{max}}}.

(24)

Based on the analysis in (23), we have

$\displaystyle F(\Delta)$	$\displaystyle\geq$	$\displaystyle-\frac{1}{16}\lambda_{min}^{2}R+\frac{1}{2}\lambda_{min}R^{2}+\frac{1}{2}C_{min}R^{2}-\frac{1}{2}R^{3}$
	$\displaystyle=$	$\displaystyle\frac{R}{2}\left(-R^{2}+(\lambda_{min}+C_{min})R-\frac{1}{8}\lambda_{min}^{2}\right)$
	$\displaystyle=$	$\displaystyle\frac{R}{2}\left(R^{2}-\frac{1}{8}\lambda_{min}^{2}\right)>\frac{\lambda_{min}^{3}}{32}.$

Notice that $F(0)=0$ and recall that $\bm{w}^{n}$ minimizes $f(\bm{w})$ so we have

F(\bm{w}^{n}-\bm{w}^{*})=f(\bm{w}^{n})-f(\bm{w}^{*})\leq 0.

Then based on the convexity of $F$ we know that $\|\bm{w}^{n}-\bm{w}^{*}\|_{2}\leq\frac{1}{2}(C_{min}+\lambda_{min})$ , otherwise the values of $F$ at 0, $\bm{w}^{n}-\bm{w}^{*}$ and the intersection between the all $\|\Delta\|_{2}=\frac{1}{2}(C_{min}+\lambda_{min})$ and line from 0 to $\bm{w}^{n}-\bm{w}^{*}$ form a concave pattern, which is contradictory. Similarly, based on the assumption that $f(\hat{\bm{w}}^{n})\leq\min_{\bm{w}}f(\bm{w})+\lambda_{min}^{3}/32$ , we have

F(\hat{\bm{w}}^{n}-\bm{w}^{*})=f(\hat{\bm{w}}^{n})-f(\bm{w}^{*})\leq\frac{\lambda_{min}^{3}}{32}.

Recall that $F(0)=0$ and $F(\Delta)>\lambda_{min}^{3}/32$ , then the convexity of $F$ also ensures that $\|\hat{\bm{w}}^{n}-\bm{w}^{*}\|_{2}\leq\frac{1}{2}(C_{min}+\lambda_{min})$ .

Now we start to calculate the probability for (24) to hold. Recall that $\bm{w}^{*}$ has a prior distribution $\bm{w}^{*}\sim\mathcal{N}(\bm{m},\bm{\Sigma})$ . Then by denoting the right hand side of Eq. (24) as $M$ , we have

			$\displaystyle\text{Pr}\Big{(}D\leq\frac{1}{16}\frac{\lambda_{min}^{2}}{\sqrt{\lambda_{max}}}\Big{)}$
		$\displaystyle=$	$\displaystyle\int_{\\|\bm{\Sigma}^{-\frac{1}{2}}(\bm{a}-\bm{m})\\|\leq M}(2\pi)^{-\frac{d}{2}}\|\bm{\Sigma}\|^{-\frac{1}{2}}\exp\big{(}-\frac{1}{2}(\bm{a}-\bm{m})^{T}\bm{\Sigma}^{-1}(\bm{a}-\bm{m})\big{)}\mathrm{d}\bm{a}$
		$\displaystyle=$	$\displaystyle\int_{\\|\bm{b}\\|_{2}\leq M}(2\pi)^{-\frac{d}{2}}\exp\big{(}-\frac{1}{2}\bm{b}^{T}\bm{b}\big{)}\mathrm{d}\bm{b},$

which is the probability of a d-dimension standard normal random variable appears in the ball with radius $M$ , $P_{d}(M)$ . This completes the proof.

11 Proofs of Asymptotic Optimality

In this section, we provide detailed proofs of all the asymptotic optimal results in Section 5.3.

11.1 Proof of Proposition 5.5

We use $q(\bm{w})$ to denote the predictive distribution under the state $s$ and use $p(\bm{w}|s,\bm{x},y_{\bm{x}})$ to denote the posterior distribution after we observe the outcome of $\bm{x}$ to be $y$ . By Jensen’s inequality, we have

	$\displaystyle\mu_{\bm{x}}^{\text{KG}}(s)$	$\displaystyle=$	$\displaystyle\mathbb{E}\Big{[}\max_{\bm{x}^{\prime}}p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},y_{\bm{x}}))\Big{]}-\max_{\bm{x}^{\prime}}p(y=+1\|\bm{x}^{\prime},s)$
		$\displaystyle\geq$	$\displaystyle\max_{\bm{x}^{\prime}}\mathbb{E}\big{[}p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},y_{\bm{x}}))\big{]}-\max_{\bm{x}^{\prime}}p(y=+1\|\bm{x}^{\prime},s).$

We then show that $\mathbb{E}\big{[}p(y=+1|\bm{x}^{\prime},T(s,\bm{x},y_{\bm{x}}))\big{]}=p(y=+1|\bm{x}^{\prime},s)$ for any $\bm{x},\bm{x}^{\prime}$ and $s$ , which leads to $\mu_{\bm{x}}^{\text{KG}}(s)\geq 0$ . Since $y_{\bm{x}}$ is binomial distributed with mean $p(y_{\bm{x}}=+1|\bm{x},s)$ , we have

			$\displaystyle\mathbb{E}\big{[}p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},y_{\bm{x}}))\big{]}$
		$\displaystyle=$	$\displaystyle p(y_{\bm{x}}=+1\|\bm{x},s)p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},+1))$
		$\displaystyle+$	$\displaystyle(1-p(y_{\bm{x}}=+1\|\bm{x},s))p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},-1)).$

Recall that $p(y=+1|\bm{x},s)=\int\sigma(\bm{w}^{T}\bm{x})p(\bm{w}|s)\text{d}\bm{w}.$ By Bayes’ Theorem, the posterior distribution in the updated state $T(s,\bm{x},y_{x})$ becomes

p(\bm{w}^{\prime}|T(s,\bm{x},+1))=\frac{\sigma((\bm{w}^{\prime})^{T}\bm{x})p(\bm{w}^{\prime}|s)}{\int\sigma(\bm{w}^{T}\bm{x})p(\bm{w}|s)\text{d}\bm{w}},

and

p(\bm{w}^{\prime}|T(s,\bm{x},-1))=\frac{\left(1-\sigma((\bm{w}^{\prime})^{T}\bm{x})\right)p(\bm{w}^{\prime}|s)}{\int\left(1-\sigma(\bm{w}^{T}\bm{x})\right)p(\bm{w}|s)\text{d}\bm{w}}.

Notice that

$\displaystyle p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},+1))$	$\displaystyle=$	$\displaystyle\int\sigma((\bm{w}^{\prime})^{T}\bm{x}^{\prime})p(\bm{w}^{\prime}\|T(s,\bm{x},+1))\text{d}\bm{w}^{\prime}$
	$\displaystyle=$	$\displaystyle\int\sigma((\bm{w}^{\prime})^{T}\bm{x}^{\prime})\frac{\sigma((\bm{w}^{\prime})^{T}\bm{x})p(\bm{w}^{\prime}\|s)}{\int\sigma(\bm{w}^{T}\bm{x})p(\bm{w}\|s)\text{d}\bm{w}}\text{d}\bm{w}^{\prime}$
	$\displaystyle=$	$\displaystyle\frac{\int\sigma((\bm{w}^{\prime})^{T}\bm{x}^{\prime})\sigma((\bm{w}^{\prime})^{T}\bm{x})p(\bm{w}^{\prime}\|s)\text{d}\bm{w}^{\prime}}{\int\sigma(\bm{w}^{T}\bm{x})p(\bm{w}\|s)\text{d}\bm{w}},$

and similarly, we have

p(y=+1|\bm{x}^{\prime},T(s,\bm{x},-1))=\frac{\int\sigma((\bm{w}^{\prime})^{T}\bm{x}^{\prime})(1-\sigma((\bm{w}^{\prime})^{T}\bm{x}))p(\bm{w}^{\prime}|s)\text{d}\bm{w}^{\prime}}{\int(1-\sigma(\bm{w}^{T}\bm{x}))p(\bm{w}|s)\text{d}\bm{w}}.

Therefore,

			$\displaystyle\mathbb{E}\big{[}p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},y_{\bm{x}}))\big{]}$
		$\displaystyle=$	$\displaystyle p(y=+1\|s,\bm{x})p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},+1))+p(y=-1\|s,\bm{x})p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},-1))$
		$\displaystyle=$	$\displaystyle\int\sigma((\bm{w}^{\prime})^{T}\bm{x}^{\prime})\sigma((\bm{w}^{\prime})^{T}\bm{x})p(\bm{w}^{\prime}\|s)\text{d}\bm{w}^{\prime}+\int\sigma((\bm{w}^{\prime})^{T}\bm{x}^{\prime})(1-\sigma((\bm{w}^{\prime})^{T}\bm{x}))p(\bm{w}^{\prime}\|s)\text{d}\bm{w}^{\prime}$
		$\displaystyle=$	$\displaystyle\int\sigma((\bm{w}^{\prime})^{T}\bm{x}^{\prime})p(\bm{w}^{\prime}\|s)\text{d}\bm{w}^{\prime}$
		$\displaystyle=$	$\displaystyle p(y=+1\|s,\bm{x}^{\prime}),$

and thus we obtain

\mathbb{E}\big{[}p(y=+1|\bm{x}^{\prime},T(s,\bm{x},y_{\bm{x}}))\big{]}=p(y=+1|s,\bm{x}^{\prime}).

11.2 Proof of Proposition 5.7

The proof is similar to that by [15] with additional tricks for Bernoulli distributed random variables. Let $\mathcal{G}$ be the sigma-algebra by the collection $\{\hat{y}^{n+1}\bm{1}_{\{\bm{x}^{n}=\bm{x}\}}\}$ . Since if the policy $\pi$ measures alternative $\bm{x}$ infinitely often, this collection is an infinite sequence of independent random variables with common Bernoulli distribution with mean $\sigma(\bm{w}^{T}\bm{x})$ , the strong law of large numbers implies $\sigma(\bm{w}^{T}\bm{x})\in\mathcal{G}$ . Since $\mathcal{G}\in\mathcal{F}^{\infty}$ , we have $\sigma(\bm{w}^{T}\bm{x})\in\mathcal{F}^{\infty}$ . Let $U$ be a uniform random variable in $[0,1]$ . Then the Bernoulli random variable $y_{\bm{x}}$ can be rewritten as $\bm{1}_{U\leq\sigma(\bm{w}^{T}\bm{x})}$ . Since $U$ is independent with $\mathcal{F}^{\infty}$ and the $\sigma$ -algebra generated by $\sigma(\bm{w}^{T}\bm{x})$ , it can be shown that by properties of conditional expectations,

\mathbb{E}\big{[}\sigma(\bm{w}^{T}\bm{x}^{\prime})|\mathcal{F}^{\infty},\bm{1}_{U\leq\sigma(\bm{w}^{T}\bm{x})}\big{]}=\mathbb{E}\big{[}\sigma(\bm{w}^{T}\bm{x}^{\prime})|\mathcal{F}^{\infty}\big{]}.

We next show that the knowledge gradient value of measuring alternative $\bm{x}$ is zero by substituting this relation into the definition of the knowledge gradient. We have

$\displaystyle\nu_{\bm{x}}(\mathcal{F}^{\infty})$	$\displaystyle=$	$\displaystyle\mathbb{E}\bigg{[}\max_{\bm{x}^{\prime}}\mathbb{E}\big{[}\sigma(\bm{w}^{T}\bm{x}^{\prime})\|\mathcal{F}^{\infty},\bm{1}_{U\leq\sigma(\bm{w}^{T}\bm{x})}\big{]}\|\mathcal{F}^{\infty}\bigg{]}-\max_{\bm{x}^{\prime}}\mathbb{E}\big{[}\sigma(\bm{w}^{T}\bm{x}^{\prime})\|\mathcal{F}^{\infty}\big{]}$
	$\displaystyle=$	$\displaystyle\mathbb{E}\bigg{[}\max_{\bm{x}^{\prime}}\mathbb{E}\big{[}\sigma(\bm{w}^{T}\bm{x}^{\prime})\|\mathcal{F}^{\infty}\big{]}\|\mathcal{F}^{\infty}\bigg{]}-\max_{\bm{x}^{\prime}}\mathbb{E}\big{[}\sigma(\bm{w}^{T}\bm{x}^{\prime})\|\mathcal{F}^{\infty}\big{]}$
	$\displaystyle=$	$\displaystyle\max_{\bm{x}^{\prime}}\mathbb{E}\big{[}\sigma(\bm{w}^{T}\bm{x}^{\prime})\|\mathcal{F}^{\infty}\big{]}-\max_{\bm{x}^{\prime}}\mathbb{E}\big{[}\sigma(\bm{w}^{T}\bm{x}^{\prime})\|\mathcal{F}^{\infty}\big{]}$
	$\displaystyle=$	$\displaystyle 0.$

11.3 Proof of the Theorem 5.9: Consistency of the KG Policy

It has been established that almost surely the knowledge gradient policy will achieve a state (at time T) that the KG values for all the alternatives are smaller than $\epsilon$ , after which the probability of selecting each alternative is $1/M$ where $M$ is the number of alternatives. Notice that in each round, KG policy first picks one out of $M$ alternatives, then a feedback of either 1 or -1 is observed. Equivalently, we can interpret the two procedures as one of the $2M$ possible outcomes from the set $\tilde{\mathcal{Y}}:=\{(0,\dots,+1/-1,\dots 0)\}$ , where each $\tilde{y}\in\tilde{\mathcal{Y}}$ is a $M$ -dimensional vector with only element being +1 or -1. It should be noted that by changing the feedback schema in this way will not affect the Bayesian update equations because the likelihood function and the normalization factor in the posterior will both multiply by a factor of $1/M$ . On the other hand, this combined feedback schema makes it possible to treat each measurement $(\bm{x}^{n},y^{n})$ as i.i.d. samples in $\tilde{\mathcal{Y}}$ .

Define the K-L neiborhood as $K_{\epsilon}(u)=\{v:KL(u,v)<\epsilon\}$ , where the K-L divergence is defined as $KL(u,v):=\int v\log(v/u)$ . Since the prior distribution is Gaussian with positive definite covariance matrix, and the likelihood function is the sigmoid function which only takes positive values, then after time $T$ , the posterior probability in the K-L neighborhood of $w^{*}$ is positive. Based on standard results on the consistency of Bayes’ estimates [20, 19, 54, 17], the posterior is weakly consistent at $w^{*}$ in the sense that for any neighborhood $U$ of $w^{*}$ , the probability that $\mu(w)$ lies in $U$ converges to 1.

\mathbb{P}[U|\tilde{y}^{1},\tilde{y}^{2},\dots,\tilde{y}^{n}]\rightarrow 1.

Without loss of generality, assume that the alternative $\bm{x}^{*}$ with the largest probability of +1 is unique, which means $\sigma((\bm{x}^{*})^{T}\bm{w}^{*})>\sigma(\bm{x}^{T}\bm{w}^{*})$ for any alternative $\bm{x}$ other than $\bm{x}^{*}$ . Then we can pick $U$ to be the neighborhood of $\bm{w}^{*}$ such that $\sigma((\bm{x}^{*})^{T}\bm{w})>\sigma(\bm{x}^{T}\bm{w})$ holds for any $\bm{w}\in U$ . The neighborhood $U$ exists because we only have finite number of altervatives. From the consistency results, the probability that the best arm under posterior estimation is the true best alternative goes to 1 as the measurement budget goes to infinity.

$\displaystyle\|J_{i}(\bm{w},\eta)-J_{i}(\bm{w}^{*},\eta)\|$	$\displaystyle=$	$\displaystyle\|\eta(\bm{w}-\bm{w}^{*})^{T}\bm{x}_{i}J_{i}^{\prime}\|$
	$\displaystyle\leq$	$\displaystyle\|\sigma(1-\sigma)(1-2\sigma)\|\|(\bm{w}-\bm{w}^{*})^{T}\bm{x}_{i}\|$
	$\displaystyle\leq$	$\displaystyle\\|\bm{w}-\bm{w}^{*}\\|_{2},$

			$\displaystyle\mathbb{E}\big{[}p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},y_{\bm{x}}))\big{]}$
		$\displaystyle=$	$\displaystyle p(y_{\bm{x}}=+1\|\bm{x},s)p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},+1))$
		$\displaystyle+$	$\displaystyle(1-p(y_{\bm{x}}=+1\|\bm{x},s))p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},-1)).$

$\displaystyle p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},+1))$	$\displaystyle=$	$\displaystyle\int\sigma((\bm{w}^{\prime})^{T}\bm{x}^{\prime})p(\bm{w}^{\prime}\|T(s,\bm{x},+1))\text{d}\bm{w}^{\prime}$
	$\displaystyle=$	$\displaystyle\int\sigma((\bm{w}^{\prime})^{T}\bm{x}^{\prime})\frac{\sigma((\bm{w}^{\prime})^{T}\bm{x})p(\bm{w}^{\prime}\|s)}{\int\sigma(\bm{w}^{T}\bm{x})p(\bm{w}\|s)\text{d}\bm{w}}\text{d}\bm{w}^{\prime}$
	$\displaystyle=$	$\displaystyle\frac{\int\sigma((\bm{w}^{\prime})^{T}\bm{x}^{\prime})\sigma((\bm{w}^{\prime})^{T}\bm{x})p(\bm{w}^{\prime}\|s)\text{d}\bm{w}^{\prime}}{\int\sigma(\bm{w}^{T}\bm{x})p(\bm{w}\|s)\text{d}\bm{w}},$

			$\displaystyle\mathbb{E}\big{[}p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},y_{\bm{x}}))\big{]}$
		$\displaystyle=$	$\displaystyle p(y=+1\|s,\bm{x})p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},+1))+p(y=-1\|s,\bm{x})p(y=+1\|\bm{x}^{\prime},T(s,\bm{x},-1))$
		$\displaystyle=$	$\displaystyle\int\sigma((\bm{w}^{\prime})^{T}\bm{x}^{\prime})\sigma((\bm{w}^{\prime})^{T}\bm{x})p(\bm{w}^{\prime}\|s)\text{d}\bm{w}^{\prime}+\int\sigma((\bm{w}^{\prime})^{T}\bm{x}^{\prime})(1-\sigma((\bm{w}^{\prime})^{T}\bm{x}))p(\bm{w}^{\prime}\|s)\text{d}\bm{w}^{\prime}$
		$\displaystyle=$	$\displaystyle\int\sigma((\bm{w}^{\prime})^{T}\bm{x}^{\prime})p(\bm{w}^{\prime}\|s)\text{d}\bm{w}^{\prime}$
		$\displaystyle=$	$\displaystyle p(y=+1\|s,\bm{x}^{\prime}),$