Budgeted Batch Bayesian Optimization With Unknown Batch Sizes

Parameter settings profoundly impact the performance of machine learning algorithms and laboratory experiments. The classical grid search or trial-error methods are exponentially expensive in large parameter spaces, and Bayesian optimization (BO) offers an elegant alternative for global optimization of black box functions. In situations where the black box function can be evaluated at multiple points simultaneously, batch Bayesian optimization is used. Current batch BO approaches are restrictive in that they fix the number of evaluations per batch, and this can be wasteful when the number of specified evaluations is larger than the number of real maxima in the underlying acquisition function. We present the Budgeted Batch Bayesian Optimization (B3O) for hyper-parameter tuning and experimental design - we identify the appropriate batch size for each iteration in an elegant way. To set the batch size flexible, we use the infinite Gaussian mixture model (IGMM) for automatically identifying the number of peaks in the underlying acquisition functions. We solve the intractability of estimating the IGMM directly from the acquisition function by formulating the batch generalized slice sampling to efficiently draw samples from the acquisition function. We perform extensive experiments for both synthetic functions and two real world applications - machine learning hyper-parameter tuning and experimental design for alloy hardening. We show empirically that the proposed B3O outperforms the existing fixed batch BO approaches in finding the optimum whilst requiring a fewer number of evaluations, thus saving cost and time.

Introduction
Global optimization is fundamental to diverse real-world problems where parameter settings and design choices are pivotal -as an example, in algorithm performance (deep learning networks (Bengio, 2009)) or quality of the products (chemical processes or engineering design (Wang & Shan, 2007)). This requires us to find the global maximum of a non-concave objective function using sequential, and often, noisy observations. Critically, the objective functions are unknown and expensive to evaluate. Therefore, the challenge is to find the maximum of such expensive objective functions in few sequential queries, thus minimizing time and cost.
Bayesian optimization (BO) is an approach to find the global optimum of such expensive, blackbox objective functions f using limited evaluations (Snoek, Larochelle, & Adams, 2012;Shahriari, Swersky, Wang, Adams, & de Freitas, 2016;Wang, Hutter, Zoghi, Matheson, & de Feitas, 2016). Instead of optimizing the real function, BO utilizes a cheaper to evaluate surrogate function, called the acquisition function. This function is used to suggest the next point by balancing exploitation (knowledge of what has been observed) and exploration (where a function has not been investigated). Then, after evaluating the objective function at the suggested point, a Gaussian process (Rasc 2017 AI Access Foundation. All rights reserved. mussen, 2006) is updated and thus gradually a profile of the mean and uncertainty of the exploration space is built up. Bayesian optimization has been demonstrated to outperform other state-of-the-art black-box optimization techniques when function evaluations are expensive and the number of allowed function evaluations is low . Thus, BO has received increasing attention in machine learning community (Thornton, Hutter, Hoos, & Leyton-Brown, 2013;Khajah, Roads, Lindsey, Liu, & Mozer, 2016).
Most selection strategies in BO are sequential, wherein only one experiment is tested at a timethe experiment selection at each iteration is optimized by using the available observed information. However, such methods are inefficient when parallel evaluations are possible. For examples, many alloy samples can be placed in an oven simultaneously (for testing the alloy quality), or several machine learning algorithms can be run in parallel at the same time using multiple cores. This motivates batch Bayesian optimization algorithms that select multiple experiments, a batch of q experiments, at each iteration for evaluation.
Existing batch BO approaches are mostly greedy, sequentially visiting all the maxima of the acquisition function. After the first maximum is found, such methods modify the acquisition function by suppressing the current maximum and then move on to find the next best maximum. This repeats till a batch of maxima is collected. Different algorithms implement this philosophy - (Ginsbourger, Le Riche, & Carraro, 2008) modifies the found maximum point by replacing it with a "fake" or constant value, thus the algorithm is called "constant liar". Other approaches, including GP-BUCB (Desautels, Krause, & Burdick, 2014) and the GP-UCB-PE (Contal, Buffoni, Robicquet, & Vayatis, 2013), exploit the predictive variance of GPs that only depend on the features x, but not the outcome values y. Local Penalization (González, Dai, Hennig, & Lawrence, 2016), on the other hand, using the estimation of Lipschitz constant to penalize the peaks. All these methods update the posterior variance sequentially to modify the acquisition function and thereby derive a batch. Most importantly, these batch BO methods are essentially greedy, choosing individual points until the batch is filled. This is not optimal as after the "real" maxima in the function are found, noisy points will get added simply to complete the batch. An alternate and non-greedy approach is to select a batch of points that reduce the uncertainty of the global maximizer (Shah & Ghahramani, 2015) by extending the Predictive Entropy Search (PES) algorithm (Hernández-Lobato, Hoffman, & Ghahramani, 2014) for a batch setting. But again, a fixed batch size is used. Fixed batch size approaches can restrict the flexibility in finding the global solution and require unnecessary evaluations or miss evaluation of important points. If we over-specify the batch size, we waste resources in evaluating redundant points. In contrast, if we under-specify the batch size, we may miss important points that could potentially be the optimal solution of f .
The best solution is to fill each batch flexibly, but the challenge is to know exactly how many maxima are there in each batch. Our key contribution is to solve this problem, by estimating the number of maxima in the acquisition function, we let the algorithm adjust the batch size at each iteration. This flexible and non-greedy approach has the advantage of being efficient in that it suggests only as many maxima locations as needed. The solution comes from our intuition that the multiple peaks in the acquisition function can be approximated by a mixture of Gaussians -the means of the Gaussian correspond to the underlying peaks. However, since the number of underlying peaks is unknown, we use the infinite Gaussian mixture model (IGMM) (Rasmussen, 1999) so that the number of peaks in the acquisition function can be estimated automatically. Because fitting the IGMM directly to the acquisition function is intractable, we present an efficient batch generalized slice sampler approach to draw samples from acquisition function which are then used to learn the unknown number of peaks. Although the IGMM and the slice sampler have existed for more than a decade, the idea of making use of these techniques for batch Bayesian optimization is novel.
In the experiments, we first validate our methods in 8 benchmark functions and compare it against 9 baselines. We use several criteria for comparison: a) in best-found-value, our method achieves better results in 5 out of 8 cases; and b) our method matches the baseline performance but requires far fewer evaluation points. We further demonstrate the algorithm on machine learning hyper-parameter tuning for support vector regression, multi-label classification and deep learning. We show that our model outperforms the baselines in finding the best values (RMSE, F1 score and accuracy) whilst we require fewer experimental evaluations to reach targets. In addition, we perform a real world experimental design in which we formulate the heat-treatment process required for producing an Aluminum-scandium alloy. We show that our method produces the highest hardness using the fewest number of evaluations. Our result is significant as searching for the best heattreatment process is costly and making this process efficient results in cost and time savings.

Preliminary
We first review the Gaussian process (GP), Bayesian optimization (BO) and the acquisition functions. Then, we summarize the batch Bayesian optimization setting and existing approaches.

Gaussian process
Gaussian processes (GP) (Rasmussen, 2006) extends a multivariate Gaussian distribution to infinite dimensionality. Formally, Gaussian process generates data located throughout some domains such that any finite subset of the range follows a multivariate Gaussian distribution. Given N observations Y = {y 1 , y 2 , ...y N } which can always be imagined as a single point sampled from some multivariate Gaussian distributions.
The mean of GP is assumed to be zero everywhere. What relates one observation to another in such cases is just the covariance function, k (x, x ′ ). From the assumption of GP, we have y ∼ N (0, K) where the covariance matrix is defined as follows: A popular choice for the covariance function k is the squared exponential function: The length parameter l will control this separation when x is not closed to x ′ .
For prediction on a new data point y * , we can update the covariance matrix with k * * = k (x * , x * ) and k * = k (x * , x 1 ) k (x * , x 2 ) · · · k (x * , x N ) . Hence, we can write The conditional probability is followed Gaussian distribution as p (y * | y) ∼ N µ (x * ) , σ 2 (x * ) (Rasmussen, 2006;Ebden, 2008) where its mean and variance are given by GPs provide a full probabilistic model of the data, and allow us to compute not only the model's prediction at input points but also to quantify the uncertainty in the predictions. Therefore, GP is flexible as a nonparametric prior for Bayesian optimization.

Bayesian optimization
We assume that f is a black-box function, that is, its form is unknown and further it is expensive to evaluate. Perturbed evaluations of the type y i = f (x i ) + ε i , where ε i ∼ N 0, σ 2 , are available. Example of the Branin function is in Fig. 1. Bayesian optimization makes a series of evaluations x 1 , ..., x T of f such that the maximum of f is found in the fewest iterations (Snoek et al., 2012;Shahriari et al., 2016;Dai Nguyen, Gupta, Rana, Nguyen, Venkatesh, Deane, & Sanders, 2016;Nguyen, Rana, Gupta, Li, & Venkatesh, 2016c). Formally, let f : X → R be a well behaved function defined on a compact subset X ⊆ R D . Our goal is solving the global optimization problem Bayesian optimization reasons about f by building a Gaussian process through evaluations (Rasmussen, 2006). This flexible distribution allows us to associate a normally distributed random variable at every point in the continuous input space.

ACQUISITION FUNCTIONS
As the original function is expensive to evaluate, the acquisition function acts as a surrogate that determines which point should be selected next. Therefore, instead of maximizing the original function f , we maximize the acquisition function α to select the next point to evaluate In this auxiliary maximization problem, the objective is known and easy to evaluate and can be easily carried out with standard numerical techniques such as multi-start or DIRECT (Jones, Perttunen, & Stuckman, 1993). We consider the case of computing the acquisition function α (x) from the posterior distribution of GP. These acquisition functions are carefully designed to trade off exploration of the search space and exploitation of current promising regions. Among many existing acquisition functions in literature (Hennig & Schuler, 2012;Hernández-Lobato et al., 2014;Srinivas, Krause, Kakade, & Seeger, 2010;Mockus, Tiesis, & Zilinskas, 1978;Jones, 2001;Freitas, Zoghi, & Smola, 2012;Nguyen, Gupta, Rana, Li, & Venkatesh, 2016b), we briefly describe three common acquisition functions including probability of improvement, expected improvement and upper confidence bound. The early work of (Kushner, 1964) suggested maximizing the probability of improvement (PI) over the incumbent α PI ( , where Φ is the standard normal cumulative distribution function (cdf) and the incumbent y + = max x i ∈D t f (x i ). However, the PI exploits quite aggressively. Thus, one could instead measure the expected improvement (EI) (Mockus et al., 1978;Jones, Schonlau, & Welch, 1998). The expected improvement chooses the next point with highest expected improvement over the current best outcome, i.e. the maxi- where φ is the standard normal probability distribution function (pdf), τ is the current best value, u (x) = µ(x)−τ σ (x) , for σ > 0 and zero otherwise. The Gaussian process upper confidence bound (GP-UCB) acquisition function is given by α GP-UCB (x) = µ(x) + β σ (x), where β is a domain-specific positive parameter which trades exploration and exploitation. There are theoretically motivated guidelines for setting and scheduling the hyper-parameter β to achieve sublinear regret (Srinivas et al., 2010).

Batch Bayesian optimization
Bayesian optimization is conventionally posed as a sequential problem where each experiment is completed before taking a new one. In practice it may be advantageous to run multiple function evaluations in parallel. Therefore, we consider evaluating f using a batch of points. An example of batch Bayesian optimization versus sequential Bayesian optimization is illustrated in Fig. 2. As discussed earlier, such scenarios appear, for instance, in the optimization of computer models where several machines (or cores) are available to run experiments in parallel, or in wet-lab experiments when the time of testing one experimental design is the same as testing a batch. Formally, we maximize the acquisition function by finding a collection of points as where n t is the batch size at iteration t. Existing approaches often fix the batch size n t to a constant value for all iterations.

Existing batch Bayesian optimization methods
Since the optimization in Eq. 2 is intractable, batch BO techniques avoid this computational burden by resorting to different strategies. The first BO parallelization was used in the context of the learning Bayesian networks (Očenášek & Schwarz, 2000), to the best of our knowledge. Then batch BO simulation matching (Azimi, Fern, & Fern, 2010) aims to select a batch of size q, which includes points 'close to' the best point. Developing from the expected improvement acquisition function (EI), the multi-point expected improvement (q-EI) (Ginsbourger, Le Riche, & Carraro, 2007) treats the acquisition function as the conditional expectation of the improvement obtained by q points. Constant Liar (CL) is a heuristic q-EI algorithm (Ginsbourger et al., 2008), which uses a greedy approach to iteratively construct a batch of q points. At each iteration, the heuristic uses the sequential algorithm to find a point that maximizes the expected improvement as follows: first, the maximum of the acquisition is found and to move to the next maximum by suppressing this point. This is done by inserting the outcome at this point as a constant value. The GP posterior is updated through the "fake" outcome which then is used to compute the acquisition function. This process is repeated until the batch is filled. Further parallel extension for q-EI are proposed in (Frazier & Clark, 2012;Wang, Clark, Liu, & Frazier, 2015).
Another direction (Contal et al., 2013;Desautels et al., 2014) in batch Bayesian optimization exploits an interesting fact about GPs: the predictive variance of GPs depends only on the feature x, but not the outcome values y. The GP-BUCB algorithm (Desautels et al., 2014) and GP-UCB-PE (Contal et al., 2013) extend the sequential UCB to a batch setting by first selecting the next point, updating the predictive variance which in turn alters the acquisition function, and then selecting the next point. This is repeated till the batch is filled. In particular, the GP-UCB-PE (Contal et al., 2013) chooses the first point of the batch via the UCB score and then defines a "relevance region" and selects the remaining points from this region greedily to maximize the information gain, to focus on pure exploration (PE).
Batch BO can also be developed using information-based policies (Hernández-Lobato et al., 2014). PES aim to select the point which maximizes the information gain of the model. This is solved by expressing the expected reduction in the differential entropy of the predicted distribution. Parallel Predictive Entropy Search (PPES) (Shah & Ghahramani, 2015) extends the Predictive Entropy Search (PES) algorithm of (Hernández-Lobato et al., 2014) to a batch setting.
More recently, Local Penalization (LP) (González et al., 2016) presents a heuristic approach for batch BO by iteratively penalizing the current peak in the acquisition function to find the next peak. LP depends on the estimation of Lipschitz constant to flexibly penalize the peaks. However, in general scenarios, the Lipschitz constant L is unknown. For ease of computation, (González et al., 2016) estimates the Lipschitz constant of the GP predictive mean instead of the actual acquisition function. Furthermore, the use of a unique value of L assumes that the function is Lipschitz homoscedastic. Gonzalez et al (González et al., 2016) has demonstrated that LP outperforms a wide range of baselines in batch Bayesian optimization. We thus consider LP as the most competitive baseline in batch BO.
Although the batch algorithms can speedup the selection of experiments, there are two main drawbacks. First, most of the proposed batch BO (Ginsbourger et al., 2008;Azimi et al., 2010;Azimi, Jalali, & Zhang-fern, 2012;Contal et al., 2013;Desautels et al., 2014;González et al., 2016) involve a greedy algorithm, which chooses individual points until the batch is filled. This is often detrimental as the requirement to fill a batch enforces the selection of not only the true maxima but also noisy peaks. Second, most approaches use a predefined and fixed batch size q for all iterations. Fixed batch size may be inefficient because the number of peaks in the acquisition function is unknown and constantly changing when the GP posterior is updated. For example, if a function has two real maxima, a batch size of more than two will force the selection of noisy points. Therefore, time and resources are wasted for evaluating these noisy points that contravening the goal of BO is to save the number of evaluations. In contrast, if we under-specify the number of peaks, we could miss important points that may be the optimal solution of f . Hence, it is necessary to identify the appropriate, if not exact, number of peaks for evaluation while preserving the ability of finding the optimal value of f .

Budgeted Batch Bayesian Optimization
We propose the novel batch Bayesian optimization method that learns the suitable batch size for each iteration. We term our approach budgeted batch Bayesian optimization (B3O). The proposed method is economic in terms of number of optimal evaluations whilst preserving the performance.
We first describe our approach to approximate the acquisition function using the infinite Gaussian mixture model (IGMM). Then, we present the batch generalized slice sampler to efficiently draw samples under the acquisition function. Next, we briefly describe the existing variational inference technique for IGMM. Finally, we summarize our algorithm.

Batch Bayesian optimization
We consider the batch Bayesian optimization. As a Bayesian optimization task, our ultimate goal is solving the global optimization problem of finding x * = argmax f (x) x∈X by making a series of batch such that the maximum of f is found (Snoek et al., 2012;Shahriari et al., 2016) where n t is the batch size.

Approximating the acquisition function with infinite Gaussian mixture models
The acquisition function α is often multi-modal with the unknown number of peaks. It is intuitive to assume that the acquisition function can be approximated by a mixture of Gaussians (cf. Fig. 3) wherein the peaks in the acquisition function are equivalent to the mean locations of the underlying Gaussians. Because the number of peaks are unknown, we borrow the elegance of Bayesian nonparametrics (Hjort, Holmes, Müller, & Walker, 2010) to identify the unknown number of Gaussian components. We use the infinite Gaussian mixture model (IGMM) (Rasmussen, 1999) which induces the Dirichlet Process prior over possibly infinite number of Gaussian components. In IGMM, each Gaussian component representing a peak k is parameterized by the mean µ k and the covariance Σ k . It is interesting that the estimated mean µ k is the location of the peaks while the estimated covariance Σ k will capture the shape of the peaks in the acquisition function. In contrast, previous work (González et al., 2016) relies on an unique Lipschitz constant to represent the shape of the peak that is not reasonable for the heteroscedastic setting when the peaks have different shapes. Formally, the probability density function of the IGMM is defined as follows ∑ ∞ k=1 π k × N (s | µ k , Σ k ) where s is an observation, π k is the mixing proportion in IGMM and N (s | µ k , Σ k ) is the Gaussian density function of the component k.
Our primary goal in utilizing IGMM is to approximately find the mean µ k (s) of Gaussian distributions as the unknown peaks from the acquisition function. However, directly estimating IGMM from acquisition function is intractable. Therefore, we use the intermediate step to draw samples under the acquisition function and then fit the IGMM to learn the mean locations as the batch of points to evaluate.
Input: #MaxIter, acquisition function α, x min , x max 1: s 0 ∼ uniform (x min , x max ) 2: α min = min x∈X α (x) 3: for i = 1 to #MaxIter do 4: while (notAccept) do 6: if α(s i ) > u i #Accept then We now present the sampling technique to draw samples under the acquisition function for fitting to the IGMM. We note that the sampling process, by its nature, generates more samples from high probability regions. This implies that most of the samples come from the peaks. Thus, even with small number of samples, the sampling process approximates the peaks and is then efficient. There are many existing sampling methods to draw samples from the probability distributions, for example, Metropolis-Hastings (Metropolis, Rosenbluth, Rosenbluth, Teller, & Teller, 1953). However, to sample under the D-dimensional acquisition function which is not strictly probability distribution, we utilize the accept-reject sampling (Casella, Robert, & Wells, 2004) where we keep samples in the region under the density function and ignore samples if they are outside the curve.
In particular, we select to use the slice sampling (Neal, 2003) because it is easily implemented for univariate distributions, and can be used to sample from a multivariate distribution by updating each variable in turn. Suppose we wish to sample from a distribution for a variable, x, taking values in some subset of R D , whose density is proportional to some functions g(x). We can do this by sampling uniformly from the (D + 1)-dimensional region that lies under the plot of g(x). This idea can be formalized by introducing an auxiliary real variable, u, and defining a joint distribution over x and u that is uniform over the region R = {(x, u) : 0 < u < g(x)} below the curve or surface defined by g(x). The joint density for (x, u) is Z . To sample for x, we can sample jointly for (x, u), and then ignore u.
We extend the standard slice sampler to the generalized slice sampler that draws samples from D-dimensional acquisition function α which is not a proper distribution and can be negative. To draw u uniformly over the region R below curve of the acquisition function, we get the α min = min x∈X α (x) obtained by one of the non-convex optimization toolbox (e.g., DIRECT (Jones et al., . Accept Figure 4: The illustration of the proposed generalized slice sampling to draw samples under the acquisition functions. We consider an example as follows. The first point s 1 is uniformly drawn from Uniform (x min , x max ), then u 1 ∼ Uniform (α min , α(s 1 )). Next, the second point is sampled s 2 ∼ Uniform (x min , x max ). Due to α(s 2 ) < u 1 , we reject s 2 . We continue to draw s 3 and accept it as α(s 3 ) ≥ u 1 . We repeat the algorithm by again sampling u 3 ∼ Uniform (α min , α(s 3 )).

1993)
). Therefore, we can define the joint density To sample for x, we can sample jointly for (x, u), and then ignore u following the Algorithm 1 and Fig. 4 where uniform (x min , x max ) denotes for drawing randomly from the uniform distribution within the bound of [x min , x max ], which is predefined.
Slice sampling approach (Neal, 2003) also brings a challenge for high-dimensional data that it is hard to find accepted area R such that g (R) > u. There are two main reasons. First, the set of accepted samples could be very small with respect to the hyper-rectangle R. Second, R usually lies in a high-dimensional space, implying the curse of dimensionality. Therefore, it is a big drawback of slice sampling to deal with high dimensions (Tibbits, Haran, & Liechty, 2011;Pietrabissa & Rusconi, 2014).
To overcome the problem in high-dimensional functions, we utilize the generalized slice sampler in batch setting where a bunch of samples are drawn i.i.d at different places. Since our goal is to draw a collection of i.i.d. samples x 1,2,...N from the acquisition function α (to later estimate the IGMM). We define a collection of joint density distribution p (x 1 , u 1 ) , ...p (x M , u) and perform the generalized slice sampling for each joint distribution independently. In the experiment, we set M as 200. In Sec 4.5, we further study the computational efficiency of the batch generalized slice sampling (BGSS) for drawing samples from high dimensional acquisition functions (e.g., D = 50).

Variational inference for infinite Gaussian mixture model
IGMM is the nonparametric mixture model where the prior distribution over the mixing proportion is a Dirichlet process (Ferguson, 1973). There are few existing approaches to learn a IGMM, such as collapsed Gibbs sampler (Kamper, 2013) and variational inference (Blei, Jordan, et al., 2006). In this paper, we follow (Blei et al., 2006) to derive the variational inference for IGMM (Rasmussen, 1999) since the variational approach is generally faster than the Gibbs sampler. After fitting the IGMM using variational inference, we obtain the mean locations µ 1....K as the selected points for evaluations in the batch BO setting. We note that K is unknown and identified automatically in this Bayesian nonparametric setting.
Using the collection of samples (s i ) N i=1 drawn from our batch generalized slice sampler, we consider the IGMM generated as follows. We first draw a mixing proportion from the stick-breaking process given the concentration parameter γ as π where v k ∼ Beta (1, γ). Next, for each component k in the model, we sample the mean µ k ∼ N (µ 0 , I) and the covariance matrix Σ k ∼ Wish (τ k , I). Finally, we generate the assignment z i ∼ Stick (π) and the observation For posterior inference using variational inference, we first write the bound on the likelihood of the data as follows: To define the variational distribution q in the Eq.
(3), we need to construct a distribution on an infinite set of random variables v k ,µ k and Σ k . For this approach be tractable, we truncate the variational distribution at some value K by setting q(v K = 1) = 1 and we can ignore µ k , Σ k for k > K (Blei et al., 2006). Then, we factorize the variational distribution using mean field assumption as the following Particularly, the variational distributions for these variables are defined as q After factorizing the variational distribution using mean field, we maximize the lower bound in Eq. (3). We first take the partial derivative w.r.t. each variational variables. We then equate it to zero and solve the optimization for each variables. Due to the space restriction, we shall refer the detailed inference of the infinite Gaussian mixture model to (Rasmussen, 1999;Blei et al., 2006).

Algorithm
Although the technique of variational inference for IGMM and the slice sampling are existed for years, the idea of connecting these things at the right place for batch BO is novel. We summarize the steps for B3O in Algorithm 2. At an iteration t, we will find a batch including n t points where the number of point n t varies at each iteration.
We note that steps 2,3,7 and 8 are standard in Bayesian optimization techniques. Our proposed method is highlighted in steps 4,5 and 6 to find a batch of points X t from the acquisition function. Particularly, Fig. 3 illustrates and summarizes our steps 4,5 and 6 in 1D and 2D, respectively. We summarize the computation time of steps 4, 5, and 6 in Sec 4.4 and further analyze step 4 by simulation in Sec 4.5 which requires more computation than other steps.

Algorithm 2 Algorithm for Budgeted Batch Bayesian Optimization (B3O).
Input: D 0 = {x i , y i } n 0 i=1 , #iter T , acquisition function α 1: for t = 1 to T do 2: Fit a GP from the data D t .

5:
Fit the Infinite GMM using s.

Experiments
We evaluate the B3O algorithm using both synthetic functions and real-world experimental settings. For the synthetic function evaluation, we utilize 8 functions across dimensions (1-10). In addition, we perform hyper-parameter tuning for three machine learning algorithms -support vector regression , Bayesian nonparametric multi-label classification (Nguyen, Gupta, Rana, Li, & Venkatesh, 2016a) and multi layer perceptron (Ruck, Rogers, Kabrisky, Oxley, & Suter, 1990). We further consider real-world experimental design for Aluminum-Scandium hardening. This heat treatment design involves searching for the best combination of temperatures and times in two stages to achieve the requisite alloy properties -in our case, hardness. We compare our results with baselines on best-found-value and the total number of evaluations, given a fixed number of iterations. We show that we outperform baselines in the best-found-value whilst requiring fewer evaluations. We compare the performance of our proposed approach with the state-of-the-art methods for batch BO.
• Gaussian process -Batch Upper Confidence Bound (GP-BUCB) (Desautels et al., 2014): this is a batch BO utilizing the variance of GP for finding the next point until the batch is filled.
• Rand (EI and UCB): this optimizes the acquisition function to find the first element in the batch (this step is identical to the sequential BO), then chooses a random sample within the bounds until the batch is filled.

Function, Dimension Objective function f (x)
Optimum f (x * ) Table 1: Benchmark functions and dimensions used ranging from 1 to 10.
• Local Penalization (LP) (González et al., 2016): This is currently the state-of-the-art method for batch BO which has been demonstrated to outperform most of other baselines (González et al., 2016). The source code is available in GPyOpt toolbox 1 .

EXPERIMENTAL SETTINGS
We denote the number of iterations by T , and the number of evaluations by N. For the sequential setting of BO, the number of iterations is identical to the number of evaluations, that is N = T while N > T for batch setting. The number of initial points n 0 for all methods is set as 3× D where D is the dimension. For fixed-batch approaches, the batch size at each iteration is set as n t = 3 for functions with D < 5 or n t = D for functions with D ≥ 5. In contrast, our method automatically determines the batch size n t which could differ at each iteration. Thus, we do not need to set the batch size n t for B3O in advance. The performance of the algorithms is compared for a fixed number of iterations T , given as 10 × D, and the total number of evaluated points is N = ∑ T t=0 n t . In all the experiments, the squared exponential (SE) kernel given as k(x, x ′ ) = exp −γ||x − x ′ || 2 is used in the underlying GP, where γ = 0.1 × D and D is the dimension. For methods using the UCB, β is fixed to 2 (following the setting used in (González et al., 2016)), which allows us to compare the different batch designs using the same acquisition function. The results are taken over 20 replicates with different initial values. We always maximize the objective function, maximizing − f for cases in which the goal is to find the minimum. All implementations are in Python. All simulations are done on Windows machine Core i7 Ram 24GB.
Although B3O is designed to work with any kind of acquisition function, in this paper we use B3O with UCB since we empirically notice that the UCB generally works better than EI for B3O. All the source codes and data are available for reproducibility at the link 2 .

TEST FUNCTIONS AND EVALUATION CRITERIA
Test functions are important to validate the performance of optimization algorithms. There have been many benchmark functions reported in the literature (Jamil & Yang, 2013). We select 5 popular 1. https://github.com/SheffieldML/GPyOpt 2. https://github.com/ntienvu/ICDM2016 B3O  benchmark functions, see Table 1. We consider two major criteria for evaluation including the bestfound-value f (x best ) and the total number of evaluations ∑ T t=0 n t .

Comparison -best found value
We examine the performance of B3O in finding the optimum of the chosen benchmark functions. We demonstrate that our method can find better optimal values (minimum) for a fixed amount of iterations T . We report the best-found-value w.r.t. iterations in Fig. 5. Our method achieves significantly better values than the baselines in 5 over 8 functions. In particular, B3O outperforms all baselines in Forrester 1D, Alpine2 5D, gSobol 5D, Alpine2 10D, gSobol 10D. In other three cases, B3O still provides relatively good values. The fixed batch approaches (e.g., LP) may suffer the under-specification (the number of specified batch size n t is smaller than the number of real peaks) at some iterations and thus negatively affect the final performance. As expected, the batch approaches are always better than the sequential approaches using EI and UCB because the number of evaluated points for batch methods is much higher than the sequential. This fact also highlights the advantage of batch BO over the standard (sequential) BO.

Comparison -number of evaluations
The next criteria for comparison is the total number of evaluations, N = ∑ T t=0 n t . The requirement of Bayesian optimization is to keep the number of evaluations to find the optimal value as low as possible, as these evaluations can be costly. For example, it can very expensive to perform a single experiment in material science (to cast an alloy testing), or it takes a few days to train a deep network on a large-scale dataset.
It is not natural to fix the number of evaluations per batch n t , since the number of peaks is unknown, and importantly, these peaks will be changed after new evaluations are done. Therefore, we may waste time and resources if we over-evaluate the number of points than we need. In particular, after all the actual peaks are detected, if we set the batch size n t large, we will get the While the existing approaches fix the batch size n t to a constant for all iterations, our B3O is more efficient by flexibly creating a suitable batch size per iterations. One the one hand, at the iteration 2 in gSobol function, the fixed-batch-size approaches suggest 5 points while we have 6 peaks. Because of underspecifying the peaks, the baselines may miss the optimum. On the other hand, at the iterations 8, 9 and 10, the baselines suggest 5 points while we have only a single peak. This overspecifying the peaks will waste time and resources to evaluate at redundant points.
noisy points which are possibly close to the already detected ones due to the effect of penalizing the peaks (González et al., 2016). Hence, these noisy points are not useful for evaluation. In contrast, under-evaluating the number of necessary points also brings detrimental effects of losing the optimal points.
B3O automatically identifies the suitable number of peaks from the acquisition function in each iteration, and thus does not have to resort to a fixed batch sizes. This is efficient without suffering any performance loss. Given a fixed number of iterations T , we summarize the total number of evaluations in Fig. 6   compared to the fixed-batch baselines, especially for high dimensional settings, e.g., Alpine2 10D and gSobol 10D.

Analysis of B3O -number of points per iteration
We demonstrate that B3O is flexible in identifying the required number of evaluations per batch. We study the number of estimated points in B3O at each iteration. In particular, we record the number of points n t at each iteration, see Fig. 7 for two functions -Dropwave 2D and gSobol 5D. The estimated points per batch in our approach is flexible and determined automatically in this Bayesian nonparametric setting. For some iterations, B3O only recommends one point per batch (see gSobol 5D in Fig. 7). This fact is essential for the homogeneous acquisition function which contains a single peak, whereas the fixed-batch approaches may be wasted to evaluate at some unnecessary points. Moreover, given the fixed batch size for all iterations, other approaches may suffer the effect of over-specification at some iterations and under-specification at other iterations.

Analysis of computational time
Next we compare the computational time at each iteration for the different batch approaches. Specifically, we are interested in the CPU time for finding the batch between different approaches. The sequential methods and the random batch approaches are the cheapest to compute than the other batch approaches. Thus, we do not compare with the sequential and the random batch approaches. Constant Liar (Ginsbourger et al., 2008) and BUCB (Desautels et al., 2014) take time for reestimation of the GP when the fake observations are added and noticeably slower when the number of data points N is large. The Local Penalization approach (González et al., 2016) consumes a considerable amount of time to estimate the Lipschitz constant. In addition, LP computes and maintains the penalized cost around the visited points. The cost for penalizing and optimizing will grow for high dimensions and/or large number of observed points.
B3O is generally competitive to LP (González et al., 2016) in terms of computation. We run faster for low dimensional functions (e.g., less than 6 dimensions). However, LP tends to run faster than B3O for 10D functions. We note that in high dimensional functions, the bottle neck in LP is estimating the Lipschitz constant whilst for B3O the batch generalized slice sampling is the bottle neck due to the nature of the sampling algorithm.

Analysis of batch generalized slice sampling (BGSS)
We investigate the efficiency of the BGSS presented in Section 3.3 for drawing samples under the acquisition function. We vary the observation dimension D -5, 10, 20, 40, 50 -and build the acquisition function α UCB . Then, we design two BGSS settings of size 100 and 200 and compare with the generalized slice sampling which is presented in Fig. 4. We record two important factors including computational time and the number of accepted samples under the curve. In Fig. 8, we present the simulation results that BGSS significantly outperforms the standard slice sampling in terms of computation (running faster) and efficiency (higher number of accepted points). For example, BGSS of the size 200 takes 32 seconds to get 2485 accepted data points at D = 10 dimension. These accepted points are lying under the curve and generally surrounding the peak. Therefore, these samples are beneficial to fit the IGMM since we are particularly interested in finding the peaks.

Tuning hyper-parameters for machine learning algorithms
Hyper-parameter settings greatly impact prediction accuracy of machine learning algorithms. The traditional way of performing hyper-parameter tuning has been to perform a grid search. In practice, however, Bayesian optimization has been shown (Bergstra, Bardenet, Bengio, & Kégl, 2011;Snoek et al., 2012;Thornton et al., 2013) to obtain better results in fewer experiments than required by a full grid search. We employ the batch Bayesian optimization to find the optimal hyper-parameters for support vector regression (SVR), Bayesian nonparametric multi-label classification (BNMC), and multi layer perceptron (MLP).   : Performance comparison on real-world experiments. We consider tuning hyperparameters for three machine learning algorithm: Support vector regression (SVR), Bayesian nonparametric multi-label classification (BNMC), and multi layer perceptron (MLP). We also conduct the Alluminium-scandium heat-treatment design.
Support Vector Regression. The support vector machines (Cortes & Vapnik, 1995) are wellknown for classification problem. It is also applicable for regression task as the support vector regression (SVR) model (Drucker, Burges, Kaufman, Smola, Vapnik, et al., 1997). Using Abalone dataset, we consider tuning the hyper-parameters for SVR with three parameters: C (regularizer parameter), epsilon (ε-insensitive loss) for regression and γ (RBF kernel function). We utilize the Python code in sklearn package.
Bayesian Nonparametric Multi-label Classification. We select to tune the hyper-parameters for the multi-label classification machine learning algorithm of BNMC (Nguyen et al., 2016a) on Scene dataset using the released source code, already constructed as a black-box function 3 . BNMC uses stochastic variational inference (SVI) and stochastic gradient descent (SGD) for learning the hidden correlation of label and feature for multi-label prediction. In particular, we optimize 6 hyperparameters: Dirichlet symmetric for feature and for label, learning rate for SVI and for SGD, truncation threshold and stick-breaking parameter. We aim to maximize the F1 score.
Multi Layer Perceptron. For the deep learning model, we use three layers in the model and evaluate it on the MNIST dataset. Our MLP model has 7 parameters including 4 parameters that include the number of hidden units and dropout rates for the first two layers and 3 parameters of learning rate, decay, and momentum for SGD which is used to learn the MLP. We utilize the Python code in sklearn package.
We compare B3O with baselines in Table 9a (columns 2-3). Our model performs better than all baselines in terms of RMSE for SVR, F1 score for BNMC and has the highest accuracy of 98.47% for MLP. We run 10 iterations for all methods in which the number of initial points is 9 and the batch size n t is 3 for the fixed batch BO. The number of evaluations in B3O is the lowest compared to other methods in this setting. Our model identifies the unknown number of peaks in the acquisition function, then evaluate at these points. Other fixed-batch approaches may evaluate points which are unnecessary. Thus, B3O requires less number of evaluation than the baselines (cf. Fig. 9b). Although we can run these machine learning algorithms in parallel, each evaluation takes a significant amount of time.

Bayesian optimization for experimental design
We consider the alloy hardening process of Aluminum-scandium (Wagner, Kampmann, & Voorhees, 1991) consisting of two stages: nucleation and precipitation coarsening. We aim to maximize the hardness for Aluminum-scandium alloys by designing the appropriate times and temperatures for the two stages in the process. The traditional approach is to use iterative trial and error, in which we choose some material design (temperature and time) based on intuition, past experience, or theoretical knowledge; and test the material in physical experiments; and then use what we learn from these experiments in choosing the material design to try next. This iterative process is repeated until some combination of success or exhaustion is reached! Since the parameter space grows exponentially, trial and error approach is ill-affordable in terms of both time and cost. Therefor, Bayesian optimization is a necessary choice for this experimental design.
We use B3O for the experimental design problem above and show that our method achieves the highest hardness using the smallest number of experiments (see Table 9a). The heat-treatment profile is shown in Fig. 10 using the best estimated parameters (times and temperatures at two steps) 3. https://github.com/ntienvu/ACML2016 BNMC Hour Celsius degree 0 Nucleation Precipitate coarsening Hardness = 88.3 Figure 10: Aluminum-scandium hardening process. We optimize the hardness y w.r.t four features x as temperatures T 1 , T 2 and times t 1 ,t 2 . and the best hardness achieved is 88.3. Given the times and temperatures setting, we evaluate the hardness using the standard kinematic KWN model used by metallurgists (Kampmann & Wagner, 1983) .

Conclusion
We have introduced a novel approach for batch Bayesian optimization. The proposed approach of B3O can identify the suitable batch size at each iteration that most existing approaches in batch BO are unable to do. We have presented the batch generalized slice sampling for drawing samples under the acquisition function. We perform extensive experiments on finding the optimal solution for 8 synthetics functions and evaluate the performance further on 4 real-world tasks. The experimental results highlight the ability of B3O in finding the optimum whilst requiring fewer evaluations than the baselines.

Acknowledgment
We thank Dr Paul G Sanders and Kyle J Deane for the real-world case study of Aluminum-scandium hardening process.


extfieldsofstudy:
Mathematics
Computer Science

text:
Online Learning with Local Permutations and Delayed Feedback

We propose an Online Learning with Local Permutations (OLLP) setting, in which the learner is allowed to slightly permute the \emph{order} of the loss functions generated by an adversary. On one hand, this models natural situations where the exact order of the learner's responses is not crucial, and on the other hand, might allow better learning and regret performance, by mitigating highly adversarial loss sequences. Also, with random permutations, this can be seen as a setting interpolating between adversarial and stochastic losses. In this paper, we consider the applicability of this setting to convex online learning with delayed feedback, in which the feedback on the prediction made in round $t$ arrives with some delay $\tau$. With such delayed feedback, the best possible regret bound is well-known to be $O(\sqrt{\tau T})$. We prove that by being able to permute losses by a distance of at most $M$ (for $M\geq \tau$), the regret can be improved to $O(\sqrt{T}(1+\sqrt{\tau^2/M}))$, using a Mirror-Descent based algorithm which can be applied for both Euclidean and non-Euclidean geometries. We also prove a lower bound, showing that for $M<\tau/3$, it is impossible to improve the standard $O(\sqrt{\tau T})$ regret bound by more than constant factors. Finally, we provide some experiments validating the performance of our algorithm.

Introduction
Online learning is traditionally posed as a repeated game where the learner has to provide predictions on an arbitrary sequence of loss functions, possibly even generated adversarially. Although it is often possible to devise algorithms with non-trivial regret guarantees, these have to cope with arbitrary loss sequences, which makes them conservative and in some cases inferior to algorithms not tailored to cope with worstcase behavior. Indeed, an emerging line of work considers how better online learning can be obtained on "easy" data, which satisfies some additional assumptions. Some examples include losses which are sampled i.i.d. from some distribution, change slowly in time, have a consistently best-performing predictor across time, have some predictable structure, mix adversarial and stochastic losses, etc. (e.g. Sani et al. [2014], Karnin and Anava [2016], Bubeck and Slivkins [2012], Seldin and Slivkins [2014], Hazan and Kale [2010], Chiang et al. [2012], Steinhardt and Liang [2014], Hazan and Kale [2011], Rakhlin and Sridharan [2013], Seldin and Slivkins [2014]).
In this paper, we take a related but different direction: Rather than explicitly excluding highly adversarial loss sequences, we consider how slightly perturbing them can mitigate their worst-case behavior, and lead to improved performance. Conceptually, this resembles smoothed analysis Spielman and Teng [2004], in which one considers the worst-case performance of some algorithm, after performing some perturbation to their input. The idea is that if the worst-case instances are isolated and brittle, then a perturbation will lead to easier instances, and better reflect the attainable performance in practice.
Specifically, we propose a setting, in which the learner is allowed to slightly reorder the sequence of losses generated by an adversary: Assuming the adversary chooses losses h 1 , . . . , h T , and before any losses are revealed, the learner may choose a permutation σ on {1, . . . , T }, satisfying max t |t − σ(t)| ≤ M for some parameter M , and then play a standard online learning game on losses h σ(1) , . . . , h σ(T ) . We denote this as the Online Learning with Local Permutations (OLLP) setting. Here, M controls the amount of power given to the learner: M = 0 means that no reordering is performed, and the setting is equivalent to standard adversarial online learning. At the other extreme, M = T means that the learner can reorder the losses arbitrarily. For example, the learner may choose to order the losses uniformly at random, making it a quasistochastic setting (the only difference compared to i.i.d. losses is that they are sampled without-replacement rather than with-replacement).
We argue that allowing the learner some flexibility in the order of responses is a natural assumption. For example, when the learner needs to provide rapid predictions on a high-frequency stream of examples, it is often immaterial if the predictions are not provided in the exact same order at which the examples arrived. Indeed, by buffering examples for a few rounds before being answered, one can simulate the local permutations discussed earlier.
We believe that this setting can be useful in various online learning problems, where it is natural to change a bit the order of the loss functions. In this paper, we focus on one well-known problem, namely online learning with delayed feedback. In this case, rather than being provided with the loss function immediately after prediction is made, the learner only receives the loss function after a certain number τ ≥ 1 of rounds. This naturally models situations where the feedback comes much more slowly than the required frequency of predictions: To give a concrete example, consider a web advertisement problem, where an algorithm picks an ad to display, and then receives a feedback from the user in the form of a click. It is likely that the algorithm will be required to choose ads for new users while still waiting for the feedback from the previous user.
For convex online learning with delayed feedback, in a standard adversarial setting, it is known that the attainable regret is on the order of O( √ τ T ), and this is also the best possible in the worst case Weinberger and Ordentlich [2002], Mesterharm [2005], Langford et al. [2009], Joulani et al. [2013], Quanrud and Khashabi [2015]. On the other hand, in a stochastic setting where the losses are sampled i.i.d. from some distribution, Agarwal and Duchi [2011] show that the attainable regret is much better, on the order of O( √ T + τ ). This gap between the worst-case adversarial setting, and the milder i.i.d. setting, hints that this problem is a good fit for our OLLP framework.
Thus, in this paper, we focus on online learning with feedback delayed up to τ rounds, in the OLLP framework where the learner is allowed to locally permute the loss functions (up to a distance of M ). First, we devise an algorithm, denoted as Delayed Permuted Mirror Descent, and prove that it achieves an expected regret bound of order O( T (τ 2 /M + 1)) assuming M ≥ τ . As M increases compared to τ , this regret bound interpolates between the standard adversarial √ τ T regret, and a milder √ T regret, typical of i.i.d. losses. As its name implies, the algorithm is based on the well-known online mirror descent (OMD) algorithm (see Hazan et al. [2016], Shalev-Shwartz et al. [2012), and works in the same generality, involving both Euclidean and non-Euclidean geometries. The algorithm is based on dividing the entire sequence of functions into blocks of size M and performing a random permutation within each block. Then, two copies of OMD are ran on different parts of each block, with appropriate parameter settings. A careful analysis, mixing adversarial and stochastic elements, leads to the regret bound.
In addition, we provide a lower bound complementing our upper bound analysis, showing that when M is significantly smaller than τ (specifically, τ /3), then even with local permutations, it is impossible to obtain a worse-case regret better than Ω( √ τ T ), matching (up to constants) the attainable regret in the standard adversarial setting where no permutations are allowed. Finally, we provide some experiments validating the performance of our algorithm. The rest of the paper is organized as follows: in section 2 we formally define the Online Learning with Local Permutation setting, section 3 describes the Delayed Permuted Mirror Descent algorithm and outlines its regret analysis, section 4 discusses a lower bound for the delayed setting with limited permutation power, section 5 shows experiments, and finally section 6 provides concluding remarks, discussion, and open questions. Appendix A contains most of the proofs.

Setting and Notation
Convex Online Learning. Convex online learning is posed as a repeated game between a learner and an adversary (assumed to be oblivious in this paper). First, the adversary chooses T convex losses h 1 , . . . , h T which are functions from a convex set W to R. At each iteration t ∈ {1, 2, . . . , T }, the learner makes a prediction w t , and suffers a loss of h t (w t ). To simplify the presentation, we use the same notation ∇h t (w) to denote either a gradient of h t at w (if the loss is differentiable) or a subgradient at w otherwise, and refer to it in both cases as a gradient. We assume that both w ∈ W and the gradients of any function h t in any point w ∈ W are bounded w.r.t. some norm: Given a norm · with a dual norm · * , we assume that the diameter of the space W is bounded by B 2 and that ∀w ∈ W, ∀h ∈ {h 1 , h 2 , ..., h T } : ∇h (w) * ≤ G. The purpose of the learner is to minimize her (expected) regret, i.e.
where the expectation is with respect to the possible randomness of the algorithm. Learning with Local Permutations. In this paper, we introduce and study a variant of this standard setting, which gives the learner a bit more power, by allowing her to slightly modify the order in which the losses are processed, thus potentially avoiding highly adversarial but brittle loss constructions. We denote this setting as the Online Learning with Local Permutations (OLLP) setting. Formally, letting M be a permutation window parameter, the learner is allowed (at the beginning of the game, and before any losses are revealed) to permute h 1 , . . . , h T to h σ −1 (1) , . . . , h σ −1 (T ) , where σ is a permutation from the set P erm := {σ : ∀t, |σ (t) − t| ≤ M }. After this permutation is performed, the learner is presented with the permuted sequence as in the standard online learning setting, with the same regret as before. To simplify notation, we let f t = h σ −1 (t) , so the learner is presented with the loss sequence f 1 , . . . , f T , and the regret is the same as the standard regret, i.e.
. Note that if M = 0 then we are in the fully adversarial setting (no permutation is allowed). At the other extreme, if M = T and σ is chosen uniformly at random, then we are in a stochastic setting, with a uniform distribution over the set of functions chosen by the adversary (note that this is close but differs a bit from a setting of i.i.d. losses). In between, as M varies, we get an interpolation between these two settings.
Learning with Delayed Feedback. The OLLP setting can be useful in many applications, and can potentially lead to improved regret bounds for various tasks, compared to the standard adversarial online learning. In this paper, we focus on studying its applicability to the task of learning from delayed feedback.
Whereas in standard online learning, the learner gets to observe the loss f t immediately at the end of iteration t, here we assume that at round t, she only gets to observe f t−τ for some delay parameter τ < T (and if t < τ , no feedback is received). For simplicity, we focus on the case where τ is fixed, independent of t, although our results can be easily generalized (as discussed in subsection 3.3). We emphasize that this is distinct from another delayed feedback scenario sometimes studied in the literature (Agarwal and Duchi [2011], Langford et al. [2009]), where rather than receiving f t−τ the learner only receives a (sub)gradient of f t−τ at w t−τ . This is a more difficult setting, which is relevant for instance when the delay is due to the time it takes to compute the gradient.

Algorithm and Analysis
Our algorithmic approach builds on the well-established online mirror descent framework. Thus, we begin with a short reminder of the Online Mirror Descent algorithm (see e.g. Hazan et al. [2016] for more details). Readers who are familiar with the algorithm are invited to skip to Subsection 3.1.
The online mirror descent algorithm is a generalization of online gradient descent, which can handle non-Euclidean geometries. The general idea is the following: we start with some point w t ∈ W, where W is our primal space. We then map this point to the dual space using a (striclty convex and continuously differentiable) mirror map ψ, i.e. ∇ψ (w t ) ∈ W * , then perform the gradient update in the dual space, and finally map the resulting new point back to our primal space W again, i.e. we want to find a point where ψ * is the dual function of ψ. This point, w t+ 1 2 , might lie outside our hypothesis class W, and thus we might need to project it back to our space W. We use the Bregman divergence associated to ψ to do this: Specific choices of the mirror map ψ leads to specific instantiations of the algorithms for various geometries. Perhaps the simplest example is ψ (x) = 1 2 x 2 2 , with associated Bregman divergence ψ (x, y) = 1 2 · x − y 2 . This leads us to the standard and well-known online gradient descent algorithm, where w t+1 is the Euclidean projection on the set W of w t − η · g t . Another example is the negative entropy mirror map ψ (x) = n i=1 x i · log (x i ), which is 1-strongly convex with respect to the 1-norm on the simplex W = x ∈ R n + : n i=1 x i = 1 . In that case, the resulting algorithm is the well-known multiplicative updates algorithm, where w t+1,i = w t,i ·exp(−ηg t,i )/ n j=1 w t,i · exp(−ηg t,j ). Instead of the 1-norm on the simplex, one can also consider arbitrary p-norms, and take ψ(x) = 1 2 · x 2 q , where q is the dual norm (satisfying 1/p + 1/q = 1).

The Delayed Permuted Mirror Descent Algorithm
Before describing the algorithm, we note that we will focus here on the case where the permutation window parameter M is larger than the delay parameter τ . If M < τ , then our regret bound is generally no better than the O( √ τ T ) obtainable by a standard algorithm without any permutations, and for M < τ /3, this is actually tight as shown in Section 4.
We now turn to present our algorithm, denoted as The Delayed Permuted Mirror Descent algorithm (see algorithm 1 below as well as figure 1 for a graphical illustration). First, the algorithm splits the time horizon T into M consecutive blocks, and performs a uniformly random permutation on the loss functions within each block. Then, it runs two online mirror descent algorithms in parallel, and uses the delayed gradients in  Figure 1: Scheme of predictions and updates of both parallel algorithms (best viewed in color; see text for details). Top color bars mark which iterations are in T 1 (purple lines) and which are in T 2 (green lines). Top timeline shows which predictor, w f or w s , is used to predict in each iteration. Middle timeline shows where gradients for updating w f come from (first τ iterations of the previous block), and lower timeline shows where gradients for updating w s come from (first M − τ iterations of the same block, each gradient from exactly τ rounds back). order to update two separate predictors -w f and w s , where w f is used for prediction in the first τ rounds of each block, and w s is used for prediction in the remaining M − τ rounds (here, f stands for "first" and s stands for "second"). The algorithm maintaining w s crucially relies on the fact that the gradient of any two functions in a block (at some point w) is equal, in expectation over the random permutation within each block. This allows us to avoid most of the cost incurred by delays within each block, since the expected gradient of a delayed function and the current function are equal. A complicating factor is that at the first τ rounds of each block, no losses from the current block has been revealed so far. To tackle this, we use another algorithm (maintaining w f ), specifically to deal with the losses at the beginning of each block. This algorithm does not benefit from the random permutation, and its regret scales the same as standard adversarial online learning with delayed feedback. However, as the block size M increases, the proportion of losses handled by w f decreases, and hence its influence on the overall regret diminishes.
The above refers to how the blocks are divided for purposes of prediction. For purposes of updating the predictor of each algorithm, we need to use the blocks a bit differently. Specifically, we let T 1 and T 2 be two sets of indices. T 1 includes all indices from the first τ time points of every block, and is used to update w f . T 2 includes the first M − τ indices of every block, and is used to update w s (see figure 1). Perhaps surprisingly, note that T 1 and T 2 are not disjoint, and their union does not cover all of {1, . . . , T }. The reason is that due to the random permutation in each block, the second algorithm only needs to update on some of the loss functions in each block, in order to obtain an expected regret bound on all the losses it predicts on.

Algorithm 1 Delayed Permuted Mirror Descent
Input: M , η f , η s Init: w f 1 = 0, w s 1 = 0, j f = j s = 1 Divide T to consecutive blocks of size M , and permute the losses uniformly at random within each block. Let f 1 , . . . , f T denote the resulting permuted losses.
Receive a loss function from τ places back: . If none exists (in the first τ iterations), take the 0 function.
Predict using w s js Receive a loss function from τ places back:

Analysis
The regret analysis of the Delayed Permuted Mirror Descent algorithm is based on a separate analysis of each of the two mirror descent sub-algorithms, where in the first sub-algorithm the delay parameter τ enters multiplicatively, but doesn't play a significant role in the regret of the second sub-algorithm (which utilizes the stochastic nature of the permutations). Combining the regret bound of the two sub-algorithms, and using the fact that the portion of losses predicted by the second algorithm increases with M , leads to an overall regret bound improving in M .
In the proof, to analyze the effect of delay, we need a bound on the distance between any two consequent predictors w t , w t+1 generated by the sub-algorithm. This depends on the mirror map and Bregman divergence used for the update, and we currently do not have a bound holding in full generality. Instead, we let Ψ (η f ,G) be some upper bound on w t+1 − w t , where the update is using step-size η f and gradients of norm ≤ G. Using Ψ (ηf ,G) we prove a general bound for all mirror maps. In Lemmas 3 and 4 in Appendix A.1, we show that for two common mirror maps (corresponding to online gradient descent and multiplicative weights), Ψ (η f ,G) ≤ c · η f G for some numerical constant c, leading to a regret bound of O( T (τ 2 /M + 1)). Also, we prove theorem 1 for 1-strongly convex mirror maps, although it can be generalized to any λ-strongly convex mirror map by scaling.
Theorem 1. Given a norm · , suppose that we run the Delayed Permuted Mirror Descent algorithm using a mirror map ψ which is 1-strongly convex w.r.t. · , over a domain W with diameter B 2 w.r.t the bregman divergence of ψ: ∀w, v ∈ W : ψ (w, v) ≤ B 2 , and such that the (sub)-gradient g of each loss function on any w ∈ W satisfies g * ≤ G (where · * is the dual norm of · ). Then the expected regret, given a delay parameter τ and step sizes η f , η s satisfies: , the regret is bounded by . similar to the standard adversarial learning case. However, as M increases, the regret gradually improves to O( √ T + τ ), which is the regret attainable in a purely stochastic setting with i.i.d. losses. The full proof can be found in appendix A.1.1, and we sketch below the main ideas.
First, using the definition of regret, we show that it is enough to upper-bound the regret of each of the two sub-algorithms separately. Then, by a standard convexity argument, we reduce this to bounding sums of terms of the form E for the second sub-algorithm (where w * f and w * s are the best fixed points in hindsight for the losses predicted on by the first and second sub-algorithms, respectively, and where for simplicity we assume the losses are differentiable). In contrast, we can use the standard analysis of mirror descent, using delayed gradients, to get a bound for the somewhat different terms E[ w f t − w * f , ∇f t−τ (w f t−τ ) ] for the first sub-algorithm, and E [ w s t − w * s , ∇f t−τ (w s our assumptions on the gradients of the losses and the distance between consecutive predictors produced by the first sub-algorithm. Overall, we get two regret bounds, one for each sub-algorithm. The regret of the first sub-algorithm scales with τ , similar to the no-permutation setting, but the sub-algorithm handles only a small fraction of the iterations (the first τ in every block of size M ). In the rest of the iterations, where we use the second sub-algorithm, we get a bound that resembles more the stochastic case, without such dependence on τ . Combining the two, the result stated in Theorem 1 follows.

Handling Variable Delay Size
So far, we discussed a setting where the feedback arrives with a fixed delay of size τ . However, in many situations the feedback might arrive with a variable delay size τ t at any iteration t, which may raise a few issues. First, feedback might arrive in an asynchronous fashion, causing us to update our predictor using gradients from time points further in past after already using more recent gradients. This complicates the analysis of the algorithm. A second, algorithmic problem, is that we could also possibly receive multiple feedbacks simultaneously, or no feedback at all, in certain iterations, since the delay is of variable size. One simple solution is to use buffering and reduce the problem to a constant delay setting. Specifically, we assume that all delays are bounded by some maximal delay size τ . We would like to use one gradient to update our predictor at every iteration (this is mainly for ease of analysis, practically one could update the predictor with multiple loss functions in a single iteration). In order to achieve this, we can use a buffer to store loss functions that were received but have not been used to update the predictors yet. We define Grad f and Grad s , two buffers that will contain gradients from time points in T 1 or T 2 , correspondingly. Each buffer is of size τ . If we denote by F t the set of function that have arrived in time t, we can simply store loss functions that have arrived asynchronously in the buffers defined above, sort them in ascending order, and take the delayed loss function from exactly τ iterations back in the update step. This loss function must be in the appropriate buffer since the maximal delay size is τ . From this moment on, the algorithm can proceed as usual and its analysis still applies.

Lower Bound
In this section, we give a lower bound in the setting where M < τ 3 with all feedback having delay of exactly τ . We will show that for this case, the regret bound cannot be improved by more than a constant factor over the bound of the adversarial online learning problem with a fixed delay of size τ , namely Ω √ τ T for a sequence of length T . We hypothesize that this regret bound also cannot be significantly improved for any M = O(τ ) (and not just τ /3). However, proving this remains an open problem.
Theorem 2. For every (possible randomized) algorithm A with a permutation window of size M ≤ τ 3 , there exists a choice of linear, 1-Lipschitz functions over [−1, 1] ⊂ R, such that the expected regret of A after T rounds (with respect to the algorithm's randomness), is For completeness, we we also provide in appendix A.2 a proof that when M = 0 (i.e. no permutations allowed), then the worst-case regret is no better than Ω( √ τ T ). This is of course a special case of Theorem 2, but applies to the standard adversarial online setting (without any local permutations), and the proof is simpler. The proof sketch for the setting where no permutation is allowed was already provided in Langford et al. [2009], and our contribution is in providing a full formal proof. The proof in the case where M = 0 is based on linear losses of the form f t = α t · w t over [−1, +1], where α t ∈ {−1, +1}. Without permutations, it is possible to prove a Ω( √ τ T ) lower bound by dividing the T iterations into blocks of size τ , where the α values of all losses at each block is the same and randomly chosen to equal either +1 or −1. Since the learner does not obtain any information about this value until the block is over, this reduces to adversarial online learning over T /τ rounds, where the regret at each round scales linearly with τ , and overall regret at least Ω(τ T /τ ) = Ω( √ τ T ). In the proof of theorem 2, we show that by using a similar construction, even with permutations, having a permutation window less than τ /3 still means that the α values would still be unknown until all loss functions of the block are processed, leading to the same lower bound up to constants.
The formal proof appears in the appendix, but can be sketched as follows: first, we divide the T iterations into blocks of size τ /3. Loss functions within each block are identical, of the form f t = α t · w t , and the value of α per block is chosen uniformly at random from {−1, +1}, as before. Since here, the permutation window M is smaller than τ /3, then even after permutation, the time difference between the first and last time we encounter an α that originated from a single block is less than τ . This means that by the time we get any information on the α in a given block, the algorithm already had to process all the losses in the block, which leads to the same difficulty as the no-permutation setting. Specifically, since the predictors chosen by the algorithm when handling the losses of the block do not depend on the α value in that block, and that α is chosen randomly, we get that the expected loss of the algorithm at any time point t equals 0. Thus, the cumulative loss across the entire loss sequence is also 0. In contrast, for w * , the optimal predictor in hindsight over the entire sequence, we can prove an expected accumulated loss of −Ω( √ τ T ) after T iterations, using Khintchine inequality and the fact that the α's were randomly chosen per block. This leads us to a lower bound of expected regret of order √ τ T , for any algorithm with a local permutation window of size M < τ /3.

Experiments
We consider the adversarial setting described in section 4, where an adversary chooses a sequence of functions such that every τ functions are identical, creating blocks of size τ of identical loss functions, of the form f t (w t ) = α t · w t where α t is chosen randomly in {−1, +1} for each block. In all experiments we use T = 10 5 rounds, a delay parameter of τ = 200, set our step sizes according to the theoretical analysis, and report the mean regret value over 1000 repetitions of the experiments.
In our first experiment, we considered the behavior of our Delayed Permuted Mirror Descent algorithm, for window sizes M > τ , ranging from τ + 1 to T . In this experiment, we chose the α values randomly, while ensuring a gap of 200 between the number of blocks with +1 values and the number of blocks with −1 values (this ensures that the optimal w * is a sufficiently strong competitor, since otherwise the setting is too "easy" and the algorithm can attain negative regret in some situations). The results are shown in Figures 2  and 3, where the first figure presents the accumulated regret of our algorithm over time, whereas the second figure presents the overall regret after T rounds, as a function of the window size M .
When applying our algorithm in this setting with different values of M > τ , ranging from M = τ + 1 and up to M = T , we get a regret that scales from the order of the adversarial bound to the order of the stochastic bound depending on the window size, as expected by our analysis. For all window sizes greater than 5 · τ , we get a regret that is in the order of the stochastic bound -this is not surprising, since after the permutation we get a sequence of functions that is very close to an i.i.d. sequence, in which case any algorithm can be shown to achieve O( √ T ) regret in expectation. Note that this performance is better than that predicted by our theoretical analysis, which implies an O( √ T ) behavior only when M ≥ Ω(τ 2 ). It is an open and interesting question whether it means that our analysis can be improved, or whether there is a harder construction leading to a tighter lower bound.
In our second experiment, we demonstrate the brittleness of the lower bound construction for standard online learning with delayed feedback, focusing on the M < τ regime. Specifically, we create loss functions with blocks as before (where following the lower bound construction, the α values in each block of size τ = 200 is chosen uniformly at random). Then, we perform a random permutation over consecutive windows of size M (ranging from M = 0 up to M = 9 10 τ in intervals of 1 10 τ ). Finally, we run standard Online Gradient Descent with delayed gradients (and fixed step size 1/ √ T ), on the permuted losses. The results are presented in Figure 4.
For window sizes M < τ 2 we see that the regret is close to the adversarial bound, whereas as we increase the window size the regret decreases towards the stochastic bound. This experiment evidently shows that this hardness construction is indeed brittle, and easily breaks in the face of local permutations, even for window sizes M < τ .

Discussion
We presented the OLLP setting, where a learner can locally permute the sequence of examples from which she learns. This setting can potentially allow for improved learning in many problems, where the worst-case regret is based on highly adversarial yet brittle constructions. In this paper, we focused on the problem of learning from delayed feedback in the OLLP setting, and showed how it is possible to improve the regret by allowing local permutations. Also, we proved a lower bound in the situation where the permutation window is significantly smaller than the feedback delay, and showed that in this case, permutations cannot allow for a better regret bound than the standard adversarial setting. We also provided some experiments, demonstrating the power of the setting as well as the feasibility of the proposed algorithm. An interesting open question is what minimal permutation size allows non-trivial regret improvement, and whether our upper bound in Theorem 1 is tight. As suggested by our empirical experiments, it is possible that even small local permutations are enough to break highly adversarial sequences and improve performance in otherwise worst-case scenarios. Another interesting direction is to extend our results to a partial feedback (i.e. bandit) setting. Finally, it would be interesting to study other cases where local permutations allow us to interpolate between fully adversarial and more benign online learning scenarios.

A.1 Analysis Of The Delayed Permuted Mirror Descent Algorithm
We will use throughout the proofs the well known Pythagorean Theorem for Bregman divergences, and the 'projection' lemma that considers the projection step in the algorithm.

Lemma 1. Pythagorean Theorem for Bregman divergences
Let v be the projection of w onto a convex set W w.r.t Bregman divergence ψ : v = argmin u∈W ψ (u, w), Lemma 2. Projection Lemma Let W be a closed convex set and let v be the projection of w onto W, namely, v = argmin x∈W x − w 2 . Then, for every u ∈ W, w − u 2 − v − u 2 ≥ 0 The following lemma gives a bound on the distance between two consequent predictions when using the Euclidean mirror map: Lemma 3. Let g ∈ R n s.t. g 2 < G, W a convex set, and η > 0 be fixed. Let w ∈ W and w 2 = w − η · g. Then, for w = argmin u∈W w 2 − u 2 2 , we have that w − w ≤ η · G Proof. From the projection lemma: w 2 − w 2 2 ≥ w − w 2 2 and so: w 2 − w 2 ≥ w − w 2 . From definition: w 2 − w 2 = η · g 2 ≤ η · G. and so we get: We prove a modification of Lemma 2 given in Menache et al. [2014] in order to bound the distance between two consequent predictions when using the negative entropy mirror map: Lemma 4. Let g ∈ R n s.t. g 1 ≤ G for some G > 0 and let η > 0 be fixed, with η < 1 √ 2·G . For any distribution vector w in the n − simplex, if we define w to be the new distribution vector ∀i ∈ {1, ..., n} , w i = w i · exp (−η · g i ) n j=1 w j · exp (−η · g j ) Proof. Since g ∞ < G and η < 1 √ 2·G we get that ∀i : |η · g i | < 1. We have that: Since w 1 = 1, we can apply Holder's inequality, and upper bound the above by Using the inequality 1 − x ≤ exp (−x) ≤ 1 1+x for all |x| ≤ 1, we know that 1 − η · g i ≤ exp (−η · g i ) ≤ 1 1 + η · g i and since −ηG ≤ η · g i ≤ ηG we have that 1 − η · g i ≤ exp (−η · g i ) ≤ 1 1 + η · g i ≤ 1 1 − ηG and so we get: Using again the fact that −ηG ≤ η · g i ≤ ηG, we have Now, since ηG < 1, we get that: and so we can conclude that Since η < 1 √ 2G , we get (η · G) 2 < 1 2 . Thus we get: which gives us our desired bound.
With the above two lemmas in hand, we bound the distance between consequent predictors by cηG, where c is a different constant in each mirror map: c = 1 for the euclidean case, and c = 3 for the negative entropy mirror map. Note that both mapping are 1-strongly convex with respect to their respective norms. For other mappings with a different strong convexity constant, one would need to scale the step sizes according to the strong convexity parameter in order to get the bound.

A.1.1 Proof of Theorem 1
We provide an upper bound on the regret of the algorithm, by competing against the best fixed action in each one of the sets of iterations-the first τ iterations and the last M − τ iterations in each block. This is an upper bound on competing against the best fixed predictor in hindsight for the entire sequence. Formally, we bound: f t (w) and w * s = argmin where expectation is taken over the randomness of the algorithm. The diameter of the domain W is bounded by B 2 , and so ψ w * f , w f 0 ≤ B 2 and ψ (w * s , w s 0 ) ≤ B 2 . We start with a general derivation that will apply both for w s and for w f simultaneously. For the following derivation we use the notation w j , w j+1 omitting the f, s superscript, for denoting subsequent updates of the predictor vector, whether it is w s or w f . Denote by g j the gradient used to update w j , i.e., ∇ψ w j+ 1 2 = ∇ψ (w j ) − η · g j , and w j+1 = argmin w∈W ψ w, w j+ 1 2 .
Then the expected sum of losses is: Eᾱ (1→i,j−τ −1) w i,j · 1 2 · 1 + 1 2 · (−1) = 0 The last equality is true because every first α in any block has probability 1 2 to be either +1 or −1. We now continue to the expected sum of losses for the optimal choice of w * = argmin w∈W


extfieldsofstudy:
Computer Science
Mathematics

text:
Slow mixing for Latent Dirichlet allocation

Markov chain Monte Carlo (MCMC) algorithms are ubiquitous in probability theory in general and in machine learning in particular. A Markov chain is devised so that its stationary distribution is some probability distribution of interest. Then one samples from the given distribution by running the Markov chain for a"long time"until it appears to be stationary and then collects the sample. However these chains are often very complex and there are no theoretical guarantees that stationarity is actually reached. In this paper we study the Gibbs sampler of the posterior distribution of a very simple case of Latent Dirichlet Allocation, an attractive Bayesian unsupervised learning model for text generation and text classification. It turns out that in some situations, the mixing time of the Gibbs sampler is exponential in the length of documents and so it is practically impossible to properly sample from the posterior when documents are sufficiently long.

Introduction
Markov chain Monte Carlo (MCMC) is a powerful tool for sampling from a given probability distribution on a very large state space, where direct sampling is difficult, in part because of the size of the state space and in part because of normalizing constants that are difficult to compute.
In machine learning in particular, MCMC algorithms are extremely common for sampling from posterior distributions of Bayesian probabilistic models. The posterior distribution given observed data is then difficult to sample from for the reasons above. One then designs an (irreducible aperiodic) Markov chain whose stationary distribution is precisely the targeted posterior. This is usually fairly easy since the posterior is usually easy to compute up to the normalizing constant (the denominator in Bayes formula). One usually uses Gibbs sampling or the related Metropolis-Hastings algorithm.
Gibbs sampling in these situations can be described as follows. Our state space is a finite set of random variables X = {X a } a∈A , where X a ∈ T for some measurable space T , so that X ∈ T A and the targeted distribution is a given probability measure P on T A . In order to sample from P one starts a Markov chain on T A whose updates are given by first choosing an index a ∈ A at random and then choosing a new value of X a according to the conditional distribution of X a given all X b , b ∈ A \ {a}. Under mild conditions, this chain converges in distribution to P. The Markov chain is then run for a "long time" and then a sample, hopefully approximately from P, is collected. A key question here is for how long the chain actually has to be run, in order for the distribution after that time to be a good approximation of P. Since A is usually large, the number of steps needed should at least be no more than polynomial in the size of A for Gibbs sampling to be feasible. In almost all practical cases, the structure of the sample space and the probability measure P is so complex that is virtually impossible to make a rigorous analysis of the mixing rate. However it may be possible to consider some very simplified special cases. In this paper, we will analyse a special case of Latent Dirichlet allocation, henceforth LDA for short, and demonstrate for such a simple special case that mixing can indeed be a problem.
LDA is model used to classify documents according to their topics. We have a large corpus of documents and we want to determine for each word in each document which topic it comes from. Knowing this we can also classify the documents according to the proportion of words of the different topics it contains. The setup in LDA is that one has a fixed number D of documents of lengths N d a fixed set of topics t 1 , t 2 , . . . , t s and a fixed set of words w 1 , w 2 , . . . , w v . These are specified in advance. The number of topics is usually not large, whereas the number of words is. Next, for each document d = 1, . . . , D, a multinomial distribution θ d = (θ d (1), . . . , θ d (s)) over topics is chosen according to a Dirichlet prior with a known parameter α = (α 1 , . . . , α s ). For each topic t i a multinomial distribution φ i = (φ i (1), . . . , φ i (v)) according to a Dirichlet prior with parameter β = (β 1 , . . . , β v ) independently of each other and of the θ d :s. Given these, the corpus is then generated by for each position p = 1, . . . , N d in each document d picking a topic z d,p according to θ d and then picking the word at that position according to φ z d,p , doing this independently for all positions. Note that the LDA model is a so called "bag of words" model, i.e. it is invariant under permutations within each document.
In this paper, we will consider the very simplified case with D = s = v = 2, N 1 = N 2 = m and α = β = (1, 1) and study the mixing time asymptotics as m → ∞. To simplify the notation, denote the two topics by A and B and the two words by 1 and 2. Define n ij as the number of occurrences of the word j in document i, i, j = 1, 2 and write n i. = n i1 + n i2 (which by assumption equals m) and n .j = n 1j + n 2j and n .. = i,j n ij = 2m. We consider the mixing time for Gibbs sampling of the posterior in a seemingly typical case, namely that the number of 1:s in the first document is 3m/10 and in the second document 6m/10. Our result is the following. Theorem 1.1 Consider the case n 11 = 3m/10 and n 21 = 6m/10. Then there exists a λ > 0 such that for each 0 < κ < 3/4, Remark. The point of this paper is to demonstrate that the mixing time issue can be a real problem for Bayesian inference in machine learning in general and not for LDA in particular. There are also other methods for estimation that seem to work well for LDA, in particular variational inference. Furthermore, experimental results seem to be fairly well in line with what one would expect from a topic classification algorithm.
Before moving on to the proof of Theorem 1.1, we formally introduce the concept of mixing time. Let {X t } ∞ t=0 be a discrete time Markov chain on the finite state space S and for s ∈ S, let P s be the underlying probability measure under X 0 = s. Definition 1.2 Let µ and ν be two probability measures on S. The the total variation distance between µ and ν is given by The essence of Theorem 1.1 is that even after an exponentially long time, the distribution of the state of the Gibbs sampler is concentrated on a set whose probability mass according to the targeted posterior is at most 1/4.

Proof of Theorem 1.1
The following two combinatorial lemmas will be needed. Lemma 2.1 Let X be a standard uniform random variable and 0 ≤ k ≤ n. Then Proof. First recall that for any m = 1, 2, . . ., Next observe that the result holds true for any n with k = n. We want to prove that the claimed result holds for (n, k) = (p, r) for some arbitrary 0 ≤ r ≤ p. We do this by induction. We may then assume that the result is true with (n, k) = (m, l) for any m and l with either m < p or m = p and l > r. Then, taking where the fourth equality follows form the induction hypothesis and the last equality holds because Proof. Both sides of the equality are equal to the multinomial coefficient

✷
Let W = (w 11 , w 12 , . . . , w 1m , w 21 , . . . , w 2m ) be the words in our corpus, let Z = (z 11 , z 12 , . . . , z 2m ) be the latent topics and Z d be the latent topics in document d. Let also θ i be the probability that z d1 = A, d = 1, 2 and let φ t be the conditional probability that w dj = 1 given that z dj = t, t = A, B. In the case under study, these four quantities are all independent standard uniform random variables. We begin by determining the posterior distribution.
. Here we have used the notation k dj for the number words in document d with the topic being A and the word being j and the same dot notation for the k:s as for the n:s. The last equality is Lemma 2.1. For the second factor we have analogously, again using Lemma 2.1, Hence, ignoring factors that do not depend on the k:s and using Lemma 2.2 for the second equality This expression only depends on z via k = k(z) := (k 11 , k 12 , k 21 , k 22 ). Identifying all z having the same k(z), we have (regarding ,with some abuse of notation, a k also as the equivalence class consisting of all z:s having that particular k dj :s) that for any k 1 and k 2 , all z ∈ k 1 have the same probability of transitioning into k 2 . Hence the process where we only record the k:s is a lumped Markov chain and it suffices to regard this chain, whose state space is where of course K = k(Z) and C 0 is the normalizing constant.
For a = 0, the right hand side equals − log(11/7) which is smaller than −9/20. Analogously, the partial derivatives of log g with respect to b, c and d are less than −1/4 (we do not get exactly the same numbers here). It follows that for a sufficiently small ǫ > 0, For such a k, let V k be the number of visits to k of the lumped Markov chain starting at k 0 until it returns to k 0 . Then by an intuitive and well known fact in Markov theory (see e.g. [1], Lemma 6 of Chapter 2) Hence By a Bonferroni bound we get for sufficiently large m (ǫm > 170 is sufficient). Hence, starting at k 0 , the probability of getting more than distance ǫm from k 0 in e Taking the logarithm and differentiating with respect to a gives ∂ ∂a log g(a, 0, 0, 0) = log 1 2 − a 2 )(1 − 10a 3 )( 2 3 + 10a 9 )a ( 1 2 + a 2 ) 10a 3 ( 1 3 − 10a 9 )(1 − a) and taking a = 0, this becomes − log(5/3) < −1/2. Repeating the above argument, it follows for some ν > 0 that starting from k 1 , the probability of getting more than distance νm from k 1 in e ν 9 m steps is at most e − ν 72 m . Combining this with the analogous result for k 0 , this proves Theorem 1.1.