CGO rbfSolve description
This page is part of the CGO Manual. See CGO Manual. |
{{#switch: | left =
{{#switch:{{#if: | {{{smallimage}}} | }} | none =| #default =
}} {{#if:{{#if: | {{{smallimageright}}} | }} | {{#ifeq:{{#if: | {{{smallimageright}}} | }}|none | | }} }}{{
#switch:left | left =| #default = }} {{#if:{{#if: | {{{smallimage}}} | }} | {{#if: | {{{smallimage}}} | }} | [[File:{{#switch:caution | critical = Ambox speedy deletion.png | important = Ambox deletion.png | warning = Ambox content.png | caution = Cleanup.png | move = Ambox move.png | protection = Ambox protection.png | notice | #default = Cleanup.png }} | {{#switch:left | left = 20x20px | #default = 40x40px }} |link=|alt=]] }}{{#switch:left | left =| #default = | {{#if:
| {{{smalltext}}} | Cleanup is needed}} | {{#switch:left
| left = {{#if: | {{{smallimageright}}} | }}| #default = {{#if:
}}| {{{smallimageright}}} |}} |
| #default =
{{#switch: | none =| #default =
}}{{#if: | {{#ifeq:|none
|| }} }}
{{
#switch: | left =| #default = }} {{#if: | | [[File:{{#switch:caution | critical = Ambox speedy deletion.png | important = Ambox deletion.png | warning = Ambox content.png | caution = Cleanup.png | move = Ambox move.png | protection = Ambox protection.png | notice | #default = Cleanup.png }} | {{#switch: | left = 20x20px | #default = 40x40px }} |link=|alt=]] }}{{#switch: | left =| #default = | Cleanup is needed Cross-referencing of the equations. |
{{#switch:
| left =| #default = |
}} Following is a detailed description of the rbfSolve algorithm.
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.
Introduction
The task of global optimization is to find the set of parameters x in the feasible region for which the objective function f(x) obtains its smallest value. In other words, a point is a global optimizer to f (x) on , if . On the other hand, a point 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 has been made, the new constrained DIRECT and the Efficient Global Optimization (EGO) algorithms. The implementations are part of the TOMLAB optimization environment. The implementation of the DIRECT algorithm is further discussed and analyzed in Bjorkman, Holmström. 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. The idea of the RBF algorithm is to use radial basis function interpolation to define a utility function (Powell). 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.
The RBF Algorithm
Our RBF algorithm is based on the ideas presented by Gutmann, 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 and . 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.
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 is chosen that is an estimate of the global minimum of f . For each there exists a radial basis function that satisfies the interpolation conditions
The next point x_{n+1} is calculated as the value of y in the feasible region that minimizes . It turns out that the function is much cheaper to compute than the original function.
Here, the radial basis function interpolant s_{n} has the form
with , , and is either cubic with or the thin plate spline . Gutmann considers other choices of and of the additional polynomial, but later concludes that the situation in the multiquadric and Gaussian cases is disappointing.
The unknown parameters , b and a are obtained as the solution of the system of linear equations
where is the n × n matrix with and
could be obtained accordingly, but there is no need to do that as one is only interested in . 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.
For s_{n} it is
Further, it is shown that is
.
Thus minimizing subject to constraints is equivalent to minimizing g_{n} defined as
where is the coefficient corresponding to y of the Lagrangian function L that satisfies , and L(y) = 1. It can be computed as follows. F is extended to
where , and P is extended to
Then is the (n + 1)-th component of that solves the system
We use the notation 0_{n} and for column vectors with all entries equal to zero and with dimension n and (d + 1), respectively. The computation of is done for many different y when minimizing . This requires operations if not exploiting the structure of and . 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 operations is described in #Factorizations and Updates.
When there are large differences between function values, the interpolant has a tendency to oscillate strongly. It might also happen that is much lower than the best known function value, which leads to a choice of 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 and are not defined at and
This will cause problems when is evaluated at a point close to one of the known points. The function defined by
is differentiable everywhere on , and is thus a better choice as objective function. Instead of minimizing one may minimize . In our implementation we instead minimize . By this we avoid a flat minimum and numerical trouble when is very small.
The Choice of f *
For the value of it should hold that
The case is only admissible if , <mah>i=1, \ldots, n</math>. There are two special cases for the choice of . In the case when , then minimizing is equivalent to
In the case when , then minimizing is equivalent to
So how should be chosen? If , 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 , 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 close to and far away from . Gutmann describes two different strategies for the choice of .
The first strategy, denoted idea 1, is to perform a cycle of length N + 1 and choose as
with
where is the number of initial points. Here, N = 5 is fixed and is not taken over all points, except for the first step of the cycle. In each of the subsequent steps the points with largest function value are removed (not considered) when taking the maximum. Hence the quantity is decreasing until the cycle is over. Then all points are considered again and the cycle starts from the beginning. More formally, if , , otherwise
The second strategy, denoted idea 2, is to consider as the optimal value of
and then perform a cycle of length N + 1 on the choice of an . Here, N = 3 is fixed and
where is set to n at the beginning of each cycle. For this strategy, is taken over all points in all parts of the cycle.
Note that for a fixed y the optimal is the one for which
Substituting this equality constraint into the objective simplifies the problem to the minimization of
Denoting the minimizer by , and choosing , it is evident that minimizes and hence .
For both strategies (idea 1 and idea 2), a check is performed when . This is the stage when a purely local search is performed, so it is important to make sure that the minimizer of sn is not one of the interpolation points or too close to one. The test used is
where is the best function value found so far, i.e. , . For the first strategy (idea 1), then
otherwise is set to 0. For the second strategy (idea 2), then an (or more correctly ) is set
otherwise is set to 0.
Factorizations and Updates
In Powell, a factorization algorithm is presented for the solution. The algorithm makes use of the conditional definiteness of , i.e. and . If
is the QR decomposition of P , then the columns of Z span the null space of , and every with can be expressed as for some vector z. Thus the conditional positive definiteness implies that
This shows that is positive definite, and thus its Cholesky factorization
exists. This property can be used to solve (4) as follows. Consider the interpolation condition in (4). Multiply from left by and replace by . Because , the interpolation condition simplifies to
Solving this system using the Cholesky factorization gives z. Then compute 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 . In analogy to the discussion above, if the extended matrices and in (10) and (11), respectively, are given, and if
and
is the Cholesky factorization, then the vector
yields , where z(y) solves
The Cholesky factorization is the most expensive part of this procedure. It requires operations. As 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 FZ , the QR factorization of P_{y} and the new Cholesky factor L_{y} can be computed in operations.
The new is
where , . The new P (y) is
Compute the QR factorization of P_{y}, defined in (10). Given , the QR factorization of may be written as
where is an orthogonal matrix obtained by d + 1 Givens rotations and for the i-th column of H is the i-th unit vector. Denote . Using 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 is already known. The new Cholesky factor is found by solving the lower triangular system for l, computing , and setting
It is easily seen that because
Note that in practice we do the following: First compute the factorization of P , i.e. , 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 ( is symmetric), gives v and by multiplying the last (n - d) columns in Q_{y} by , 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 = xn+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 . 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 is chosen.
A Compact Algorithm Description
Section #Description of the Algorithm-#Factorizations and Updates 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 evalua- tions reached, quit. | |
4: | Compute the radial basis function interpolant sn by solving the system of linear equations (4). | |
5: | Solve the minimization problem min sn (y). y?? | |
6: | Compute f * in (18) corresponding to the current position in the cycle. | Compute an in (22) corresponding to the current position in the cycle. |
7: | New point xn+1 is the value of y that mini- mizes gn (y) in (9). | New point x_{n+1} is the value of y that min- imizes 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. |
Some Implementation Details
The first question that arise is how to choose the points to include in the initial set. We only consider box constrained problems, and choose the corners of the box as initial points, i.e. . 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. , , 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 sn and perform the local optimization using ucSolve in TOMLAB running the inverse BFGS algorithm.
To minimize for the first strategy, or 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 (when a purely local search is performed) and the minimizer of sn is not too close to any of the interpolation points, i.e. (25) is not true, glbSolve is not used to minimize or . 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.