Adaptive initial step size selection for Simultaneous Perturbation Stochastic Approximation

A difficulty in using Simultaneous Perturbation Stochastics Approximation (SPSA) is its performance sensitivity to the step sizes chosen at the initial stage of the iteration. If the step size is too large, the solution estimate may fail to converge. The proposed adaptive stepping method automatically reduces the initial step size of the SPSA so that reduction of the objective function value occurs more reliably. Ten mathematical functions each with three different noise levels were used to empirically show the effectiveness of the proposed idea. A parameter estimation example of a nonlinear dynamical system is also included.

evaluations of SPSA is D-fold smaller compared to FDSA (Spall 1998b). An extension to this method exists to include second-order (Hessian) effects to accelerate convergence (Spall 2000(Spall , 2009Zhu and Spall 2002). However, we will not treat this enhancement here.
The problem solved by SPSA in this paper can be formulated as following.
where f (θ) is the objective function and θ is a D-dimensional vector of parameters. We assume that each element in the vector θ is a real number and has upper and lower bounds that defines the Cartesian product domain Θ. The SPSA and FDSA procedures are in the general recursive form: where ĝ k (θ k ) is the estimate of the gradient vector g(θ) at iteration k based on the measurements of the objective function. The a k is the step size at iteration k. Equation (2) is analogous to the gradient descent algorithm in nonlinear programming, in which g k is the gradient of the objective function ∇f (θ k ). The difference is that in Eq.
(2), ĝ k represent gradients stochastically and the effect of the noise or deviation from the true gradient is expected to cancel out as the iteration count k increases. The step sizes a k are normally prescribed in SPSA and FDSA as a function of k just like the Simulated-Annealing's (Kirkpatrick et al. 1983) cooling schedule. This is because these methods do not assume deterministic responses in the measurements of the objective function values. Thus, unlike the nonlinear programming counterparts, adaptation of step sizes based on gradients and amount of descent achieved (such as in the line search) is usually not done in the stochastic approximation optimization methods. The rationale of the Eq.
(2) is intuitively depicted in Fig. 1 for one variable case. Under appropriate conditions, the iteration in Eq.
(2) will converge to the optimum θ * in some stochastic sense. The hat symbol indicates an "estimate". Thus, θ k denotes the estimate of the optimum θ * at iteration k. Let y(·) denote a measurement of the objective (2) θ k+1 =θ k − a kĝk (θ k ), Fig. 1 Objective value minimization using gradient descent (one variable): if gradient g is positive at θ k then move to θ k+1 < θ k , if gradient g is negative then move to θ k+1 > θ k function f (·) at parameter value denoted by " · " and c k be some small positive number. The measurements are assumed to contain some noise, i.e. y(·) = f (·) + noise. In SPSA, the ith component ĝ ki (θ k ) of the gradient vector ĝ k (θ k ) is formed from a ratio involving the individual components in the perturbation vector and the difference in the two corresponding measurements. For two-sided simultaneous perturbations, we have where the D-dimensional random perturbation vector follows a specific statistical distribution criterion. Here, i is the parameter index. A simple choice for each component of k is to use Bernoulli ±1 distribution, which is essentially a random switching between +1 and − 1. The Bernoulli distribution is proven to be an optimal distribution for the simultaneous perturbation (Sadegh and Spall 1997). Note also that in the Eq. (3), we do not evaluate y(θ k ). The recursive equation (2) proceeds with only the responses from the two perturbed inputs y(θ k + c k k ) and y(θ k − c k k ).
The choice of a k and c k is critical to the performance of SPSA and suggested values can be found in Spall (1998a). At given iteration k: where α = 0.602, γ = 0.101, c ≃ standard deviation of measurement noise , A ≤ 10 % of maximum number of iterations , a = δθ 0 min (A+1) α |ĝ 0i (θ 0 )| , k = iteration index starting with 0, δθ 0 min = smallest initial change desired in a parameter.
The setting for α and γ above are not optimal in the asymptotic sense, but are adapted to finite iteration settings. In practice, one of the drawbacks of SPSA is that one has to find good values for a and c, as both affect the performance of the algorithm Spall (2003, pp. 165-166) (Altaf et al. 2011;Shen et al. 2012;Radac et al. 2011;Easterling et al. 2014;Taflanidis and Beck 2008). However, for c, we have a tangible measure, which is the output measurement error (Spall 1998a), to select a proper value up front. If the function response is noiseless, c is usually not a critical parameter. On the other hand, a is more problematic, because no clear measure exists. It is possible to work with δθ 0 min instead of a, but a priori assignment of its value is still non-trivial if little is known about the function that we are trying to optimize.
A larger value of a generally produces better results compared to a smaller value of a. However, this also increases the chance that the optimization diverges to a worse solution than the starting point. Very often, the user of SPSA has to find as big a as possible that would not cause divergence.
To avoid divergence, an adaptation called "blocking" exists (Easterling et al. 2014;Spall 1998a) in which the objective values at θ k is evaluated in addition to the two perturbations. If the new objective function value is "significantly worse" than the current objective function value, the updating of θ k does not happen. The extra function evaluation at each iteration increases the cost of iteration by 33 %. In addition, a problem dependent threshold parameter to block the θ k update needs to be set up by the user.
Another way to mitigate divergence is to modify the gradient approximation ĝ k by "scaling" and "averaging" (Andradóttir 1996;Xu and Wu 2013). However, the methods proposed in the literature require set up of additional threshold parameters critical to their performance. Furthermore, their methods require additional gradient estimations per iteration.
Stochastic Gradient Descent (SGD) methods use noisy information of the gradient of the objective functions. On the other hand, Stochastic Approximation methods such as FDSA and SPSA only uses measurement of noisy objective values. Therefore, adaptive determination of step sizes based on (approximate) gradients and inverse Hessians in SGD literature (such as in Zeiler (2012), Bottou (2010)) may not be directly applicable to or feasible in SPSA. Convergence conditions also differ between the two. Although this does not exclude the possibility of successful import of ideas from SGD literature, in this paper, we will not delve into this direction.
This paper provides a solution to determine the appropriate values of a by introducing an adaptive scheme as discussed in "Adaptive initial step sizes" section. It does not require any additional objective function evaluations per iteration nor extra problem dependent parameters to set up.

Adaptive initial step sizes
To remedy the sensitivity to a, we propose an adaptive stepping algorithm. At the end of each iteration k, we perform the adjustment described in Algorithm 1.

Algorithm 1 Adaptive Initial
Step The condition requires that at least one of the two parameter perturbations produce a better (smaller) measurement of the objective function than that of initial guess of parameters θ 0 to proceed without modifying a. Therefore, at each iteration k, the smaller of the two measurements of the objective function values of perturbed parameters is compared to that of the initial value at iteration k = 0. If the measurements of the objective values of the perturbed parameters are larger, θ k is reset to θ b , which is the point that gave the minimum in the history of iteration and a is reduced to half of its previous value. A pseudocode of the proposed SPSA with the adaptive initial step is shown in Algorithm 2. The difference between the standard SPSA and our SPSA is in line 10.

Comments on convergence
Currently available theories of stochastic algorithms are almost all based on asymptotic properties with k → ∞, and SPSA is no exception. For given conditions Spall (2003, p. 183), SPSA is proven to converge to a local optima almost surely. However, under limited function evaluation budget, we frequently encounter situations in which SPSA returns worse solution than the initial i.e. divergence. The method we propose is a practical remedy conceived in a finite k setting. We will show, in the next section, its effectiveness empirically via numerical experiments with k in the order of 10 3 .
For θ k to converge to the optimal solution θ * in infinite steps, the following conditions are required for a k and c k (Spall 1992 a k c k 2 < ∞. With Algorithm 1, ∞ k=0 a k = ∞ is not guaranteed. For example, if the reduction of a happens in every iteration k, the sum is convergent. In practice, the numbers of function evaluations are finite, and reductions of a are expected to happen only a limited number of times. Therefore, this violation is expected to pose little problem. The intention of the proposed method is not to modify the asymptotic convergence rate of the original SPSA algorithm Spall (2003, pp. 186-188). The adaptive step takes place only if it is suspected that the objective value has become larger than at the starting point θ 0 . The probability of Algorithm 1 taking place is expected to go to zero under reasonable signal-to-noise ratio as f (θ k ) decreases. The worst situation that can happen is that the every perturbation c k k produces worsening moves and no improvement is obtained compared to the starting point θ 0 . In "Computational results" section, we will confirm empirically what we have described about the convergence in finite k settings (k ∼ 10 3 ).
Another reason to take the objective value at the starting point as the threshold value to judge divergence is that if we update this value with y(θ k ), where k > 0, we may risk picking a point that is too low due to the noise incurred in the measurement y. This in turn inhibits further improvement of θ k for lower objective values.
In the following section, the smallest output of mathematical functions will be sought using the standard SPSA and our adaptive initial stepping SPSA. This will show the sensitivity of the function value in the final iteration to the initial step size δθ 0 min and so the sensitivity to a, and how the adaptive initial stepping substantially mitigates the difficulty to find the proper initial perturbation magnitude.

Computational results
In this section, we will compare the original SPSA and our modified SPSA as described in Algorithm 2 using 10 analytical test functions and a parameter estimation example of a nonlinear dynamic system.

Test functions
To see the effect of the new adaptive stepping algorithm in SPSA, the minimum points of ten different mathematical test functions were sought. Except for Griewank function, the following conditions were applied. The functions' responses were minimized from arbitrary starting points θ 0 ∈ [−2, 2] D (D-dimensional product space with lower bound -2 and upper bound 2). If θ k = [θ k0 ,θ k1 , · · · ,θ ki , · · · ,θ k(D−1) ] T exceeded [−10, 10] in any of its D dimensions, that parameter was replaced by -10 if it was less than -10 or was replaced by 10 if it was larger than 10. For Griewank function, it was randomly started from θ 0 ∈ [−120, 120] D . If θ k exceeded [−600, 600] in any of its D dimensions, that parameter was replaced by −600 if it was less than −600 or was replaced by 600 if it was larger than 600. For all ten functions, the iteration was stopped when 2000 evaluations of the objective function were reached. For convenience, we will label our proposed algorithm as "A_SPSA" and the standard SPSA as "SPSA".
The optimizations for each of the ten objective functions were started from 20 different starting points. After the 2000 iterations, the distributions of objective values were plotted with respect to δθ 0 min . Eleven different values of δθ 0 min between 1.0 × 10 −4 and 1.0 × 10 1 (up to 1.0 × 10 2 for Griewank) were used to make the plot. The dimensions of the functions were set to be 20, i.e. D = 20.
The definitions of the ten functions are given in the following. The Rosenbrock function is described as The Sphere function is described as The Schwefel function is described as (8) The Rastrigin function is described as The Skewed Quartic function Spall (2003, ex. 6.6) is described as where the matrix B in the Skewed Quartic function is a square matrix with upper triangular elements set to 1 and the lower triangular elements set to zero. The Griewank function is described as The Ackley function is described as The Manevich function is described as The Ellipsoid function is described as (11) (12) (13) (14) The Rotated Ellipsoid function is described as Each of Figs. 2, 3, 4, 5, 6, 7, 8, 9, 10 and 11 show three different cases of noisy measurements of the outputs. The subfigures (a) have no noise added, subfigures (b) and (c) have Gaussian noise added to the true output with standard deviation σ of 0.1 and 1.0 respectively. In all the three noise levels of the ten functions, c = 0.2 was used. A general trend observed from the figures is that when the initial step size is large, the original SPSA tends to diverge to big objective values. The SPSA with the proposed initial step size reduction, on the other hand, effectively mitigates this divergence problem producing smaller objective values in general as the (a priori) initial step size is increased. This is because if the two function evaluations in the iteration are not smaller than the starting point value f (θ 0 ), the algorithm will reduce the step size (by halving a) and restart at θ b , which is the point that gave the smallest output in the history of iterations. However, note that the iteration index k in a k and c k is not reinitialized. For the ten functions tested, A_SPSA achieved its best performance when δθ 0 min was close to 10 or 100 for Griewank function. This indicates that one can simply set the minimum perturbation δθ 0 min close to the magnitude of the difference between upper and lower bound As mentioned earlier, the value for c is important when the measurements of y contain noise. Figure 12  narrower than that of the standard SPSA. On the other hand, the choice of δθ 0 min (and therefore a) is much easier for A_SPSA. We can, for example, let δθ 0 min ≃ min(U − L), where min(U − L) is the minimum difference between upper and lower bounds of the domain of parameter vector θ. In practice, it is better to scale all the input dimensions to fall in similar or equal intervals. Figure 13 shows the results of optimizing the Rosenbrock and Rastrigin functions using three different values of multiplication factor of a: 0.1, 0.5, and 0.9. The difference in multiplication factor does not change the general trend that larger δθ 0 min produces better results and that divergence does not occur. One could tune the value of the multiplication factor, but the default value of 0.5 that we showed in the Algorithm 1 generally produces satisfactory results compared to other values of multiplication factors between 0 and 1. The Fig. 13 (b) also shows that δθ 0 min ≃ min(U − L) may not be an optimal setting since smaller value δθ 0 min ≃ 10 −1.5 is shown to produce better optimization results when the reduction rate is slow at 0.9. This implies that in a bumpy (highly multimodal) function like Rastrigin, the slow decrease in a can adversely affect the minimization of the objective value by a large number of resets to θ b . The opposite is true with Rosenbrock function in (a), in which the slow reduction factor 0.9 gave the best result at δθ 0 min ≃ 10 1 .
For all the mathematical functions tested in this paper, optimization using SPSA diverges almost surely if the δθ 0 min is large. However, A_SPSA and SPSA give closely matching results when the initial step sizes are relatively small (i.e., the left hand side of the plots in Figs. 2,3,4,5,6,7,8,9,10 and 11). This is because, in cases that divergence does not happen, the adaptation of a does not take place in A_SPSA and therefore SPSA and A_SPSA have identical behavior. This is a confirmation that Algorithm 1 does not alter, in any significant way, the finite sample convergence characteristics of the original SPSA when the divergence does not manifest.

Nonlinear dynamics example
We consider a parameter estimation problem with Lorenz attractor. Its nonlinear dynamics is described as We seek to identify the system parameters θ = [s, r, b] by minimizing the one-time-stepahead prediction error L k of the state x k+1 given the current state We use fourth-order Runge-Kutta method to obtain x k+1 . Let us denote x k+1 as one-time-step-ahead prediction given by the estimated system with parameters θ k . Then, we can define the prediction error as Thus, the optimization to be solved is The index k above is the same as the index k in the SPSA algorithms. So the SPSA iteration proceeds along with the time steps of the dynamic system to compute L k .
For this problem, we set the parameter space as three-dimensional product space = [0, 500] 3 . The initial state is x 0 = [2, 3, 4] T . The initial guess (starting point) of the parameter set θ 0 is a random pick from . Figure 14 show the box plots of final L k when started from different values of δθ 0 min . The smallest median of final L k is obtained at δθ 0 min = 10 for SPSA and δθ 0 min = 100 and 1000 for A_SPSA. The best medians of final L k obtained for A_SPSA (5.62 × 10 −15 ) is smaller compared to that of SPSA (3.10 × 10 −13 ). However, both SPSA and A_SPSA had some runs that did not converge to the above mentioned near-zero L k values even at these δθ 0 min .
Again, for A_SPSA, the best setting were obtained when δθ 0 min was set to large values near the order of magnitude of the distance between upper and lower bound of the domain, while for SPSA, the best δθ 0 min was at an interior value between 10 −3 and 10 3 . Figure 15 shows the trajectory of the reference Lorenz attractor and the simulation of the Lorenz attractor whose system parameters s, r, and b were successfully identified by A_SPSA. The time t is run from 0 to 20 starting from the same initial condition used in the identification. The figure shows excellent match. Figure 16 shows the box plots of parameters estimated by A_SPSA and SPSA starting at their best δθ 0 min settings. The corresponding statistics are shown in Tables 1 and 2. The boxes appear collapsed as single horizontal lines at medians since the spaces between first quartiles and third quartiles are very narrow. Some non-converging cases are visible as dots on the figure. The figure and the tables show that the parameter estimates are more consistent from run to run in A_SPSA than that of SPSA as A_SPSA has narrower first and third quartile differences.

Conclusion
With the adaptive initial step algorithm, one can avoid divergence in SPSA iterations. Moreover, with a large initial step size, the SPSA algorithm with the adaptive initial step algorithm was able to find equal or better solutions compared to the original SPSA for all the ten mathematical function minimization problems that we have tested. In the nonlinear dynamics example, the new algorithm was able to find system parameters more precisely. The proposed method may not eliminate the need of tuning the parameters of SPSA algorithms, but it facilitates the process by eliminating the risk of solution divergence and reducing the trial-and-error effort. Further testing of the algorithm with different test functions, noise distributions, and industrial use-cases would be beneficial. The improvement proposed in this paper is expected to be valuable when the objective functions are costly to evaluate or if the algorithm is employed inside another algorithm such as machine learning or target tracking, for manual tuning of the parameters would be cumbersome in such cases. As a future work, it would be beneficial to investigate under what conditions the probability of the proposed adaptation (i.e. going into ifbranch in Algorithm 1) happening tends to zero as iteration k tends to infinity.
Authors' contributions KI has conceived the Algorithm 1, conducted the numerical experiments, and has written this manuscript. TD has provided guidance by supervising the progress of the research, and has reviewed the manuscript for intellectual feedback and editing.