Correlational Dueling Bandits with Application to Clinical Treatment in Large Decision Spaces

1 Introduction

In many online learning settings, particularly those that involve human feedback, reliable feedback is often limited to pairwise preferences instead of real valued feedback. Examples include implicit or subjective feedback for information retrieval and recommender systems, such as clicks on search results, and subjective feedback on the quality of recommended care [Chapelle et al., 2012; Sui and Burdick, 2014]. This setup motivates the dueling bandits problem [Yue and Joachims, 2009], which formalizes the problem of online regret minimization via preference feedback (e.g., choosing a pair of arms to be compared at each time step). Many dueling bandits algorithms [Yue and Joachims, 2009; 2011; Zoghi et al., 2014; Ailon et al., 2014; Komiyama et al., 2015;

Wu and Liu, 2016] have been developed for efficiently computing this problem with independent arms. However, these algorithms are not efficient in situations involving a large number of dependent arms. Specifically, when the time horizon T is smaller than the number of arms K, it is hopeless to achieve low regret without leveraging structure among arms.

Our problem is motivated by clinical research for recovering motor function after severe spinal cord injury. Previous research [Harkema et al., 2011] has shown that electrical stimulation applied to the spinal cord via electrode arrays implanted in the epidural space over the lumbosacral area enables paralyzed patients to achieve full weight-bearing standing, improvements in stepping, and partial recovery of lost autonomic functions. Stimulation consists of electrical pulse trains applied to selected electrodes. The challenge is that the optimal stimulus pattern (the choice of active electrodes and their polarities, the pulse amplitude and width, and the pulse train frequency) varies significantly across patients. And even for the same patient, the response to the same stimulus has some variation across trials. Hence, clinicians must determine the optimal stimulus for each patient under noisy conditions, which currently a laborious and ad-hoc approach.

Figure 1 shows the clinical treatment procedure for standtraining of paraplegics. During a treatment/optimization session, new stimuli are recommended by the algorithm to be applied to the electrode implanted in the patient. The patient then attempts to stand using the given stimuli, and the observing clinicians compare the patient's standing performance. The total number of different stimulating configurations is ∼ 4.3 × 107 due to the complexity of electrodes, and so it is not feasible to search through the whole space. The goal is to develop a algorithm that can automatically select stimuli in order to quickly converge to good treatments.

Motivated by this application, we consider the problem of finding optimal stimuli based on the general setting of the multi-armed bandit problem. The classical bandit problem trades off between exploration and exploitation among a number of different arms, each having a quantifiable but stochastic reward with an initially unknown distribution. In contrast, for our clinical problem, the patient's motor response to stimulation is hard to quantify. Neither video motion capture nor electromyographic (EMG) recordings of muscle activity can yet provide a consistent and satisfactory measure of motor skill under stimulation. One reasonably reliable measure is that of pairwise comparisons, e.g., whether one stimulus more effective than another. While the patient's performance under a specific stimulus is hard to quantify in the clinical setting, we can obtain comparisons of stimuli which are tested within the short time period of one training session.

Our Contributions. In this paper, we show how to cast the problem of online learning of personalized clinical treatment as a dueling bandits problem with a correlated action space, which we call correlational dueling bandits . We present an algorithm which meets the demands of such clinical settings, and can effectively model such correlation dependencies to achieve good performance. Our algorithm takes advantage of the correlations among different arms to update the whole active set of arms instead of only updating the two dueling arms. This approach achieves fast convergence to the (near) optimal decisions regardless of the large decision space. We deployed CORRDUELas the first algorithmic approach to the control of spinal cord stimulation in clinical experiments. We find that CORRDUELcan identify a group of optimal stimuli and help paraplegic human patients to achieve full-weight standing.

2 RelatedWork

2.1 Multi-Armed Bandits

The stochastic multi-armed bandits problem [Robbins, 1952] refers to an iterative decision making problem in which one repeatedly chooses among K options, such as pulling one of K arms of a bandit machine. In each round, we receive a reward that depends on the arm being selected. Without loss of generality, assume that every reward is bounded between [0, 1]. The goal then is to minimize the cumulative regret compared to the best arm.

Popular algorithms for the stochastic setting include UCB (upper confidence bound) algorithms [Auer et al., 2002a; Bubeck and Cesa-Bianchi, 2012], and Thompson Sampling [Chapelle and Li, 2011; Russo and Van Roy, 2014].

In the adversarial setting, the rewards are chosen in an adversarial fashion, rather than sampled independently from some underlying distribution. In this case, regret is rephrased as the difference in the sum of rewards. The predominant algorithm for the adversarial setting is EXP3 [Auer et al., 2002b].

2.2 Correlated Bandits

The set of candidate actions is very large (or even infinite) in many applications. When that is the case, one must exploit dependencies between payoffs of different decisions in order to arrive at an efficient algorithm.

In some applications, the underlying problem comes equipped with a correlational structure. Various methods of introducing dependence include bandits on trees [Kocsis and Szepesv´ari, 2006], bandits with linear correlations [Dani et al., 2008; Abernethy et al., 2008; Abbasi-Yadkori et al., 2011; Gentile et al., 2014] or Lipschitz continuous payoffs [Kleinberg et al., 2008; Bubeck et al., 2008], and Gaussian payoffs [Srinivas et al., 2010].

2.3 Dueling Bandits

Dueling bandits problem [Yue et al., 2012], as a variant of the multi-armed bandits, takes (noisy) comparative feedback instead of real-valued feedback. It is under the general framework of preference learning (learning with preferential feedback). The dueling bandits problem can also be viewed as a special case of partial monitoring problems [Cesa-Bianchi et al., 2006]. Its problem setting naturally fits in with many applications such as information retrievals and recommender systems. The stochastic dueling bandits problem has been extensively studied in [Yue et al., 2012; Ailon et al., 2014; Zoghi et al., 2014; Komiyama et al., 2015; Wu and Liu, 2016].

Beyond the stochastic K-armed dueling bandits setting, other dueling bandit settings include multi-way preference feedback [Sui and Burdick, 2014], continuous-armed convex dueling bandits [Yue and Joachims, 2009], contextual dueling bandits which also introduces the von Neumann winner solution concept [Dud´ık et al., 2015], sparse dueling bandits that focus on the Borda winner solution concept [Jamieson et al., 2015], Copeland dueling bandits that focus on the Copeland winner solution concept [Zoghi et al., 2015], and adversarial dueling bandits [Gajane et al., 2015]. It would be interesting to study how to extend our analysis to these other settings as well.

3 Problem Statement

In the classical dueling bandits problem, at each iteration t, the following happens:

• The algorithm chooses a pair of actions b(1) (t) and b(2) (t) from a set of K possible actions.
• The algorithm duels b(1) (t) and b(2) (t) and receives (noisy) feedback corresponding to the winner.

Our procedure can be described as follows. There is a set of arms B = {b1, · · · , bK}, and a total number of T tests to be performed. At each time step, a pair of arms are chosen from the set B and a (noisy) comparison of them is observed. T is determined before we run the algorithm. The set of arms are correlated and T ≤ |B| = K in general.

We follow the original notation of the dueling bandit problem. For two arms bi and bj sampled from B, we write the comparison factor as

where P(bi ≻ bj) is the probability that bi dominates bj and ϵ(bi, bj) ∈ [-1/2, 1/2] represents the priority between bi and bj. We define bi ≻ bj ⇔ ϵ(bi, bj) > 0. We use the notation ϵi,j ≡ ϵ(bi, bj) for convenience. Note that ϵ(bi, bj) = -ϵ(bj, bi) and ϵ(bi, bi) = 0. We assume the distribution of reward for each arm is stationary so that all comparison factors converge in [-1/2,1/2]. We also assume w.l.o.g. that the arms are indexed in preferential order b1 ≻ b2 ≻ · · · ≻ bK so that there is one preferred arm.

Our goal is to minimize the total regret:

The total regret RT = 0 if we constantly choose b(t) = b1 during the experiment. RT = Θ(T) is linear w.r.t. T if we constantly choose b(t) ∈ B.

We also inherit two properties of the comparison factors from the original dueling bandit problem:

Strong Stochastic Transitivity. For any triplet of arms bi ≻ bj ≻ bk, we assume ϵi,k ≥ max{ϵi,j , ϵj,k}.

Stochastic Triangle Inequality. For any triplet of arms bi ≻ bj ≻ bk, we assume ϵi,k ≤ ϵi,j + ϵj,k. This can be viewed as a diminishing returns property.

Correlational Dueling Bandits. When the size of the decision set, K, is large, it is unavoidable to carry out a very large number of tests before the algorithm converges to its optimal solution. In some applications like our clinical example, each test is expensive and time-consuming. The number of tests - the time horizon of an algorithm - is often predetermined by clinical conditions. We thus augment the dueling bandits problem into correlational dueling bandits, which takes the correlations among arms into consideration. For any pair of arms (bi, bj) ∈ B2 , we consider the dependence between them are captured by some similarity function rij ∈ [0, 1], and it satisfies:

• rij = rji;
• rij = 0 ⇐⇒ bi and bj are not correlated;
• rij = 1 ⇐⇒ bi = bj.

For all tuples (bi; bj , bk) ∈ B3 , if we play pair (bi, bj) once and observe bi ≻ bj, we define κ(bk; bi, bj) to be the update of wins of bk and τ(bk; bi, bj) to be the update of plays of bk. κ(·; ·, ·) and τ(·; ·, ·) represent the dependent structure of the tuple arms. They could be functions of rij, rik, and rjk.

In our synthetic experiments, we assume the input space (set of arms) B has dependent structure and there exists an underlying utility function f(b) : B → R over input space which we cannot observe directly. Our observations are the noisy comparisons between pairs of arms (e.g., bi and bj) which can be viewed as the noisy comparison of utility values (e.g., f(bi) and f(bj)). The properties of strong stochastic transitivity, stochastic triangle inequality, and the dependency assumptions on rij generally hold for a wide range of applications. In the clinical experiments, we extract comparisons from physician's online judgment.

4 Algorithm

Our algorithm, CORRDUEL as shown in Algorithm 1, is a correlational dueling bandits algorithm based on the Beat-theMean algorithm [Yue and Joachims, 2011]. It uses observational feedback and the correlational structure to successively remove suboptimal arms, while keeping the optimal one(s) in the sample space with high probability. The inputs to CORRDUEL are the set of arms B, the total number of iterations T , and the correlational structure (κ, τ).

Parameters-Initialization (Algorithm 2) defines the set of active arms Wℓ, whose size shrinks as more tests are completed. For each arm b, let nb be the total number of comparisons between b and other arms, and let wb be the total number of wins against all other arms. Let ˆ Pb be the empirical average of P(b ≻ b′ ) for all b′ in Wℓ, and let ˆ Pb,n be the value of ˆ Pb after n comparisons between arm b and any other arms. Set the confidence interval of P(b ≻ b′ ) as:

where cδ(n) = √(1/n)log(1/δ), and δ is the confidence that P(b ≻ b′ ) lies in ˆ Cb,n. The function cδ(n) decreases as the number of comparisons n increases. By properly setting parameter δ, the optimal reward can be reached within the fixed time horizon as shown in Proposition 1.

Active-Elimination (Algorithm 3) is the key part of CORRDUEL. For each pair of tests, two arms are randomly chosen fromWℓ. The randomized selection method enjoys lowvariance total regret in general. For each arm b, the values of wb, nb and ˆ Pb are updated, as is the corresponding confidence radius c∗ . An arm b dominates another arm b′ , if their confidence intervals do not overlap, and the inferior arm is eliminated fromWℓ. The algorithm runs until the time horizon T is reached, or only one active arm remains.

CORRUPDATE (Algorithm 4) is the subroutine of ActiveElimination (Algorithm 3) which updates the weights of bk by rules κ(·; ·, ·) and τ(·; ·, ·). In the classical dueling bandits setting, we assume arms are independent. For independent arms, if we have one comparison between bi and bj and gets bi ≻ bj , we only update the weights for arm bi and bj :

For a large decision space, existing dueling bandits algorithms are extremely slow if one does not exploit dependencies among arms, even if they can achieve provably optimal cumulative regret (w.r.t. independent arms). When the arms are correlated and the correlation between any pair of arms bi and bj is measured properly by rij, we can update all active arms at each iteration.

As shown in Algorithm 4, we update every arm bk after comparing arms bi and bj (w.l.o.g. assume bi ≻ bj) via:

Algorithm 1 CORRDUEL

• 1: Input: B, T, (κ, τ)
• 2: Input: cδ(n) = √(1/n)log(1/δ)

3: Run: [Parameters-Initialization]

• 4: Run: [Active-Elimination]
• 5: return b∗ // Optimal arm

Algorithm 2 Parameters-Initialization

• 1: W1 ←B // set of active arms
• 2: ℓ ← 1 // rounds
• 3: ∀b ∈Wℓ, nb ← 0 // comparisons
• 4: ∀b ∈Wℓ, wb ← 0 // priorities
• 5: ∀b ∈Wℓ, ˆ Pb ≡ wb/nb, or 1/2 if nb = 0
• 6: n∗ ≡minb∈Wℓ nb
• 7: c∗ ≡ cδ(n∗ ), or 1 if n∗ = 0 // confidence radius
• 8: t ← 0 // total number of iterations
• 9: return all new parameters

where κ(·; ·, ·) and τ(·; ·, ·) represent the correlational structure, which is assumed to satisfy:

• 0 ≤ κ(bk; bi, bj) ≤ τ(bk; bi, bj) ≤ 1;
• if bk = bi, κ(bk; bi, bj) = τ(bk; bi, bj) = 1;
• if bk = bj, κ(bk; bi, bj) = 0, τ(bk; bi, bj) = 1.

These updates are based on the assumption that κ(·; ·, ·) τ(·; ·, ·) is an unbiased estimation of the dependent structure. The CORRUPDATE subroutine (Algorithm 4) can efficiently update all arms at each iteration. So CORRDUEL enjoys fast convergence towards the near optimal arms.

Definition 1. ε-optimal arm. If arm b satisfies ϵ(b1, b) ≤ ε, then b is an ε-optimal arm.

Proposition 1. If ∃µ > 0 such that τ(bk; bi, bj) ≥ µ for every (bi, bj, bk) ∈ B3 , then with high probability, the cumulative time to achieve purely ε-optimal arms T(ε) is bounded by:

Proof. Proposition 1 holds based on the Theorem 1 of [Yue and Joachims, 2011]. After t iterations, since τ(bk; bi, bj) ≥ µ, we have n∗ ≥ µt. Then c∗ = cδ(n∗ ) = √(1/n∗ )log(1/δ) ≤ √(1/µt)log(1/δ). Notice c∗ is a function of time step t.

For any arm b which is not ε-optimal (satisfies ϵ(b1, b) > ε), with probability 1 - δ, ˆ Pb1 - ˆ Pb > εCδ holds for some fixed concentration parameter Cδ. Suppose arm b has not been eliminated at iteration t. Then from elimination criterion Line 7 of Algorithm 3 we have εCδ < ˆ Pb1 - ˆ Pb < 2c∗ ≤ 2√(1/µt)log(1/δ). The inequality breaks when t ≥ 4 µε2C2 δ log 1δ = O ( 1 µε2 log 1δ ).

Notice, the iteration time T(ε) in Propositions 1 does not depend on |B| = K, which suggests the fast convergence of CORRDUEL in large decision spaces.

Algorithm 3 Active-Elimination

• 1: while |Wℓ| > 1 and t ≤ T do
• 2: select bi, bj ∈Wℓ at random
• 3: compare selected arms (assume bi ≻ bj)
• 4: for all bk ∈ Wℓ do
• 5: update wk, nk by CORRUPDATE
• 6: end for
• 7: if minb′∈Wℓ ˆ Pb′ + c∗ ≤maxb∈Wℓ ˆ Pb - c∗ then
• 8: b′ ← argminb∈Wℓ ˆ Pb
• 9: ∀b ∈ Wℓ, delete comparisons with b′ from wb, nb
• 10: Wℓ+1 ←Wℓ\{b′ } // update working set
• 11: ℓ ← ℓ +1 // new round
• 12: end if
• 13: end while
• 14: return b∗ = argmaxb∈Wℓ ˆ Pb

Algorithm 4 CORRUPDATE

• 1: Input: bk, bi ≻ bj
• 2: wk ←wk +κ(bk; bi, bj)
• 3: nk ←nk + τ(bk; bi, bj)
• 4: return wk, nk

In our application of CORRDUEL to selection of optimal multi-electrode stimulating parameters for paraplegic, we define the similarity of different configurations to be the correlation coefficient of electrical potential fields generated by the two different electrode stimulating configurations. Since the correlation coefficient function r(·, ·) has support on [-1, 1], we only update with the CORRUPDATE rule when r(·, ·) ≥ 0. The existence of negative r values is based on clinical observations. The correlational property arises from analysis of electric fields applied by the array as shown in Figure ??.

The standard notion of correlation coefficient, rXY = E[XY -E[X]E[Y ]]/√V ar[X]V ar[Y ], is used in our experiments. However, one can use any measure as a basis for rXY as long as rXY ∈ [0, 1], rXY = 1 when X = Y , and rXY = 0 when X has an "irrelevant" relation to Y . The coefficient r can take negative values, but the algorithm doesn't use negative values for its updates.

For correlated arms, we perform an update for every arm k for which rik, rjk > 0 as follows:

Proposition 2. If ∃µ > 0 such that rij ≥ µ for every pair (bi, bj) ∈ B2 , then with high probability, the cumulative time to achieve purely ε-optimal arms T (ε) satisfies:

Proof. If rij ≥ µ for every pair (bi, bj) ∈ B2 , since rij ≤ 1, τ(bk; bi, bj) = rik+rjk 1+rij ≥ 2µ 2 = µ for every tuple (bi, bj, bk). The result follows from substituting it into Proposition 1.

The CORRUPDATE subroutine above updates the dueling pair bi, bj in the same way as if they are independent since (5) and (6) will collapse to (1) and (2) for bi and bj. For extreme cases, if bi ≻ bj and arm bk is very close to bj. We have rik ≃ 1 and rjk ≃ rij, the updating rules for arm bk will be close to the updates of arm bj. If bk is far from both bi and bj, (5) and (6) guarantees that the update for bk is very small since we acquire little information about bk from far away comparisons. Also, if bi and bj are less dependent (with smaller rij), we would expect to acquire larger updates for the points in between.

One can also consider a Bayesian version, e.g., by using Gaussian processes. In this paper, we focus on a frequentist approach, which is a better model of the clinical application.

5 Experiments

We evaluated our approach in two settings, synthetic simulations and a real clinical application of online optimization for spinal cord stimulation therapy. In our controlled synthetic experiments, we seek to address the following questions:

• How does the algorithm compare against standard dueling bandit algorithms?
• How effective is it in terms of convergence?

We compare the algorithm against Beat-the-Mean, RUCB, and Sparring algorithm with UCB1. These three algorithms are the representative dueling bandits algorithms designed for independent arms, which do not, however, leverage the correlations between arms.

5.1 Simulation Experiments

Setup. We first evaluate the algorithm with simulation experiments. The purpose of this experiment is to validate our algorithm, and demonstrate its quick convergence when the arms are dependent. To generate the underlying utility function over correlated arms, we sampled random functions from a zero-mean Gaussian Process with squared exponential kernel over the sample space B = [0, 1] × [0, 1], uniformly discretized into 50 × 50 points (set of arms) and used this function as the mean function for the 2500 arms. We chose σ = 0.5 as the standard deviation of arms. One evaluation of the mean functions is shown in Figure 3. The utility function is not necessarily convex or simple. Within each iteration, we sample 2 points in the active set and compare their (noisy) sampling values to get the {0, 1} feedback of the duel. We run the duel for T = 100 iterations for 10000 trials for each of the 4 comparing algorithms.

Figure 4: Regret versus iteration. The dashed lines represent one standard deviation.

Results. We report a notion of regret as the stepwise regret instead of the cumulative regret. It converges to zero as iteration number goes to infinity for every no-regret algorithm. As seen in Figure 4, CORRDUEL converges much faster than the other three algorithms since it takes the advantage of the dependent arms. The independent-armed dueling bandits algorithms require an exhaustive searching period which is significantly larger than the time horizon we use here before concentrating on the (near) optimal arms.

5.2 Human Experiments

Background. As depicted in Figure 1 from before, our human clinical experiments involve optimizing a system for stand training under spinal cord stimulation with spinal cord injury patients. The subject practices standing under spinal stimulation using a stand frame for assistant in balance. The training processes largely follow the procedures in [Rejc et al., 2015]. Two trainers on the subject's left and right protect and assist the subject. Within each experiment, a specific stimulating pattern (a combination of active electrode selections, the polarity of the actively selected electrodes, and the stimulation amplitude and frequency) is applied through the implanted electrode array and its controlling circuitry. An anonymous short video 1 shows the standing quality under different stimuli. The first part shows a low quality bipedal standing and the second part shows a better standing, both with electrical spinal cord stimulation. Different standings could look similar for the non-specialist.

The participants are under stable medical condition and have no musculoskeletal dysfunction that might interfere with stand training. They have no motor response present in leg muscles during transcranial magnetic stimulation, indicating that there are no strongly active neural pathways connecting cortex and lower limb muscles. No volitional control can be achieved during voluntary movement attempts in leg muscles as measured by EMG activity.

Setup. We use clinical knowledge to restrict the decision space from around 4.3×107 to be on the order of 103 ∼ 104 . It is still a very large decision space considering the number of trials, or arm pulls, are on the order of 102 .

A total of 414 experimental comparisons were done with two patients under the CORRDUEL algorithm. Each trial lasted for about 5 minutes. Within each trial, one stimulating pattern was generated by the 16-channel electrode. The patterns were unchanged within each trial. For a fixed electrode configuration, the stimulation frequency and amplitude were modulated synergistically in order to find the best values for effective weight-bearing standing. We optimized the electrode patterns with CORRDUEL and performed exhaustive search for stimulation frequency and amplitude over a narrow range.

Stimulation began while the patient was seated. Then the participant initiated the sit to stand transition by positioning his feet shoulder width apart and shifting his weight forward to begin loading the legs.

Results. For the clinical experiments, we cannot create a direct plot for regrets since the ground truth optimal stimulation is unknown. In the experiments, we observed the convergence of CORRDUEL, which is not possible for independentarmed dueling bandits algorithms. The set of (near) optimal configurations found by CORRDUEL is shown in Figure 5. We compared the performance of CORRDUEL to the optimal selections found heuristically for each patient by clinicians, which are shown in Figure 6. We found that the manual selection is a subset of the algorithm's selection, and there exist high performing configurations (e.g., the 2nd in Figure 5) found by the algorithm which are not in the manual selection. This shows that CORRDUEL is performing no worse than specialized physicians.

6 Conclusion and Discussion

Our analysis and simulation demonstrate that CORRDUEL indeed exhibits fast convergence properties compared to independent-armed dueling bandits algorithms when correlation information is available. We deployed this algorithm in clinical experiments for the control of spinal cord stimulation and showed that CORRDUEL performs no worse than specialized physicians. We believe that our result provides an important step towards employing machine learning algorithms in many problems with a large volume of parameter selection and sequential decision making. These problems could be facilitated by our algorithm, which simultaneously delivers effective decisions and explores the decision space based on comparative feedback.

1 https://youtu.be/loJLtbcUBDM

The CORRUPDATE subroutine is easy to incorporate with Beat-the-Mean algorithm to achieve efficient CORRDUEL. Although we developed CORRDUEL specifically based on Beat-the-Mean, CORRUPDATE is a more general approach which has potential to incorporate with the existing dueling bandits algorithms. For instance, it can incorporate with RUCB to realize a variant of RUCB for dependent arms by updating the wins wij with CORRUPDATE.

To our knowledge, our work is the first to apply an algorithmic approach towards spinal cord injury treatments. The algorithm could find a proper set of optimal stimulating configurations within the test time horizon. We achieved good performance in both simulations and human experiments. The paraplegic human patients could achieve full-weight standing under the stimulation provided by our algorithm.

Acknowledgements

Yisong Yue (yyue@caltech.edu) is a coauthor of this paper. This research was supported in part by Caltech/JPL PDF IAMS100224, NIH-U01-EB007615-08, NIHU01-EB015521-05, and a gift from Northrop Grumman.