« Previous « Start » Next »

6  rbfSolve description

Following is a detailed description of the rbfSolve algorithm.

6.1  Summary

The manual considers global optimization of costly objective functions, i.e. the problem of finding the global minimum when there are several local minima and each function value takes considerable CPU time to compute. Such problems often arise in industrial and financial applications, where a function value could be a result of a time-consuming computer simulation or optimization. Derivatives are most often hard to obtain, and the algorithms presented make no use of such information.

The emphasis is on a new method by Gutmann and Powell, A radial basis function method for global optimization. This method is a response surface method, similar to the Efficient Global Optimization (EGO) method of Jones. The TOMLAB implementation of the Radial Basis Function (RBF) method is described in detail.

6.2  Introduction

The task of global optimization is to find the set of parameters x in the feasible region Ω ⊂ R^d for which the objective function f(x) obtains its smallest value. In other words, a point x^* is a global optimizer to f(x) on Ω, if f(x^*) ≤ f(x) for all x Ω . On the other hand, a point x is a local optimizer to f(x), if f( x ) ≤ f(x) for all x in some neighborhood around x. Obviously, when the objective function has several local minima, there could be solutions that are locally optimal but not globally optimal and standard local optimization techniques are likely to get stuck before the global minimum is reached. Therefore, some kind of global search is needed to find the global minimum with some reliability.

Previously a Matlab implementations of the DIRECT [3] has been made, the new constrained DIRECT [12] and the Efficient Global Optimization (EGO) [4] algorithms. The implementations are part of the TOMLAB optimization environment. The implementation of the DIRECT algorithm is further discussed and analyzed in Björkman, Holmstrom [1]. Since the objective functions in our applications often are expensive to compute, we have to focus on very efficient methods. At the IFIP TC7 Conference on System Modelling and Optimization in Cambridge 1999, Hans-Martin Gutmann presented his work on the RBF algorithm [5]. The idea of the RBF algorithm is to use radial basis function interpolation to define a utility function (Powell [15]). The next point, where the original objective function should be evaluated, is determined by optimizing on this utility function. The combination of our need for efficient global optimization software and the interesting ideas of Powell and Gutmann led to the development of an improved RBF algorithm implemented in Matlab.

6.3  The RBF Algorithm

Our RBF algorithm is based on the ideas presented by Gutmann [5, 8], with some extensions and further development. The algorithm is implemented in the Matlab routine rbfSolve.

The RBF algorithm deals with problems of the form where f R and x, x_L, x_U R^d. We assume that no derivative information is available and that each function evaluation is very expensive. For example, the function value could be the result of a time-consuming experiment or computer simulation.

6.3.1  Description of the Algorithm

We now consider the question of choosing the next point where the objective function should be evaluated. The idea of the RBF algorithm is to use radial basis function interpolation and a measure of `bumpiness' of a radial function, σ say. A target value f_n^* is chosen that is an estimate of the global minimum of f. For each y {x₁, … ,x_n} there exists a radial basis function s_y that satisfies the interpolation conditions The next point x_n+1 is calculated as the value of y in the feasible region that minimizes σ(s_y). It turns out that the function y σ(s_y) is much cheaper to compute than the original function.

Here, the radial basis function interpolant s_n has the form with λ₁,…,λ_n R, b R^d, a R and is either cubic with (r)=r³ or the thin plate spline (r)=r² logr. Gutmann considers other choices of and of the additional polynomial in [6], but later in [7] concludes that the situation in the multiquadric and Gaussian cases is disappointing.

The unknown parameters λ_i, b and a are obtained as the solution of the system of linear equations where Φ is the n× n matrix with Φ_ij=(x_i−x_j2) and s_y could be obtained accordingly, but there is no need to do that as one is only interested in σ(s_y). In [13] Powell shows that if the rank of P is d+1, then the matrix is nonsingular and the linear system (4) has a unique solution.

Gutmann defines σ in [8]. For s_n in (3) it is Further, it is shown that σ(s_y) is Thus minimizing σ(s_y) subject to constraints is equivalent to minimizing g_n defined as where μ_n(y) is the coefficient corresponding to y of the Lagrangian function L that satisfies L(x_i)=0, i=1, …, n and L(y)=1. It can be computed as follows. Φ is extended to where (_y)_i=(y−x_i2), i=1, … ,n, and P is extended to Then μ_n(y) is the (n+1)-th component of v R^n+d+2 that solves the system We use the notation 0_n and 0_d+1 for column vectors with all entries equal to zero and with dimension n and (d+1), respectively. The computation of μ_n(y) is done for many different y when minimizing g_n(y). This requires O(n³) operations if not exploiting the structure of Φ_y and P_y. Hence it does not make sense to solve the full system each time. A better alternative is to factorize the interpolation matrix and update the factorization for each y. An algorithm that requires O(n²) operations is described in Section 6.3.3.

When there are large differences between function values, the interpolant has a tendency to oscillate strongly. It might also happen that mins_n(y) is much lower than the best known function value, which leads to a choice of f_n^* that overemphasizes global search. To handle these problems, large function values are in each iteration replaced by the median of all computed function values.

Note that μ_n and g_n are not defined at x₁, … ,x_n and This will cause problems when μ_n is evaluated at a point close to one of the known points. The function h_n(x) defined by is differentiable everywhere on Ω, and is thus a better choice as objective function. Instead of minimizing g_n(y) in (9) one may minimize −h_n(y). In our implementation we instead minimize −log(h_n(y)). By this we avoid a flat minimum and numerical trouble when h_n(y) is very small.

6.3.2  The Choice of f_n^*

For the value of f_n^* it should hold that The case f_n^*=min_{y Ω} s_n(y) is only admissible if min_{y Ω} s_n(y)<s_n(x_i), i=1, …, n. There are two special cases for the choice of f_n^*. In the case when f_n^*=min_{y Ω} s_n(y), then minimizing (9) is equivalent to In the case when f_n^*=−∞, then minimizing (9) is equivalent to So how should f_n^* be chosen? If f_n^*=−∞, then the algorithm will choose the new point in an unexplored region, which is good from a global search point of view, but the objective function will not be exploited at all. If f_n^*=min_{y Ω} s_n(y), the algorithm will show good local behaviour, but the global minimum might be missed. Therefore, there is a need for a mixture of values for f_n^* close to and far away from min_{y Ω} s_n(y). Gutmann describes two different strategies for the choice of f_n^* in [8].

The first strategy, denoted idea 1, is to perform a cycle of length N+1 and choose f_n^* as with where n_init is the number of initial points. Here, N=5 is fixed and max_i f(x_i) is not taken over all points, except for the first step of the cycle. In each of the subsequent steps the n−n_max points with largest function value are removed (not considered) when taking the maximum. Hence the quantity max_i f(x_i) is decreasing until the cycle is over. Then all points are considered again and the cycle starts from the beginning. More formally, if (n−n_init)mod(N+1)=0, n_max=n, otherwise

The second strategy, denoted idea 2, is to consider f_n^* as the optimal value of and then perform a cycle of length N+1 on the choice of α_n. Here, N=3 is fixed and where n₀ is set to n at the beginning of each cycle. For this strategy, max_i f(x_i) is taken over all points in all parts of the cycle.

Consider equation (21). Note that for a fixed y the optimal f^*(y) is the one for which Substituting this equality constraint into the objective of (21) simplifies the problem to the minimization of Denoting the minimizer of (24) by y^*, and choosing f_n^*=f^*(y^*), it is evident that y^* minimizes μ_n(y)[s_n(y)−f_n^*]² and hence g_n(y) in (9).

For both strategies (idea 1 and idea 2), a check is performed when (n−n_init)mod(N+1)=N. This is the stage when a purely local search is performed, so it is important to make sure that the minimizer of s_n is not one of the interpolation points or too close to one. The test used is where f_min is the best function value found so far, i.e. min_i f(x_i), i=1, … ,n. For the first strategy (idea 1), if (25) is true, then otherwise f_n^* is set to 0. For the second strategy (idea 2), if (25) is true, then α_n (or more correctly α_n0+3) is set to otherwise α_n0+3 is set to 0.

6.3.3  Factorizations and Updates

In Powell [14] a factorization algorithm is presented for the solution of (4). The algorithm makes use of the conditional definiteness of Φ, i.e. λ^T Φ λ > 0,  λ ≠ 0 and P^T λ=0. If is the QR decomposition of P, then the columns of Z span the null space of P^T, and every λ with P^T λ=0 can be expressed as λ=Z z for some vector z. Thus the conditional positive definiteness implies that This shows that Z^T Φ Z is positive definite, and thus its Cholesky factorization exists. This property can be used to solve (4) as follows. Consider the interpolation condition Φ λ+Pc=F in (4). Multiply from left by Z^T and replace λ by Z z. Because Z^T P=0, the interpolation condition simplifies to Solving this system using the Cholesky factorization gives z. Then compute λ=Z z and solve for c using the QR decomposition of P as

The same principle can be applied to solve (12) for a given y to get μ_n(y). In analogy to the discussion above, if the extended matrices Φ_y and P_y in (10) and (11), respectively, are given, and if and is the Cholesky factorization, then the vector yields μ_n(y) = v_n+1, where z(y) solves The Cholesky factorization is the most expensive part of this procedure. It requires O(n³) operations. As μ_n(y) must be computed for many different y this is inacceptable. However, if one knows the QR factors of P and the Cholesky factor of Z^T Φ Z, the QR factorization of P_y and the new Cholesky factor L_y can be computed in O(n²) operations.

The new Φ(y) is where (_y)_i=(y−x_i2), i=1, … ,n. The new P(y) is Compute the QR factorization of P_y, defined in (10). Given P=QR, the QR factorization of P_y may be written as where H is an orthogonal matrix obtained by d+1 Givens rotations and for i=d+2, …, n the i-th column of H is the i-th unit vector. Denote B=Q^T Φ Q. Using Φ_y as defined in (10) consider the expanded B matrix Multiplications from the right and left with H affects only the first (d+1) rows and columns and the last row and the last columns of the matrix in the middle. (Remember, d is the dimension of the problem). Hence where * denotes entries not important for the moment. From the form of B_y it follows that holds. The Cholesky factorization of Z^T Φ Z is already known. The new Cholesky L_y factor is found by solving the lower triangular system L l = v for l, computing β=√γ−l^Tl, and setting It is easily seen that L_y L_y^T = Z_y^T Φ_y Z_y because

Note that in practice we do the following: First compute the factorization of P, i.e. P_y=Q_yR_y, using Givens rotations. Then, since we are only interested in v and γ in (42), it is not necessary to compute the matrix B_y in (41). Setting v to the last column in Q_y and computing v=Φ_y^T v=Φ_y v (Φ_y is symmetric), gives v and γ by multiplying the last (n−d) columns in Q_y by v, i.e. Using this algorithm, v and γ are computed using ((n+1) + (n−d)) inner products instead of the two matrix multiplications in (41).

Note that the factorization algorithm is a normal `null-space' method for solving an optimization problem involving linear equality constraints. The system of linear equations in (4) defines the necessary conditions for a stationary point to the unconstrained quadratic programming (QP) problem Viewing c as Lagrange multipliers for the linear equalities in (4), (47) is equivalent to the QP problem in λ defined as The first condition in the conditional positive definiteness definition is the same as saying that the reduced Hessian must be positive definite at the solution of the QP problem if that solution is to be unique.

The type of update procedure described above is suitable each time an optimal point y = x_n+1 is added. However, when evaluating all candidates y an even more efficient algorithm can be formulated. What is needed is a black-box procedure to solve linear systems with a general right-hand side: Using the QR-factorization in (28) the steps simplify when r=0 as in (4), but all steps are useful for solving the extended system (49); see next.

For each of many vectors y, the extended system takes the form where p = ( y^T 1 ). This permutes to which may be solved by block-LU factorization (also known as the Schur-complement method). It helps that most of the right-hand side is zero. The solution is given by the steps Thus, each y requires little more than solving for (,) using the current factorizations (two operations each with Q, R and L). This is cheaper than updating the factors for each y, and should be reliable unless the matrix in (4) is nearly singular. The updating procedure is best numerically, and it is still needed once when the final y = x_n+1 is chosen.

6.3.4  A Compact Algorithm Description

Section 6.3.1-6.3.3 described all the elements of the RBF algorithm as implemented in our Matlab routine rbfSolve, but our discussion has covered several pages. We now summarize everything in a compact step-by-step description. Steps 2, 6 and 7 are different in idea 1 and idea 2.

	idea 1		idea 2
1:	Choose n initial points x1, … ,xn (normally the 2d box corner points defined by the variable bounds). Compute Fi=f(xi), i=1,2, … ,n and set ninit=n.
2:	Start a cycle of length 6.		Start a cycle of length 4.
3:	If the maximum number of function evaluations reached, quit.
4:	Compute the radial basis function interpolant sn by solving the system of linear equations (4).
5:	Solve the minimization problem miny Ω sn(y).
6:	Compute fn* in (18) corresponding to the current position in the cycle.		Compute αn in (22) corresponding to the current position in the cycle.
7:	New point xn+1 is the value of y that minimizes gn(y) in (9).		New point xn+1 is the value of y that minimizes f*(y) in (24).
8:	Compute Fn+1=f(xn+1) and set n=n+1.
9:	If end of cycle, go to 2. Otherwise go to 4.

6.3.5  Some Implementation Details

The first question that arise is how to choose the points x₁, … ,x_ninit to include in the initial set. We only consider box constrained problems, and choose the corners of the box as initial points, i.e. n_init=2^d. Starting with other points is likely to lead to the corners during the iterations anyway. But as Gutmann suggests, having a "good" point beforehand, one can include it in the initial set.

The subproblem is itself a problem which could have more than one local minima. To solve (51) (at least approximately), we start from the interpolation point with the least function value, i.e. argmin f(x_i), i=1, … ,n, and perform a local search. In many cases this leads to the minimum of s_n. Of course, there is no guarantee that it does. We use analytical expressions for the derivatives of s_n and perform the local optimization using ucSolve in TOMLAB [9, 10] running the inverse BFGS algorithm [11].

To minimize g_n(y) for the first strategy, or f^*(y) for the second strategy, we use our Matlab routine glbSolve implementing the DIRECT algorithm (see the TOMLAB manual). We run glbSolve for 500 function evaluations and choose x_n+1 as the best point found by glbSolve. When (n−n_init)mod(N+1)=N (when a purely local search is performed) and the minimizer of s_n is not too close to any of the interpolation points, i.e. (25) is not true, glbSolve is not used to minimize g_n(y) or f^*(y). Instead, we choose the minimizer of (51) as the new point x_n+1. The TOMLAB routine AppRowQR is used to update the QR decomposition.

Our experience so far with the RBF algorithm shows that for the second strategy (idea2), the minimum of (24) is very sensitive for the scaling of the box constraints. To overcome this problem we transform the search space to the unit hypercube. This algorithm improvement is necessary to avoid rank deficiency in the interpolation matrix for the train design problem.

In our implementation it is possible to restart the optimization with the final status of all parameters from the previous run.

« Previous « Start » Next »