Simultaneous perturbation stochastic approximation (SPSA) is an algorithmic method for optimizing systems with multiple unknown parameters. It is a type of stochastic approximation algorithm. As an optimization method, it is appropriately suited to large-scale population models, adaptive modeling, simulation optimization, and atmospheric modeling. Many examples are presented at the SPSA website http://www.jhuapl.edu/SPSA. A comprehensive recent book on the subject is Bhatnagar et al. (2013). An early paper on the subject is Spall (1987) and the foundational paper providing the key theory and justification is Spall (1992).
SPSA is a descent method capable of finding global minima, sharing this property with other methods as simulated annealing. Its main feature is the gradient approximation that requires only two measurements of the objective function, regardless of the dimension of the optimization problem. Recall that we want to find the optimal control
u
∗
with loss function
J
(
u
)
:
u
∗
=
arg
min
u
∈
U
J
(
u
)
.
Both Finite Differences Stochastic Approximation (FDSA) and SPSA use the same iterative process:
u
n
+
1
=
u
n
−
a
n
g
^
n
(
u
n
)
,
where
u
n
=
(
(
u
n
)
1
,
(
u
n
)
2
,
…
,
(
u
n
)
p
)
T
represents the
n
t
h
iterate,
g
^
n
(
u
n
)
is the estimate of the gradient of the objective function
g
(
u
)
=
∂
∂
u
J
(
u
)
evaluated at
u
n
, and
{
a
n
}
is a positive number sequence converging to 0. If
u
n
is a p-dimensional vector, the
i
t
h
component of the symmetric finite difference gradient estimator is:
FD:
(
g
n
^
(
u
n
)
)
i
=
J
(
u
n
+
c
n
e
i
)
−
J
(
u
n
−
c
n
e
i
)
2
c
n
,
1 ≤i ≤p, where
e
i
is the unit vector with a 1 in the
i
t
h
place, and
c
n
is a small positive number that decreases with n. With this method, 2p evaluations of J for each
g
n
are needed. Clearly, when p is large, this estimator loses efficiency.
Let now
Δ
n
be a random perturbation vector. The
i
t
h
component of the stochastic perturbation gradient estimator is:
SP:
(
g
n
^
(
u
n
)
)
i
=
J
(
u
n
+
c
n
Δ
n
)
−
J
(
u
n
−
c
n
Δ
n
)
2
c
n
(
Δ
n
)
i
.
Remark that FD perturbs only one direction at a time, while the SP estimator disturbs all directions at the same time (the numerator is identical in all p components). The number of loss function measurements needed in the SPSA method for each
g
n
is always 2, independent of the dimension p. Thus, SPSA uses p times fewer function evaluations than FDSA, which makes it a lot more efficient.
Simple experiments with p=2 showed that SPSA converges in the same number of iterations as FDSA. The latter follows approximately the steepest descent direction, behaving like the gradient method. On the other hand, SPSA, with the random search direction, does not follow exactly the gradient path. In average though, it tracks it nearly because the gradient approximation is an almost unbiased estimator of the gradient, as shown in the following lemma.
Denote by
b
n
=
E
[
g
^
n
|
u
n
]
−
∇
J
(
u
n
)
the bias in the estimator
g
^
n
. Assume that
{
(
Δ
n
)
i
}
are all mutually independent with zero-mean, bounded second moments, and
E
(
|
(
Δ
n
)
i
|
−
1
)
uniformly bounded. Then
b
n
→0 w.p. 1.
The main idea is to use conditioning on
Δ
n
to express
E
[
(
g
^
n
)
i
]
and then to use a second order Taylor expansion of
J
(
u
n
+
c
n
Δ
n
)
i
and
J
(
u
n
−
c
n
Δ
n
)
i
. After algebraic manipulations using the zero mean and the independence of
{
(
Δ
n
)
i
}
, we get
E
[
(
g
^
n
)
i
]
=
(
g
n
)
i
+
O
(
c
n
2
)
The result follows from the hypothesis that
c
n
→0.
Next we resume some of the hypotheses under which
u
t
converges in probability to the set of global minima of
J
(
u
)
. The efficiency of the method depends on the shape of
J
(
u
)
, the values of the parameters
a
k
and
c
k
and the distribution of the perturbation terms
Δ
k
i
. First, the algorithm parameters must satisfy the following conditions:
a
t
>0,
a
t
→0 when t→∝ and
∑
t
=
1
∞
a
t
=
∞
. A good choice would be
a
t
=
a
t
;
a>0;
c
t
=
c
t
γ
, where c>0,
γ
∈
[
1
6
,
1
2
]
;
∑
t
=
1
∞
(
a
t
c
t
)
2
<
∞
Δ
t
i
must be mutually independent zero-mean random variables, symmetrically distributed about zero, with
Δ
k
i
<
a
1
<
∞
. The inverse first and second moments of the
Δ
t
i
must be finite.
A good choice for
Δ
k
i
is the Rademacher distribution, i.e. Bernoulli +-1 with probability 0.5. Other choices are possible too, but note that the uniform and normal distributions cannot be used because they do not satisfy the finite inverse moment conditions.
The loss function J(u) must be thrice continuously differentiable and the individual elements of the third derivative must be bounded:
|
J
(
3
)
(
u
)
|
<
a
3
<
∞
. Also, |J(u)|→∝ as u→∝.
In addition,
∇
J
must be Lipschitz continuous, bounded and the ODE
u
˙
=
g
(
u
)
must have a unique solution for each initial condition. Under these conditions and a few others,
u
k
converges in probability to the set of global minima of J(u) (see Maryak and Chin, 2008).
