 # A Nonparametric Conjugate Prior Distribution for the Maximizing Argument of a Noisy Function

We propose a novel Bayesian approach to solve stochastic optimization problems that involve finding extrema of noisy, nonlinear functions. Previous work has focused on representing possible functions explicitly, which leads to a two-step procedure of first, doing inference over the function space and second, finding the extrema of these functions. Here we skip the representation step and directly model the distribution over extrema. To this end, we devise a non-parametric conjugate prior based on a kernel regressor. The resulting posterior distribution directly captures the uncertainty over the maximum of the unknown function. We illustrate the effectiveness of our model by optimizing a noisy, high-dimensional, non-convex objective function.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Historically, the fields of statistical inference and stochastic optimization have often developed their own specific methods and approaches. Recently, however, there has been a growing interest in applying inference-based methods to optimization problems and vice versa [1, 2, 3, 4]. Here we consider stochastic optimization problems where we observe noise-contaminated values from an unknown nonlinear function and we want to find the input that maximizes the expected value of this function.

The problem statement is as follows. Let be a metric space. Consider a stochastic function mapping a test point to real values characterized by the conditional pdf . Consider the mean function

 ¯f(x):=E[y|x]=∫yP(y|x)dy. (1)

The goal consists in modeling the optimal test point

 x∗:=argmaxx{¯f(x)}. (2)

Classic approaches to solve this problem are often based on stochastic approximation methods . Within the context of statistical inference, Bayesian optimization methods have been developed where a prior distribution over the space of functions is assumed and uncertainty is tracked during the entire optimization process [6, 7]. In particular, non-parametric Bayesian approaches such as Gaussian Processes have been applied for derivative-free optimization [8, 9], also within the context of the continuum-armed bandit problem . Typically, these Bayesian approaches aim to explicitly represent the unknown objective function of (1) by entertaining a posterior distribution over the space of objective functions. In contrast, we aim to model directly the distribution of the maximum of (2) conditioned on observations.

The paper is structured as follows. Section 2 gives a brief description of the model suitable for direct implementation. The model is then derived in Section 3. Section 4 presents experimental results. Section 4 concludes.

## 2 Description of the Model Figure 1: a) Given an estimate h of the mean function ¯f(left), a simple probability density function over the location of the maximum x∗ is obtained using the transformation P(x∗)∝exp{αh(x∗)}, where α>0 plays the role of the precision (right). b) Illustration of the Gramian matrix for different test locations. Locations thar are close to each other produce large off-diagonal entries.

Our model is intuitively straightforward and easy to implement111Implementations can be downloaded from http://www.adaptiveagents.org/argmaxprior. Let be an estimate of the mean constructed from data (Figure 1a, left). This estimate can easily be converted into a posterior pdf over the location of the maximum by first multiplying it with a precision parameter and then taking the normalized exponential (Figure 1a, right)

 P(x∗|Dt)∝exp{α⋅h(x∗)}.

In this transformation, the precision parameter controls the certainty we have over our estimate of the maximizing argument: expresses almost no certainty, while expresses certainty. The rationale for the precision is: the more distinct inputs we test, the higher the precision—testing the same (or similar) inputs only provides local information and therefore should not increase our knowledge about the global maximum. A simple and effective way of implementing this idea is given by

 (3)

where , , , and are parameters of the estimator: is the precision we gain for each new distinct observation; is the number of prior points; is a finite, symmetric kernel function; is a prior precision function; and is a prior estimate of .

In (3), the mean function is estimated with a kernel regressor , and the total effective number of locations is calculated as the sum of the prior locations and the number of distinct locations in the data . The latter is estimated by multiplying the number of data points with the coefficient

 ∑iK(xi,xi)∑i∑jK(xi,xj)∈(0,1],

i.e. the ratio between the trace of the Gramian matrix and the sum of its entries. Inputs that are very close to each other will have overlapping kernels, resulting in large off-diagonal entries of the Gramian matrix—hence decreasing the number of distinct locations (Figure 1b).

The expression for the posterior can be calculated, up to a constant factor, in quadratic time in the number of observations. It can therefore be easily combined with Markov chain Monte Carlo methods (MCMC) to implement stochastic optimizers as illustrated in Section

4.

## 3 Derivation of the Model

### 3.1 Function-Based, Indirect Model

Our first task is to derive an indirect Bayesian model for the optimal test point that builds its estimate via the underlying function space. Let be the set of hypotheses, and assume that each hypothesis corresponds to a stochastic mapping . Let be the prior222For the sake of simplicity, we neglect issues of measurability of . over and let the likelihood be . Then, the posterior of is given by

 P(g|{yt},{xt})=P(g)P({yt}|g,{xt})P({yt}|{xt})=P(g)∏tP(yt|g,xt)P({yt}|{xt}). (4)

For each , let be the subset of functions such that for all , 333Note that we assume that the mean function is bounded and that it has a unique maximizing test point.. Then, the posterior over the optimal test point is given by

 P(x∗|{yt},{xt})=∫G(x∗)P(g|{yt},{xt})dg, (5)

This model has two important drawbacks: (a) it relies on modeling the entire function space , which is potentially much more complex than necessary; (b) it requires calculating the integral (5), which is intractable for virtually all real-world problems.

### 3.2 Domain-Based, Direct Model

We want to arrive at a Bayesian model that bypasses the integration step suggested by (5) and directly models the location of optimal test point . The following theorem explains how this direct model relates to the previous model.

###### Theorem 1.

The Bayesian model for the optimal test point is given by

 P(x∗) =∫G(x∗)P(g)dg (prior) P(yt|x∗,xt,Dt−1) =∫G(x∗)P(yt|g,xt)P(g)∏t−1k=1P(yk|g,xk)dg∫G(x∗)P(g)∏t−1k=1P(yk|g,xk)dg, (likelihood)

where is the set of past tests.

###### Proof.

Using Bayes’ rule, the posterior distribution can be rewritten as

 P(x∗)∏tP(yt|x∗,xt,Dt−1)P({yt}|{xt}). (6)

Since this posterior is equal to (5), one concludes (using (4)) that

 P(x∗)∏tP(yt|x∗,xt,Dt−1)=∫G(x∗)P(g)∏tP(yt|g,xt)dg.

Note that this expression corresponds to the joint . The prior is obtained by setting . The likelihood is obtained as the fraction

 P(yt|x∗,xt,Dt−1)=P(x∗,{yk}tk=1|{xk}tk=1)P(x∗,{yk}t−1k=1|{xk}t−1k=1),

where it shall be noted that the denominator doesn’t change if we add the condition . ∎

From Theorem 1 it is seen that although the likelihood model for the indirect model is i.i.d. at each test point, the likelihood model for the direct model depends on the past tests , that is, it is adaptive. More critically though, the likelihood function’s internal structure of the direct model corresponds to an integration over function space as well—thus inheriting all the difficulties of the indirect model.

### 3.3 Abstract Properties of the Likelihood Function

There is a way to bypass modeling the function space explicitly if we make a few additional assumptions. We assume that for any , the mean function is continuous and has a unique maximum. Then, the crucial insight consists in realizing that the value of the mean function inside a sufficiently small neighborhood of is larger than the value outside of it (see Figure 2a).

We assume that, for any and any , let denote the open -ball centered on . The functions in fulfill the following properties:

1. Continuous: Every function is such that its mean is continuous and bounded.

2. Maximum: For any , the functions are such that for all and all , . Figure 2: Illustration of assumptions. a) Three functions from G(x∗). They all have their maximum located at x∗∈X. b) Schematic representation of the likelihood function of x∗∈X conditioned on a few observations. The curve corresponds to the mean and the shaded area to the confidence bounds. The density inside of the neighborhood is unique to the hypothesis x∗, while the density outside is shared amongst all the hypotheses. c) The log-likelihood ratio of the hypotheses x∗1 and x∗2 as a function of the test point x. The kernel used in the plot is Gaussian.

Furthermore, we impose a symmetry condition on the likelihood function. Let and be in , and consider their associated equivalence classes and . There is no reason for them to be very different: in fact, they should virtually be indistinguishable outside of the neighborhoods of  and . It is only inside of the neighborhood of when becomes distinguishable from the other equivalence classes because the functions in systematically predict higher values than the rest. This assumption is illustrated in Figure 2b. In fact, taking the log-likelihood ratio of two competing hypotheses

 logP(yt|x∗1,xt,Dt−1)P(yt|x∗2,xt,Dt−1)

for a given test location should give a value equal to zero unless is inside of the vicinity of or (see Figure 2c). In other words, the amount of evidence a hypothesis gets when the test point is outside of its neighborhood is essentially zero (i.e. it is the same as the amount of evidence that most of the other hypotheses get).

### 3.4 Likelihood and Conjugate Prior

Following our previous discussion, we propose the following likelihood model. Given the previous data and a test point , the likelihood of the observation is

 P(yt|x∗,xt,Dt−1)=1Z(xt,Dt−1)λ(yt|xt,Dt−1)exp{αt⋅ht(x∗)−αt−1⋅ht−1(x∗)}, (7)

where: is a normalizing constant;

is a posterior probability over

given and the data ; is a precision measuring the knowledge we have about the whole function given by

 α0:=ρ⋅ξandαt:=ρ⋅(ξ+∑iK(xi,xi)∑i∑jK(xi,xj))

where is a precision scaling parameter; is a parameter representing the number prior locations tested; and is an estimate of the mean function given by

 h0(x∗):=y0(x∗)andht(x∗):=∑ti=1K(xi,x∗)yi+K0(x∗)y0(x∗)∑ti=1K(xi,x∗)+K0(x∗).

In the last expression, corresponds to a prior estimate of with prior precision . Inspecting (7), we see that the likelihood model favours positive changes to the estimated mean function from new, unseen test locations. The pdf does not need to be explicitly defined, as it will later drop out when computing the posterior. The only formal requirement is that it should be independent of the hypothesis .

We propose the conjugate prior

 P(x∗)=1Z0exp{α0⋅g0(x∗)}=1Z0exp{ξ⋅y0(x∗)}. (8)

The conjugate prior just encodes a prior estimate of the mean function. In a practical optimization application, it serves the purpose of guiding the exploration of the domain, as locations with high prior value are more likely to contain the maximizing argument.

Given a set of data points , the prior (8) and the likelihood (7) lead to a posterior given by

 P(x∗|Dt) =P(x∗)∏tk=1P(yk|x∗,xk,Dk−1)∫XP(x′)∏tk=1P(yk|x′,xk,Dk−1)dx′ =exp{∑tk=1αk⋅hk(x∗)−αk−1⋅hk−1(x∗)}Z−10∏tk=1Z(xk,Dk−1)−1∫Xexp{∑tk=1αk⋅hk(x′)−αk−1⋅hk−1(x′)}Z−10∏tk=1Z(xk,Dk−1)−1dx′ =exp{αt⋅ht(x∗)}∫Xexp{αt⋅ht(x′)}dx′. (9)

Thus, the particular choice of the likelihood function guarantees an analytically compact posterior expression. In general, the normalizing constant in (9) is intractable, which is why the expression is only practical for relative comparisons of test locations. Substituting the precision  and the mean function estimate  yields

 P(x∗|Dt)∝exp{ρ⋅(ξ+t⋅∑iK(xi,xi)∑i∑jK(xi,xj))⋅∑iK(xi,x∗)yi+K0(x∗)y0(x∗)∑iK(xi,x∗)+K0(x∗)}.

## 4 Experimental Results

### 4.1 Parameters.

We have investigated the influence of the parameters on the resulting posterior probability distribution. Figure

3 shows how the choice of the precision  and the kernel width  affect the shape of the posterior probability density. We have used the Gaussian kernel

 K(x,x∗)=exp{−12σ2(x−x∗)2}. (10)

In this figure, 7 data points are shown, which were drawn as , where the mean function is

 f(x)=cos(2x+32π)+sin(6x+32π). (11)

The functions and were chosen as

 K0(x)=1andy0(x)=−12σ20(x−μ0)2, (12)

where the latter corresponds to the logarithm of a Gaussian with mean

and variance

. Choosing a higher value for  leads to sharper updates, while higher values for the kernel width  produce smoother posterior densities. Figure 3: Effect of the change of parameters on the posterior density over the location of the maximizing test point. Panel (a) shows the 7 data points drawn from the noisy function (solid curve). Panel (b) shows the effect of diminishing the precision on the posterior, where solid and shaded curves correspond to ρ=0.2 and ρ=0.1 respectively. Panel (c) shows the effect of increasing the width of the kernel (here, Gaussian). The solid and dotted curves correspond to σ=0.01 and σ=0.1 respectively.

### 4.2 Application to Optimization.

#### Comparison to Gaussian Process UCB.

We have used the model to optimize the same function (11) as in our preliminary tests but with higher additive noise equal to one. This is done by sampling the next test point directly from the posterior density over the optimum location , and then using the resulting pair

to recursively update the model. Essentially, this procedure corresponds to Bayesian control rule/Thompson sampling

[12, 13].

We compared our method against a Gaussian Process optimization method using an upper confidence bound (UCB) criterion . The parameters for the GP-UCB were set to the following values: observation noise and length scale . For the constant that trades off exploration and exploitation we followed Theorem in  which states with . We have implemented our proposed method with a Gaussian kernel as in (10) with width . The prior sufficient statistics are exactly as in (12). The precision parameter was set to .

Simulation results over ten independent runs are summarized in Figure 4. We show the time-averaged observation values of the noisy function evaluated at test locations sampled from the posterior. Qualitatively, both methods show very similar convergence (on average), however our method converges faster and with a slightly higher variance. Figure 4: Observation values obtained by sampling from the posterior over the maximizing argument (left panel) and according to GP-UCB (right panel). The solid blue curve corresponds to the time-averaged function value, averaged over ten runs. The gray area corresponds to the error bounds (1 standard deviation), and the dashed curve in red shows the time-average of a single run.

#### High-Dimensional Problem.

To test our proposed method on a challenging problem, we have designed a non-convex, high-dimensional noisy function with multiple local optima. This Noisy Ripples function is defined as

 f(x)=−11000∥x−μ∥2+cos(23π∥x−μ∥)

where is the location of the global maximum, and where observations have additive Gaussian noise with zero mean and variance . The advantage of this function is that it generalizes well to any number of dimensions of the domain. Figure 5a illustrates the function for the 2-dimensional input domain. This function is difficult to optimize because it requires averaging the noisy observations and smoothing the ridged landscape in order to detect the underlying quadratic form.

We optimized the 50-dimensional version of this function using a Metropolis-Hastings scheme to sample the next test locations from the posterior over the maximizing argument. The Markov chain was started at , executing 120 isotropic Gaussian steps of variance before the point was used as an actual test location. For the arg-max prior, we used a Gaussian kernel with lengthscale , precision factor , prior precision and prior mean estimate . The goal was located at the origin.

The result of one run is presented in Figure 5b. It can be seen that the optimizer manages to quickly ( samples) reach near-optimal performance, overcoming the difficulties associated with the high-dimensionality of the input space and the numerous local optima. Crucial for this success was the choice of a kernel that is wide enough to accurately estimate the mean function. The authors are not aware of any method capable of solving a problem is similar characteristics. Figure 5: a) The Noisy Ripples objective function in 2 dimensions. b) The time-averaged value and the regret obtained by the optimization algorithm on a 50-dimensional version of the Noisy Ripples function.

## 5 Discussion & Conclusions

We have proposed a novel Bayesian approach to model the location of the maximizing test point of a noisy, nonlinear function. This has been achieved by directly constructing a probabilistic model over the input space, thereby bypassing having to model the underlying function space—a much harder problem. In particular, we derived a likelihood function that belongs to the exponential family by assuming a form of symmetry in function space. This in turn, enabled us to state a conjugate prior distribution over the optimal test point.

Our proposed model is computationally very efficient when compared to Gaussian process-based (cubic) or UCB-based models (expensive computation of ). The evaluation time of the posterior density scales quadratically in the size of the data. This is due to the calculation of the effective number of previously seen test locations—the kernel regressor requires linear compuation time. However, during MCMC steps, the effective number of test locations does not need to be updated as long as no new observations arrive.

In practice, one of the main difficulties associated with our proposed method is the choice of the parameters. As in any kernel-based estimation method, choosing the appropriate kernel bandwidth can significantly change the estimate and affect the performance of optimizers that rely on the model. There is no clear rule on how to choose a good bandwidth.

In a future research, it will be interesting to investigate the theoretical properties of the proposed nonparametric model, such as the convergence speed of the estimator and its relation to the extensive literature on active learning and bandits.

## References

• Brochu et al.  E. Brochu, V. Cora, and N. de Freitas.

A tutorial on bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning.

Technical Report TR-2009-023, University of British Columbia, Department of Computer Science, 2009.
• Rawlik et al.  K. Rawlik, M. Toussaint, and S. Vijayakumar. Approximate inference and stochastic optimal control. arXiv:1009.3958, 2010.
• Shapiro  A. Shapiro. Probabilistic Constrained Optimization: Methodology and Applications, chapter Statistical Inference of Stochastic Optimization Problems, pages 282–304. Kluwer Academic Publishers, 2000.
• Kappen et al.  H.J. Kappen, V. Gómez, and M. Opper. Optimal control as a graphical model inference problem. Machine Learning, 87(2):159–182, 2012.
• Kushner and Yin  H.J. Kushner and G.G. Yin. Stochastic Approximation Algorithms and Applications. Springer-Verlag, 1997.
• Mockus  J. Mockus. Application of bayesian approach to numerical methods of global and stochastic optimization. Journal of Global Optimization, 4(4):347–365, 1994.
• Lizotte  D. Lizotte. Practical Bayesian Optimization. Phd thesis, University of Alberta, 2008.
• Jones et al.  D.R. Jones, M. Schonlau, and W.J. Welch. Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 13(4):455–492, 1998.
• Osborne et al.  M.A. Osborne, R. Garnett, and S.J. Roberts. Gaussian processes for global optimization. In 3rd International Conference on Learning and Intelligent Optimization (LION3), 2009.
• Srinivas et al.  N. Srinivas, A. Krause, S. Kakade, and M. Seeger. Gaussian process optimization in the bandit setting: No regret and experimental design. In International Conference on Machine Learning, 2010.
• Hastie et al.  T. Hastie, R. Tbshirani, and J. Friedman. The Elements of Statistical Learning. Springer, second edition, 2009.
• May and Leslie  B.C. May and D.S. Leslie. Simulation studies in optimistic Bayesian sampling in contextual-bandit problems. Technical Report 11:02, Statistics Group, Department of Mathematics, University of Bristol, 2011.
• Ortega and Braun  P.A. Ortega and D.A. Braun. A minimum relative entropy principle for learning and acting.

Journal of Artificial Intelligence Research

, 38:475–511, 2010.

## References

• Brochu et al.  E. Brochu, V. Cora, and N. de Freitas.

A tutorial on bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning.

Technical Report TR-2009-023, University of British Columbia, Department of Computer Science, 2009.
• Rawlik et al.  K. Rawlik, M. Toussaint, and S. Vijayakumar. Approximate inference and stochastic optimal control. arXiv:1009.3958, 2010.
• Shapiro  A. Shapiro. Probabilistic Constrained Optimization: Methodology and Applications, chapter Statistical Inference of Stochastic Optimization Problems, pages 282–304. Kluwer Academic Publishers, 2000.
• Kappen et al.  H.J. Kappen, V. Gómez, and M. Opper. Optimal control as a graphical model inference problem. Machine Learning, 87(2):159–182, 2012.
• Kushner and Yin  H.J. Kushner and G.G. Yin. Stochastic Approximation Algorithms and Applications. Springer-Verlag, 1997.
• Mockus  J. Mockus. Application of bayesian approach to numerical methods of global and stochastic optimization. Journal of Global Optimization, 4(4):347–365, 1994.
• Lizotte  D. Lizotte. Practical Bayesian Optimization. Phd thesis, University of Alberta, 2008.
• Jones et al.  D.R. Jones, M. Schonlau, and W.J. Welch. Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 13(4):455–492, 1998.
• Osborne et al.  M.A. Osborne, R. Garnett, and S.J. Roberts. Gaussian processes for global optimization. In 3rd International Conference on Learning and Intelligent Optimization (LION3), 2009.
• Srinivas et al.  N. Srinivas, A. Krause, S. Kakade, and M. Seeger. Gaussian process optimization in the bandit setting: No regret and experimental design. In International Conference on Machine Learning, 2010.
• Hastie et al.  T. Hastie, R. Tbshirani, and J. Friedman. The Elements of Statistical Learning. Springer, second edition, 2009.
• May and Leslie  B.C. May and D.S. Leslie. Simulation studies in optimistic Bayesian sampling in contextual-bandit problems. Technical Report 11:02, Statistics Group, Department of Mathematics, University of Bristol, 2011.
• Ortega and Braun  P.A. Ortega and D.A. Braun. A minimum relative entropy principle for learning and acting.

Journal of Artificial Intelligence Research

, 38:475–511, 2010.