Simultaneous perturbation stochastic approximation

Simultaneous perturbation stochastic approximation (SPSA) is an <a href="/facts/Algorithm/fnl5NmRt">algorithmic</a> method for optimizing systems with multiple unknown <a href="/facts/Parameters/zFvVFMv9">parameters</a>. It is a type of <a href="/facts/Stochastic_approximation/3C3BeMy3">stochastic approximation</a> algorithm. As an optimization method, it is appropriately suited to large-scale population models, adaptive modeling, simulation <a href="/facts/Optimization/oRn8Iv5I">optimization</a>, and <a href="/facts/Atmospheric_model/j7epJxSw">atmospheric modeling</a>. Many examples are presented at the SPSA website <a href="http://www.jhuapl.edu/SPSA">http://www.jhuapl.edu/SPSA</a>. A comprehensive 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 such as <a href="/facts/Simulated_annealing/JCmHqd3t">simulated annealing</a>. 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
 
 ∗
 
 
 
 
 {\displaystyle u^{*}}
 
 with loss
function 
 
 
 
 J
 (
 u
 )
 
 
 {\displaystyle J(u)}
 
:

u
          
            ∗
          
        
        =
        arg
        ⁡
        
          min
          
            u
            ∈
            U
          
        
        J
        (
        u
        )
        .
      
    
    {\displaystyle u^{*}=\arg \min _{u\in 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
          
        
        )
        ,
      
    
    {\displaystyle u_{n+1}=u_{n}-a_{n}{\hat {g}}_{n}(u_{n}),}

where 
 
 
 
 
 u
 
 n
 
 
 =
 (
 (
 
 u
 
 n
 
 
 
 )
 
 1
 
 
 ,
 (
 
 u
 
 n
 
 
 
 )
 
 2
 
 
 ,
 …
 ,
 (
 
 u
 
 n
 
 
 
 )
 
 p
 
 
 
 )
 
 T
 
 
 
 
 {\displaystyle u_{n}=((u_{n})_{1},(u_{n})_{2},\ldots ,(u_{n})_{p})^{T}}

represents the 
 
 
 
 
 n
 
 t
 h
 
 
 
 
 {\displaystyle n^{th}}
 
 iterate, 
 
 
 
 
 
 
 
 g
 ^
 
 
 
 
 n
 
 
 (
 
 u
 
 n
 
 
 )
 
 
 {\displaystyle {\hat {g}}_{n}(u_{n})}
 
 is the estimate of the gradient of the objective function 
 
 
 
 g
 (
 u
 )
 =
 
 
 ∂
 
 ∂
 u
 
 
 
 J
 (
 u
 )
 
 
 {\displaystyle g(u)={\frac {\partial }{\partial u}}J(u)}
 
 evaluated at 
 
 
 
 
 
 u
 
 n
 
 
 
 
 
 {\displaystyle {u_{n}}}
 
, and 
 
 
 
 {
 
 a
 
 n
 
 
 }
 
 
 {\displaystyle \{a_{n}\}}
 
 is a positive number sequence converging to 0. If 
 
 
 
 
 u
 
 n
 
 
 
 
 {\displaystyle u_{n}}
 
 is a p-dimensional vector, the 
 
 
 
 
 i
 
 t
 h
 
 
 
 
 {\displaystyle i^{th}}
 
 component of the <a href="/facts/Symmetric/hpu7DK8p">symmetric</a> 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
 
 
 
 
 
 ,
 
 
 {\displaystyle ({\hat {g_{n}}}(u_{n}))_{i}={\frac {J(u_{n}+c_{n}e_{i})-J(u_{n}-c_{n}e_{i})}{2c_{n}}},}

1 ≤i ≤p, where 
 
 
 
 
 e
 
 i
 
 
 
 
 {\displaystyle e_{i}}
 
 is the unit vector with a 1 in the 
 
 
 
 
 i
 
 t
 h
 
 
 
 
 {\displaystyle i^{th}}

place, and 
 
 
 
 
 c
 
 n
 
 
 
 
 {\displaystyle c_{n}}
 
 is a small positive number that decreases with n. With this method, 2p evaluations of J for each 
 
 
 
 
 g
 
 n
 
 
 
 
 {\displaystyle g_{n}}
 
 are needed. When p is large, this estimator loses efficiency.
Let now 
 
 
 
 
 Δ
 
 n
 
 
 
 
 {\displaystyle \Delta _{n}}
 
 be a random perturbation vector. The 
 
 
 
 
 i
 
 t
 h
 
 
 
 
 {\displaystyle i^{th}}
 
 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
 
 
 
 
 
 .
 
 
 {\displaystyle ({\hat {g_{n}}}(u_{n}))_{i}={\frac {J(u_{n}+c_{n}\Delta _{n})-J(u_{n}-c_{n}\Delta _{n})}{2c_{n}(\Delta _{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
 
 
 
 
 {\displaystyle g_{n}}
 
 is always 2, independent of the <a href="/facts/Dimension/0ktmUyac">dimension</a> 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 <a href="/facts/Approximation/LvYHsXzq">approximately</a> the <a href="/facts/Steepest/a2L3k0j5">steepest</a> 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 <a href="/facts/Unbiased/geuh5KrE">unbiased</a>
estimator of the gradient, as shown in the following lemma.

Simultaneous perturbation stochastic approximation open-in-new

Simultaneous perturbation stochastic approximation