Best Arm Identification in Generalized Linear Bandits

TL;DR

This paper introduces the first algorithm for best-arm identification in generalized linear bandits, providing theoretical guarantees and demonstrating its effectiveness through simulations, with applications like drug design.

Contribution

It develops a novel algorithm for best-arm identification in generalized linear bandits, extending previous linear bandit methods with theoretical analysis.

Findings

Algorithm achieves near-optimal sample complexity

Theoretical guarantees on accuracy and efficiency

Effective performance demonstrated in simulations

Abstract

Motivated by drug design, we consider the best-arm identification problem in generalized linear bandits. More specifically, we assume each arm has a vector of covariates, there is an unknown vector of parameters that is common across the arms, and a generalized linear model captures the dependence of rewards on the covariate and parameter vectors. The problem is to minimize the number of arm pulls required to identify an arm that is sufficiently close to optimal with a sufficiently high probability. Building on recent progress in best-arm identification for linear bandits (Xu et al. 2018), we propose the first algorithm for best-arm identification for generalized linear bandits, provide theoretical guarantees on its accuracy and sampling efficiency, and evaluate its performance in various scenarios via simulation.

Equations

$${\mu }_a=\mu ({\theta }^{⊺}{x}^a),$$

$$\mathbb{E}[r_t|a_t]={\mu }_{a_t}=\mu ({\theta }^{⊺}{x}^{a_t}).$$

$$\mathbb{P}\left[\mu ({\theta }^{⊺}{x}^{a^{\ast }})-\mu ({\theta }^{⊺}{x}^{{\hat{a}}_{\tau }})\ge ϵ\right]\le \delta .$$

$${\Delta }_i=\{\begin{array}{rlrl}& \mu ({\theta }^{\top }{x}^{a^{\ast }})-\mu ({\theta }^{\top }{x}^i)& & \text{mbox}ifi\ne a^{\ast }\\ & \mu ({\theta }^{\top }{x}^i)-\underset{j\ne i}{\max }\mu ({\theta }^{\top }{x}^j)& & \text{mbox}ifi=a^{\ast }.\end{array}$$

$$\sum _{l=1}^{t-1}(r_l-\mu ({\theta }_t^{⊺}{x}_l))x_l=0.$$

$${\Delta }_t(i,j)=\mu ({\theta }_t^{\top }{x}^i)-\mu ({\theta }_t^{\top }{x}^j)$$

$$C_{t}=\alpha\sqrt{2d\log t\log\Bigl{(}\frac{\pi^{2}dt^{2}}{6\delta}\Bigr{)}},$$

$$\alpha =\frac{2\kappa R}{c_{\mu }}$$

$${\mathcal{C}}_t(i,j)=\left\{\Delta \in \mathbb{R}:|\Delta -{\Delta }_t(i,j)|\le C_t\underset{c,c^{\prime }\in [c_{\mu },k_{\mu }]}{\max }\Vert c{x}^i-c^{\prime }{x}^j{\Vert }_{M_t^{-1}}\right\}.$$

$${\beta }_t(i,j)=C_t\underset{c,c^{\prime }\in [c_{\mu },k_{\mu }]}{\max }\Vert c{x}^i-c^{\prime }{x}^j{\Vert }_{M_t^{-1}}$$

$$p_a(y)=\frac{|w_a^{\ast }(y)|}{{\sum }_{k=1}^K|w_k^{\ast }(y)|},$$

$$\underset{w\in {\mathbb{R}}^K}{\min }\Vert w{\Vert }_1\text{mbox}s.t.\sum _{k=1}^K{w}_k{x}^k=y.$$

$$a_{t+1}=\underset{a\in [K]:p_a(y_t)>0}{\mathrm{argmin}}\frac{T_a(t)}{p_a(y_t)},$$

$$\forall x,x^{\prime }\in \mathcal{X}:|\left[\mu ({\theta }_t^{⊺}x)-\mu ({\theta }_t^{⊺}{x}^{\prime })\right]-\left[\mu ({\theta }^{⊺}x)-\mu ({\theta }^{⊺}{x}^{\prime })\right]|\le D_t(\delta )\underset{c,c^{\prime }\in [c_{\mu },k_{\mu }]}{\max }\Vert cx-c^{\prime }{x}^{\prime }{\Vert }_{M_t^{-1}},$$

$$D_t(\delta )=\frac{2\kappa R}{c_{\mu }}\sqrt{2d\log t\log (d/\delta )}.$$

$$|\left[\mu ({\theta }_t^{⊺}x)-\mu ({\theta }_t^{⊺}{x}^{\prime })\right]-\left[\mu ({\theta }^{⊺}x)-\mu ({\theta }^{⊺}{x}^{\prime })\right]|$$

$$=\left|\left[\dot{\mu }(z)({\theta }_t-\theta {)}^{\top }x-\dot{\mu }(z^{\prime })({\theta }_t-\theta {)}^{\top }{x}^{\prime }\right]\right|,$$

$$=\left|({\theta }_t-\theta {)}^{\top }(cx-c^{\prime }{x}^{\prime })\right|.$$

$$g_t({\theta }_t)-g_t(\theta )=G_t({\theta }_t-\theta ),$$

$$G_t={\int }_0^1\nabla g_t(s\theta +(1-s){\theta }_t)ds.$$

$$\nabla g_t(\beta )=\sum _{l=1}^{t-1}{x}_l{x}_l^{\top }\dot{\mu }({\beta }^{\top }{x}_l).$$

$$\left|\left[\mu ({\theta }_t^{⊺}x)-\mu ({\theta }_t^{⊺}{x}^{\prime })\right]-\left[\mu ({\theta }^{⊺}x)-\mu ({\theta }^{⊺}{x}^{\prime })\right]\right|$$

$$\le \Vert cx-c^{\prime }{x}^{\prime }{\Vert }_{G_t^{-1}}\Vert g_t({\theta }_t)-g_t(\theta ){\Vert }_{G_t^{-1}}.$$

$$\left|\left[\mu ({\theta }_t^{⊺}x)-\mu ({\theta }_t^{⊺}{x}^{\prime })\right]-\left[\mu ({\theta }^{⊺}x)-\mu ({\theta }^{⊺}{x}^{\prime })\right]\right|\le \frac{2}{c_{\mu }}\Vert cx-c^{\prime }{x}^{\prime }{\Vert }_{M_t^{-1}}\Vert g_t({\theta }_t)-g_t(\theta ){\Vert }_{M_t^{-1}}.$$

$$\Vert g_t({\theta }_t)-g_t(\theta ){\Vert }_{M_t^{-1}}\le \kappa R\sqrt{2d\log t\log (d/\delta )}$$

$$\left|\left[\mu ({\theta }_t^{⊺}x)-\mu ({\theta }_t^{⊺}{x}^{\prime })\right]-\left[\mu ({\theta }^{⊺}x)-\mu ({\theta }^{⊺}{x}^{\prime })\right]\right|\le D_t(\delta )\Vert cx-c^{\prime }{x}^{\prime }{\Vert }_{M_t^{-1}}$$

$$\mathcal{E}=\left\{\forall t>\max (2,d),\forall i,j\in [K]:\Delta (i,j)\in {\mathcal{C}}_t(i,j)\right\},$$

$${\mathcal{E}}_t=\left\{\forall x,x^{\prime }\in \mathcal{X}:|\left[\mu ({\theta }_t^{⊺}x)-\mu ({\theta }_t^{⊺}{x}^{\prime })\right]-\left[\mu ({\theta }^{⊺}x)-\mu ({\theta }^{⊺}{x}^{\prime })\right]|\le D_t({\delta }_t)\underset{c,c^{\prime }\in [c_{\mu },k_{\mu }]}{\max }\Vert cx-c^{\prime }{x}^{\prime }{\Vert }_{M_t^{-1}}\right\}.$$

$$\mathbb{P}[{\cap }_{t=1}^{\infty }{\mathcal{E}}_t]\ge 1-\sum _{t=1}^{\infty }{\delta }_t=1-\frac{6\delta }{{\pi }^2}\sum _{t=1}^{\infty }\frac{1}{t^2}=1-\delta .$$

$$\mathbb{P}[\Delta (a^{\ast },{\hat{a}}_{\tau })>ϵ]<\delta .$$

$$\Delta (a^{\ast },{\hat{a}}_{\tau })>ϵ\ge B(\tau )\ge {\Delta }_{\tau }(a^{\ast },{\hat{a}}_{\tau })+{\beta }_{\tau }(a^{\ast },{\hat{a}}_{\tau }).$$

$$\Delta (a^{\ast },{\hat{a}}_{\tau })-{\Delta }_{\tau }(a^{\ast },{\hat{a}}_{\tau })>{\beta }_{\tau }(a^{\ast },{\hat{a}}_{\tau }),$$

Full text

\DoubleSpacedXI\TheoremsNumberedThrough\ECRepeatTheorems\EquationsNumberedThrough

\RUNAUTHOR

Kazerouni and Wein

\RUNTITLE

Best Arm Identification in Generalized Linear Bandits

\TITLE

Best Arm Identification in Generalized Linear Bandits

\ARTICLEAUTHORS\AUTHOR

Abbas Kazerouni \AFFDepartment of Electrical Engineering, Stanford University, Stanford, CA 94305, \EMAILabbask@stanford.edu \AUTHORLawrence M. Wein \AFFGraduate School of Business, Stanford University, Stanford, CA 94305, \EMAILlwein@stanford.edu

\ABSTRACT

Motivated by drug design, we consider the best-arm identification problem in generalized linear bandits. More specifically, we assume each arm has a vector of covariates, there is an unknown vector of parameters that is common across the arms, and a generalized linear model captures the dependence of rewards on the covariate and parameter vectors. The problem is to minimize the number of arm pulls required to identify an arm that is sufficiently close to optimal with a sufficiently high probability. Building on recent progress in best-arm identification for linear bandits (Xu et al. 2018), we propose the first algorithm for best-arm identification for generalized linear bandits, provide theoretical guarantees on its accuracy and sampling efficiency, and evaluate its performance in various scenarios via simulation.

\KEYWORDS

best arm identification, generalized linear bandits, sequential clinical trial

1 Introduction

The multi-armed bandit problem is a prototypical model for optimizing the tradeoff between exploration and exploitation. We consider a pure-exploration version of the bandit problem known as the best-arm identification problem, where the goal is to minimize the number of arm pulls required to select an arm that is - with sufficiently high probability – sufficiently close to the best arm. We assume that each arm has an observable vector of covariates or features, and there is an unknown vector of parameters (of the same dimension as the vector of features) that is common across arms. Whereas in a linear bandit the mean reward of an arm is the linear predictor (i.e., the inner product of the parameter vector and the feature vector), in our generalized linear model the mean reward is related to the linear predictor via a link function, which allows for mean rewards that are nonlinear in the linear predictor, as well as binary or integer rewards (via, e.g., logistic or Poisson regression). Hence, with every pull of an arm, the decision maker refines the estimate of the unknown parameter vector and learns simultaneously about all arms.

Our motivation for studying this version of the bandit problem comes from drug design. The field of drug design is immense (Martin 2010). There are different types of drugs (e.g., organic small molecules or biopolymer-based, i.e. biopharmaceutical, drugs), several phases of development, and a variety of experimental and computational methods used to select compounds with favorable properties. Three characteristics of our bandit problem make it ideally suited to certain subproblems within the drug design process. The first characteristic is the desire to run as few experiments as possible with the goal of selecting a particular compound (i.e., an arm in our model) for the next stage of analysis or testing. Hence, the actual performance of the selected arms during testing is of secondary importance (these pre-clinical tests do not involve human subjects), giving rise to a pure exploration problem. Second, drugs can often be described by a vector of features, which may include calculated properties of molecules, or the three-dimensional structure of the molecule and how it relates to the three-dimensional structure of the biological target. Finally, in many of these experimental settings, a generalized linear model provides a better fit than the linear model. This is true when the outcome is binary; e.g., in animal experiments, there is either survival or death, or either the presence or absence of disease, and some in vitro assays are qualitative (i.e., binary output). It is also true when the outcome is counting data (e.g., the number of cells or animals that survived or died), or when the relationship between the outcome and the best linear predictor is nonlinear.

There may also be other potential applications of this model, e.g., where the aim is to choose the best (non-personalized) ad for an advertising campaign, and the observed output is binary (e.g., the ad led to a purchase or a click-through).

Literature Review. A recent review (Bubeck and Cesa-Bianchi 2012) categorizes bandit problems into three types: Markovian (maximizing discounted rewards, often in a Bayesian setting, and exemplified by the celebrated result in Gittins (1979)), stochastic (typically minimizing regret in a frequentist setting, with independent and identically distributed rewards (Lai and Robbins 1985)), and adversarial (a worst-case setting, where an adversary chooses the rewards (Auer et al. 2002)). We study the stochastic problem, but consider best-arm identification rather than regret minimization. In passing, we note that the ranking-and-selection problem (Kim and Nelson 2006) in the simulation literature has a similar goal to best-arm identification, although evaluating an alternative requires running simulation experiments.

There is a vast literature on best-arm identification in the multi-arm bandit setting with independent arms, where pulling one arm does not reveal any information about the reward of other arms (Even-Dar et al. 2006, Bubeck et al. 2009, Audibert et al. 2010, Gabillon et al. 2012, Kalyanakrishnan et al. 2012, Karnin et al. 2013, Chen et al. 2014). Different algorithms have been developed for variants of this setting, most of which use gap-based exploration where arms are played to reduce the uncertainty about the gaps between the rewards of pairs of arms. Because playing an arm reveals information only about that arm, each arm needs to be played several times to reduce the uncertainty about its reward. As such, these algorithms can be practically implemented only when the number of available arms is relatively small.

Rather than assume independent arms, our formulation considers parametric arms, where each arm has a covariate vector and there is an unknown parameter vector that is common across arms. There has been considerable work on the linear parametric bandits (i.e., the mean reward of an arm is the inner product of its covariate vector and the parameter vector) under the minimum-regret objective (e.g., Auer (2002), Rusmevichientong and Tsitsiklis (2010) and references therein) as well as alternative probabilistic models of arm dependence (e.g., Russo and Van Roy (2014) and references therein). In addition, contextual bandit models allow additional side information in each round, which can model, e.g., patient information in clinical trials or consumer information in online advertising (e.g., Wang et al. (2005), Seldin et al. (2011), Goldenshluger and Zeevi (2013)). Relevant for our purposes, the generalized linear parametric bandit, which uses an inverse link function to relate the linear predictor and the mean reward, has been studied under regret minimization (Filippi et al. 2010, Li et al. 2017).

However, relatively little work exists on best-arm identification in parametric bandits. The first analysis of best-arm identification in linear bandits (Soare et al. 2014) proposed a static exploration algorithm that was inspired by the transductive experiment design principle (Yu et al. 2006). This algorithm determines the sequence of to-be-played actions before making any observations and fails to adapt to the observed rewards. Recently (Xu et al. 2018), major progress has been made via the first adaptive algorithm for best-arm identification in linear bandits. These authors design a gap-based exploration algorithm by employing the standard confidence sets constructed in the literature for linear bandits under regret minimization.

Our Contribution. We adapt the gap-based exploration algorithm of Xu et al. (2018) from the linear setting to the generalized linear case, which requires us to derive confidence sets for reward gaps between different pairs of arms. In the regret-minimization setting, the typical approach to the linear bandit is to develop a confidence set for the unknown parameter vector that governs the rewards of all arms, whereas in the best-arm setting, a confidence set on the reward gaps is needed. In the best-arm identification for the linear bandit (Xu et al. 2018), the authors were able to convert the confidence set for the parameter vector into efficient confidence sets for the reward gaps. However, this approach breaks down in the generalized linear bandit; i.e., naively converting the confidence set for the parameter vector in Filippi et al. (2010) into confidence sets for reward gaps between arms leads to extremely loose confidence sets, which strongly degrades the performance of the gap-based exploration algorithm. Rather than use this indirect method, we build gap confidence sets directly from the data.

The remainder of this paper is organized as follows. In Section 2, we formulate the best-arm identification problem for generalized linear bandits. We describe our algorithm in Section 3 and establish theoretical guarantees in Section 4. We provide simulation results in Section 5 and offer concluding remarks in Section 6.

2 Problem Formulation

Consider a decision maker who is seeking to find the best among a set of $K$ available arms. We let $[K]=\{1,2,\cdots,K\}$ denote the set of possible arms. There is a feature vector $x^{a}\in\mathbb{R}^{d}$ associated with arm $a$, for $a\in[K]$. These feature vectors are known to the decision maker and each summarizes the available information about the corresponding arm. We employ a generalized linear model (McCullagh and Nelder 1989) and assume that the reward of each arm $a$ has a particular distribution in an exponential family with mean

[TABLE]

where $\theta\in\mathbb{R}^{d}$ is an unknown parameter that governs the reward of all arms and $\mu:\mathbb{R}\to\mathbb{R}$ is a strictly increasing function known as the inverse link function. Different choices for the function $\mu(z)$ in (1) result in modeling different reward structures inside the exponential family. For example, choosing $\mu(z)=e^{z}$ and $\mu(z)=1/(1+e^{-z})$ correspond to a Poisson regression model and a logistic regression model, respectively.

The decision maker chooses an arm to play in each round $t$. If arm $a_{t}$ is played in round $t$, a stochastic reward $r_{t}$ is observed, which satisfies

[TABLE]

Let $a^{*}=\argmax_{a\in[K]}~{}\mu_{a}$ be the optimal arm; i.e., the arm with the highest expected reward. By exploring different arms, the decision maker is trying to find the optimal arm as soon as possible based on the noisy observations. Let $\tau$ be a stopping time that dictates whether enough evidence has been gathered to declare the optimal arm. The declared optimal arm is denoted by $\hat{a}_{\tau}$. An exploration strategy can be represented by $\mathcal{S}=(\mathcal{A},\tau)$, where, at any time $t$, $\mathcal{A}$ is a function mapping from the previous observations $\{(a_{i},r_{i})\}_{i=1}^{t-1}$ to the arm $a_{t}$ to be played next, and $\tau$ determines whether enough information has been gathered to declare the optimal arm, $\hat{a}_{\tau}$. Because finding the exact optimal arm may require a prohibitively large amount of exploration, the performance of an exploration strategy is evaluated via the following relaxed criterion.

Definition 2.1

Given $\epsilon>0$ and $\delta\in(0,1)$, an exploration strategy $\mathcal{S}=(\mathcal{A},\tau)$ is said to be $(\epsilon,\delta)-$optimal if

[TABLE]

In this definition, $\epsilon$ denotes an acceptable region around the optimal arm and $1-\delta$ represents the confidence in identifying an arm within this region. This criterion relaxes the notion of optimality by allowing the exploration strategy to return a sufficiently good – but not necessarily optimal – arm. With this definition in place, the decision maker’s goal is to design an $(\epsilon,\delta)-$optimal exploration strategy with the smallest possible stopping time.

Before proceeding to the algorithm, we introduce additional notation and state a set of regularity assumptions. We let $\mathcal{X}=\{x^{a}\}_{a\in[K]}$ denote the set of feature vectors and assume that feature vectors, the unknown reward parameter and the rewards are bounded; i.e., there exist $S,L,R>0$ such that $\|\theta\|_{2}\leq S$, $\|x\|_{2}\leq L~{}\forall x\in\mathcal{X}$, and $r_{t}\leq R$ almost surely for all $t$. We also assume that $\mu$ is continuously differentiable, Lipschitz continuous with constant $k_{\mu}$ and satisfies $c_{\mu}=\inf_{\theta:\|\theta\|\leq S,~{}x\in\mathcal{X}}\dot{\mu}(\theta^{\intercal}x)>0$. For example, in the case of logistic regression, $R=1,k_{\mu}=1/4$ and $c_{\mu}$ depends on $\sup_{\theta:\|\theta\|\leq S,~{}x\in\mathcal{X}}|\theta^{\intercal}x|$. We define the gap between any two arms $i,j\in[K]$ to be $\Delta(i,j)=\mu_{i}-\mu_{j}=\mu(\theta^{\top}x^{i})-\mu(\theta^{\top}x^{j})$ and define the optimal gap associated to an arm $i\in[K]$ as

[TABLE]

Finally, for any positive semi-definite matrix $A$, we let $\|x\|_{A}=\sqrt{x^{\top}Ax}$.

3 The Proposed Algorithm

In this section, we propose an exploration strategy for the problem formulated in Section 2. Following Xu et al. (2018), our algorithm consists of the following steps:

Build confidence sets for the pairwise gaps between arms, 2. 2.

Identify the potential best arm and an alternative arm that has the most ambiguous gap with the best arm, 3. 3.

Play an arm to reduce this ambiguity.

These steps are repeated sequentially until the ambiguity in step 2 drops below a certain threshold.

The confidence sets are derived in Subsection 3.1 and the algorithm is presented in Subsection 3.2.

3.1 Confidence Sets for Gaps

To build the confidence sets for reward gaps, we follow ideas in Filippi et al. (2010) but develop confidence sets directly for gaps instead of arm rewards.

Let $(a_{1},r_{1}),(a_{2},r_{2}),\cdots,(a_{t-1},r_{t-1})$ be the history of actions played and random rewards observed prior to period $t$, and let $x_{l}=x^{a_{l}}$ be shorthand notation for the feature vector associated with the arm played in period $l$. For any $t>0$, let $M_{t}=\sum_{l=1}^{t-1}x_{l}x_{l}^{\top}$ be the empirical covariance matrix and assume it is nonsingular for any $t>E$ for some fixed value $E>d$. We let $\lambda_{0}>0$ be the minimum eigenvalue of $M_{E+1}$ and define $\kappa=\sqrt{3+2\log(1+2L^{2}/\lambda_{0})}$. Given the observations by the start of period $t$, the Maximum Likelihood (ML) estimate of the reward parameter, $\theta_{t}$, solves the equation

[TABLE]

Based on the estimated reward parameter, we can take

[TABLE]

as an estimate for the gap between arms $i,j\in[K]$, which is a function of the observations made prior to period $t$.

Given $\delta\in(0,1)$, let

[TABLE]

where $\alpha$ is a tunable parameter and $C_{t}$ is a time-varying quantity that scales the width of the confidence sets for all pairs of arms. We set

[TABLE]

in the theoretical analysis in Section 4, and set $\alpha$ so as to achieve robust performance across a variety of scenarios in the computational study in Section 5. We consider the following confidence set for the gap between any two arms $i,j\in[K]$ based on the observations made prior to period $t$:

[TABLE]

The confidence set defined in (8) is centered around the estimated gap between the two arms in (5) and, as we will show in Section 4, contains the true gap with high probability. To further simplify the representation, we let

[TABLE]

represent the width of the confidence set in (8).

Although we follow the basic approach in Filippi et al. (2010) in deriving these confidence sets, it is worth noting that they cannot be deduced from the results presented in Filippi et al. (2010). More precisely, an immediate use of the results in that paper gives rise to a very loose confidence set that is similar to (8) except that $\max_{c,c^{\prime}\in[c_{\mu},k_{\mu}]}\|cx^{i}-c^{\prime}x^{j}\|_{M_{t}^{-1}}$ is replaced by the looser factor $\left(\|x^{i}\|_{M_{t}^{-1}}+\|x^{j}\|_{M_{t}^{-1}}\right)$.

\TableSpaced

\DoubleSpacedXI

3.2 The Algorithm

With the confidence sets established, we are now ready to describe our proposed algorithm. Details and successive steps of the proposed algorithm are presented in Algorithms 1-3. Following the gap-based exploration scheme, the proposed algorithm consists of two major components that are described below.

Selecting a Gap to Explore:

The algorithm starts by playing $E>d$ random arms such that the empirical covariance matrix $M_{E+1}$ is nonsingular. At any subsequent period $t>E$, the algorithm first finds the empirically best arm $i_{t}=\argmax_{i\in[K]}\mu(\theta_{t}^{\top}x^{i})$. Then, to check whether this arm is within $\epsilon$ distance of the true optimal arm, the algorithm takes a pessimistic approach. In particular, it finds another arm that is the most advantageous over arm $i_{t}$ within the gap confidence sets; i.e., $j_{t}=\argmax_{j\in[K],j\neq i_{t}}\Delta_{t}(j,i_{t})+\beta_{t}(j,i_{t})$. Note that this pessimistic gap consists of two components: the estimated gap and the uncertainty in the gap.

If this pessimistic gap is less than $\epsilon$, the algorithm stops and declares the empirically best arm as the optimal arm. Otherwise, it selects an arm to reduce the uncertainty component in the identified gap. According to (9), this uncertainty is governed by $\|y_{t}\|_{M_{t}^{-1}}$ where $y_{t}=c_{1}x^{i_{t}}-c_{2}x^{j_{t}}$ and $c_{1},c_{2}=\argmax_{c,c^{\prime}\in[c_{\mu},k_{\mu}]}\|cx^{i_{t}}-c^{\prime}x^{j_{t}}\|_{M_{t}^{-1}}$.

Selecting an Arm:

While there are different ways to reduce $\|y_{t}\|_{M_{t}^{-1}}$, we follow the approach in Xu et al. (2018). With a slight abuse of notation, let us define $M_{\mathbf{z}_{n}}=\sum_{i=1}^{n}z_{i}z_{i}^{\top}$ for any sequence of feature vectors $\mathbf{z}_{n}=(z_{1},z_{2},\ldots,z_{n})\in\mathcal{X}^{n}$, where $n$ represents a generic time period. Let $\mathbf{z}_{n}^{*}(y)=\argmin_{\mathbf{z}_{n}\in\mathcal{X}^{n}}\|y\|_{M^{-1}_{\mathbf{z}_{n}}}$ be the sequence of feature vectors that would have minimized the uncertainty in the direction of $y$. For each arm $a\in[K]$, let $p_{a}(y)$ denote the relative frequency of $x^{a}$ appearing in the sequence $\mathbf{z}_{n}$ when $n\to\infty$. As has been shown in Section 5.1 of Xu et al. (2018),

[TABLE]

where $w^{*}(y)\in\mathbb{R}^{K}$ is the solution of the linear program

[TABLE]

To minimize the uncertainty in the direction of $y_{t}$, the algorithm plays the arm

[TABLE]

where $T_{a}(t)$ is the number of times arm $a$ has been played prior to period $t$.

4 Theoretical Analysis

In this section, we provide theoretical guarantees for the performance of the proposed algorithm. We prove that the algorithm indeed finds an $(\epsilon,\delta)-$optimal arm in Subsection 4.1 and provide an upper bound on the stopping time of the algorithm in Subsection 4.2.

4.1 $(\epsilon,\delta$)-Optimality of the Proposed Algorithm

We start by proving that the confidence sets constructed in Section 3.1 hold with high probability at all times.

Proposition 4.1

Fix $\delta$ and $t$ such that $0<\delta<\min(1,d/e)$ and $t>\max(2,d)$, where $d$ is the feature dimension. Then the following holds with probability at least $1-\delta$:

[TABLE]

where

[TABLE]

The proof is a slight modification of the proof of Proposition 1 in Filippi et al. (2010).

Proof 4.2

Proof. Let $g_{t}(\beta)=\sum_{l=1}^{t-1}\mu(\beta^{\top}x_{l})x_{l}$. According to (4), the ML estimate $\theta_{t}$ satisfies $g_{t}(\theta_{t})=0$. By the mean value theorem, there exist points $z,z^{\prime}$ with $c=\dot{\mu}(z),c^{\prime}=\dot{\mu}(z^{\prime})$ such that

[TABLE]

On the other hand, by the fundamental theorem of calculus, we have

[TABLE]

where

[TABLE]

The definition of $g_{t}(\beta)$ implies that for any $\beta$,

[TABLE]

It follows that $G_{t}\succeq c_{\mu}M_{t}\succeq c_{\mu}M_{d}\succeq\lambda_{0}I\succ 0$. Hence, $G_{t}$ and $G_{t}^{-1}$ are positive definite and nonsingular. Therefore, from (15) and (16), it follows that

[TABLE]

The inequality $G_{t}\succeq c_{\mu}$ implies that $G_{t}^{-1}\preceq c_{\mu}^{-1}M_{t}^{-1}$, and hence $\|y\|_{G_{t}^{-1}}\leq\frac{1}{\sqrt{c_{\mu}}}\|y\|_{M_{t}^{-1}}$ for arbitrary $y\in\mathbb{R}^{d}$. Thus, by (17) we get

[TABLE]

As has been shown in the proof of Proposition 1 in Filippi et al. (2010),

[TABLE]

holds with probability at least $1-\delta$. Combining (18) and (19) shows that

[TABLE]

holds with probability at least $1-\delta$. Finally, because $c,c^{\prime}\in[c_{\mu},k_{\mu}]$, taking the maximum over $c,c^{\prime}$ on the right side of (20) completes the proof. \Halmos

The following theorem, which is a direct consequence of Proposition 4.1, shows that the confidence sets in (8) contain the true gaps at all times with high probability. We define event $\mathcal{E}$ to be

[TABLE]

and throughout this section we set $\alpha$ according to (7).

Theorem 4.3

*Let $\delta$ be such that $0<\delta<\min(1,d/e)$. Then event $\mathcal{E}$ occurs with probability at least $1-\delta$. *

Proof 4.4

Proof. For any $t>\max(2,d)$, let $\delta_{t}=\frac{6\delta}{\pi^{2}t^{2}}$ and define the event $\mathcal{E}_{t}$ as

[TABLE]

Applying Proposition 4.1 for $\delta_{t}$ implies that $\mathbb{P}[\mathcal{E}_{t}]\geq 1-\delta_{t}$. The union bound then gives

[TABLE]

The proof is completed by noting that $\mathcal{E}=\cap_{t=1}^{\infty}\mathcal{E}_{t}$.\Halmos

With the confidence sets established, the next theorem proves $(\epsilon,\delta)-$optimality of the proposed algorithm.

Theorem 4.5

Let $\epsilon>0$ and $0<\delta<\min(1,d/e)$ be arbitrary. Then the proposed algorithm is $(\epsilon,\delta)-$optimal; i.e., at its stopping time, the algorithm returns an arm $\hat{a}_{\tau}$ such that

[TABLE]

Proof 4.6

Proof. Let $\tau$ be the stopping time of the algorithm and let $\hat{a}_{\tau}$ be the returned arm. Suppose that $\Delta(a^{*},\hat{a}_{\tau})>\epsilon$. Then, according to line 5 of Algorithm 1, we have

[TABLE]

From this, we get that

[TABLE]

which means that event $\mathcal{E}$ does not occur. According to Theorem 4.3, this can happen with probability at most $\delta$. Thus, $\Delta(a^{*},\hat{a}_{\tau})>\epsilon$ happens with probability at most $\delta$. \Halmos

4.2 An Upper Bound on the Stopping Time

In this subsection, we study the sample complexity of the proposed algorithm. The following theorem provides an upper bound on the number of experiments the proposed algorithm needs to carry out before identifying the optimal arm. Before stating the result, let us take $\Delta_{\min}=\min_{i\in[K]}\Delta_{i}$ to be the smallest gap and for any $\epsilon>0$, define

[TABLE]

which represents the complexity of the exploration problem in terms of the problem parameters.

Theorem 4.7

Let $\tau$ be the stopping time of the proposed algorithm. Then

[TABLE]

is satisfied with probability at least $1-\delta$. In asymptotic notation, (23) can be expressed as

[TABLE]

Theorem 4.7 provides an upper bound on the stopping time of the proposed algorithm in terms of the parameters of the exploration problem. As expected, the number of experiments required by the proposed algorithm before declaring a near-optimal arm decreases in the reward tolerance ($\epsilon$) and the error probability ($\delta$), and increases in the number of features ($d$) and arms ($K$). In terms of the dependence on dimension ($d$), reward bound ($R$) and the complexity parameter ($H_{\epsilon}$), the sample complexity in (24) is similar to that derived in Xu et al. (2018)** for linear bandits. The main difference is the appearance of the factor $\kappa/c_{\mu}$ in the complexity bound (24), which encodes the difficulty of learning the inverse link function $\mu$. As the inverse link function $\mu$ becomes flatter on the boundaries of the input domain (i.e., has smaller $c_{\mu}$), more samples are required to distinguish between different pairs of arms.

The remainder of this subsection is devoted to proving Theorem 4.7, which requires a set of preliminaries. We start by introducing some additional notation. For any $i,j\in[K]$ and any $t$, let

[TABLE]

and define $y^{ij}_{t}=c_{1}^{ij}(t)x^{i}-c_{2}^{ij}(t)x^{j}$. Note that with this notation, $y_{t}$ defined in Algorithm 2 can be represented as $y_{t}=c_{1}^{i_{t}j_{t}}(t)x^{i_{t}}-c_{2}^{i_{t}j_{t}}(t)x^{j_{t}}$. For any $y\in\mathbb{R}^{d}$, let $\rho(y)$ be the optimal value of the linear program in (11); i.e.,

[TABLE]

For any two real numbers $a,b$, we use the shorthand notation $a\vee b=\max(a,b)$.

Our proof for Theorem 4.7 relies on a number of results from the literature, which we state here for completeness. The following two lemmas are proved in Xu et al. (2018)**.

Lemma 4.8** **(Lemma 1 of Xu et al. (2018))

For any $y\in\mathbb{R}^{d}$, we have

[TABLE]

where

[TABLE]

Lemma 4.9** **(Lemma 4 of Xu et al. (2018))

When event $\mathcal{E}$ holds, $B(t)$ satisfies the following bounds:

If either $i_{t}$ or $j_{t}$ is the best arm:

[TABLE] 2. 2.

If neither $i_{t}$ nor $j_{t}$ is the best arm:

[TABLE]

The following lemma is proved in Antos et al. (2010)**.

Lemma 4.10** **(Proposition 6 of Antos et al. (2010))

For any $t>0$, let $q(t)=a\sqrt{t}+b$ and $l(t)=\log(t)$, for some $a>0$. Define t^{*}=\frac{4}{a^{2}}\left[\log\Bigl{(}\frac{4}{a^{2}}\Bigr{)}-b\right]^{2}. Then, for any positive $t$ such that $t\geq t^{*}$, we have $q(t)>l(t)$.

The following lemma establishes an upper bound on the solution of the linear program in (11).

Lemma 4.11

Let $y_{t}=c_{1}^{i_{t}j_{t}}(t)x^{i_{t}}-c_{2}^{i_{t}j_{t}}(t)x^{j_{t}}$ and let $w^{*}(y_{t})$ be the solution of (11). Then we have

[TABLE]

Proof 4.12

Proof. Let $w^{\prime}\in\mathbb{R}^{K}$ be such that $w^{\prime}_{i_{t}}=c_{1}^{i_{t}j_{t}}(t),w^{\prime}_{j_{t}}=-c_{2}^{i_{t}j_{t}}(t)$ and all other elements of $w^{\prime}$ are zero. Clearly, $w^{\prime}$ satisfies the constraint in (11). Therefore, we have

[TABLE]

\Halmos

With the above lemmas in place, we are ready to prove the following theorem.

Theorem 4.13

If event $\mathcal{E}$ occurs, then the stopping time of the proposed algorithm satisfies

[TABLE]

Proof 4.14

Proof. Suppose that event $\mathcal{E}$ occurs. Let $k\in[K]$ be an arbitrary arm, and let $t_{k}<\tau$ be the last round in which arm $k$ was pulled. Because $B(t_{k})>\epsilon$, Lemma 4.9 implies that

[TABLE]

which in turn gives rise to the following three inequalities:

[TABLE]

Rearranging (26) yields

[TABLE]

From Lemma 4.8, we have

[TABLE]

which by substituting $\beta_{t_{k}}(i_{t_{k}},j_{t_{k}})=C_{t_{k}}\|y_{t_{k}}\|_{M_{t_{k}}^{-1}}$ from (9) gives

[TABLE]

Combining (27) and (28) gives

[TABLE]

Because arm $k$ was selected by Algorithm 3 in round $t_{k}$, it follows that

[TABLE]

Then, from (29), we get

[TABLE]

Note that $\tau=1+\sum_{k=1}^{K}T_{k}(\tau)$. Hence, it follows from (22) and (31) that

[TABLE]

which completes the proof. \Halmos

We are now in a position to prove Theorem 4.7.

Proof 4.15

Proof of Theorem 4.7. Suppose that event $\mathcal{E}$ holds. Then

[TABLE]

where (32) follows by Jensen’s inequality and the logarithmic arithmetic-geometric mean inequality (i.e., \sqrt{\log x\log y}\leq\frac{1}{2}(\log x+\log y)\leq\log\Bigl{(}\frac{x+y}{2}\Bigr{)}). Applying the inequality $x+y\leq(\sqrt{x}+\sqrt{y})^{2}$ for any $x,y>0$ to (32) yields

[TABLE]

Let $c_{1}=\frac{4\kappa R\sqrt{d}}{c_{\mu}}\sqrt{H_{\epsilon}},~{}c_{2}=\pi\sqrt{\frac{d}{6\delta}}+1$ and define $a=\frac{1}{c_{1}\sqrt{c_{2}}},~{}b=-\frac{\sqrt{K+1}}{c_{1}}$. With a change of variable $z=c_{2}\tau$, (33) can be written as

[TABLE]

According to Lemma 4.10, this is possible only if

[TABLE]

Substituting for $a,b$ and $z$ shows that (33) holds only if

[TABLE]

On the other hand, according to Theorem 4.3, event $\mathcal{E}$ holds with probability at least $1-\delta$. This completes the proof.\Halmos

5 Numerical Results

In this section, we test the proposed algorithm in various scenarios using synthetic data in Subsection 5.2 and drug design data in Subsection 5.3. The experimental details are described in Subsection 5.1.

5.1 Experimental Details

All scenarios considered here involve a logistic inverse link function with binary observations. Recall that the proposed algorithm starts with an exploratory phase in which $E$ random arms are played to ensure that $M_{E+1}$ is nonsingular. Although we have not optimized this parameter in our numerical studies, we use $E=\min(K,3d)$ throughout all simulations here.

The confidence sets defined in (6) and (8) rely on the free parameter $\alpha$. Our theoretical analysis in Section 4 suggests the use of $\alpha=\frac{2\kappa R}{c_{\mu}}$ from (7). However, the theoretical confidence sets derived for generalized linear settings in Filippi et al. (2010)* are loose and the performance can be improved by scaling these sets. Although Section 4.2 of Filippi et al. (2010)** uses an asymptotic analysis to shrink these sets in the regret-minimization setting (a finite-sample result in Theorem 1 of Li et al. (2017)** provides a rigorous justification of this analysis), we simply note that all the rewards and reward gaps are in $[0,1]$ in the case of a logistic reward function while the width of the confidence set defined in (8) may be larger than one. Consequently, we set $\alpha$ so as to constrain the confidence set widths in (8) to the interval $[0,1]$, which can be achieved by setting $\frac{2\kappa R}{c_{\mu}}$ times $\max_{i,j\in[K]}\beta_{E(i,j)}$ equal to 1, yielding via (6) and (9)*

[TABLE]

which is determined at the end of the exploratory phase of length $E$.

Throughout our simulations, we fix the confidence parameter $\delta=0.05$, which theoretically guarantees the accuracy of the proposed algorithm with $95\%$ confidence. The results in Figs. 1-5 are generated by averaging over tens of random realizations of the problem. In all cases, the proposed algorithm achieved the desired accuracy (as dictated by $\epsilon$) with probability more than 0.95.

5.2 Synthetic Data Simulations

*To create a synthetic data set, we generate the reward parameter $\theta$ by sampling from a *$d-$dimensional standard multivariate Gaussian distribution. The feature vector associated with each of the $K$ arms is randomly generated such that each of its components is uniformly distributed in $[-1,1]$. Acting as if we do not know the reward parameter $\theta$, we run the proposed algorithm to identify the best among the available $K$ arms under different scenarios. Recalling that the parameters $L$, $R$ and $S$ are defined via $\|x\|_{2}\leq L$, $r_{t}\leq R$ and $\|\theta\|_{2}\leq S$, we have $L\leq d$ and $R=1$. Note that the parameter $S$ is used only in determining $c_{\mu}$, which we treat as given in our simulations; hence, there is no need to specify $S$.

We first fix the tolerance parameter $\epsilon=0.1$ and the feature dimension $d=10$, and plot the stopping time $\tau$ as a function of the number of arms $K$ (Fig. 1). We then fix $\epsilon=0.1$ and the number of arms $K=100$, and compute $\tau$ as a function of the feature dimension $d$ (Fig. 2). These plots reveal that the amount of exploration undertaken by the proposed algorithm depends more strongly on the feature dimension $d$ than on the number of arms $K$. This is due to the fact that the proposed algorithm builds confidence sets for the reward parameter, which is of dimension $d$, rather than for the reward of each arm separately.

Fixing the feature dimension $d=10$ and the number of arms $K=100$, we plot the stopping time $\tau$ versus the tolerance parameter $\epsilon$ (Fig. 3). As expected, a smaller tolerance requires a larger amount of exploration in order to guarantee that the identified arm is within the acceptable distance of the optimal arm.

Recall that exploration algorithms for independent arms require playing each arm many times. To assess the improvement over this approach, we implement the GapE algorithm (Gabillon et al. 2012), which treats arms as independent with binary rewards, to an instance of this synthetic scenario with $K=50$ arms, $d=10$ features and $\epsilon=0.1$. On average, the GapE algorithm requires 60k experiments before identifying a near-optimal arm, while our proposed algorithm identifies a near-optimal arm after only 436 experiments.

We know of no other algorithms that are natural comparisons to our proposed algorithm: the algorithms in Xu et al. (2018)* and Soare et al. (2014)** use a linear model that does not accommodate binary outcomes, and regret-minimization algorithms are not typically tested against best-arm identification algorithms because the former do not have any stopping criteria, which precludes us from comparing our algorithm to the one in Filippi et al. (2010)**.*

5.3 Drug Design Simulations

In this subsection, we consider the problem of identifying the most effective among a large set of designed drugs for a certain disease. We use the ChEMBL data repository (Gaulton et al. 2011) and download the data set describing 982 small molecules designed for targeting the Cyclin-Dependent Kinase 2/Cyclin A protein complex. The covariates for each small molecule are associated with its three-dimensional structure, which can be represented as a 1024-dimensional binary vector. To construct such a representation, we compute the Morgan Fingerprint (Riniker and Landrum 2013) of each molecule with radius 4. Then, by applying a Principal Components Analysis (PCA) to these sparse high-dimensional binary vectors, we generate a feature vector with the desired dimension.

The generalized linear model allows for the random experimental output to have values that are continuous, integer or binary, depending on the link function that is used. If we want to use a model with binary data (i.e., the logit inverse link function), we would need binomial or Bernoulli data for each arm. However, we have been unable to find data like these. Rather, the output data in the data set from Gaulton et al. (2011)* are the affinity of each small molecule, as measured by the IC50 (the concentration of an inhibitor where the biological function, typically the binding to a protein, is reduced by half). The IC50 is a deterministic number (there are no data in Gaulton et al. (2011)** on the possible measurement errors) that is an indicator of the effectiveness of each small molecule.*

*To adapt these outcome data into our modeling framework, we convert these deterministic IC50 values into random binary outcomes as follows. First, we assume the drug produces a cure if the IC50 value is less than *$10^{6}$*nM (nanomolar), and does not produce a cure if the IC50 is greater than or equal to *$10^{6}$nM. Next, given the feature vectors and these artificial cure outcomes, we fit a logistic regression model and estimate the cure rate of each drug. Finally, by assuming that these estimated cure rates are the unknown expected rewards of the drugs, we apply the proposed algorithm to find the best drug; i.e., in our simulations, each time a drug is tested, we observe a Bernoulli random variable with its mean equal to the drug’s expected reward.

We first fix the tolerance parameter $\epsilon=0.1$ and the number of arms $K=400$ by randomly choosing 400 of the 982 molecules in Gaulton et al. (2011)**. The stopping time $\tau$ is plotted against the feature dimension $d$, which is varied using PCA (Fig. 4). As in the synthetic data simulations, this relationship is increasing and – when the number of features is smaller than 50 – the proposed algorithm finishes exploration before all the arms are explored.

*To roughly assess the effectiveness of the proposed algorithm in its choice of best arm, we note that the 982 IC50 values in the data set range from 0.3nM to *$1.5\times 10^{6}$nM, and that the best (i.e., lowest-valued) among the 400 molecules in the experiment has an IC${}_{50}=0.4$nM. The algorithm’s declared optimal arm has an IC${}_{50}=11.99$*nM. However, because considerable information is lost when we binarize the data (e.g., two drugs with IC50 values of 0.5nM and *$10^{5}$nM would both be assigned a binary value of one for a cure), we should not expect our algorithm to identify the drug with the lowest IC50.

We next fix $\epsilon=0.1$ and the feature dimension $d=20$, and plot the stopping time $\tau$ versus the number of arms $K$ (Fig. 5), where for each value of $K$ we randomly choose $K$ arms out of the 982 small molecules in Gaulton et al. (2011)**. When there are more than 200 arms in Fig. 5, the number of experiments carried out by the proposed algorithm is smaller than the total number of arms.

6 Conclusion

Motivated by drug design, we consider the best-arm identification problem for a generalized linear bandit, where arms correspond to drugs, the covariate vectors describe biological properties of the drugs, and the inverse link function in the generalized linear model captures the nonlinear and perhaps binary relationship between the covariates and the experimental outcomes. We use ideas from Filippi et al. (2010)* on building confidence sets, but apply these directly to confidence sets for gaps between pairs of arms rather than confidence sets for rewards of arms, which allows us to generalize the results in Xu et al. (2018)** from the linear bandit model to the generalized linear bandit model. After performing our analysis, we became aware of recent work (Li et al. 2017) that improves the regret rate in Filippi et al. (2010)** by a factor of $\sqrt{d}$. We leave it to future work to investigate whether the results in Li et al. (2017)** can be used to improve our theoretical bounds.*

Numerical results show that – at least with $\delta=0.05$ and $\epsilon=0.1$ (i.e., being 95% confident of selecting an arm with a cure rate within 10% of optimal) – can be identified very quickly: when there are a large number of arms, the number of experiments performed can be much less than the total number of arms (Fig. 5). This is achieved by learning about all arms whenever an arm is played. Our algorithm, coupled with an approach to dimensionality reduction of the feature dimension such as PCA or machine learning tools, has the potential to streamline some nonlinear problems in drug design and other experimental settings.

\ACKNOWLEDGMENT

This research was supported by the Graduate School of Business, Stanford University (L.M.W.) and a Stanford Graduate Fellowship (A.K.).