National University of Kaohsiung Repository System:Item 310360000Q/10518

全文

(1)國立高雄大學統計學研究所碩士論文. A Bayesian Approach for Optimization Problems in Computer Experiments 利用一個貝氏方法來解決電腦模擬實驗的極值問題. 研究生：謝鐘毅撰指導教授：陳瑞彬. 中華民國九十五年七月.

(2) A Bayesian Approach for Optimization Problems in Computer Experiments. by Jhong-Yi Sie Advisor Ray-Bing Chen. Institute of Statistics National University of Kaohsiung Kaohsiung, Taiwan 811, R.O.C. July 2006.

(3) Contents. Z` zZ`. ii iii. 1 Introduction. 1. 2 Metropolized Matching Pursuit Algorithm. 2. 3 A Noise Basis-based Response Surface Methodology. 4. 4 Numerical Experiments. 7. 5 Some Comparisons with Design and Analysis of Computer Experiments 15 5.1. DACE on grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 16. 5.2. Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 18. 6 Conclusions. 25. References. 27. i.

(4) ¿à×ÍP]°¼Xé\ÿa@ÝÁÂ®Þ ¼0>0Wü}ÿ »ñ{..Ù.@~X .ß Ö. »ñ{..Ù.@~X ` ÍS¡Z&ÆèÝ×ÍÝÕ°à¼Xé\ÿa@ (computer experiments) îÝÁÂ®Þ (optimization problem)3&ÆÝÕ°Ý×Í¥F-Î¿àï ÝÃ9 (basis) ÝaPà)¼¿ÎáÝË@DT`«¬2à×Íh]° (Bayesian approach) ¼Oÿ8ETÝ;ó3h&Æ¿à×°ÿa|C@jÝ@¼Ó¨&ÆÝ± Õ°Ý[t¡f´&ÆÝÕ°CÍÙîOÁÂÝ]° n"ÞMetropolized matching pursuit Õ°ÄÃ9DT`«à°. ii.

(5) A Bayesian Approach for Optimization Problems in Computer Experiments Advisor: Dr. Ray-Bing Chen Institute of Statistics National University of Kaohsiung. Student: Jhong-Yi Sie Institute of Statistics National University of Kaohsiung. ABSTRACT In this thesis, we propose an iterative algorithm for solving the optimization problems in computer experiments. The key point of our algorithm is to approximate the true yet unknown response surface by a linear combination of predefined bases, and a Bayesian approach is applied here for inferring the corresponding coefficients. Several simulations and real experiments are displayed to show the performances of our new algorithm. Finally we compare our algorithm with other statistical optimization methods. Keywords: Metropolized matching pursuit algorithm, overcomplete basis, response surface, surrogate approach. iii.

(6) 1. Introduction The optimization problem in computer experiment is considered here, and this. optimization problem can be represented as maximize. y = f (x),. subject to. x ∈ X,. (1). where y ∈ R is the response variable, f : Rk → R is the response function, x =. (x1 , · · · , xk ) is the vector of k factors, and X ⊂ Rk is the compact experimental region. Typically, here the response function f is assumed to be unknown or very complicated. An intuitive idea for finding optimums is that given a grid of X first, we compute the re-. sponses for all grid points, and then the optimum is found by maximizing these response values. However, the charge for computing all grid points is costly, especially when the experimentation is very expensive or time-consuming. There is another possible approach for solving this kind of optimization problem, which is based on the surrogate models.. In this surrogate approach, the response function f is approximated by an inexpensive surrogate fˆ, and then we optimize fˆ instead of f . One popular statistical optimization method with surrogate models is Response Surface Methodology (RSM), proposed by Box and Wilson (1951). In RSM, the surrogate models are the low-order polynomial models, and the central composite designs (CCD) are used to collecting the data for fitting these surrogate models. The first step of RSM is to use the first-order polynomial model for finding the direction of the maximum improvement in the response. If surface curvature exists, then the first-order polynomial model would be modified by adding higher-order terms into the model. At this time, usually the second-order polynomial model is employed to approximate the response to locate the optimum. Besides RSM, Design and Analysis of Computer Experiment (DACE, Sacks et al., 1989b) is another useful statistical method for optimization by using surrogate approach. In DACE, the surrogate approximation is constructed by an interpolation scheme, kriging method (Matheron, 1963), and the space filling designs, for example, the Latin hypercube designs, are chosen for sampling experimental points. Here kriging method can be understood as linear prediction. Hence the kriging predictor for the unobserved response is a weighted linear combination of all the observed experiment points, and these weight coefficients are estimated by minimizing the Mean Square Error (MSE) of the predictor. Wang and Chen (2004) and Chen et al. (2006) proposed the Basis-based Response Surface Method (BRSM), which is also an optimization algorithm by using surrogate 1.

(7) models. The key idea of BRSM is to consider response surfaces as images and then the P surrogate surfaces are linear combinations of pre-defined bases, i.e. fˆ = ci φi (x), where φi is an element of the basis dictionary defined on the experimental region.. Due to the noiseless model of response, y = f (x), in BRSM, the coefficients of P this surrogate surface, fˆ = ci φi (x), are inferred by a practical image representation technique, matching pursuit algorithm (MP, Mallat and Zhang, 1993, and Bergeaud and P Mallat, 1996). Basically, MP is to directly minimize ky − ci φi (x)k. In fact, we can consider the difference between f and the surrogate response fˆ as a noise ε. Hence in P this thesis, the response y is modelled as a noise model, y = ci φi + ε. Indeed, in image. representation (Lewicki and Olshausen, 1999, Lewicki and Sejnowski, 2000 and Wu et al., 2003) and signal processing (Wolfe et al, 2004), the models of the images (or signals) also can be represented as noise model. Thus, the goal of this thesis is to modify BRSM for this noise model, and the key point of this modification is that we need to have a new method for inferring the unknown coefficients instead of MP. Here, based on the noise model, the Bayesian inference approach is considered, and a Markov Chain Monte Carlo algorithm, Metropolized matching pursuit algorithm (Wu et al., 2003), is adopted for inferring the coefficients. This thesie is organized as follows. In Section 2, Metropolized matching pursuit algorithm is introduced. Then the modification of BRSM for the noise model is proposed in Section 3. In Section 4, the results of experiments for exploring the performances of our new algorithm are presented. In Section 5, our novel algorithm is compared with DACE numerically. Finally, a conclusion is given in Section 6.. 2. Metropolized Matching Pursuit Algorithm. In image (or signal) representation, the model of an image I is usually represented P as a linear superposition of basis dictionary. i.e. I = ci bi , where bi is a basis in the. dictoinary and ci is the corresponding basis coefficient. To infer the coefficients ci , several. iterative algorithms have been proposed, and the matching pursuit algorithm is the most prominent one. Basically, the MP algorithm is an iterative algorithm to infer ci ’s by minimizing the ℓ2 norm of the residual of the image I and the weighted bases. However, P the difference between the image I and the approximation ci bi can be considered as a. noise. Thus, the image I can also be modelled by a noise model. In fact, the noise models for images have been applied in image representation and signal processing. For detail, see Lewicki and Olshausen (1999), Lewicki and Sejnowski (2000), Wu et al. (2003), and. 2.

(8) Wolfe et al. (2004). As mentioned above, the image I is modelled as a noise model: X I= c(x,y,l,θ,e) b(x,y,l,θ,e) + ε,. (2). (x,y,l,θ,e). where ε is a vector of white noise with variance σ 2 , the bases b(x,y,l,θ,e) are selected. from an overcomplete base dictionary, and c(x,y,l,θ,e) is the corresponding base coefficient of b(x,y,l,θ,e) . For each base, (x, y) is the central position of the base on the image domain, and (l, θ, e) is the type of the base. l is the scale or the length of the base, θ is the orientation, and e is the indicator for even/odd bases. Based on this noise model (2), to infer the coefficients ci ’s can be accomplished by a Bayesian framework. Specifically, these coefficients ci ’s are inferred by sampling a posterior probability, P (c(x,y,l,θ,e) | I) ∝ P (I |c(x,y,l,θ,e) )P (c(x,y,l,θ,e) ),. (3). where the prior of the coefficients c(x,y,l,θ,e) is assumed to be independently distributed according to a distribution p(c), which is the mixture of a point mass at 0 and a Gaussian distribution with a large variance, i.e., c(x,y,l,θ,e) ∼ p(c) = ρδ0 + (1 − ρ)N (0, τ 2 ), independently.. (4). In this section, a Bayesian approach for inferring ci ’s, the Metropolized matching pursuit algorithm, is introduced. For simplicity, we use i and j to index (x, y, l, θ, e) and we define zi = 1 if bj is active, i.e., cj 6= 0, and zj = 0 otherwise. Here we also assume. that the parameters (ρ, τ 2 , σ 2 ) are known for our model. Then the Metropolized matching pursuit algorithm iterate the following steps: 1. Randomly select a window W on the image domain, let A be the number of active bases within W , and let B be the number of inactive bases within W . With probability pbirth (a pre-designed number), go to [2]. With probability pdeath = 1 − pbirth , go to [4].. 2. For each inactive base j = (x, y, l, θ, e) with (x, y) ∈ W and zj = 0, compute the P residual image, R = I − i6=j ci bi and γj =. rj2 p p(I |zj = 1; {ci , ∀i 6= j}, I) = exp{ 2 } σ⋆2 /τ 2 , p(I |zj = 0; {ci , ∀i 6= j}, I) 2σ⋆. (5). where rj =< R, bj > /(1+σ 2 /τ 2 ), and σ⋆2 = 1/(1/σ 2 +1/τ 2 ). Then with probability P pdeath j:zj =0; (x,y)∈W (1 − ρ)γj /ρ , paccept = pbirth × (A + 1) go to [3], with probability 1 − paccept go back to [1]. 3.

(9) Figure 1: the observed natural image (left) and the reconstructed image (right).. 3. Among all the inactive bases j with (x, y) ∈ W and zj = 0, sample a base j with. probability proportional to γj , then let zj = 1 and sample cj ∼ N (rj , σ⋆2 ). Go back. to [1].. 4. If A > 0, then randomly select an active base bk with zk = 1 and (x, y) ∈ W . Then. temporarily turn off bk , i.e., set ck = 0, and A → A − 1 temporarily. Then for all. the inactive bases j with zj = 0 and (x, y) ∈ W , including base k, do the same computation as [2] (as if ck = 0), and compute paccept as in [2].. 5. With probability 1/paccept , accept the proposal of deleting the base k, i.e., set ck = 0. Go back to [1]. With probability 1 − 1/paccept , reject the proposal of deleting base k, i.e., recover the original ck . Go back to [1].. Figure 1 shows the result of natural image with size 70 × 70 by Metropolized MP. algorithm. In this example, the dictionary is composed of two-dimensional Gabor func-. tions, and the details of Gabor bases will be described in section 4. Here the parameters (ρ, τ 2 ) of the prior are assumed to be (0.5, 1) , the variance σ 2 of the noise is assumed to be 1 and the pre-designed number pbirth is assumed to be 0.7.. 3. A Noise Basis-based Response Surface Methodology Here, the noise model of response is considered, and a modification of the Basis-. based Response Surface Methodology is proposed for the noise model. This novel algo4.

(10) rithm is called as Noise Basis-based Response Surface Methodology (NBRSM). The response surface y is modelled as y=. X. ci φi + ε.. (6). i. As mentioned in Section 2, the Metropolized matching pursuit algorithm is used here for inferring coefficients in the surrogate model. Just like BRSM, NBRSM is an iterative algorithm for numerical optimization, and the outline of NBRSM is shown in Figure 2. (1). Generate grid over the experimental region. (2). Choose initial experiment points. (3). Generate basis functions. (4). Repeat until the optimal points are found (4.1). Evaluate the response variables. (4.2). Construct the surrogate surface. (4.3). Predict possible optimal points and choose the next experiment points according to the constructed surrogate surface. Figure 2: Schema of the Noise Basis-based Response Surface Method.. The details of NBRSM are introduced as follows. Grid generation. The first step of NBRSM is to choose a grid in the experimental region. Specifically, the k-dimensional experimental region X is discretized into a. grid P = (p1 × p2 × · · · × pN ). These grid points are denoted as x1 , · · · , xN , where N = p1 p2 · · · pk . Therefore the response surface on the grid is defined as SP = {(xi , f (xi ))|xi ∈ P},. (7). and the discretized response surface can be treated as an (multi-dimension) image with finite number of pixels. Initial design. Since the response surface is unknown, a set of grid points is chosen as the initial experimental points, Pinit ⊂ P. Then we compute the corresponding response for the experimental points to construct the surrogate surface. Here, we. choose a space filling designs, the uniform design (Fang, 1980, and Fang et al., 2000), to sample the initial experimental points. Note that we denote the set containing all the experiment points as Pexp and assume Pexp = Pinit initially. 5.

(11) Surrogate approximation. At the tth iteration, the response values f (x)’s are evaluated for all x ∈ Pexp . Thus, we only have the response values over the grid P with respect to Pexp , and let VPexp be the corresponding response vector, i.e., VPexp = (f (x1 ), . . . , f (xN ))T where f (xi ) =. (. f (xi ), if xi ∈ Pexp ,. 0,. if xi ∈ P\Pexp .. As mentioned before, the response vector VPexp is treated as an image and is modelled as (2), i.e.,. VPexp =. M X. ctj φj + ε,. (8). j=1. where ε is a vector of white noises with variance σ 2 , φj ’s are elements of basis dictionary over the whole grid P, and ctj ’s are the corresponding coefficients at the. tth iteration. In addition, in this thesis, the basis dictionaries are considered as overcomplete dictionaries. That is, the number of the elements of the dictionary, M , is larger than the dimensionality of the gird P, i.e., M > N . However, as mentioned in Chen et al. (2006), it may circumvent the possible pitfall. that force the respone variables in P\Pexp to be a constant. Thus, following the suggestion of Chen et al. (2006), we only collect the response values of VPexp that. have been evaluated on the grid P, i.e.,. V˜Pexp = (f (x1 ), . . . , f (xp ))T .. (9). Let exi is the N ×1 unit vector whose values are all zero expect the one corresponding. to the point xi in which the value is assigned to be one, and Ip be the p × N identification matrix that the ith row of Ip is (ex )T . Then V˜Pexp can be represented i. as. V˜Pexp = Ip VPexp , and from (8), we then have V˜Pexp = =. M X. ctj φ˜j + ε˜. j=1. M X. ctj. φ˜j ˜ kφj k + ε˜ kφ˜j k. (11). c˜tj. φ˜j + ε˜, kφ˜j k. (12). j=1. =. (10). M X j=1. 6.

(12) where ε˜ is still a vector of white noises with variance σ 2 , and φ˜j = Ip φj is an element. of basis dictionary corresponding to the experiment points. Here, the coefficients P c˜tj of the current surrogate surface M ˜tj φ˜j /kφ˜j k can be inferred by using the j=1 c. Metropolized matching pursuit algorithm. After obtaining the coefficient c˜tj , we construct the surrogate surface defined over the whole experiment region P as M X j=1. c˜tj φj (x). kIp φj k2. (13). From the above discussion, the algorithm of NBRSM is summerized in Figure 3. (1) Generate a grid P containing N points on the experimental region X .. (2) Choose Ninit initial experiment points by a uniform design and define Pexp = Pinit .. (3) Generate a basis dictionary {φj , j = 1, . . . , M }.. (4) Repeat until the optimal points are found.. (4.1) Evaluate f (x)’s for all x ∈ Pexp . (4.2) Construct V˜Pexp and φ˜j = Ip φj accordingly.. Estimate the coefficients c˜j by the “Metropolized matching pur-. ‘. suit algorithm”. Construct the surrogate surface. c˜j j=1 kIp φj k φi (x).. PM. (4.3) Find new possible optimal point xnew by minimizing the surrogate surface from P\Pexp . Check if f (xnew ) is optimal. Update Pexp = Pexp ∪ xnew .. Figure 3: The Noise Basis-based Response Surface Method Algorithm.. 4. Numerical Experiments In this section, several numerical experiments for illustrating the performances. of NBRSM on the cases of two factors, i.e. k = 2, are described. Before showing these numerical experiments, we introduce the basis dictionary which is used here. In this thesis, our basis dictionary is an overcomplete dictionary, which is widely used in image representation, because the wider range of generating elements allows more flexibility in the image representation and may avoid some numerical difficulties while approximating an image by limited bases. For more details about the necessary of overcomplete dictionary, please see Donoho et al. (2004). For the overcomplete dictionary, 7.

(13) 1. 0.5. 0. 10 8. −0.5 10. 6 8. 6. 4 4. 2. 2. Figure 4: The diagram of 2D Gabor function (upper) and the 8 orientations of the Gabor function for σu = 1 (lower).. a popular dictionary is a set of Gabor functions. The two-dimensional Gabor functions are defined by 1 2πu 1 2 2 exp − (σu u + σv v ) cos +ϕ , g(u, v) = Z 2 λ u = u0 + x1 cos θ − x2 sin θ v = v0 + x1 sin θ − x2 cos θ. (14) (15) (16). where Z is the normalizing constant, (x1 , x2 ) are coordinates of X , u0 , v0 , σu , and σv. are user chosen parameters of a two-dimensional Gaussian window satisfying relations √ √ σv = 2σu and λ = 2πσu , λ and ϕ are parameters of a sinusoidal grating, and θ is the angle between the x1 -axis of the image and the u-axis of the Gabor functions. In Figure 4, the surface of a 2D Gabor functions for σu = 1, θ = 0 and ϕ = 0, and the 8 orientations of the Gabor dictionary for σu = 1 are shown. To illustrate the performance of NBRSM, five experiments are demonstrated below, and the parameters of Gabor dictionary are: the angle θ = {0, π8 , 2π , 3π , 4π , 5π , 6π , 7π }, σu = 8 8 8 8 8 8 10 √ , t 2π. for t = { 81 , 14 , 12 , 1, 2, 4, 8}, and phase ϕ = 0. Thus, our Gabor dictionary contains. M = 8 × 7 × N bases for N grid points. Meanwhile, the parameters in the Metropolized. MP algorithm is set to be (φ, τ 2 , σ 2 ) = (0.5, 4, 0.04) and pbirth = 0.5. 8.

(14) Figure 5: The response surface and its corresponding contours of the objective function defined in problem PB .. Experiment 1. We first consider a problem which is defined as Problem PB : max fB (x, y) = s.t.. 100 , 10((x+1)2 −(y+1)2 )+x2 +4. − 2.5 ≤ x ≤ 0.5 and − 1.5 ≤ y ≤ 1.5.. (17). This response surface is a banana-shaped function with unique maximum, and it has been studied in Balkin and Lin (2000) and Chen et al. (2006). The real response surface of PB and its corresponding contour are shown in (a) and (b) of Figure 5, respectively. To apply NBRSM to solve the problem, we choose the grid set , − 13 , . . . , 0.5} and y ∈ {−1.5, − 68 , − 76 , . . . , 1.5}}. In PB = {(x, y)|x ∈ {−2.5, − 14 6 6. other words, the grid contains 361 (19 × 19) points in the two-dimentional domain. [−2.5, 0.5] × [−1.5, 1.5]. Then the initial experimental points are chosen from the uniform design with 19 levels.. Figure 6 shows the evolution of NBRSM for solving problem PB , and these four sub-figures illustrate the surrogate response surface that contains 19, 23, 27, and 35 experimental points, respectively. Note that, in the figures, the extreme point is marked by △, the initial uniform design points are marked by •, the chosen. experimental points are marked by ◦, and the next experimental point determined by the corresponding surrogate surface is marked by ×.. From Figure 6, it seems that the trend of the true response surface is quickly detected by NBRSM. In fact, at the first iteration, the shape of the true response surface is soon sketched in rough as shown in part (a) and the “hot spot”, or the area which contains the extreme point probably is also identified by NBRSM efficiently. Moreover, as shown in part (b)-(d), most of the new experimental points are chosen 9.

(15) Figure 6: Evolution of the surrogate response surfaces for solving Problem PB by NBRSM. The four surrogate response surfaces containing 19, 23, 27, and 35 experimental points are shown in each sub-figures, respectively.. in the hot spot area, and at the 17th iteration (35 experimental points), the extreme point of the response surface with 361 grid points is located by NBRSM. Experiment 2. In this experiment, we considered a problem which is called Weber’s problem. Weber’s problem is to locate a central facility so that the total cost associated with distribution to several demand centers is minimized. For more details about Weber’s problem, please see (Francis et al., 1992). The cost function in this problem is fW (x, y) =. X i. wi. p. (x − z1i )2 + (y − z2i )2. + sin(0.0035(x2 + y 2 )) + 5 sin(0.03((x + 20)2 + y 2 )).. where w = (2, −4, 2, 1), z1 = (−10, 0, 5, 25), and z2 = (−10, 0, 8, 30). Here the Weber’s problem is transformed as the following optimization problem Problem PW : max fW (x, y), s.t.. − 20 ≤ x ≤ 60, and − 20 ≤ y ≤ 60.. (18). The real response surface of fW and its corresponding contour are shown in (a) and 10.

(16) Figure 7: The response surface and its corresponding contours of the objective function defined in problem PW .. (b) of Figure 7, respectively. To apply NBRSM to solve the problem, we choose the , 40 . . . , 60} and y ∈ {−20, 50 , 40 . . . , 60}}. In grid set PW = {(x, y)|x ∈ {−20, 50 3 3 3 3. other words, the grid contains 625 (25 × 25) points in the two-dimentional domain. [−20, 60] × [−20, 60]. Then the initial experimental points are chosen from the uniform design with 25 levels.. Figure 8 shows the evolution of NBRSM for solving problem PW , and these four sub-figures illustrate the surrogate response surface that contains 25, 27, 30, and 33 experimental points, respectively. From Figure 8, NBRSM can also quickly detect the hot spot area which contains the maximum point, and finally finds the maximum point at 8 iterations (33 experimental points). Experiment 3. Here we consider a problem which is defined as Problem PLI : max fLI (x, y) = (e−x. 2 −y 2. − 1). × max{(x + c1 )2 + (y + c2 )2 , (x + d1 )2 + (y + d2 )2 }, s.t.. − 8.4 ≤ x ≤ 7, and − 8.4 ≤ y ≤ 7,. (19). where (c1 , c2 ) = −(d1 , d2 ) = (30, 40). This problem has been presented in Kolda. et al. (2003), and the true response surface fLI is locally Lipschitz and is strictly differentiable at its maximum point (0, 0). The real response surface of fLI and its corresponding contour are shown in (a) and (b) of Figure 9, respectively. In this problem, the grid is set to be PLI = {(x, y)|x ∈ {−8.4, −7.7, −7, . . . , 7} and y ∈. {−8.4, −7.7, −7, . . . , 7}}. In other words, the grid contains 529 (23 × 23) points. in the two-dimentional domain [−8.4, 7] × [−8.4, 7]. Then the initial experimental points are chosen from the uniform design with 23 levels. 11.

(17) Figure 8: Evolution of the surrogate response surfaces for solving Problem PW by NBRSM. The four surrogate response surfaces containing 25, 27, 30, and 33 experimental points are shown in each of the sub-figures, respectively.. Figure 10 demonstrates the evolution of NBRSM for solving problem PLI , and these four sub-figures illustrate the surrogate response surface that contains 23, 24, 27, and 30 experimental points, respectively. Observing the figure, we see that NBRSM quickly arrive the vicinity of the peak and finds the maximum point at (0,0) after 7 iterations. Also NBRSM captures the trend of the true response surface. Experiment 4. In this case, we consider an artificial problem with two symmetric humps, which has been studied in Chen et al. (2006), and apply NBRSM to find the two local extreme points of the response surface. The problem is defined as follows. Problem PT H : max fT H (x, y) = s.t.. −x4 +4.5x2 +2 , 2 e2y. − 2 ≤ x ≤ 2 and − 2 ≤ y ≤ 2.. (20). As shown in Figure 11, there is a saddle point at (0, 0) and two maximum points at (1.5, 0) and (−1.5, 0) in the true response surface. The two maximum values of this are all equal to 113/16. To apply NBRSM to solve problem PT H , we choose the , − 10 , . . . , 2} and y ∈ {−2, − 11 , − 10 , . . . , 2}}. grid set PT H = {(x, y)|x ∈ {−2, − 11 6 6 6 6 12.

(18) Figure 9: The response surface and its corresponding contours of the objective function defined in problem PLI .. That is, the grid contains 625 (25 × 25) points in the domain [−2, 2] × [−2, 2]. Then the initial experimental points are chosen from the uniform design with 25 levels.. The evolution of NBRSM for solving problem PT H is shown in Figure 12, and the surrogate response surface that contains 25, 26, 47, and 52 experimental points is illustrated in part (a)-(d) of Figure 12, respectively. In part (a) of Figure 12, two hot spot areas are identified at the first iteration. Then, NBRSM successfully locates one maximum at (1.5, 0) after 2 iterations as shown in part (b). In part (c)-(d), NBRSM keeps examining the area for a while and switches to the opposite area that contains another maximum point at (−1.5, 0). Finally, the second maximum is located by NBRSM at 28th iteration (52 experimental points). Observing the searching process, we know that NBRSM not only finds both the optimums, NBRSM also detects the shape of the true response surface. Experiment 5. In this case, we show the performance of NBRSM for solving an oscillatory response surface, and the problem is also adopted in Wang and Chen (2004) and Chen et al. (2006). Here in this problem, a dynamical system modelling an absorptive bistable laser diodes with an electronic-controlled external drive is considered, and in this dynamical system, the rate equations are dNe1 2 2 = Sp1 + mc sin(2π · mf · T ) − α1 Ne1 − (α2 Ne1 + Ne1 + α3 )Np dT −Ne1 , (21) dNe2 2 2 = Sp2 − α1 Ne2 − (α2 Ne2 + Ne2 + α3 )Np − α4 Ne2 , (22) η dT dNp 2 2 = Np [γ1 (α2 Ne1 + Ne1 + α3 ) + γ2 (α2 Ne2 + Ne2 + α3 )] − Np dT η. 13.

(19) Figure 10: Evolution of the surrogate response surfaces for solving Problem PLI by NBRSM. The four surrogate response surfaces containing 23, 24, 27, and 30 experimental points are shown in each sub-figures, respectively.. 2 +ε(γ1 Ne1 + γ2 Ne2 ).. (23). One essential indicator in characterizing a dynamical system is Lyapunov exponents. A positive Lyapunov exponent implies that the system is chaotic for the corresponding parameter seeting. More details regarding the definition and computation of Lyapunov exponents can be found in Parker and Chua (1989). Here we intend to find at least one set of the parameters Sp1 and mc in Eq. (21), such that the associated Lyapunov exponent is positive. However, there are two difficulties for finding suitable these parameters. First, it is extremely time consuming for computing Lyapunov exponents. Secondly, the relationship between the adjustable parameters and the resulting Lyapunov exponents are exceeding complicated. Therefore, an efficient numerical scheme for solving the problem is needed. The problem for searching positive Lyapunov exponents is considered as an optimization problem, and as mentioned in Chen et al. (2006), the goal is the find at least on effective point such that the corresponding Lyapunov exponent is positive. Now the problem is shown as follows. Problem PL : max fL (Sp1 , mc ), 14.

(20) Figure 11: The response surface and its corresponding contours of the objective function defined in problem PT H .. s.t. 20 ≤ Sp1 ≤ 30, and 5 ≤ mc ≤ 15,. (24). where fL (Sp1 , mc ) is the value of Lyapunov exponent at (Sp1 , mc ). To apply NBRSM to solve problem PL , we first discretize the experimental region into the grid PL =. {(Sp1 , mc )|Sp1 ∈ {20, 20.5, 21, . . . , 30} and y ∈ {5, 5.5, 6, . . . , 15}}, which contains. 441 (21 × 21) points in the domain [20, 30] × [5, 15]. Then the initial experimental. points are chosen from the uniform design with 21 levels. Figure 13 shows the true response surface and its projection over the 21 × 21 grid of problem PL .. Figure 14 shows the surrogate response surfaces containing 21, 24, 59, and 101 experimental points. Observing the searching process, we see that NBRSM successfully found all 21 positive Lyapunov exponents. To be precise, expect the three positive Lyapunov are found by the initial points, NBRSM identifies the other 18 experimental points by using 24, 30, 33, 36, 38, 42, 43, 47, 49, 59, 60, 75, 76, 77, 97, 99, 100, and 101 experimental points. Besides, NBRSM also detects the trend of the true response surface, and keeps searching in the right quarter of the domain which contains most positive Lyapunov exponents possibly.. 5. Some Comparisons with Design and Analysis of Computer Experiments In this section, first we briefly introduce a popular statistical optimization method,. Design and Analysis of Computer Experiments (DACE), and then we compare the performance between DACE and NBRSM by using above five problems. 15.

(21) Figure 12: Evolution of the surrogate response surfaces for solving Problem PT H by NBRSM. The four surrogate response surfaces containing 25, 26, 47, and 52 experimental points are shown in each sub-figures, respectively.. 5.1. DACE on grid In DACE, the response model is typically expressed as y(x) = β T F (x) + Z(x),. (25). where F (x) is a q × 1 vector of the q chosen regression functions, β is a q × 1 vector of. the regression coefficients, and Z(x) is a Gaussian stochastic process with zero mean and covariance function of the form c(s, t) = σ 2 rθ (s, t). Note that σ 2 > 0 is unknown and rθ is an correlation function, for example, Gauss correlation model, i.e., rθ (s, t) = exp(−θks − tk2 ). Further, in DACE, a space filling designs, for example, the Latin hypercube designs, is chosen for sampling a set of experimental points {xi |i = 1, · · · , p}. Let Ys = (y1 , · · · , yp )T. be the corresponding response vector, and Fs = (F (x1 ), · · · , F (xp )), where yi is the. response for xi . Here we also define R as the matrix of stochastic process correlations. between the experimental points, i.e., R = [rθ (xsi , xsj )], and define r(x) as the vector of correlations between experimental points and an untried point x, i.e., r(x) = [rθ (xsi , xu )], for n = 1, · · · , k. Then the surrogate response surfaces is constructed using the kriging 16.

(22) (a). (b) 5. 50 0. 40. 40. 20. 20. 0 −50 −20 5 0 5 −5 −5. 0. 0 0. −20. −40. −40. −60. −60. −80. −5 −5. −80 0. 5. Figure 13: The three-dimensional true response surfaces and two-dimensional projection over the 21 × 21 grid of the objective functions defined in problem PL are shown in each. sub-figure, respectively.. method. Here the unknown parameters (β, σ 2 , θ) is estimated by the method of maximum likelihood, and the surrogate response surface is constructed by computing ˆ yˆ(x) = βˆT F (x) + r(x)R−1 (YsT − FsT β).. (26). Hence, the kriging method can be considered as an exact interpolation. That is, the estimated values at experimental points are exactly equal to the response values. Schonlau et al. (1998) noted that the general framework is proceeded as follows: 1. Choose a small initial experimental design spread over the entire input space. Run the computer code at these points. 2. Use all previous function evaluations to fit a statistical model for the objective function. 3. Based on the fitted model, find the “most promising” point in the input space for the next run. 4. Compute a stopping criterion. If it is met, then stop. 5. Run the computer code at the selected point in the input space. Go to Step 2. According to above framework, an optimization algorithm on grid is proposed by applying DACE for surrogate surface construction, and this modified algorithm is called G-DACE. The first step of G-DACE is to discretize the experimental region as a grid. A uniform design is chosen for drawing the initial experimental points, and then the corresponding responses are evaluated. Based on these responses, DACE could be used for constructing the surrogate surfaces. Therefore, exclude the model assumption, the difference between 17.

(23) Figure 14: Evolution of the surrogate response surfaces for solving Problem PL by NBRSM. The four surrogate response surfaces containing 21, 24, 59, and 101 experimental points are shown in each sub-figures, respectively.. G-DACE and NBRSM is the surrogate surface construction. The algorithm of G-DACE for optimization problem is shown in Figure 15. (1) Generate a grid P containing N points on the experimental region X .. (2) Choose Ninit initial experiment points by a uniform design and define Pexp = Pinit .. (3) Repeat until the optimal points are found. (3.1) Evaluate f (x)’s for all x ∈ Pexp .. (3.2) Construct the surrogate surface VˆPexp by the kriging method.. (3.3) Find new possible optimal point xnew by minimizing the surrogate surface from P\Pexp . Check if f (xnew ) is optimal. Update Pexp = Pexp ∪ xnew .. Figure 15: The G-DACE Algorithm.. 5.2. Comparison In this subsection, we would show the results of G-DACE for solving the five. problems in Section 4, and then compare these results with those of NBRSM. Here in 18.

(24) Figure 16: Evolution of the surrogate response surfaces for solving Problem PB by GDACE. The four surrogate response surfaces containing 19, 20, 22, and 23 experimental points are shown in each sub-figures, respectively.. G-DACE, the function of the kriging method is from a Matlab toolbox which is available at the web-site http://www2.imm.dtu.dk/∼hbn/dace. For more details about this Matlab toolbox, please see Lophaven et al. (2002). In our experiments, the regression function F (·) is chosen as the second-order polynomial model, and the correlation function r(·, ·) is Gaussian function. The results of these five problem are shown in the following:. Problem PB : Figure 16 shows the evolution of DACE for solving problem PB , and these four sub-figures illustrate the surrogate surface that contains 19, 20, 22 and 23 experimental points, respectively. In Figure 16, G-DACE finds the possible area which contains optimal point and searches around the area. After 5 iterations (23 experimental points), G-DACE successfully identify the optimal point at (0, 0). By comparing Figure 6 and 16, it seems that the trend of surface detects by NBRSM are is more similar to true surface than G-DACE does. Problem PW : The evolution of G-DACE for solving problem PW is shown in Figure 17, and the surrogate surface that contains 25, 54, 60, and 63 experimental points is shown in each sub-figure, respectively. From Figure 17, though G-DACE detects 19.

(25) Figure 17: Evolution of the surrogate response surfaces for solving Problem PW by GDACE. The four surrogate response surfaces containing 25, 54, 60, and 63 experimental points are shown in each sub-figures, respectively.. the trend of the response surface roughly, G-DACE keeps searching in the left side of the surface. After 39 iterations, G-DACE finds the optimal point. Problem PLI : In Figure 18, the searching process of G-DACE for solving problem PLI shows that G-DACE detects the trend quickly, and uses only 26 experimental points to identify the optimal point at (0, 0). The surrogate surface that contains 23, 24, 25 and 26 experimental points is illustrated in part (a)-(d) of Figure 18, respectively. Problem PT H : From Figure 19, we see that the trend of true response surface is detected by G-DACE at the first iteration, and then one optimal point at (1.5, 0) is successfully identified after 2 iterations. Then, G-DACE keeps searching around the first optimal point, and G-DACE finally finds the second optimal point at (−1.5, 0) after 11 iterations. Problem PL : Figure 20 shows the evolution of G-DACE for solving problem PL , and these four sub-figures illustrate the surrogate surface that contains 21, 50, 150 and 439 experimental points, respectively. In Figure 20, G-DACE first searches in the right quarter of the surface, but after 30 iterations, G-DACE identifies 11 optimal points and starts to search in the median of the surface as shown in part (b) and 20.

(26) Figure 18: Evolution of the surrogate response surfaces for solving Problem PLI by GDACE. The four surrogate response surfaces containing 23, 24, 25, and 26 experimental points are shown in each sub-figures, respectively.. (c). In part (d), we see that almost all the grid points are evaluated, and then the last optimal point is just identified after 419 iterations. To be precise, besides three points are found by 21 initial points, G-DACE identified the other 18 positive Lyapunov exponents when 26, 28, 29, 32, 34, 40, 42, 47, 60, 61, 261, 340, 341, 355, 413, 436, 437, and 440 experimental points are used. By observing above results, except the problem PL , it seems that G-DACE also have well performances for solving these problems. Here we compare the performances of NBRSM and G-DACE from these five problems. Since NBRSM is a Bayesian approach, we repeat 1000 times for each problem, and let S be a 1 × 1000 vector of the number. of experimental points for identifying the optimal point(s). Then from S, we obtain the sample mean, sample variance, and sample standard deviation. Beside, we compute the 2.5th and 97.5th sample percentile for NBRSM. Here the (100p)th sample percentile of n samples is the (np)th order statistic. In addition, we define the efficient probability of G-DACE via NBRSM as p = P (X ≤ c),. where X is the number of experimental points for identifying the optimal point(s) by 21.

(27) Figure 19: Evolution of the surrogate response surfaces for solving Problem PT H by GDACE. The four surrogate response surfaces containing 25, 26, 34, and 35 experimental points are shown in each sub-figures, respectively.. NBRSM, c is the number of experimental points for identifying the optimal point(s) by G-DACE. Here, the efficient probability p of G-DACE via NBRSM can be considered as that the NBRSM explores less points to identify the extremes than those used in GDACE with probability p. The performances for these problems are summarized in the following tables. In Table 1, besides the performances of NBRSM and G-DACE, the BRSM algorithm in Chen et al. (2006) is also applied here to solve these problems, but for problem PLI , we add a step of basis selection into BRSM algorithm for getting the better performance. In NBRSM, since the number of explored experimental points for problem PLI is always 31, there is no sample variance in this problem. Hence, we denote the sample variance, standard deviation, and sample percentile for PI to be N/A. Moreover, for PT H , there are two extreme points in the response surface. Thus we denote PT H : 1 Max. to be the problem for identifying one optimal point, and PT H : 2 Max. to be the problem for identifying two optimal points in problem PT H . From Table 1, we see that all these three methods can efficiently find the optimal point for each problem. For NBRSM and BRSM, we find that the results of these two algorithms are similar. The reason might be that the difference of these two algorithms is only the method for inferring coefficients in the 22.

(28) Figure 20: Evolution of the surrogate response surfaces for solving Problem PL by GDACE. The four surrogate response surfaces containing 21, 50, 150, and 439 experimental points are shown in each sub-figures, respectively.. Table 1: Performance summary of NBRSM, G-DACE and BRSM for PB , PW , PLI and PT H .. NBRSM. Problem. PB. PW. PLI. Grid points. 361. 625. 529. 625. 625. sample mean. 30.0820. 42.1820. 31. 26.158. 52.561. sample standard deviation. 6.4078. 7.3851. N/A. 0.8437. 4.3863. 22. 31. N/A. 26. 43. 44. 60. N/A. 29. 60. 24. 64. 27. 27. 36. efficient probability. 0.2. 0.982. N/A. 0.96. 0.006. explored point. 30. 39. 29. 29. 46. 0.64. 0.344. N/A. 0.98. 0.068. th. 2.5. sample percentile. 97.5th sample percentile G-DACE explored point BRSM. efficient probability. 23. PT H : 1 Max. PT H : 2 Max..

(29) surrogate models. Furthermore, according to efficient probabilities, except the problems PW and PT H : 1 Max., it seems that G-DACE is quicker than NBRSM and BRSM to identify the optimal point. For example, in problem PB , the efficient probability is equal to 0.2, and it means that G-DACE use less points to find the optimal point than those of NBRSM with probability 0.8. However, in the problem PW , we see that the efficient probability is equal to 0.982, and it can be considered as that NBRSM use less points to identify the optimal point than G-DACE with probability 0.982. On the other hand, by observing the sample mean of NBRSM and the explored point of G-DACE in problem PB , we find that the efficient probability of G-DACE is smaller than NBRSM, but the difference of explored point between NBRSM and G-DACE is equal to 6. It is because that in the 1000 repeated experiments, the results of explored points may contain 29 or 30 mostly. However, in problem PW , the sample mean of NBRSM is only 42.1820, but the explored point of G-DACE is 64, and the difference is about 22. Thus, it seems that the performances of NBRSM is better than G-DACE. Here we show the response surface of fB and fW in Figure 21. From this figure, it seems that fW is an oscillatory surface, and fB is more smoother than fW . From Table 1, we think that the reason for such kinds of results is caused by the surrogate construction and the smoothness of response surfaces. Since the kriging method is an exact interpolation, then the estimated values at experimental points are exactly equal to the response values. Thus, while the response surfaces are oscillatory or not smooth and the experimental points are fewer, the interpolation may be improper. In the real problem PL , the corresponding response surface is not smooth. In our experiments of G-DACE in Section 5, we find that G-DACE uses a lot of experimental points to search the median region and explores almost all grid points to find all the optimal points. Table 2 compares the performance of NBRSM and G-DACE for problem PL . In this table, we see that the efficient probabilities of G-DACE via NBRSM are small,. Figure 21: The response surfaces of fB and fW .. 24.

(30) Table 2: Performance summary of NBRSM, G-DACE and BRSM for PL . positive L.E.. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. NBRSM. 21. 21. 21. 25. 29.5264. 32.336. 40.7615. 46.0046. 49.6468. 54.2993. 58.3463. G-DACE. 21. 21. 21. 26. 28. 29. 32. 34. 40. 42. 47. efficient prob.. N/A. N/A. N/A. N/A. 0.2924. 0.2844. 0.0034. 0. 0.0378. 0.0333. 0.0528. BRSM. 21. 21. 21. 22. 31. 36. 37. 42. 46. 49. 56. efficient prob.. N/A. N/A. N/A. N/A. 1. 0.9920. 0.3326. 0.3647. 0.3739. 0.3268. 0.4186. positive L.E.. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. NBRSM. 62.0505. 66.0642. 70.4381. 75.1663. 80.0252. 86.2741. 91.9839. 99.2317. 108.914. 116.9633 440. G-DACE. 60. 61. 261. 340. 341. 355. 413. 436. 437. efficient prob.. 0.4232. 0.3073. 1. 1. 1. 1. 1. 1. 1. 1. BRSM. 65. 69. 70. 76. 84. 93. 100. 120. 122. 139. efficient prob.. 0.6560. 0.6468. 0.5103. 0.5723. 0.7076. 0.8119. 0.8291. 0.9702. 0.8578. 0.9461. when we only need to identify few positive Lyapunov exponents. It means at this time, G-DACE explores less points than NBRSM with higher probability. However, in order to find more than 14 positive Lyapunov exponents, NBRSM seems have better performance than G-DACE. Thus, we believe that NBRSM is enable to deal the cases of unsmooth surface.. 6. Conclusions In this thesis, we propose a Bayesian algorithm, Noise Basis-based Response. Surface Methodology (NBRSM), for optimization problems in computer experiments. Basically, NBRSM is the modification of Basis-based Response Surface Methodology (BRSM) for noise model. Like BRSM, the key point of NBRSM is to consider the response surface as an image, and to approximate the true yet unknown response surface by a linear combination of overcomplete Gabor bases. To infer the coefficients in the surrogate model, we apply a Bayesian approach, Metropolized matching pursuit algorithm. Then based on the is surrogate surface, we find the possible optimal points. Five optimization problems are illustrated here for showing the performances of NBRSM. For the smooth surfaces, we see that NBRSM quickly sketches the trend of the unknown response surfaces, and efficiently identifies the hot spot areas, which may contains the possible optimal point. For the unsmooth surface, NBRSM also quickly finds the hot spots, and then NBRSM still successfully identify all extremes. For the real problem PL , NBRSM explores 101 points of 441 grid points to find the all the positive Lyapunov exponents. In Section 5, we also show the results of these five problems with G-DACE. It seems that G-DACE does not work well on unsmooth surfaces. According 25.

(31) to these experiments, we believe that NBRSM is enable efficiently to find the multiple extremes, no matter the true surface is smooth or not. A weakness of NBRSM is that the optimal points identified by NBRSM may not be the “exact” optimal point, because NBRSM need to discretize the whole experimental region first. But due to the results of our experiments, it seems that NBRSM can quickly identify the trend of the unknown response surface and the hot spot areas. Thus, in order to search the exact optimal point, we would suggest to run our NBRSM first for identifying the hot spots for possible extremes, and then the iterative optimization method, for example, Response Surface Methodology (RSM) or Pattern Search (Torczon, 1997), can be applied to search the exact optimal point.. 26.

(32) References [1] Balkin, S. D. and Lin, D. K. J. (2000). Aneural Network Approach to Response Surface Methodology. Commun. Statist.-thoery Meth., 29, 2215-2227. [2] Bergeaud, F. and Mallat, S. (1995). Matching Pursuit of Images. Proc. IEEE ICIP, 1: 53 - 56. [3] Box, G. E. P. and Wilson, K. B. (1951). On the Experimental Attainment of Optimum Condition. Journal of the Royal Statistical Socienty, Ser. B., 13: 1 - 45. [4] Chen, R.-B., Wang, W. and Tsai, F. (2006). Basis-based Response Surface Method in Computer Experiments Optimization. [5] Donoho, D. L. Elad, M. and Temlyakov, V. (2004). Stable Recovery of Sparse Overcomplete Representations in the Presence of Noise. to appear in the IEEE Trans. On Information Theory. [6] Fang, K. T. (1980). The Uniform Design: Application of Number-theoretic Methods in Experimental Design. Acta Math. Appl. Sinica, 3: 363-372. [7] Fang, K. T., Lin, D. K. J., Winker, P., and Zhang, Y. (2000). Uniform Design: Theory and Applicatios. Technometrica, 42: 237-248. [8] Francis, R. L., McGinnis Jr. L. F. and White J. A. (1992). Facility Layout and Location: An Analytical Approach. Prentice-Hall. [9] Kolda, T. G., Lewis, R. M. and Torczon, V. (2003). Optimization by Direct Search: New Perspectives on Some Classical and Modern Methods. SIAM Review, 45(3): 385-482. [10] Lewicki, M. S. and Olshausen, B. A. (1999). Probabilistic Framework for the Adaptation and Comparison of Image Codes. Journal of the Optical Society of America, 16(7): 1587 - 1601. [11] Lewicki, M. S. and Sejnowski, T. J. (2000) Learning Overcomplete Representations. Neural Computation, 12,37 - 365. [12] Lophaven, S. N., Nielsen, H. B. and Søndergaard, J. (2002). DACE - A Matlab Kriging Toolbox. Technical Report IMM-REP-2002-12.. 27.

(33) [13] Mallat, S. and Zhang, Z. (1993). Matching Pursuit with Time-frequency Dictionaries. IEEE Transaction on Signal Processing, 41, 3397 - 3415. [14] Matheron, G. (1963). Principles of Geostatistics. Economic Geology 58: 1246 - 1266. [15] Parker T. S. and Chua L. O. (1989). Practical Numerical Algorithms for Chaotic Systems. Springer-Verlag. [16] Sacks, J., Welch, W. J., Mitchell, T. J. and Wynn, H. P. (1989b). Design and Analysis of Computer Experiments. Statistical Science, 4(4): 409 - 435. [17] Schonlau, M., Welch, W. and Jones, D. (1998). Global Versus Local Search in Constrained Optimization of Computer Models, in New Developments and Applications in Experimental Design, N. Flournoy, W. F. Rosenberger, and W. K. Wong (editors), Institute of Mathematical Statistics, Hayward, California, 34, 11-25. [18] Torczon, V. (1997). On the Convergence of Pattern Search Algorithms. SIAM J. Optimization, 7: 1-25. [19] Wang, W. and Chen, R.-B. (2004). Basis Representation Methodology for Response Surface, Technical report, Department of Applied Mathematics and Institute of Statistics, National University of Kaohsiung. [20] Wolfe, P. J., Godsill, S. J. and Ng, W. J. (2004). Bayesian Variable Selection and Regularisation for Time-frequency Surface Estimation. Journal of the Royal Statistical Society, Series B, 66 part3, 575 - 589. [21] Wu, Y. N., Zhu, S. C. and Guo, C. (2003). Statitical Modeling of Texture Sketch. Proc. of 7th European Conf. on Computer Vision, 240-254.. 28.

(34)