Global Optimization via SPSA

John L. Maryak and Daniel C. Chin

The Johns Hopkins University, Applied Physics Laboratory

11100 Johns Hopkins Road, Laurel, Maryland 20723-6099

john.maryak@jhuapl.edu and daniel.chin@jhuapl.edu

Phone (240) 228-4959 and (240) 228-4956

FAX (240) 228-1093

ABSTRACT

A desire with iterative optimization techniques is that the algorithm reaches the global optimum rather

than get stranded at a local optimum value. In this paper, we examine the theoretical and numerical

global convergence properties of a certain “gradient free” stochastic approximation algorithm called

“SPSA,” that has performed well in complex optimization problems. We establish two theorems on the

global convergence of SPSA. The first provides conditions under which SPSA will converge in

probability to a global optimum using the well-known method of injected noise. The injected noise

prevents the algorithm from converging prematurely to a local optimum point. In the second theorem,

we show that, under different conditions, “basic” SPSA without injected noise can achieve

convergence in probability to a global optimum. This occurs because of the noise effectively (and

automatically) introduced into the algorithm by the special form of the SPSA gradient approximation.

This global convergence without injected noise can have important benefits in the setup (tuning) and

performance (rate of convergence) of the algorithm. The discussion is supported by numerical studies

showing favorable comparisons of SPSA to simulated annealing and genetic algorithms.

KEYWORDS: Stochastic Optimization, Global Convergence, Stochastic Approximation,

Simultaneous Perturbation Stochastic Approximation (SPSA), Recursive Annealing

1. INTRODUCTION

A problem of great practical importance is the problem of stochastic optimization, which may be stated

as the problem of finding a minimum point, p ∈ R

*

q , of a real-valued function L(q ) , called the “loss

function,” that is observed in the presence of noise. Many approaches have been devised for numerous

applications over the long history of this problem. A common desire in many applications is that the

algorithm reaches the global minimum rather than get stranded at a local minimum value. In this paper,

we consider the popular stochastic optimization technique of stochastic approximation (SA), in

particular, the form that may be called “gradient-free” SA. This refers to the case where the gradient,

g(q ) = ¶L(q ) /¶q , of the loss function is not readily available or not directly measured (even with noise).

This is a common occurrence, for example, in complex systems where the exact functional relationship

between the loss function value and the parameters, q , is not known and the loss function is evaluated

by measurements on the system (or by other means, such as simulation). In such cases, one uses

instead an approximation to g(q ) (the well-known form of SA called the Kiefer-Wolfowitz type is an

example).
The usual form of this type of SA recursion is:

)

ˆ ˆ (ˆ ˆ qk +1 = qk − ak gk qk

, (1)

where ˆ (q ) gk is an approximation (at the k

th step of the recursion) of the gradient g(q ) , and { } ak

 is a

sequence of positive scalars that decreases to zero (in the standard implementation) and satisfies other

properties. This form of SA has been extensively studied, and is known to converge to a local minimum

of the loss function under various conditions.

Several authors (e.g., Chin (1994), Gelfand and Mitter (1991), Kushner (1987), and Styblinski

and Tang (1990)) have examined the problem of global optimization using various forms of gradientfree SA. The usual version of this algorithm is based on using the standard “finite difference” gradient

approximation for ˆ (q ) gk. It is known that carefully injecting noise into the recursion based on this

standard gradient can result in an algorithm that converges (in some sense) to the global minimum. For a

discussion of the conditions, results, and proofs, see, e.g., Fang et al. (1997), Gelfand and Mitter

(1991), and Kushner (1987). The amplitude of the injected noise is decreased over time (a process

called “annealing”), so that the algorithm can finally converge when it reaches the neighborhood of the

global minimum point.

A somewhat different version of SA is obtained by using a “simultaneous perturbation” gradient

approximation, as described in Spall (1992) for multivariable ( p >1 ) problems. The gradient

approximation in simultaneous-perturbation SA (SPSA) is much faster to compute than the finitedifference approximation in multivariable problems. More significantly, using SPSA often results in a

recursion that is much more economical, in terms of loss-function evaluations, than the standard version

of SA. The loss function evaluations can be the most expensive part of an optimization, especially if

computing the loss function requires making measurements on the physical system. Several studies

(e.g., Spall (1992), Chin (1997)) have shown SPSA to be very effective in complex optimization

problems. A considerable body of theory has been developed for SPSA (Spall (1992), Chin (1997),

Dippon and Renz (1997), Spall (2000), and the references therein), but, because of the special form of

its gradient approximation, existing theory on global convergence of standard SA algorithms is not

directly applicable to SPSA. In Section 2 of this paper, we present a theorem showing that SPSA can

achieve global convergence (in probability) by the technique of injecting noise. The “convergence in

probability” results of our Theorem 1 (Section 2) and Theorem 2 (Section 3) are standard types of

global convergence results. Several authors have shown or discussed global convergence in probability

or in distribution (Chiang et al. (1987), Gelfand and Mitter (1991), Gelfand and Mitter (1993), Geman

and Geman (1984), Fang et al. (1997), Hajek (1988), Kushner (1987), Yakowitz et al. (2000), and

Yin (1999)). Stronger “almost sure” global convergence results seem only to be available by using

generally infeasible exhaustive search (Dippon and Fabian (1994)) or random search methods

(Yakowitz (1993)), or for cases of optimization in a discrete (q -) space (Alrefaei and Andradottir

(1999)).

The method of injection of noise into the recursions has proven useful, but naturally results in a

relative slowing of the rate of convergence of the algorithm (e.g., Yin (1999)) due to the continued

injection of noise when the recursion is near a global solution. In addition, the implementation of the

extra noise terms adds to the complexity of setting up the algorithm. In Section 3, we present a theorem

showing that, under different (more demanding) conditions, the basic version of SPSA can perform as a
global optimizer without the need for injected noise. Section 4 contains numerical studies

demonstrating SPSA’s performance compared to two other popular strategies for global optimization,

namely, simulated annealing and genetic algorithms; and Section 5 is a summary. The Appendix

provides some technical details.

2. SPSA WITH INJECTED NOISE AS A GLOBAL OPTIMIZER

Our first theorem applies to the following algorithm, which is the basic SPSA recursion indicated in

equation (1), modified by the addition of extra noise terms:

q k+ =q k − ak gk q k + qkwk)

ˆ ˆ (ˆ ˆ

1

, (2)

where p wk ∈R is i.i.d. N (0,I ) injected noise, a a k k = / , q q k k

k

/ loglog 2= , a > 0 , q > 0 , and ˆ (•) gk is the

“simultaneous perturbation” gradient defined as follows: (3)

ˆ ( ) (2 ) [ ( ) ( ) ]

−1 (+) ( −)

≡ Δ + Δ − − Δ + − k k k k k k k k kg q c L q c L q c e e ,

where ( ),

±

k k

c e are scalars, p Δk ∈R , and the inverse of a vector is defined to be the vector of inverses.

This gradient definition follows that given in Spall (1992). The ke terms represent (unknown) additive

noise that may contaminate the loss function observation, kc and Δkl are parameters of the algorithm,

the kc sequence decreases to zero, and the Δkl components are chosen randomly according to the

conditions in Spall (1992), usually (but not necessarily) from the Bernoulli ( ±1 ) distribution. (Uniformly

or normally distributed perturbations are not allowed by the regularity conditions.)

In this Section, we will refer to Gelfand and Mitter (1991) as GM91. Our theorem on global

convergence of SPSA using injected noise is based on a result in GM91. In order to state the theorem,

we need to develop some notation, starting with the definition of a key probability measure, hp , used in

hypothesis H8 below. Define for any h > 0 : h hdp (q )/ dq exp( 2L(q )/h )/Z

2

= − , where

∫

= − p RZ L q h dq

h

exp( 2 ( ) / )

2

. Next, define an important constant, C0, for convergence theory as follows

(GM91). For t ∈R and p u1,u2 ∈ R , let

I t d s ds g s ds

t 2

2 0

1

1 2

( ,u ,u ) inf | f ( )/ (f ( )) |

f

= +

∫

,

where the inf is taken over all absolutely continuous functions pf : R → R such that 1 f(0) =u and 2 f(t) =u ,

and |•| is the Euclidean norm. Let ( , ) lim ( , , ) u1 u2 u1 u2 V I t

t→∞

= , and { | ( ) 0} S0 = q g q = . Then

sup ( ( , ) 2 ( )) 1 2 2

,

2

3

0

1 2 0

u u u

u u

C V L

S

≡ −

∈

.

We will also need the following definition of tightness. If K is a compact subset of p R and {Xk } is a

sequence of random p -dimensional vectors, then {Xk } is tight in K if X0 ∈K and for any e > 0 , there

exists a compact subset p Ke ⊂ R such that P(Xk ∈ Ke) >1−e , ∀k > 0 . Finally, let )

ˆ

) (

ˆ

ˆ (

*

k gk qk g qk

z ≡ − and let

superscript prime (′) denote transpose.

The following are the hypotheses used in Theorem 1.

H1.. Let p Δk ∈R be a vector of p mutually independent mean-zero random variables

{ , ,..., } 1 2 Δ Δ Δ ′

k k kp such that { } Δk

 is a mutually independent sequence that is also independent of
the sequences }

ˆ

,...,

ˆ {q1 qk −1

, { ,..., }

( )

1

( )

1

±

−

±

k

e e , and { ,..., } w1 wk−1, and such that Δki is symmetrically

distributed about zero, | Δki |≤ a1 < ∞ a.s. and Δ ≤ < ∞

−

2

2

| | a ki E , a.s. ∀ i,k .

H2.. Let ( +)

k

e and (−)

k

e represent random measurement noise terms that satisfy ( ) 0

( ) ( )

− =

+ −

k k k

E e e

a.s. ∀k , where Ek denotes the conditional expectation given Ik ≡ the sigma algebra induced by

, ,..., , ,..., }

ˆ {

*

1

*

0 1 k−1 1 k− q w w z z . The { }

(±)

k

e sequences are not assumed independent. Assume that

≤ < ∞

±

3

( ) 2

[(e ) ] a k kE a.s. ∀k .

H3. L(q ) is a thrice continuously differentiable map from p R into 1 R ; L(q ) attains the minimum value

of zero; as |q |→ ∞ , we have L(q ) → ∞ and | g(q ) |→ ∞ ; inf(| ( ) | − ( ( ))) > −∞

2

g q Lap L q (Lap here is

the Laplacian, i.e., the sum of the second derivatives of L(q ) with respect to each of its

components); (q ) ≡ ∂ (q ) / ∂q ′ ∂q ′∂q ′

(3) 3 L L exists continuously with individual elements satisfying

≤ 5 < ∞

(3)

| ( )|

1 2 3

q a i i iL .

H4. The algorithm parameters have the form ga a k c c k k k = / , = / , for k = 1,2,... , where a, c > 0 ,

0 q/ a > C , and g ∈[1/ 6, 1/ 2) .

H5. [(4 4) /(4 3)] lim inf( ( ) (| ( ) | | |))

| |

1/ 2

q q q q

q

p − p − < g ′ g

→∞

.

H6. ± Δ ≤ 4 < ∞

2

)) ˆ Ek(L(qkck k a a.s. ∀k .

H7.Let wk be an i.i.d. N (0,I ) sequence, independent of the sequences }

ˆ

,...,

ˆ {q1 qk −1

, { ,..., }

( )

1

( )

1

±

−

±

k

e e ,

and { ,..., } Δ1 Δk−1.

H8.For any h h h > 0, Z < ∞; p has a unique weak limit p as h → 0 .

H9.There exists a compact subset K of p R such that }

ˆ {qk

 is tight in K .

Comments:

 (a) Assumptions H3, H5, and H8 correspond to assumptions (A1) through (A3) of GM91;

assumptions H4 and H9 supply the hypotheses stated in GM91’s Theorem 2; and the definitions of

ak and qk given in equation (2) correspond to those used in GM91. Since we will show that

assumption (A4) of GM91 is satisfied by our algorithm, this allows us to use the conclusion of their

Theorem 2.

(b) The domain of g given in H4 is one commonly assumed for convergence results (e.g., Spall

(1992)).

We can now state our first theorem as follows:

Theorem 1: Under hypotheses H1 through H9, qk

ˆ

 converges in probability to the set of global minima

of L(q ) .

Proof: See Maryak and Chin (1999), and the remark on convergence in probability in GM91, p.

1003.

3. SPSA WITHOUT INJECTED NOISE AS A GLOBAL OPTIMIZER

As indicated in the introduction above, the injection of noise into an algorithm, while providing for global

optimization, introduces some difficulties such as the need for more “tuning” of the extra terms and
retarded convergence in the vicinity of the solution, due to the continued addition of noise. This effect on

the rate of convergence of an algorithm using injected noise is technically subtle, but may have an

important influence on the algorithm’s performance. In particular, Yin (1999) shows that an algorithm of

the form (2) converges at a rate proportional to log log (k + const) , while the nominal local convergence

rate for an algorithm without injected noise is 1/3k , i.e., )

ˆ

(

1/3 *

k q k −q converges in distribution (Spall

(1992)). These rates indicate a significant difference in performance between the two algorithms.

A certain characteristic of the SPSA gradient approximation led us to question whether SPSA

needed to use injected noise for global convergence. Although this gradient approximation tends to

work very well in an SA recursion, the SPSA gradient, evaluated at any single point in q -space, tends

to be less accurate than the standard finite-difference gradient approximation evaluated at q . So, one is

led to consider whether the effective noise introduced (automatically) into the recursion by this

inaccuracy is sufficient to provide for global convergence without a further injection of additive noise. It

turns out that basic SPSA (i.e., without injected noise) does indeed achieve the same type of global

convergence as in Theorem 1, but under a different, and more difficult to check, set of conditions.

In this Section, we designate Kushner (1987) as K87, and Kushner and Yin (1997) as KY97.

Here we are working with the basic SPSA algorithm having the same form as equation (1):

)

ˆ ˆ (ˆ ˆ q k+1 = qk − ak gk qk

, (4)

where ˆ (•) gk is the simultaneous-perturbation approximate gradient defined in Section 2, and now

(obviously) no extra noise is injected into the algorithm. For use in the subsequent discussion, it will be

convenient to define

)| )

ˆ

) (

ˆ

) (ˆ (

ˆ

( bk q k ≡ E gk qk − g qk אk, and ) | )

ˆ

) ( ˆ (

ˆ

) ˆ (

ˆ

( k k gk k E gk k ke q ≡ q − q א ,

where אk denotes the s -algebra generated by }

ˆ

,...,

ˆ

,

ˆ {q1 q 2 q k

, which allows us to write equation (4) as

)] ˆ) (

ˆ

) (

ˆ

[ (

ˆ ˆ

k 1 k k k k k bk k q + =q − a g q + e q + q . (5)

Another key element in the subsequent discussion is the ordinary differential equation (ODE):

q = g(q )

&

, (6)

which, in Lemma 1 of the Appendix is shown to be the “limit mean” ODE for algorithm (4).

Now we can state our assumptions for Theorem 2, as follows:

J.1 . Let p Δk ∈R be a vector of p mutually independent mean-zero random variables

{ , ,..., } 1 2 Δ Δ Δ ′

k k kp such that { } Δk

 is a mutually independent sequence and Δkis independent of

the sequence }

ˆ

,...,

ˆ {q1 qk −1

, and such that Δki is ∀ i,k symmetrically distributed about zero,

| Δki |≤ a1 < ∞ a.s. and Δ ≤ < ∞

−

2

2

| | a ki E .

J.2 Let (+)

k

e and (−)

k

e represent random measurement noise terms that satisfy (( ) | ) 0

( ) ( )

− א =

+ −

E k k ke e

a.s. ∀k . The { }

(±)

k

e sequences need not be assumed independent. Assume that

א ≤ < ∞

±

3

( ) 2

((e ) | k) a kE a.s. ∀k .

J.3 (a). L(q ) is thrice continuously differentiable and the individual elements of the third derivative

satisfy ≤ 5 < ∞

(3)

| ( )|

1 2 3

q a i i iL .

(b). L(q ) → ∞ as q → ∞ .
J.4 The algorithm parameters satisfy the following: the gains ak > 0 , ak →0 as k → ∞ , and

∑

∞

=

= ∞ k 1ak. The sequence { }kc is of form gc c k k = / , where c > 0 and g ∈[1/ 6, 1/ 2) , and

∑

∞

=

< ∞ 0

2

( / )

k

k k a c .

J.5 The gradient g(q ) is bounded and Lipschitz continuous.

J.6 The ODE (6) has a unique solution for each initial condition.

J.7 For the ODE (6), suppose that there exists a finite collection of disjoint compact stable invariant

sets (see K87) K1,K2,...,Km , such that UiKi contains all the limit sets for (6). These sets are

interpreted as closed sets containing all local (including global) minima of the loss function.

J.8 For any h h h > 0, Z < ∞; p has a unique weak limit p as h → 0 ( hZ and hp are defined in Section

2).

J.9 ∑ =< ∞ ∀

k

i

i i E e k

1

) |

ˆ

| (q .

J.10 For any asymptotically stable (in the sense of Liapunov) point, q , of the ODE (6), there exists

a neighborhood of the origin in p R such that the closure, Q2, of that neighborhood satisfies

q +Q2 ≡{q + y : y ∈Q2} ⊂ Θ , where p Θ ⊂ R denotes the allowable q -region. There is a

neighborhood, Q1, of the origin in p R and a real-valued function ( , ) H1 y1 y2, continuous in

Q1 ×Q2, whose y1-derivative is continuous on Q1 for each fixed y2 ∈Q2, and such that the

following limit holds. For any c,Δ > 0 , with c being an integral multiple of Δ , and any functions

( ), ( ) 1 2 y • y • taking values in Q1 ×Q2 and being constant on the intervals [iΔ,iΔ+ Δ), iΔ < c , we have

[∑ ∑ ] ∫

Δ −

=

+ −

Δ = + +

′

Δ

=

( / ) 1

0

1

1

,

2

0

1 1

)

ˆ

( ( ), ( )) lim sup log exp ( ) (

c c

y y y q i

im m

j im n j n j

m n

E i b

m

H s s ds . (7)

Also, there is a function ( ) H2 y 3 that is continuous and differentiable in a small neighborhood of

the origin, and such that

= ∫H ( (s))ds

0

2 1

c

y 







′ Δ

Δ ∑ ∑Δ −

=

+ −

=

+ +

( / ) 1

0

1

1

,

)

ˆ

lim sup log exp ( ) (

c

y q

i

im m

j im n j n j

m n

E i e

m

. (8)

A bit more notation is needed. Let T > 0 be interpreted such that [0,T] is the total time period

under consideration in ODE (6). Let

( , ) sup[ ( ( )) ( , )

( , ) 0.5[ (2 , ) (2 )

2 1 2 1 2

1 2 1 1 2 2 1

1

b y y b y y y

y y y y y

y

L g H

H H H

= ′ − −

= +

, (9)

and, for 1f(0) = x∈ R , define the function

∫

=

T

S T L s s ds

0

( ,f ) (f ( ),f ( )) &,

if f(•) is a real-valued absolutely-continuous function on [0,T] and to take the value ∞

otherwise. S(T,f ) is the usual action functional of the theory of large deviations (adapted to our

context). Define ∑

−

=

≡

1

0

n

i

n ai

t , and ∑

−

=

= +

1

0

k

i

n i

n

k

t a . Define }

ˆ {

n

k

q and (•)

n

q by

= x ∈Θ

n

0

ˆq , )ˆ

ˆ (

ˆ ˆ

1

n

n k n k k

n

k

n

k

q q a +g + q

+

= − , and n

k

n

q t q

ˆ

( ) = for [ , )

1

n

k

n

k

t t t

+

∈ .

Now we can state the last two assumptions for Theorem 2:
J.11 For each d > 0 and i =1,2,..., m , there is a r -neighborhood of Ki

, denoted ( ) Nr Ki, and

d r > 0, Tr < ∞ such that, for each , ( ) N Kix y ∈ r, there is a path, f(•) , with x T y f(0) = , f( y) = ,

where Ty ≤ Tr and f d S(Tr, ) ≤ .

J.12 There is a sphere, D1, such that D1 contains UiKi in its interior, and the trajectories of (•)

n

q

stay in D1. All paths of ODE (6) starting in D1 stay in D1.

Note 1. Assumptions J1, J2, and J3(a) are from Spall (1992), and are used here to characterize the

noise terms )

ˆ

( bk q k and )

ˆ

( k ke q . Assumption J3(b) is used on page 178 of K87. Assumption J4

expresses standard conditions on the algorithm parameters (see Spall (1992)), and implies

hypothesis (A10.2) in KY97, p. 174. Assumptions J5 and J6 correspond to hypothesis

(A10.1) in KY97, p. 174. Assumption J7 is from K87, p. 175. Assumption J8 concerns the

limiting distribution of qk

ˆ

. Assumption J9 is used to establish the “mean” criterion for the

martingale sequence in Lemma 2. Assumptions J11 and J12 are the “controllability” hypothesis

A4.1 and the hypothesis A4.2, respectively, of K87, p. 176.

Note 2. Assumption J10 corresponds to hypotheses (A10.5) and (A10.6) in KY97, pp. 179-181.

Although these hypotheses are standard forms for this type of large deviation analysis, it is

important to justify their reasonableness. The first part (equation (7), involving noise terms

)

ˆ

( bk q k) of J10 is justified by the discussion in KY97, p. 174, which notes that the results of their

subsection 6.10 are valid if the noise terms (that they denote xn) are bounded. This discussion

is applicable to our algorithm since the )

ˆ

( bk q k noise terms were shown by Spall (1992) to be

( ) ( 0) O c

k

2

ck → a.s. The second part (equation (8), involving noise terms )

ˆ

( k ke q ) is justified by

the discussion in KY97, p. 174, which notes that the results in their subsection 6.10 are valid if

the noise terms they denote dMn (corresponding to our noise terms )

ˆ

( k ke q ) satisfy the

martingale difference property that we have established in Lemma 2 of the Appendix.

Now we can state our main theorem:

Theorem 2. Under assumptions J1 through J12, qk

ˆ

 converges in probability to the set of global minima

of L(q ) .

The idea of the proof is as follows (see the Appendix for the details). This theorem follows from results

(in a different context) in K87 for an algorithm ) ]

ˆ

[ (

ˆ ˆ q k+1 = qk − ak g q k +z k

, where zk is i.i.d. Gaussian

(injected) noise. In order to prove our Theorem 2, we start by writing the SPSA recursion as

) ]

ˆ

[ (

ˆ ˆ *

k 1 k ak g k k q = q − q +z +

, where )

ˆ

) (

ˆ

ˆ (

*

k gk qk g qk

z ≡ − is the “effective noise” introduced by the inaccuracy

of the SPSA gradient approximation. So, our algorithm has the same form as that in K87. However,

since *

k

z is not i.i.d. Gaussian, we cannot use K87’s result directly. Instead, we use material in Kushner

and Yin (1997) to establish a key “large deviation” result related to our algorithm (4), which allows the

result in K87 to be used with *

k

z replacing the zk in his algorithm.

4. NUMERICAL STUDIES: SPSA WITHOUT INJECTED NOISE
4.1. Two-Dimensional Problem

A study was done to compare the performance of SPSA to a recently published application of the

popular genetic algorithm (GA). The loss function is the well-known Griewank function (see Haataja

(1999)) defined for a two-dimensional ( , ) 1 2 q = t t ′ , by:

( ) cos( 100)cos[( 100) / 2] [( 100) ( 100) ]/ 4000 1

2

2

2

L q = t

1 − t2 − − t1 − + t − − ,

which has thousands of local minima in the vicinity of a single global minimum at q = (100,100)′ at which

L(q ) = 0 . Haataja (1999) describes the application of a GA to this function (actually, to find the

maximum of −L(q ) ) based on noise-free evaluations of L(q ) (i.e., ek = 0 ). This study achieved a

success rate of 66% (see Haataja’s Table 1.3, p.16) in 50 independent trials of the GA, using 300

generations and 9000 L(q ) evaluations in each run of the GA. Haataja’s definition of a successful

solution is a reported solution where the norm of the solution minus the correct value, *q , is less than

0.2, and the value of the loss function at the reported solution is within 0.01 of the correct value of zero.

We examined the performance of basic SPSA (without adding injected noise) on this problem, using

a a a /(A k ) k = + , with A = 60, a =100 and a = .602 , a slowly decreasing gain sequence of a form that has

been used in many applications (see Spall (1998)). For the gradient approximation (equation (3)), we

chose each component of Δk to be an independent sample from a Bernoulli (±1) distribution, and

g

c c k k = / , with c = 10 and g = .101 . Since we used the exact loss function, the ke noise terms were zero.

We ran SPSA, allowing 3000 function evaluations in each of 50 runs, and starting the algorithm (each

time) at a point randomly chosen in the domain [−200, 400]×[−200, 400]. Haataja’s q -domain was also

constrained to lie in a box, but the dimensions of the box were not specified. Hence we chose a domain

that is a cube centered at the global minimum, in which there are many local minima of L(q ) (as seen in

Haataja’s (1999) Figure 1.1). SPSA successfully located the global minimum in all 50 runs (100%

success rate).

4.2. Ten-Dimensional Problem

For a more ambitious test of the global performance of SPSA, we applied SPSA to a loss function

given in Example 6 of Styblinski and Tang (1990), which we will designate for convenience as ST90.

The loss function is:

( ) (2 ) 4 cos( )

1 1

1 2 ∑ ∏

= =

−

= −

p

i

i

p

i

i L q p t p t ,

where p =10 and ( ,..., ) 1 = ′ p q t t . This function has the global minimum value of −40 at the origin, and a

large number of local minima. As in the two-dimensional study above, we used the exact loss function.

Our goal is to compare the performance of SPSA without injected noise with simulated annealing and

with a GA.

For the simulated annealing algorithm, we use the results reported in ST90. They used an

advanced form of simulated annealing called fast simulated annealing (FSA). According to ST90, FSA

has proven to be much more efficient than classical simulated annealing due to using Cauchy (rather than

Gaussian) sampling and using a fast (inversely linear in time) cooling scheme. For more details on FSA,
see ST90. The results of their application of FSA to the above L(q ) are given in Table 1 below (FSA

values taken from Table 10 of ST90). Table 1 shows the results of 10 independent runs of each

algorithm. In each case (each run of each algorithm), the best value of L(q ) found by the algorithm is

shown. In their study, although FSA was allowed to use 50,000 function evaluations for each of the

runs, the algorithm showed very limited success in locating the global minimum. It should be noted that

the main purpose of the ST90 paper was to examine a relatively new algorithm, stochastic

approximation combined with convolution smoothing. This algorithm, which they call SAS, was much

more effective than FSA, yielding results between those shown in Table 1 for GA and SPSA.

For the genetic algorithm (GA), we implemented a GA using the popular features of elitism (elite

members of the old population pass unchanged into the new population), tournament selection

(tournament size = 2), and real-number encoding (see Mitchell (1996), pp. 168, 170, and 157,

respectively). After considerable experimentation, we found the following settings for the GA algorithm

to provide the best performance on this problem. The population size was 100, the number of elite

members (those carried forward unchanged) in each generation was 10, the crossover rate was 0.8,

and mutation was accomplished by adding a Gaussian random variable with mean zero and standard

deviation 0.01 to each component of the offspring. The original population of 100 (10-dimensional) q -

vectors was created by uniformly randomly generating points in the 10-dimensional hypercube centered

at the origin, with edges of length 6 (so, all components had absolute value less than or equal to 3

radians). We constrained all component values in subsequent generations to be less than or equal to 4.5

in absolute value. This worked a bit better than constraining them to be less than 3, since, with the

tighter constraints, the GA got stuck at the constraint boundary and could not reach local minima that

were just over the boundary. All runs of the GA algorithm reported here used 50,000 evaluations of the

loss function. The results of the 10 independent runs of GA are shown in Table 1. Although the

algorithm did reasonably well in getting close to the minimum loss value of –40, it only found the global

minimum in one of the 10 runs (run #8). In the other nine cases, a few (typically two or four) of the

components were trapped in a local minimum (around ± pi radians), while the rest of the components

(approximately) achieved the correct value of zero. Note that the nature of the loss function is such that

the value of L(q ) is very close to an integer (e.g., –39.0 or –38.0) when an even number (e.g., 2 or 4)

of components of q are near ± pi radians.

We examined the performance of basic SPSA (without adding injected noise), using the

algorithm parameters defined in Subsection 4.1 with A = 60, a =1 , a = .602 , c = 2 , and g = .101 . We

started q at t

i = 3 radians, i =1,..., p , resulting in an initial loss function value of −31 . This choice of

starting point was at the outer boundary of the domain in which we chose initial values for the GA

algorithm, and we did not constrain the search space for SPSA as we did for GA (the initialization and

search space for FSA were not reported in ST90). We ran 10 Monte Carlo trials (i.e., randomly

varying the choices of Δk). The results are tabulated in Table 1. The results of these numerical studies

show a strong performance of the basic SPSA algorithm in difficult global optimization problems

.
Table 1. Best Loss Function Value in Each of 10 Independent Runs of Three Algorithms

Run SPSA GA FSA

1 −40.0 −38.0 −24.9

2 −40.0 −39.0 −15.5

3 −40.0 −39.0 −29.0

4 −40.0 −38.0 −32.1

5 −40.0 −37.0 −30.2

6 −40.0 −39.0 −30.1

7 −40.0 −38.0 −27.9

8 −40.0 −40.0 −20.9

9 −40.0 −38.0 −28.5

10 −40.0 −39.0 −34.6

Average Value −40.0 −38.5 −27.4

Number of Function Evaluations 2,500 50,000 50,000

5. SUMMARY

SPSA is an efficient gradient-free SA algorithm that has performed well on a variety of complex

optimization problems. We showed in Section 2 that, as with some standard SA algorithms, adding

injected noise to the basic SPSA algorithm can result in a global optimizer. More significantly, in

Section 3, we showed that, under certain conditions, the basic SPSA recursion can achieve global

convergence without the need for injected noise. The use of basic SPSA as a global optimizer can

ease the implementation of the global optimizer (no need to tune the injected noise) and result in a

significantly faster rate of convergence (no extra noise corrupting the algorithm in the vicinity of the

solution). In the numerical studies, we found significantly better performance of SPSA as a global

optimizer than for the popular simulated annealing and genetic algorithm methods, which are often

recommended for global optimization. In particular, in the case of a 10-dimensional optimization

parameter (q ), the fast simulated annealing and genetic algorithms generally failed to find the global

solution.

APPENDIX (LEMMAS RELATED TO THEOREM 2 AND PROOF OF THEOREM 2)

In this Appendix, we designate Kushner (1987) as K87, and Kushner and Yin (1997) as KY97. Here

we are working with the basic SPSA algorithm as defined in equation (4):

)

ˆ ˆ (ˆ ˆ q k+1 = qk − ak gk qk

.

We first establish an important preliminary result that is needed in order to apply the results from K87

and KY97 in the proof of Theorem 2.

Lemma 1. The ordinary differential equation (eq. (6) above),

q = g(q )

&

,

is the “limit mean ODE” for algorithm (4).

Proof: Examining the definition of limit mean ODE given in KY97, pp. 174 & 138, it is clear that we

need to prove that )] ( )

ˆ

) (

ˆ

[ ( ) (

1 1

g q e q b q g q

m n

k n

k k k k m ∑ + + →

+ −

=

 w.p. 1 as m,n → ∞ . Since Spall
(1992) has shown that ) 0

ˆ bk

(qk → w.p. 1, we can conclude using Cesaro summability that the

contribution of the )

ˆ

( bk q k terms to the limit is zero w.p. 1. For the )

ˆ

( k ke q terms, we have by

definition that )] 0

ˆ E[ek

(qk = ; hence, by the law of large numbers, the contribution of the )

ˆ

( k ke q

terms to the limit is also zero. Q.E.D.

Our next Lemma relates to Note 2 in Section 3.

Lemma 2. Under assumptions J1, J3(a), and J9, the sequence )} ˆ { ( k ke q is a אk-martingale difference.

Proof: It is sufficient to show that ∑ =≡

k

i

k k k M e

1

)

ˆ

(q is a אk-martingale. Assumption J9 satisfies the first

requirement (see KY97, p.68) of the martingale definition, that E M k < ∞ . For the main

requirement, we have for any k:

( ) [ ( ) ] 1 1 1 1 1| ,...,

ˆ E Mk+

| Mk,...,M = E ek+ qk+ + Mk Mk M

) | ,..., ]

ˆ

[ ( 1 1 M M1 M E e k k+ k+ k= + q

)]| }

ˆ

)|

ˆ

) ( ˆ (

ˆ {[ ˆ ( = M k + E gk+1 qk+1 − E gk +1 qk +1 qk+1 M k

}

ˆ

)]| ,

ˆ

)|

ˆ

) ( ˆ (

ˆ {[ ˆ ( ˆ 1 1 1 1 1 1

1

= + + + − + + + +

+

M k E E gk k E gk k k M k k

k

q q q q q

Mk E E gk k Mk k E E E gk k Mk k Mk

k k

= + + + + − + + + =

+ +

]

ˆ

) | ,

ˆ

) [ ( ˆ (

ˆ

) | ,

ˆ

( ˆ ( ˆ 1 1 1 ˆ 1 1 1

1 1

q q q q q q

,

where 

1

ˆ

k+

E

q

 denotes expectation conditional on 1

ˆq k+

, and all equalities concerning conditional

expectations are w.p. 1. Q.E.D.

A key step in the proof of our main result (Theorem 2 below) is establishing the following “large

deviation” result (Lemma 3). Let Bx be a set of continuous functions on [0,T] taking values in Θ and

with initial value x . Let 0 Bx denote the interior of Bx, and Bx denote the closure.

Lemma 3. Under assumptions J4, J5, J6, and J10, we have

inf ( , ) lim inf log { ( ) }

0

x

n n

x

B n

S T P B

x

− ≤ • ∈

∈

f q

f

lim suplog {q ( ) } inf ( ,f )

f

P B S T

Bx

x

n n

x

n ∈

≤ • ∈ ≤ − , (9)

where n Px denotes the probability under the condition that x

n

q (0) = .

Proof: This result is adapted from Theorem 10.4 in KY97, p. 181. Note that our assumption J10 is a

modified form of their assumptions (A10.5) and (A10.6), using “equals” signs rather than

inequalities. The two-sided inequality in (9) follows from J10 by an argument analogous to the

proof of KY87’s Theorem 10.1 (p. 178), which uses an “equality” assumption ((A10.4), p.

174) to arrive at a two-sided large deviation result analogous to (9) above. Q.E.D.

We restate our main theorem:

Theorem 2: Under hypotheses J1 through J12, qk

ˆ

 converges in probability to the set of global minima

of L(q ) .

Proof: This result follows from a discussion in K87. Theorem 2 of K87, (p.177) describes

probabilities involving expected times for the SA algorithm (system (1.1) of K87) to transition
from one Ki

 to another. The SA algorithm he uses can be written in our notation as

) ]

ˆ

[ (

ˆ ˆ q k+1 = qk − ak g q k +z k

, where zk is i.i.d. Gaussian (injected) noise. The K87 Theorem 2 uses

the i.i.d. Gaussian assumption only to arrive at a large deviation result exactly analogous to our

Lemma 3. The subsequent results in K87 are based on this large deviation result. Recall that

the SPSA algorithm without injected noise can be written in the form ) ]

ˆ

[ (

ˆ ˆ *

k 1 k ak g k k q = q − q +z +

.

Since we have established Lemma 3 for SPSA, the results of K87 hold for the SPSA algorithm

with its “effective” noise { }

*

k

z replacing the { } z k sequence used in K87. In particular, K87’s

discussion (pp. 178, 179) of his Theorem 2 is applicable to our Theorem 2 context (SPSA

without injected noise), which corresponds to K87’s “potential case.” Note that our

formulation corresponds to the K87 setup where b(x,x ) = b(x) in his notation, which, by the

comment in K87, p. 179, means that his discussion is applicable to his system (1.1) and hence

to our setup. In his discussion on p. 179, K87 indicates that the difference between the

measure of Xn (which corresponds to our qk

ˆ

) and the invariant measure (which we have

denoted hp ) converges asymptotically (n,k → ∞,h →0 ) to the zero measure weakly. This means

that, in the limit as k → ∞ , qk

ˆ

is equivalent to p in the same sense as in Theorem 2 of Gelfand

and Mitter (1991), and the desired convergence in probability follows as in Theorem 1 above.

Q.E.D.

ACKNOWLEDGEMENT

This work was partially supported by the JHU/APL Independent Research and Development Program

and U.S. Navy contract N00024-98-D-8124.

REFERENCES

[1] Alrefaei, M.H. and Andradottir, S. (1999), “A Simulated Annealing Algorithm with Constant

Temperature for Discrete Stochastic Optimization,” Management Science, 45, pp. 748-764.

[2] Chaing T-S., Hwang, C-R., and Sheu, S-J. (1987), “Diffusion for Global Optimization in n R ,”

SIAM J. Control Optim., 25, pp. 737-753.

[3] Chin, D.C. (1997), “Comparative Study of Stochastic Algorithms for System Optimization Based

on Gradient Approximations,” IEEE Trans. Systems, Man, and Cybernetics – Part B:

Cybernetics, 27, pp. 244-249.

[4] Chin, D.C. (1994), “A More Efficient Global Optimization Algorithm Based on Styblinski and

Tang,” Neural Networks, 7, pp. 573-574.

[5] Dippon, J. and Fabian, V. (1994), “Stochastic Approximation of Global Minimum Points,” J.

Statistical Planning and Inference, 41, pp. 327-347.

[6] Dippon, J. and Renz, J. (1997), “Weighted Means in Stochastic Approximation of Minima,” SIAM

J. Control Optim., 35, pp. 1811-1827.

[7] Fang, H., Gong, G., and Qian, M. (1997), “Annealing of Iterative Stochastic Schemes,” SIAM J.

Control Optim., 35, pp.1886-1907.

[8] Gelfand, S.B. and Mitter, S.K. (1991), “Recursive Stochastic Algorithms for Global Optimization in

d R ,” SIAM J. Control Optim., 29, pp. 999-1018.
[9] Gelfand, S.B. and Mitter, S.K. (1993), “Metropolis-Type Annealing Algorithms for Global

Optimization in d R ,” SIAM J. Control Optim., 31, pp. 110-131.

[10] Geman S. and Geman, D. (1984), “Stochastic Relaxation, Gibbs Distributions, and the Bayesian

Restoration of Images, IEEE Trans. Pattern Analysis and Machine Intelligence, PAMI-6, pp.

721-741.

[11] Haataja, J. (1999), “Using Genetic Algorithms for Optimization: Technology Transfer in Action,”

Chapter 1, pp. 3 – 22, in Evolutionary Algorithms in Engineering and Computer Science,

Edited by K. Miettinen, M.M. Makela, P. Neittaanmaki, and J. Periaux, Wiley, Chichester.

[12] Hajek, (1988), “Cooling Schedules for Optimal Annealing,” Mathematics of Operations

Research, 13, pp. 311-329.

[13] Kushner, H.J. and Yin, G.G. (1997), Stochastic Approximation Algorithms and Applications,

Springer, New York.

[14] Kushner, H.J. (1987), “Asymptotic Global Behavior for Stochastic Approximation and Diffusions

with Slowly Decreasing Noise Effects: Global Minimization Via Monte Carlo,” SIAM J. Appl.

Math., 47, pp. 169-185.

[15] Maryak, J.L. and Chin, D.C. (1999), “Efficient Global Optimization Using SPSA,” Proc. Amer.

Control Conf., San Diego, June 2-4, pp. 890-894.

[16] Mitchell, M. (1996), An Introduction to Genetic Algorithms, MIT Press, Cambridge, Mass.

[17] Spall, J.C. (2000), “Adaptive Stochastic Approximation by the Simultaneous Perturbation

Method,” IEEE Trans. Automat. Control, 45, pp. 1839-1853.

[18] Spall, J.C. (1998), “Implementation of the Simultaneous Perturbation Algorithm for Stochastic

Optimization,” IEEE Trans. Aerospace and Electronic Systems, 34, pp. 817-823.

[19] Spall, J.C. (1992), “Multivariate Stochastic Approximation Using a Simultaneous Perturbation

Gradient Approximation,” IEEE Trans. Automat. Control, 37, pp. 332-341.

[20] Styblinski, M.A. and Tang, T.-S. (1990), “Experiments in Nonconvex Optimization: Stochastic

Approximation with Function Smoothing and Simulated Annealing,” Neural Networks, 3, pp. 467-

483.

[21] Yakowitz, S. (1993), “A Globally Convergent Stochastic Approximation,” SIAM J. Control

Optim. 31 pp. 30-40.

[22] Yakowitz, S., L’Ecuyer, P., and Vazquez-Abad, F. (2000), “Global Stochastic Optimization with

Low-Dispersion Point Sets,” Operations Research, 48, pp. 939-950.

[23] Yin, G. (1999), “Rates of Convergence for a Class of Global Stochastic Optimization Algorithms,”

SIAM J. Optim., 10, pp. 99-120.