# Equal Angle Distribution of Polling Directions in Direct-Search Methods

Authors: Benjamin Van Dyke
Journal: Journal of Optimization
 (2014)

---

## Abstract

The purpose of this paper is twofold: first, to introduce deterministic strategies for directional direct-search methods, including new instances of the mesh adaptive direct-search (Mads) and the generating set search (Gss) class of algorithms, which utilize a nice distribution of Poll directions when compared to other strategies, and second, to introduce variants of each algorithm which utilize a minimal positive basis at each step. The strategies base their Poll directions on the use of the QR decomposition to obtain an orthogonal set of directions or on using the equal angular directions from a regular simplex centered at the origin with vertices on the unit sphere. Test results are presented on a set of smooth, nonsmooth, unconstrained, and constrained problems that give comparisons between the various implementations of these directional direct-search methods.


---

## Body

## 1. Introduction

In [1], Vicente and Custódio introduced a framework for directional direct-search methods (DSM) encompassing both the Mads [2] class of algorithms (using an underlying mesh and simple decrease) and the Gss [3] class of algorithms (using sufficient decrease under a forcing function ρ) for black-box optimization problems of the form
(1)min⁡x∈Ωf(x),
where the objective function f:Ω⊆Rn→R⋃‍{∞} is nonsmooth, possibly discontinuous and Ω is the set of feasible points. The algorithms under this framework are shown to have first order convergence results in the Rockafellar sense [4] for discontinuous functions. The framework is shown in Algorithm 1.Algorithm 1:The directional direct-search method of [1].
(0)Initialization
x0∈Ω,α0>0,0<β1≤β2<1,γ≥1
letD be a (possibly infinite) set of positive spanning sets
(1)Iteration  k
Search
evaluatef on a finite set of points Sk (on a mesh if  ρ-=0)
Poll (optional if Search is successful)
ChooseDk∈D, order the set of Poll points
Pk={xk+αkd:d∈Dk}, and evaluate f on Pk
(2)Updates
if the iteration is successful (i.e. there exists xs∈Sk or xp∈Pk
such thatf(xs)<f(xk)-ρ-(αk) or f(xp)<f(xk)-ρ-(αk∥dk∥).)
xk+1←xs or xp
maintain or increase the step size:αk∈[αk,γαk]
else (the iteration fails)
xk+1←xk
decrease the step size:αk∈[β1αk,β2αk]
k←k+1
goto(1) if stopping criteria are not metWhile both algorithms fall under this framework and are shown to have the same convergence results, there are tradeoffs for either choice. The Mads class of algorithms not only requires a bit more complexity in choosing polling directions, but also requires simple decrease in order to accept an iteration as successful. Under the given framework, any instance of Gss requires sufficient decrease using a forcing function that must be chosen under unknown criteria. The sufficient decrease condition, however, can be implemented in such a way that it is mostly negligible and the method for choosing polling directions is simpler because of the absence of a mesh.Each iteration of this DSM framework is characterized by two steps, a Search and a Poll. The Poll step may be opportunistic (the Poll stops if a successful point has been found) and is governed by a Poll size parameter αk (Mads denotes this Δkp) such that liminf⁡k→∞αk=0 under certain assumptions. These assumptions include the requirement that at each iteration either the set of Search and Poll points lie on an integer mesh (Mads) or in order for the iteration to be a success the algorithm must find a point that satisfies sufficient decrease for some forcing function ρ:(0,∞)→(0,∞) that is continuous, nondecreasing with ρ(t)/t→0, when t↓0 as shown in Algorithm 1. The use of ρ¯ in the figure indicates the use of either a forcing function ρ (sufficient decrease) or the zero function (simple decrease).In [2], Audet and Dennis introduced the mesh adaptive direct-search (Mads) class of algorithms for (1), where the objective function f either has no available derivative information or the derivative information is unreliable and Ω is the set of feasible points determined by black-box constraint functions. Under certain assumptions, Mads is shown to converge to both first-order [2] and second-order [5] stationary points in the Clarke sense [6], as well as the previously mentioned first-order stationary points in the Rockafellar sense.Mads was introduced as a means of extending the Gps [7] class of algorithms by allowing for an infinite number of polling directions when exploring the space of variables. The first instance of Mads called LtMads [2] relied on a random lower triangular matrix to generate the set of Poll directions and relied on probabilistic arguments to obtain convergence results. This instance was later shown to have undesirably large angles between Poll directions at some iterations. To correct these deficiencies, a second instance of Mads called OrthoMads [8] was introduced that is deterministic and uses a maximal orthogonal basis at each iteration for the set of polling directions.Each iteration of Mads is characterized by two steps, a Search and a local Poll, performed on a mesh whose fineness approaches zero in the limit. At each iteration the mesh Mk is defined by
(2)Mk={x+ΔkmDz:x∈Uk,z∈NnD},
where Uk is the set of all evaluated points at the start of iteration k, Δkm is the mesh size parameter at iteration k, and D∈Rn×nD is a matrix with nD directions in Rn. The Poll step is executed if the Search step fails to find a lower objective function value than the current best solution. The Poll step evaluates mesh points adjacent to the current best solution in directions that form a positive spanning set for Rn. The Poll step can be characterized by the following definition taken from [2].Definition 1.
At iterationk, the set of Poll trial points for the Mads algorithm is defined to be
(3)Pk={xk+Δkmd:d∈Dk}⊂Mk,
where Dk is a positive spanning set such that 0∉Dk and for each d∈Dk, (i)
d can be written as a nonnegative integer combination of the directions of D,(ii)
the distance inl∞-norm from the Poll center xk to a Poll trial point is bounded above by the Poll size parameter Δkp, and(iii)
limits (as defined by Coope and Price [9]) of the normalized sets Dk are positive spanning sets.In [2] it is observed that the Poll and mesh size parameters produced by a Mads instance satisfy
(4)liminf⁡k→+∞Δkp=liminf⁡k→+∞Δkm=0.OrthoMads [8] uses the Halton sequence [10] and the scaled Householder transformation, H=∥q∥2(IN-2vvT) with q∈Rn and v=q/∥q∥, to construct an orthogonal basis for Rn at each iteration. It defines
(5)Δkp=2-lk,Δkm=min⁡{1,4-lk},
where lk is the mesh index, with l0=0 at the initial iteration. OrthoMads has been shown to behave well in practice, but recently a new variant of Mads, called QrMads [11], was introduced to counter the tendency of OrthoMads to generate a nonuniform set of polling directions.QrMads augments a single direction from a dense sequence on Sn-1 with a matrix of full rank and then uses the QR decomposition to construct an orthogonal basis that is then rounded to the current mesh. It defines
(6)Δkp=2-lk,Δkm=min⁡{4-lk-ln,4-ln},
where ln satisfies
(7)n2-n<22ln+1
and is chosen so that after rounding the basis does not become degenerate.Although QrMads was shown to have a more uniform set of polling directions than OrthoMads, the version tested in [11] is not deterministic. This is due to the fact that the augmented matrix with full rank was generated at random. Consequently, one outcome of this paper is to introduce a deterministic version of QrMads by deterministically generating a matrix with full rank.It has been observed, in some cases, that reducing the number of function evaluations at every iteration from a maximal positive spanning set to a minimal positive spanning set can improve the performance of these types of algorithms (see, e.g., [12–14]). For this reason, this paper introduces instances of Mads and Gss, called EadMads and EadGss, that have versions using both a maximal positive basis (2n polling directions) and a minimal positive basis (n+1 polling directions). This is accomplished by generating either an orthogonal basis or a regular simplex at each iteration. In the case of the maximal variant, the algorithmic framework for EadMads is exactly the same as QrMads and the new name supersedes the old. The new name was chosen to represent the use of a set of directions, for each Poll step, that have an “equal angle distribution.” By this it is meant that within the positive basis, the angle between any two directions (except a direction and its opposite, when applicable) is constant.Section2 gives the methods for constructing maximal and minimal polling sets by finding either an orthogonal set of directions or the vertices of a regular simplex centered at the origin on the unit sphere. Section 3 gives the method used by EadMads for constructing a polling set on the current mesh consisting of “nearly orthogonal” directions or a “nearly regular” simplex. Section 4 outlines the EadMads algorithm and shows that it is a valid Mads instance that shares the theoretical convergence results of this class of algorithms. Section 5 outlines the EadGss instance of Gss which is essentially the same as EadMads except that it requires sufficient decrease instead of restricting the Search and Poll points to a mesh. Finally, in Section 6, a method for making both EadMads and EadGss deterministic is detailed and numerical results comparing each against OrthoMads, QrMads, and one another on a set of test problems taken from the optimization literature are given.


## 2. Constructing a Basis

At each iteration of the Mads algorithm (and DSM framework) a positive spanning set Dk is required. In this section, we discuss a method for taking a single direction from a sequence that is dense on the unit sphere and forming either an orthogonal basis or the set of vertices of a regular simplex on the unit sphere centered at the origin that includes the given direction. We can then take the negatives of all the directions from the orthogonal basis to form a maximal positive spanning set or take the vertices of the regular simplex to form a minimal positive spanning set. Since one of the directions from the dense set is included in each positive spanning set, we will then have ⋃k=1∞‍Dk as dense in the unit sphere.


### 2.1. Orthogonal Basis

To generate an orthogonal basis from a given directionck from a dense set, we employ the QR decomposition of an n×n matrix Ak of full rank. This matrix Ak can be created by augmenting ck with any matrix of full rank (e.g., the identity or a rotation matrix) and taking the first n columns that form a basis. The QR decomposition then generates unique orthogonal matrix Qk and upper triangular matrix Rk such that QkRk=Ak. As this method will generate the same basis as the Gram-Schmidt process, the first direction will correspond to ck.


### 2.2. Regular Simplex

To generate a regular simplex,Vk, centered at the origin with vertices on the unit sphere from a given direction ck from a dense set, we again start from an n×n matrix Ak of full rank, created in the same manner as above, and use the following algorithm.LetAk∈Rn×n, with ai∈Sn-1, i=1,…,n, where ai is the ith column of Ak and {ai}i=1n is linearly independent. We construct a regular simplex, Vk, with vertices v1,…,vn+1 on the unit sphere by using the fact that vi·vj=-(1/n) for any two vertices of the simplex. This can be done as follows.(i)
Setv1=a1.(ii)
Forj=1,…,n-1 set
(8)λ=((n-j+1)2-jn(n-j+1))-1/2,vj+1=aj+1-λ∑i=1jvi(n-j+1)λ.(iii)
Set(9)vn+1=-∑i=1n‍vi.Proposition 2.
v
j
+
1 is the unique solution to the following equation:
(10)vj+1∈argmax⁡w∈Rnw·aj+1subjectto∥w∥2=1,w·vi=-1n,i=1,…,j.Due to rounding issues, this algorithm has the same unfortunate computational problems as the Gram-Schmidt process and may give an answer that does not closely approximate a regular simplex, possibly even containing complex directions. To ensure that a valid regular simplex is generated, an orthogonal matrix can be used in place of an arbitrary matrix of full rank. This can be accomplished by finding the QR decomposition ofAk and then applying the algorithm to the orthogonal matrix Qk. In the following proof, we will simply assume that Ak is already orthogonal.Proof.
To find the solution to (10) for each j, we employ the method of Lagrange multipliers. First we define the Lagrangian function as follows:
(11)L(w,λ)=f(w)-∑i=1j+1‍λici(w)=w·aj+1-∑i=1j‍λi(w·vi+1n)-λj+1(∥w∥2-1).
We then solve the system of the following equations:
(12)∇wL(w,λ)=aj+1-∑i=1j‍λivi-2λj+1w=0,(13)∥w∥2-1=0,(14)w·vk+1n=0.
First, we multiply (12) by vk and use (14) to obtain
(15)aj+1·vk-(∑i=1j‍λivi)·vk-2nλj+1=0.
We can then solve for λj+1 to obtain
(16)λj+1=-n2(aj+1·vk-(∑i=1j‍λivi)·vk),k=1,…,j.
From (12) and by induction, we can see that for each i, vi∈span⁡{a1,…,ai}. Consequently, aj+1·vi=0, i=1,…,j. We can then reduce (16) to
(17)λj+1=n2(∑i=1j‍λivi)·vk,k=1,…,j.
Using vk·vk=1, vi·vk=-1/n, when i≠k, and (17), we solve (n/2)(∑i=1j‍λivi)·v1=(n/2)(∑i=1j‍λivi)·vk to obtain
(18)λk=λ1,k=2,…,j.
Now we can substitute (18) into (17) and simplify to obtain
(19)λj+1=(n-j+12)λ1.
Next, we substitute (18) and (19) into (12) to obtain
(20)aj+1-∑i=1j‍λ1vi=(n-j+1)λ1w.
Then we take the norm of each side and use vk·vk=1, vi·vk=-1/n, when i≠k, (13), and the fact that ∥aj+1∥=1 to obtain the quadratic equation
(21)[(n-j+1)2-j(n-j+1)n]λ12-1=0.
We then get λ1=±1/(n-j+1)2-j(n-j+1)/n.
Finally, we observe(22)∇w2L=[-2λj+10…00-2λj+1…0⋮⋮⋱⋮00…-2λj+1].
Thus the maximum occurs when this is negative definite or when λj+1>0. We can then take λj+1=(n-j+1)λ1/2 and the corresponding value of λ1=1/(n-j+1)2-j(n-j+1)/n. We then substitute λ1, λj+1, and (18) into (12) and solve for w to obtain vj+1.Note, the above computation shows that when orthogonal directions are used, problem (10) has a unique solution. This is reflected by the fact that the above quadratic equation has two real solutions, one corresponding to the closest point to aj+1 within the constrained region and the other corresponding to the farthest point. In order to ensure that the process for generating a regular simplex is deterministic, it is important to note that this unique closest point exists and is found. We can also observe that each new vertex vj+1 of the regular simplex is linearly independent to the set {vi}i=1j, showing that the first n vertices of the simplex form a basis for Rn.


## 2.1. Orthogonal Basis

To generate an orthogonal basis from a given directionck from a dense set, we employ the QR decomposition of an n×n matrix Ak of full rank. This matrix Ak can be created by augmenting ck with any matrix of full rank (e.g., the identity or a rotation matrix) and taking the first n columns that form a basis. The QR decomposition then generates unique orthogonal matrix Qk and upper triangular matrix Rk such that QkRk=Ak. As this method will generate the same basis as the Gram-Schmidt process, the first direction will correspond to ck.


## 2.2. Regular Simplex

To generate a regular simplex,Vk, centered at the origin with vertices on the unit sphere from a given direction ck from a dense set, we again start from an n×n matrix Ak of full rank, created in the same manner as above, and use the following algorithm.LetAk∈Rn×n, with ai∈Sn-1, i=1,…,n, where ai is the ith column of Ak and {ai}i=1n is linearly independent. We construct a regular simplex, Vk, with vertices v1,…,vn+1 on the unit sphere by using the fact that vi·vj=-(1/n) for any two vertices of the simplex. This can be done as follows.(i)
Setv1=a1.(ii)
Forj=1,…,n-1 set
(8)λ=((n-j+1)2-jn(n-j+1))-1/2,vj+1=aj+1-λ∑i=1jvi(n-j+1)λ.(iii)
Set(9)vn+1=-∑i=1n‍vi.Proposition 2.
v
j
+
1 is the unique solution to the following equation:
(10)vj+1∈argmax⁡w∈Rnw·aj+1subjectto∥w∥2=1,w·vi=-1n,i=1,…,j.Due to rounding issues, this algorithm has the same unfortunate computational problems as the Gram-Schmidt process and may give an answer that does not closely approximate a regular simplex, possibly even containing complex directions. To ensure that a valid regular simplex is generated, an orthogonal matrix can be used in place of an arbitrary matrix of full rank. This can be accomplished by finding the QR decomposition ofAk and then applying the algorithm to the orthogonal matrix Qk. In the following proof, we will simply assume that Ak is already orthogonal.Proof.
To find the solution to (10) for each j, we employ the method of Lagrange multipliers. First we define the Lagrangian function as follows:
(11)L(w,λ)=f(w)-∑i=1j+1‍λici(w)=w·aj+1-∑i=1j‍λi(w·vi+1n)-λj+1(∥w∥2-1).
We then solve the system of the following equations:
(12)∇wL(w,λ)=aj+1-∑i=1j‍λivi-2λj+1w=0,(13)∥w∥2-1=0,(14)w·vk+1n=0.
First, we multiply (12) by vk and use (14) to obtain
(15)aj+1·vk-(∑i=1j‍λivi)·vk-2nλj+1=0.
We can then solve for λj+1 to obtain
(16)λj+1=-n2(aj+1·vk-(∑i=1j‍λivi)·vk),k=1,…,j.
From (12) and by induction, we can see that for each i, vi∈span⁡{a1,…,ai}. Consequently, aj+1·vi=0, i=1,…,j. We can then reduce (16) to
(17)λj+1=n2(∑i=1j‍λivi)·vk,k=1,…,j.
Using vk·vk=1, vi·vk=-1/n, when i≠k, and (17), we solve (n/2)(∑i=1j‍λivi)·v1=(n/2)(∑i=1j‍λivi)·vk to obtain
(18)λk=λ1,k=2,…,j.
Now we can substitute (18) into (17) and simplify to obtain
(19)λj+1=(n-j+12)λ1.
Next, we substitute (18) and (19) into (12) to obtain
(20)aj+1-∑i=1j‍λ1vi=(n-j+1)λ1w.
Then we take the norm of each side and use vk·vk=1, vi·vk=-1/n, when i≠k, (13), and the fact that ∥aj+1∥=1 to obtain the quadratic equation
(21)[(n-j+1)2-j(n-j+1)n]λ12-1=0.
We then get λ1=±1/(n-j+1)2-j(n-j+1)/n.
Finally, we observe(22)∇w2L=[-2λj+10…00-2λj+1…0⋮⋮⋱⋮00…-2λj+1].
Thus the maximum occurs when this is negative definite or when λj+1>0. We can then take λj+1=(n-j+1)λ1/2 and the corresponding value of λ1=1/(n-j+1)2-j(n-j+1)/n. We then substitute λ1, λj+1, and (18) into (12) and solve for w to obtain vj+1.Note, the above computation shows that when orthogonal directions are used, problem (10) has a unique solution. This is reflected by the fact that the above quadratic equation has two real solutions, one corresponding to the closest point to aj+1 within the constrained region and the other corresponding to the farthest point. In order to ensure that the process for generating a regular simplex is deterministic, it is important to note that this unique closest point exists and is found. We can also observe that each new vertex vj+1 of the regular simplex is linearly independent to the set {vi}i=1j, showing that the first n vertices of the simplex form a basis for Rn.


## 3. Rounding a Basis to the Mesh

Since each iteration of the Mads class of algorithms requires a positive spanning set to lie on a mesh, it is necessary to take any basis from Section 2 and obtain a close approximation to that basis lying on the given mesh. To obtain a “nearly orthogonal” integer basis or a “nearly regular” simplex with integer vertices we take each column vector qi of Qk or the first n columns vi of Vk, then project onto the hypercube centered at the origin with integer side length 2·2|lk|+2ln, and then round each component to the nearest integer as follows:
(23)di=round(2|lk|+2lnqi∥qi∥∞),
where round(·) refers to the usual rounding operation applied to each component of the vector. We then generate the new basis Hk with column vectors di. When Hk is formed from Qk we define Dk=[Hk-Hk] to obtain a maximal positive spanning set. When Hk is formed from Vk we set Dk=[Hk-Σi=1ndi] to get a minimal positive spanning set.Note, for small values of2|lk|, it is possible that the new set of vectors {di:1≤i≤n} will not form a basis. Consequently, we add the constant 2ln to the exponent in order to ensure that Hk will be nonsingular. In [11], it was shown that when Hk is formed from Qk and ln satisfies (7), the resulting set of directions will form a basis.To find the constantln when Hk is obtained from Vk, we use the fact that the distance from a matrix A to the nearest singular matrix B (in Frobenius norm) is bounded below by the smallest singular value of A; that is, if B is singular, then ∥B-A∥≥σn(A). Thus we need the singular values of the matrix comprised of the first n columns of Vk. Proposition 3 gives this information. (Note: in [13] the matrices V and VTV shown below are studied in a similar context.)Proposition 3.
LetV∈Rn×n be a matrix with columns corresponding to any n vertices of a regular simplex centered at the the origin with vertices on Sn-1. Then the singular values of V are σ1=(n+1)/n and σ2=(n+1)/n,…,σn=1/n.Proof.
We need to find the eigenvalues ofVTV. We use the fact that the vertices, v1,…,vn+1, of a regular simplex centered at the origin on Sn-1 have the property
(24)vi·vj={-(1n)ifi≠j1ifi=j
to obtain
(25)VTV=[1-1n…-1n-1n-1n1…-1n-1n⋮⋮⋱⋮⋮-1n-1n…1-1n-1n-1n…-1n1].
Then observe that
(26)[1-1n…-1n-1n1…-1n⋮⋮⋱⋮-1n-1n…1][11⋮1]=[1-n-1n1-n-1n⋮1-n-1n]=[1n1n⋮1n]=1n[11⋮1],
showing that VTV has eigenvalue 1/n. Next, observe that
(27)[1-1n…-1n-1n-1n1…-1n-1n⋮⋮⋱⋮⋮-1n-1n…1-1n-1n-1n…-1n1][1-10⋮0]=[1+1n-1-1n0⋮0]=[n+1n-n-1n0⋮0]=n+1n[1-10⋮0],
showing that (n+1)/n is an eigenvalue with eigenvector [1-10…0]T.
Similarly,[10-10…0]T,…,[100…0-1]T are also eigenvectors with eigenvalue (n+1)/n. Since these vectors are linearly independent, they form a basis for the corresponding eigenspace, thus accounting for all n eigenvalues. We finish by taking the square roots to obtain the singular values of V.Next, we need the fact that for any twon×n matrices B and C, σn(BC)≥σn(B)σn(C), where σn denotes the smallest singular value (follows from Theorem 3.1 of [21]). Then letting A be the matrix with columns 2|lk|+2lnvi/∥vi∥∞, that is,
(28)A=VT=2|lk|+2ln[v11∥v1∥∞v21∥v2∥∞⋯vn1∥vn∥∞v12∥v1∥∞v22∥v2∥∞⋯vn2∥vn∥∞⋮⋮⋱⋮v1n∥v1∥∞v2n∥v2∥∞⋯vnn∥vn∥∞]=[v11v21⋯vn1v12v22⋯vn2⋮⋮⋱⋮v1nv2n⋯vnn][2|lk|+2ln∥v1∥∞0⋯002|lk|+2ln∥v2∥∞⋯0⋮⋮⋱⋮00⋯2|lk|+2ln∥vn∥∞],
we obtain σn(A)=σn(VT)≥σn(V)σn(T)≥2|lk|+2ln(1/n). So for A+E to be nonsingular, we need
(29)∥A+E-A∥=∥E∥<2|lk|+2lnn,
where E represents the perturbation matrix resulting from rounding. Thus we require
(30)n2-n2<2|lk|+2lnn,
for all values of lk. Then letting lk=0, we obtain the following requirement for the constant ln determined only by the dimension n:
(31)n-1<22ln+1.Finally, we end with the following proposition from [11] which will be useful for showing EadMads as a valid instance of Mads and for establishing convergence results.Proposition 4.
Ifw∈Sn-1 is arbitrary, ε>0, and d=round(αw/∥w∥∞), then for large α>0 we have ∥d/∥d∥-w∥<ε/2.


## 4. TheEadMads Instance of Mads

The EadMads algorithm differs from the LtMads and OrthoMads algorithms only in the construction of the Poll directions Dk and the Poll set Pk. For each of the three instances, the mesh Mk is defined by the set of directions D=[In-In]. For the EadMads implementation with 2n polling directions, the choices of the mesh and Poll size parameters are
(32)Δkp=2-lk,Δkm=min⁡{4-lk-ln,4-ln},
with ln defined as in (7), so that Δkp/Δkm=2|lk|+2ln, for all k, as in (23). For the EadMads implementation with n+1 polling directions, the choices of the mesh and Poll size parameters are
(33)Δkp=n2-lk,Δkm=min⁡{4-lk-ln,4-ln},
with ln defined as in (31), so that Δkp/Δkm=n2|lk|+2ln, for all k. The use of the factor n reflects the possibility that -∑i=1n‍di may be up to n times as long as each di in l∞-norm and is the same strategy used for the n+1 implementation of LtMads [2]. Since Hk approximates a regular simplex, however, we can expect -∑i=1n‍di to remain close to the length of each di. The update rules for the mesh index lk are as in Algorithm 2, which describes the EadMads algorithm.Algorithm 2:The EadMads algorithm.
(0)Initialization
x0∈Ω,l0←0,k←0,t0←0
(1)Iteration  k
SEARCH (optional)
evaluatef on a finite subset Sk⊂Mk
POLL (optional if SEARCH is successful)
if the POLL size is the smallest one so far
tk←lk+t0
else
tk←1+max⁡{tj:j=0,1,…,k}
findctk, compute Hk and Dk=[Hk-Hk] or Dk=[Hk-∑i=1ndi]
evaluatef on Pk⊂Mk
(2)Updates
if the iteration is successful (there exists xs∈Sk or xp∈Pk
such thatf(xs)<f(xk) or f(xp)<f(xk))
xk+1←xs or xp
lk+1←lk-1
else (the iteration fails)
xk+1←xk
lk+1←lk+1
k←k+1
goto (1) if stopping criteria are not metThis algorithm is very similar to the algorithm found in [8], with tk and ctk in place of the Halton directions utk and qtk,lk. The integer tk and the update rules are chosen so that the algorithm will generate a subsequence of unsuccessful iterations with Δkm→0 and such that the directions used in this subsequence will correspond to the entire sequence {ck}k=1∞, where {ck}k=1∞ represents a sequence that is dense on Sn-1.The Poll directions Dk are determined by using any sequence {ci}i=1∞ that is dense on the unit sphere (or the hypercube [-1,1]n), the mesh index lk, a sequence of n×n matrices of full rank, and the index tk. At each iteration, a direction from the dense sequence ctk is augmented with a matrix to obtain a basis and then either QR decomposition or the algorithm from Section 2.2 is used to find a basis matrix with the direction from the first column corresponding to ctk. The columns of the resulting matrix are then projected onto the hypercube centered at zero with side length 2·2|lk|+2ln and rounded. The resulting matrix Hk will then be either “nearly orthogonal” or correspond to n vertices of a “nearly regular” simplex with the length of each column in l∞-norm bounded by Δkp/Δkm. In the case of the orthogonal basis, Hk is completed to a maximal positive basis composed of the 2n directions given by Dk=[Hk-Hk]. In the case of the regular simplex, Hk is completed to a minimal positive basis composed of the n+1 directions given by Dk=[Hk-∑i=1n‍di]. This then forces the trial point xk+ΔkmDkei to be exactly 2-lk from the Poll center in l∞-norm (except possibly in the case of the direction -∑i=1n‍di).The following proposition from [11] will be useful to show that EadMads has the same convergence properties as in [2].Proposition 5.
The set of normalized directions{di,k/∥di,k∥}i=1∞, with di,k=round(2|lk|+2lnci/∥ci∥∞) as in (23), is dense on Sn-1.Proposition6 was proven in [11] when EadMads uses 2n polling directions (Hk is generated from Qk) and so only the proof for the n+1 version (Hk is generated from Vk) is included.Proposition 6.
EadMads is a valid Mads instance with the convergence properties from [2] (and [1]).Proof.
We need to show that the Poll directions for EadMads satisfy the following five properties from [2, 22]. (i)
EachDk is a positive spanning set: this follows from the use of the constant ln to ensure that the round procedure on the first n columns of Vk results in a basis.(ii)
Any directionDkei(1≤i≤n+1) can be written as a nonnegative integer combination of the directions of D: this is true by construction since D=[In-In] and Dkei is an integer vector on the hypercube centered at the origin with side length 2·2|lk|+2ln.(iii)
The distance from the Poll center xk to a Poll trial point in l∞-norm is bounded above by Δkp: this is true by construction since we ensured ∥Dkei∥∞<n2|lk|+2ln and so ∥ΔkmDkei∥∞<n2|lk|+2lnΔkm=Δkp.(iv)
The set of normalized directions used over all failed iterations is dense on the unit sphere: this follows from the strategy chosen for updatingDk so that the entire sequence {ck}k=1∞ corresponds to the sequence of failed iterates, by Proposition 5, and (4) which ensures the existence of large |lk|.(v)
Limits (as defined by [9]) of convergent subsequences of the normalized sets D¯k={d/∥d∥:d∈Dk} are positive spanning sets: for this property we need to show det⁡(H¯k)>τ, for all k, for some constant τ>0. The result will then follow as in [22].
First, note that the setS={A∈Rn×n: columns of A correspond to n vertices of a regular simplex centered at the origin on the unit sphere} is bounded since S⊂{A∈Rn×n:∥A∥=n}⊂B(0,n+ε), where ∥·∥ denotes the Frobenius norm on Rn×n and B(0,n+ε) is the ball of radius n+ε centered around the matrix of zeros. Also, the set Wδ={B:∥A-B∥≤nδforanyA∈S} is bounded since any such B would satisfy ∥B∥≤n+nδ and so Wδ⊂B(0,n+nδ+ε). Now, since the determinant function is C1 on Rn×n, it is Lipschitz on the bounded set Wδ. That is, for any B,C∈Wδ there exists a constant K>0 such that
(34)|det⁡(C)-det⁡(B)|≤K∥C-B∥.
This then yields
(35)|det⁡(B)|≥|det⁡(C)|-K∥C-B∥.
So, in particular, if we take C∈S (since S⊂Wδ) and B such that ∥C-B∥≤nδ, then by applying Proposition 3 to compute |det⁡(C)|, we get
(36)|det⁡(B)|≥|det⁡(C)|-K∥C-B∥≥(n+1)n-1nn-Knδ.
If Ak is the matrix corresponding to the first n columns of Vk, then note that, for any k, Ak∈S and H¯k∈Wδk, for some δk>0, since H¯k=Qk+Ek for some perturbation matrix Ek resulting from rounding. We can then choose δ small enough so that ((n+1)n-1/nn)-Knδ>0. Then by Proposition 4, we can choose |lk| large enough, say |lk|>lδ, so that ||H¯k-Ak||≤nδ and conclude that
(37)|det⁡(H¯k)|≥(n+1)n-1nn-Knδ,000forallksuchthat|lk|>lδ.
Next, letHk have columns di and note that by construction det⁡(Hk)≠0. Since Hk has integer components we know that |det⁡(Hk)|≥1. Also, note that ∥di∥≤n2|lk|+2ln, 1≤i≤n, since each di is on the hypercube with side length 2·2|lk|+2ln. We then get the following estimate:
(38)|det⁡(H¯k)|=|det⁡[d11∥d1∥⋯dn1∥dn∥⋮⋱⋮d1n∥d1∥⋯dnn∥dn∥]|=|det⁡(Hk)|∏i=1n∥di∥≥|det⁡(Hk)|nn/22n(|lk|+2ln)≥n-n/22-n(|lk|+2ln).
Finally, note thatn-n/22-n(|lk|+2ln) is decreasing as |lk| increases. So |det⁡(H¯k)|≥n-n/22-n(lδ+2ln), if |lk|≤lδ, where lδ is defined above. We then let
(39)τ=min⁡{(n+1)n-1nn-Knδ,n-n/22-n(lδ+2ln)}
and conclude that
(40)det⁡(H¯k)>τ,forallk.


## 5. TheEadGss Instance of a Directional Direct-Search Method

At each iteration, the EadGss instance of Gss uses a positive spanning set whose directions have an equal angle distribution either by using the QR decomposition to generate an orthogonal basis or by generating a regular simplex centered at the origin. The integer tk and the update rules are chosen so that there will be a subsequence of unsuccessful iterations with αk→0 and such that the directions used in this subsequence will correspond to the entire sequence {ck}k=1∞, where {ck}k=1∞ is dense on Sn-1. The Poll size parameter is given by αk=2-lk. The update rules for lk are as in Algorithm 3, which describes the EadGss algorithm.Algorithm 3:The EadGss algorithm.
(0)Initialization
x0∈Ω,l0←0,k←0,t0←1
(1)Iteration  k
SEARCH (optional)
evaluatef on a finite set of points Sk
POLL (optional if SEARCH is successful)
if the POLL size is the smallest one so far
tk←lk+t0
else
tk←1+max⁡{tj:j=0,1,…,k}
findctk,
computeQk and Dk=[Qk-Qk] or Vk and Dk=[Vk-∑i=1nvi],
evaluatef on Pk={xk+2-lkd:d∈Dk}
(2)Updates
if the iteration is successful (i.e. there exists xs∈Sk or xp=xk+2lkd∈Pk
such thatf(xs)<f(xk)-ρ-(2-lk) or f(xp)<f(xk)-ρ-(2lk∥d∥).)
xk+1←xs or xp
lk+1←lk-1
else (the iteration fails)
xk+1←xk
lk+1←lk+1
k←k+1
goto (1) if stopping criteria are not metThe Poll directions Dk are determined by using any sequence {ci}i=1∞ that is dense on the unit sphere (or the hypercube [-1,1]n), the index lk, a sequence of n×n matrices of full rank, and the index tk. At each iteration, a direction from the dense sequence ctk is augmented with a matrix to obtain a basis and then either QR decomposition or the algorithm from Section 2.2 (along with the QR decomposition) is used to find a basis matrix with the direction from the first column corresponding to ctk. The resulting matrix will then be either orthogonal, Qk, or correspond to n vertices of a regular simplex, Vk, with the length of each column having l2-norm one. In the case of the orthogonal basis, Qk is completed to a maximal positive basis composed of 2n directions, Dk=[Qk-Qk]. In the case of the regular simplex, the first n columns of Vk are completed to a minimal positive basis composed of n+1 directions by finding the final vertex of the simplex, Dk=[Vk-∑i=1nvi], where {vi}i=1n are the columns of Vk.It should be noted that this algorithm is very similar to that for EadMads (Algorithm 2), with the exception that EadGss uses sufficient decrease at each iteration instead of confining the Search and Poll steps to a mesh. Consequently, there are only two main differences. The first is that in the final step of the construction of Dk, EadMads rounds and EadGss does not (hence no need for separate αk for the 2n and n+1 variants). The second difference is in the determination of a successful iterate; EadGss uses a forcing function to determine sufficient decrease and EadMads uses simple decrease.


## 6. Numerical Tests

In this section, numerical test results are given to compare various instances of Mads and Gss amongst one another. Each test was performed using a MATLAB implementation. The tests were performed on 130 problems taken from the optimization literature, with dimensions between nine and 60. The test problems were chosen to satisfy three categories. The first category consists of 33 smooth functions, the second category consists of 60 nonsmooth functions, and the third category consists of 37 constrained test problems, as listed in Table 1. Since many of the problems are variably dimensioned versions of one another, no generalization was used more than three times. The distribution of problems by dimension is given in Table 2.Table 1
List of test problems where S: smooth, N: nonsmooth, and C: constrained.
Function
n
type
Brown almost-linear [15]
15,40,55
S
Discrete boundary value [15]
10,30,60
S
Discretized variational problem [16]
15,55
S
Extended Rosenbrock [15]
20,40,60
S
Penalty 1 [15]
30,55
S
Trigonometric [15]
35,45
S
Watson [15]
10,25
S
Chained CB3 I [17]
25,40,50
N
Chained Cragg and Levy [16]
20,40,56
N
Chained crescent II [17]
26,34,60
N
Chained Mifflin 2 [17]
15,30,55
N
Generalized Broyden tridiagonal [16]
25,35,50
N
MXHILB [17]
20,40,50
N
Number of active faces [17]
20,35,50
N
Polak 2 [18]
10
N
Schwefel [19]
15,30,45
N
Sparse signomial [16]
20,40,60
N
Watson [18]
15,25,30
N
Wong 3 [18]
20
N
Brown with tridiagonal II [17]
16,34,56
C
Crescent with tridiagonal II [17]
30,44,56
C
Hypersphere problem [2]
15,40,60
C
MAD 8 [18]
20
C
Mifflin 2 with P20 II [17]
15,35,50
C
Number of active faces with P20 I [17]
20,35,55
C
Problem 7.29 [20]
30,45,55
C
Wong 3 [18]
20
C
Broyden tridiagonal [15]
20,50
S
Discrete integral equation [15]
20,35,50
S
Extended Powell singular [15]
16,32,44
S
Generalized Brown function I [16]
24,40
S
Penalty 2 [15]
10,35,50
S
Variably dimensioned [15]
24,45,60
S
Ackley [19]
25,45,60
N
Chained CB3 II [17]
10,45,60
N
Chained crescent I [17]
30,44
N
Chained LQ [17]
15,35,50
N
Generalized Broyden banded [16]
30,45,55
N
Generalized Maxq [17]
30,46,54
N
N generalized Brown [17]
26,40,60
N
Osborne 2 [18]
11
N
Polak 3 [18]
11
N
Shell dual [18]
15
N
Steiner 2 [18]
12
N
Wong 2 [18]
10
N
Yang [19]
34,56,60
N
CB3 with MAD1 I [17]
25,40,50
C
HS114 [18]
9
C
LQ with MAD1 I [17]
25,40,50
C
Maxq with MAD1 II [17]
10,30,50
C
MXHILB with P20 II [17]
20,45,60
C
Problem 7.26 [20]
16,24,60
C
Wong 2 [18]
10
CTable 2
Distribution of test problems used for the numerical results.
Dimension
9–12
15-16
20
24–26
30–32
34-35
40
44–46
50
54–56
60
Total
Smooth functions
3
3
3
3
3
3
3
3
3
3
3
33
Nonsmooth functions
6
5
5
6
6
5
5
6
5
5
6
60
Constrained functions
3
4
4
3
3
3
3
3
4
4
3
37
Total
12
12
12
12
12
11
11
12
12
12
12
130The choices forΔkm and Δkp for OrthoMads correspond to (5), for the 2n versions of EadMads they correspond to (6) and (7), and for the n+1 version of EadMads they correspond to (33) and (31). For each version of EadGss αk=2-lk. The mesh index lk is allowed to be negative for all instances of both algorithms and for each instance of Mads the initial mesh size Δ0m is set to one.The stopping criteria are satisfied when the number of function evaluations reach 3,000. For each implementation, a Poll step was required at every iteration and an opportunistic strategy was used (the Poll stops if a successful point has been found). For each test, no search was performed and constraints are handled using the extreme barrier method (f(x)=+∞ if x∉Ω).In order to create and test deterministic versions of both EadMads and EadGss, the Sobol sequence [23] is used in place of either the Halton sequence or the sequence of partition centers used for QrMads. This is done for two reasons. Although the Halton sequence is deterministic, the sequence of partition centers is not. In higher dimensions, however, it can be observed that the Halton sequence exhibits poor uniformity. Sobol sequences, on the other hand, are both deterministic and remain uniform. For each test, the Sobol sequence was generated using thesobolset command in MATLAB and used in the given order.At each iterationctk is taken to be the tkth Sobol direction. This is then augmented with the next 2n-1 Sobol directions, taken in order, and the identity matrix. This is done in order to ensure that a matrix of full rank is created. From this matrix, a new matrix is formed from the first n columns that form a basis for Rn and then the QR decomposition is applied to this new matrix. When a minimal positive basis is required, the resulting orthogonal matrix is then input into the simplex algorithm of Section 2.2. For EadMads, these directions are then rounded to the current mesh. In place of the mesh, EadGss uses a sufficient decrease requirement implemented with the forcing function ρ=t2.To compare the various instances of each algorithm, the results from the tests were used to generate data profiles, as outlined in [24]. For each test function, fL is defined to be the smallest value found for f by any of the instances of either Mads or Gss after 3,000 function evaluations. Then after each evaluation, the algorithm is declared a success if the smallest known value for f thus far satisfies the test
(41)f(x)≤fL+0.001(f(x0)-fL).The number of successes over all the test functions is then recorded as a percent of success. The results for each instance of Mads or Gss are then plotted as the number of function evaluations against the percent of problems to obtain the data profile.The first test performed was to compare both the nondeterministic and deterministic instances of EadMads with 2n polling directions against OrthoMads. The data profiles for the deterministic instance of EadMads, five runs of the nondeterministic instance of EadMads (formerly called QrMads), and OrthoMads are given in Figure 1. A data profile is given for each instance using the set of all test functions, the set of smooth test functions only, the set of nonsmooth test functions only, and the set of constrained test functions only. In this figure, we can see that the deterministic instance of EadMads performs very well compared to the other two instances when considering the set of all test functions or just the set of nonsmooth functions. When just considering the set of smooth or constrained test functions, the deterministic instance of EadMads performs roughly as well as the nondeterministic instance and better than OrthoMads.Figure 1
Data profiles for three instances of Mads, including five runs of QrMads, on 130 test functions.Figure2 illustrates the distribution and the density of directions for all three Mads instances. Each instance was run on the Rosenbrock function with dimensions n=10 and 20. In each case, there are roughly 500n directions shown and each direction is projected onto the plane defined by the first two coordinates. The resulting plots are typical of what can be expected regardless of the choices of the two coordinates. In this figure, we can see that the directions for both the deterministic and nondeterministic instances of EadMads have relatively uniform distributions with little clustering. The directions for OrthoMads, on the other hand tend to cluster within specific regions. While the directions for both instances of EadMads are very similar, with the unit sphere being covered nicely, the test results from this section indicate that the deterministic instance is expected to perform as well as or better than the nondeterministic instance (QrMads).Figure 2
First two coordinates of the normalized Poll directions for three variants of the Mads algorithm on the n-dimensional Rosenbrock function with n=10 and 20.Note 1.
In Figures2, 4, and 6, the EadMads and EadGss directions become more restricted towards the center as the dimension increases. This behavior is a reflection of the concentration of measure phenomenon and the fact that the observable diameter of Sn-1 converges to zero as n→∞ (for more on this topic see [25]). Intuitively, one can imagine that as n increases, the contribution to the norm of a vector from the first two elements gets smaller.The next test performed was to compare the deterministic instance of EadGss with 2n polling directions to a nondeterministic instance of Gss with 2n polling directions using random directions in place of the Sobol sequence (this is essentially the algorithm tested in [1]). The data profiles for the deterministic instance of EadGss and five runs of the nondeterministic instance are given in Figure 3. A data profile is given for each instance using the set of all test functions, the set of smooth test functions only, the set of nonsmooth test functions only, and the set of constrained test functions only. From this figure, we can observe that the deterministic instance performed better on the sets containing all the test functions, the nonsmooth test functions, and the constrained test functions. Both the deterministic and nondeterministic instances show similar performances on the set of smooth test functions.Figure 3
Data profiles for two variants of EadGss, including five runs of EadGss with random directions, on 130 test functions.Figure 4
First two coordinates of the normalized Poll directions for two variants of EadGss on the n-dimensional Rosenbrock function with n=10 and 20.Figure4 illustrates the distribution and the density of directions for both instances of Gss. Each instance was run on the Rosenbrock function with dimensions n=10 and 20. Similar to what is observed in Figure 2, we can see that the directions from these two instances of EadGss have relatively uniform distributions with little clustering.The last test performed was to compare the deterministic instances of EadMads and EadGss with both n+1 and 2n polling directions. The data profiles for each instance of the two algorithms are given in Figure 5. A data profile is given for each instance using the set of all test functions, the set of smooth test functions only, the set of nonsmooth test functions only, and the set of constrained test functions only. From this figure, we can observe that the n+1 instances of both EadMads and EadGss perform better than the 2n instances on both the set of all test problems and the set of smooth test problems. We can see a lot of parity amongst the four instances on the set of nonsmooth functions with the two n+1 instances doing a bit better as the number of function evaluations approaches 3,000. The results from the set of constrained test functions are a little more mixed with the n+1 instance of EadMads performing very well overall and the n+1 instance of EadGss performing well as the number of function evaluations approaches 3,000.Figure 5
Data profiles for four instances of EadMads on 130 test functions.Figure 6
First two coordinates of the normalized Poll directions for the variants of EadGss and EadMads using the Sobol sequence and n+1 polling directions on the n-dimensional Rosenbrock function with n=10 and 20.Figure6 illustrates the distribution and the density of directions for the n+1, deterministic instances of EadMads and EadGss. Each instance was run on the Rosenbrock function with dimensions n=10 and 20. In this figure, we can observe that the overall distribution of directions for the n+1 instances of EadMads and EadGss is very similar to one another and while they do exhibit some clustering, they still have relatively uniform distributions overall.


## 7. Conclusions

When choosing a direct-search method, often two desirable attributes have a deterministic algorithm and the ability to use a minimal positive spanning set. This paper has outlined strategies for making both the Mads and Gss algorithms deterministic, ensuring that the end results are repeatable, and outlined a method for using a set of n+1 directions generated from a regular simplex. These strategies were devised in a way that modifications are readily available without losing convergence results. For example, one could easily substitute any deterministic sequence that is dense on the unit sphere into the framework presented in this paper. Consequently, an area of future inquiry could include studying the behavior of these algorithms using a variety of different sequences.It was not the direct intention of this paper to make a comparison of the Mads algorithm and Gss algorithm using sufficient decrease. Although these two algorithms are similar, they each have advantages and disadvantages. The advantage to the Mads class of algorithms is that it only requires simple decrease, with the disadvantage being that it has the added complexity of requiring Poll and Search points to lie on a mesh at each iteration. Alternately, Gss has the advantage of eliminating the mesh requirement but is burdened by needing a sufficient decrease requirement before accepting a new best point. One can imagine that the advantages of the two algorithms may be more desirable for different real-world applications. Consequently, an area of future study could include a more in-depth comparison of these two algorithms.In this paper, constraints are handled using the extreme barrier method. Recently, however, there have been suggestions for better ways to handle constraints, such as the progressive barrier method [26]. Another area of future study would include analyzing the methods from this paper in conjunction with the progressive barrier method, as well as other extensions of the Mads algorithm.



---

*Source: 619249-2014-07-20.xml*