Parallel Simultaneous Perturbation Optimization

TL;DR

This paper introduces PSPO, a parallelized stochastic optimization algorithm that efficiently utilizes high-performance computing to optimize expensive, noisy black-box functions, demonstrated through epidemiological modeling.

Contribution

The paper extends the simultaneous perturbation stochastic approximation algorithm to fully leverage parallel computing, reducing convergence time and adapting step sizes automatically.

Findings

PSPO converges faster than traditional methods.

It effectively utilizes parallel computing resources.

Demonstrated on a real epidemiological model.

Abstract

Stochastic computer simulations enable users to gain new insights into complex physical systems. Optimization is a common problem in this context: users seek to find model inputs that maximize the expected value of an objective function. The objective function, however, is time-intensive to evaluate, and cannot be directly measured. Instead, the stochastic nature of the model means that individual realizations are corrupted by noise. More formally, we consider the problem of optimizing the expected value of an expensive black-box function with continuously-differentiable mean, from which observations are corrupted by Gaussian noise. We present Parallel Simultaneous Perturbation Optimization (PSPO), which extends a well-known stochastic optimization algorithm, simultaneous perturbation stochastic approximation, in several important ways. Our modifications allow the algorithm to fully take advantage of parallel computing resources, like high-performance cloud computing. The resulting PSPO algorithm takes fewer time-consuming iterations to converge, automatically chooses the step size, and can vary the error tolerance by step. Theoretical results are supported by a numerical example. To demonstrate the performance of the algorithm, we implemented the algorithm to maximize the pseudo-likelihood of a stochastic epidemiological model to data of a measles outbreak.

Equations

$${\theta }_{k+1}={\theta }_k-a_k{\left({\Pi }_{\mathcal{P}}({\bar{H}}_k)\right)}^{-1}{\hat{\mathbf{g}}}_k({\theta }_k),$$

$${\bar{H}}_k=\frac{k}{k+1}{\bar{H}}_{k-1}+\frac{1}{k+1}{\hat{H}}_k,$$

$${\hat{\mathbf{g}}}_k({\theta }_k)=\frac{y_k^{(+)}-y_k^{(-)}}{2c_k}{\Delta }_k^{-1},$$

$${\hat{\mathbf{g}}}_k({\theta }_k)=\frac{y_k^{(+)}-y_k}{c_k}{\Delta }_k^{-1}.$$

$${\hat{H}}_k=\frac{1}{2}\left[\frac{\delta G_k}{2{\tilde{c}}_k}{\tilde{\Delta }}_k^{-1}+{\left(\frac{\delta G_k}{2{\tilde{c}}_k}{\tilde{\Delta }}_k^{-1}\right)}^T\right],$$

$$\delta G_k={\hat{\mathbf{g}}}_k({\theta }_k+{\tilde{c}}_k{\tilde{\Delta }}_k)-{\hat{\mathbf{g}}}_k({\theta }_k-{\tilde{c}}_k{\tilde{\Delta }}_k).$$

$${\mathbf{g}}^i=\frac{⟨\mathbf{g}(\theta ),{\Delta }_i⟩}{\Vert {\Delta }_i{\Vert }^2}{\Delta }_i.$$

$${\hat{\mathbf{g}}}^i\approx p{\mathbf{g}}^i=⟨\mathbf{g}(\theta ),{\Delta }_i⟩{\Delta }_i\approx ⟨\hat{\mathbf{g}}(\theta ),{\Delta }_i⟩{\Delta }_i.$$

$${\hat{\mathbf{g}}}^T{\Delta }_i\approx \frac{\delta f_i}{c},$$

$$\delta f_i=L(\theta +c{\Delta }_i)-L(\theta ),$$

$${\hat{\mathbf{g}}}^T\Delta \approx \frac{1}{c}\left[\begin{array}{cccc}\delta f_1& \delta f_2& \cdots & \delta f_M\end{array}\right],$$

$$\Delta =\left[\begin{array}{cccc}{\Delta }_1& {\Delta }_2& \cdots & {\Delta }_M\end{array}\right].$$

$$\hat{\mathbf{g}}\approx \frac{1}{c}{\left(\Delta {\Delta }^T\right)}^{-1}\Delta {\left[\begin{array}{cccc}\delta f_1& \delta f_2& \dots & \delta f_M\end{array}\right]}^T.$$

$$\hat{\mathbf{g}}\approx \frac{1}{c}\Delta {\left({\Delta }^T\Delta \right)}^{-1}{\left[\begin{array}{cccc}\delta f_1& \delta f_2& \dots & \delta f_M\end{array}\right]}^T.$$

$${\Delta }_0=\left[\begin{array}{c}{a}_1\\ a_2\\ a_3\end{array}\right],$$

$${\Delta }^{(3)}=\left[\begin{array}{ccc}-a_1& a_1& a_1\\ a_2& -a_2& a_2\\ a_3& a_3& -a_3\end{array}\right].$$

$${\Delta }^{(k+1)}=\left[\begin{array}{cccc}-a_1& \dots & a_1& a_1\\ ⋮& \ddots & ⋮& ⋮\\ a_k& \dots & -a_k& a_k\\ a_{k+1}& \dots & a_{k+1}& -a_{k+1}\end{array}\right].$$

$$\mathrm{rank}({\Delta }^{\prime (k)})=\mathrm{rank}\left(\left[\begin{array}{ccc}-a_1& \dots & a_1\\ ⋮& \ddots & ⋮\\ a_k& \dots & -a_k\\ a_{k+1}& \dots & a_{k+1}\end{array}\right]\right)=k.$$

$$\sum _{i=1}^k{\alpha }_i{c}_i+{\alpha }_{k+1}{\Delta }_{k+1}=0,$$

$$\sum _{i=1}^{k+1}{\alpha }_i-2{\alpha }_j=0,j=1,\cdots ,k+1$$

$$\Delta {\Delta }^T=Q\Lambda Q^T,$$

$$Q=\left[\begin{array}{cccc}{\mathbf{v}}_1& {\mathbf{v}}_2& \dots & {\mathbf{v}}_p\end{array}\right],$$

$${\left(\Delta {\Delta }^T\right)}^{-1}=Q{\Lambda }^{-1}{Q}^T.$$

$$p$$

$$=\sum _{i=1}^M\left[\left(\sum _{j=1}^p{\mathbf{v}}_j^T{\Delta }_i{\mathbf{v}}_j^T\right)\left(\sum _{k=1}^p\frac{1}{{\lambda }_k}{\mathbf{v}}_k{\Delta }_i^T{\mathbf{v}}_k\right)\right]$$

$$=\sum _{i=1}^M\left[\sum _{j=1}^p\frac{1}{{\lambda }_j}{\mathbf{v}}_j^T{\Delta }_i{\Delta }_i^T{\mathbf{v}}_j\right]$$

$$=\sum _{j=1}^p\frac{1}{{\lambda }_j}\sum _{i=1}^M\left[{\mathbf{v}}_j^T(1-2{\mathbf{e}}_i)(1-2{\mathbf{e}}_i{)}^T{\mathbf{v}}_j\right]$$

$$=(M-4)\left(1^T{\left(\Delta {\Delta }^T\right)}^{-1}1\right)+4\mathrm{Tr}{\left(\Delta {\Delta }^T\right)}^{-1}$$

$$=\mathrm{Tr}\left(\left((M-4)11^T+4I\right){\left(\Delta {\Delta }^T\right)}^{-1}\right).$$

$$\mathrm{Tr}\left({\left(\Delta {\Delta }^T\right)}^{-1}\right)\le \frac{p}{4}.$$

$$\mathbb{E}\left[\delta f_i\right]={\mu }_i-\mu .$$

$${\mathbf{g}}^T{\Delta }_i=\frac{{\mu }_i-\mu }{c}.$$

$$\mathbb{E}\left[\delta f_i\right]=c{\mathbf{g}}^T{\Delta }_i.$$

Full text

1. Introduction

Stochastic optimization is of core practical importance in many fields of science and engineering. For instance, epidemiological systems involve individuals subject to exogenous disturbances, which emphasizes the need for models and methods capable of dealing with stochasticity. While many optimization methodologies exist for deterministic systems, these algorithms can give misleading results when applied to stochastic systems. Closed-form solutions do not generally exist for stochastic optimization problems, thus we seek an iterative algorithm which guarantees convergence to the local optimal point.

Many approximation algorithms have been developed to solve a wide variety of deterministic or stochastic problems. The steepest descent method (Wright and Nocedal, 1999) is the most prominent iterative method for optimizing a complex objective function. The gradient-based algorithms, such as Robbins-Monro (Robbins and Monro, 1951), Newton-Raphson (Fröberg, 1969), and neural network back-propagation (Rumelhart et al., 1985), rely on direct measurements of the gradient of the objective function with respect to the optimization parameter. But, in many cases the gradient of the loss function is not available. This is a common occurrence, for example, in complex systems, such as the optimization problems given in (Tsilifis et al., 2017; Alaeddini and Klein, 2017), the exact functional relationship between the loss function value and the parameters is not known, and the loss function is evaluated by measurements on the system or by running simulation. A review on the main areas of optimization via simulation can be found in (Fu, 1994; Hong and Nelson, 2009; Swisher et al., 2000).

Some algorithms have been developed specifically to optimize stochastic black-box cost functions. Of these, the Kiefer-Wolfowitz algorithm (Kiefer et al., 1952) is perhaps the most well known. This algorithm estimates the gradient from noisy measurements using finite differencing. Another well known stochastic optimization algorithm in the case of high dimensional problems is Simultaneous Perturbation Stochastic Approximation (SPSA), which is an approximation algorithm based on simultaneous perturbation (Spall, 1992). Later, J. C. Spall presented a second-order variant of SPSA (Spall, 2000). The SPSA algorithm is used extensively in many different areas, e.g.  signal timing for traffic control (Ma et al., 2013), and some large scale machine learning problems (Byrd et al., 2011). The convergence of this algorithm to the optimal value in the stochastic almost sure sense makes it suitable in many applications.

The stochastic nature of some complex models, e.g.  epidemiological models, means that each set of parameters maps to a distribution of outcomes, from which each sample (model run) can take several hours to obtain. In the quick-to-evaluate deterministic model setting, almost any classical optimization algorithm, such as steepest descent or Newton-Raphson, can be used. However, care must be taken when applying deterministic methods to stochastic objective functions, as the inherent noise causes unexpected behavior. In this paper, we introduce PSPO, an algorithm for optimization of stochastic objective functions. Despite many advantageous properties of the SPSA algorithm, mentioned earlier, it is a serial algorithm that evaluates the (stochastic) function a few points at a time, which results in low convergence rate. Researchers have looked at ways of enhancing the convergence of the SPSA algorithm, e.g.  iterate averaging is an approach aimed at achieving higher convergence rate in a stochastic setting (Polyak and Juditsky, 1992), using deterministic parameter perturbations instead of random perturbations (Bhatnagar et al., 2003), and more (Kocsis and Szepesvári, 2006; Spall, 2009).

The PSPO algorithm, introduced in this paper, takes advantage of fundamentally-parallel resources like high-performance cloud-based computing. Our method is appropriate for problems with very noisy gradients. We also derive a relationship between the number of parallel rounds of computation and error tolerance of the gradient for each iteration. The main contribution of this paper is introducing a stochastic optimization algorithm which can be easily implemented on parallel computers, and to provide the minimum number of parallel computers in order to have an upper-bound on the error for each iteration. We demonstrate that PSPO works well in practice and compares favorably to the conventional simultaneous perturbation optimization algorithm (SPSA). The rest of the paper is structured as follows. In Section 2, we give a review of the simultaneous perturbation optimization algorithm. The Parallel Simultaneous Perturbation Optimization (PSPO) algorithm is presented in Section 3. The numerical simulations are given in Section 4, and Section 5 concludes the paper.

2. Preliminaries

Consider the problem of minimizing a cost function $L(\theta):\mathbb{R}^{p}\rightarrow\mathbb{R}$. Spall (Spall, 1992, 2000, 2009) presented an efficient stochastic algorithm called Simultaneous Perturbation Stochastic Approximation (SPSA). The SPSA algorithm estimates the gradient and the Hessian matrix by finite difference. This algorithm basically consists of two parallel recursions for estimating the optimization parameter, $\theta$, and the Hessian matrix, $H(\theta)$. The first recursion is a stochastic equivalence of the Newton-Raphson algorithm, and the second one estimates the Hessian matrix. These two recursions are

[TABLE]

where $a_{k}$ is a positive scalar factor, $\mathcal{P}$ denotes the set of all positive definite matrices, and $\Pi_{\mathcal{P}}(\cdot)$ is the projection into the admissible set $\mathcal{P}$. Here, $\hat{\mathbf{g}}_{k}$ and $\hat{H}_{k}$ are the estimated gradient and Hessian at iteration $k$. The SPSA approach for estimating the $\hat{\mathbf{g}}_{k}(\theta)$ and $\hat{H}_{k}(\theta)$ follows.

Let $\Delta_{k}\in\mathbb{R}^{p}$ be vectors of $p$ mutually independent zero-mean random variables satisfying the condition of $E\{\Delta_{k}^{-1}\}$ be bounded. An admissible distribution is a Bernoulli $\pm 1$ distribution. The positive scalars $c_{k}$ are chosen such that they usually get smaller as $k$ gets larger. The two-sided estimate of the gradient at iteration $k$ is given by:

[TABLE]

where $y_{k}^{(\pm)}$ are the noisy measurements of the cost function at $\theta_{k}\pm c_{k}\Delta_{k}$. Note that in (2), $\Delta_{k}^{-1}$ is the element-wise inverse of $\Delta_{k}$. In the case of second order SPSA, it is suggested to use a one-sided gradient approximation, given by:

[TABLE]

Now let $\tilde{\Delta}_{k}\in\mathbb{R}^{p}$ be vectors of $p$ mutually independent zero-mean random variables satisfying the same condition of $\Delta_{k}$. The positive scalars $\tilde{c}_{k}$ are also chosen such that they get smaller as $k$ gets larger. The numerical value of $\tilde{c}_{k}$ is suggested to be chosen smaller than $c_{k}$. The following formula gives a per-iteration estimate of the Hessian matrix.

[TABLE]

where

[TABLE]

The practical implementation details of the algorithm is presented in (Spall, 2000).

3. Parallel Simultaneous Perturbation Optimization: PSPO

While SPSA is an efficient optimization algorithm for high dimensional problems, it requires many iterations to converge, particularly for high-noise problems. In the case of high uncertain problems, the computation of the gradient and Hessian from noisy objective function evaluations benefits from more evaluations. The idea here is using parallel computing of the gradient and Hessian with different perturbation vectors, and use the estimated gradient and Hessian in a conjugate gradient type optimization algorithm.

3.1. Gradient Estimation

The gradient vector at each iteration can be obtained by performing multiple evaluations of (3) for different values of perturbation. Let $\mathbf{g}$ be the gradient at a point $\theta$, and $\mathbf{g}^{i}$ represents the directional gradient in $\Delta_{i}$ direction at this point. Then we know that:

[TABLE]

Let $\hat{\mathbf{g}}^{i}$ be the estimated gradient using the one-sided gradient approximation (3) and perturbation vector $\Delta_{i}$. If $\Delta_{i}$ is a Bernoulli $\pm 1$ vector, we have

[TABLE]

Then,

[TABLE]

[TABLE]

and then,

[TABLE]

where,

[TABLE]

If $\{\Delta_{i}\}_{1\leq i\leq M}$ span $\mathbb{R}^{p}$, then $\Delta\Delta^{T}$ is invertible. So, the least square estimation of $\hat{\mathbf{g}}$ is given by:

[TABLE]

On the other hand, if $M<p$, then (6) is an under-determined system which has non-unique solutions. Then, the solution of the minimum Euclidean norm, $\displaystyle\|\hat{\mathbf{g}}\|_{2}$, among all solutions is given by:

[TABLE]

The Parallel Simultaneous Perturbation (PSP) algorithm for estimating the gradient at a given point is given in Algorithm 1. In order to improve the efficiency of this algorithm, we can compute $L(\mathbf{\theta})$ once, out of the loop. Doing this, we need $M+1$ function evaluations for estimating the gradient.

In Algorithm 1, we suggest an efficient strategy to choose independent $\Delta$ vectors. First, we sample $\Delta_{0}$ from a Bernoulli $\pm 1$ distribution. In round $j$, $1\leq j\leq p$, we need to switch the sign of the $j$th element in $\Delta_{0}$ to generate $\Delta_{j}$. The same procedure repeats for the next $p$ rounds of computations.

Lemma 3.1.

Given the suggested procedure in Algorithm 1 for choosing $\Delta_{i}$ vectors, the vectors $\{\Delta_{j}\}_{1\leq j\leq p}$ spans $\mathbb{R}^{p}$ for all $p\neq 2$.

Proof.

The $\{\Delta_{j}\}_{1\leq j\leq p}$ vectors span $\mathbb{R}^{p}$ if and only if $\Delta^{(p)}=\begin{bmatrix}\Delta_{1}&\Delta_{2}&\cdots&\Delta_{p}\end{bmatrix}$ is full rank. For $p=1$, it is trivial that $\operatorname{rank}(\Delta^{(p)})=p$. For $p=3$, let assume

[TABLE]

then

[TABLE]

We know that $\det{\Delta^{(3)}}=4a_{1}a_{2}a_{3}\neq 0$, for all $a_{i}\neq 0$. Since $a_{i}=\pm 1$, $\Delta^{(3)}$ is full rank. Now, assume we know that $\operatorname{rank}(\Delta^{(k)})=k$. We need to prove that $\operatorname{rank}(\Delta^{(k+1)})=k+1$. From the procedure in Algorithm 1, we have

[TABLE]

Since adding a row to a matrix does not change the rank of that matrix, it is clear that

[TABLE]

Now, we need to prove that $\Delta_{k+1}$ (the last column of $\Delta^{(k+1)}$) is not in the span of the columns of the matrix $\Delta^{\prime(k)}$. Let $\mathbf{c}_{i}$ represent the columns of $\Delta^{\prime(k)}$. If $\Delta_{k+1}$ is not independent of the columns, $\mathbf{c}_{i}$, then there should exists a nonzero vector, say $\mathbf{\alpha}\neq\mathbf{0}$, such that

[TABLE]

which results in

[TABLE]

It is easy to show that there is not any real valued $\mathbf{\alpha}$ which satisfies these conditions. Thus, $\Delta_{k+1}$ is independent of the columns of $\Delta^{(k+1)}$, and adding the $\Delta_{k+1}$ column, which generates $\Delta^{(k+1)}$, increases the rank of the matrix by one. Therefore, $\operatorname{rank}\Delta^{(k+1)}=k+1$, which proves the lemma by induction. ∎

Lemma 3.2.

Given the suggested procedure for choosing $\Delta$ vectors, if $\Delta_{0}=\mathbf{1}$, where $\mathbf{1}$ is a vector that all its elements are one, then for all $p\neq 2$ and $M\geq 4$, then $\operatorname{Tr}(\Delta\Delta^{T})^{-1}\leq\frac{p}{4}$ and $\operatorname{Tr}(\Delta^{T}\Delta)^{-1}\leq\frac{p}{4}$.

Proof.

Using $\Delta_{0}=\mathbf{1}$, we have $\Delta_{i}=\mathbf{1}-2\mathbf{e}_{i}$. Since $(\Delta\Delta^{T})$ is a symmetric full rank matrix, we can write it as its eigenvalues-decomposition, given by:

[TABLE]

where,

[TABLE]

and $\mathbf{v}_{i}$ is the eigenvector corresponding to eigenvalue $\lambda_{i}$. Then,

[TABLE]

We know

[TABLE]

Since $M-4\geq 0$ and $\mathbf{1}\mathbf{1}^{T}\succeq 0$, then $\left((M-4)\mathbf{1}\mathbf{1}^{T}+4I\right)\left(\Delta\Delta^{T}\right)^{-1}\succeq 4\left(\Delta\Delta^{T}\right)^{-1}$. Therefore,

[TABLE]

Proving $\operatorname{Tr}(\Delta^{T}\Delta)^{-1}\leq\frac{p}{4}$ is very similar to the above proof. Due to this similarity, the proof is skipped here. ∎

Now let $y=L(\theta)$ and $y^{i}=L(\theta+c\Delta_{i})$ represent the noisy measurements at $\theta$ and $\theta+c\Delta_{i}$ respectively. Assume these measurements can be expressed as normal random variables as $y\sim\mathcal{N}(\mu,\sigma^{2})$ and $y^{i}\sim\mathcal{N}(\mu_{i},\sigma^{2})$. It is assumed that the noise variance is constant over the whole state space, $\theta$.

Lemma 3.3.

Given $y\sim\mathcal{N}(\mu,\sigma^{2})$ and $y^{i}\sim\mathcal{N}(\mu_{i},\sigma^{2})$, if $M\geq p$, then the expected value of the gradient from Algorithm 1 is equal to the gradient $\mathbf{g}$.

Proof.

Based on our assumption,

[TABLE]

For small $c$, we have

[TABLE]

Then

[TABLE]

If $M\geq p$, using (7), we have

[TABLE]

∎

Theorem 3.4.

If the rounds of computation, $M$, satisfies

[TABLE]

then the error in estimated gradient from Algorithm 1 is bounded by

[TABLE]

Proof.

Using the assumptions of $y\sim\mathcal{N}(\mu,\sigma^{2})$ and $y^{i}\sim\mathcal{N}(\mu_{i},\sigma^{2})$, the function difference $\delta f_{i}(\theta)$ can be expressed as:

[TABLE]

Using (7) for estimating the gradient and Lemma 3.3, the error in gradient, $(\hat{\mathbf{g}}-\mathbf{g})$, is a normal random variable given by:

[TABLE]

where

[TABLE]

After $M$ rounds of computation of $\hat{\mathbf{g}}$, Chebyshev inequality implies that

[TABLE]

Using the results of Lemma 3.2, we have

[TABLE]

Therefore, if

[TABLE]

then

[TABLE]

Therefore, $\displaystyle\mathbb{E}\left[\|\hat{\mathbf{g}}-\mathbf{g}\|\right]\leq\epsilon\,.$ ∎

3.2. Hessian Estimation

Similar to every other Newton-based algorithm, we need to estimate the Hessian (matrix of second derivatives). In this work, we are suggesting estimation of the reduced Hessian, $H_{\mathbf{p}}$, instead of full Hessian, $H$. The reduced Hessian, $H_{\mathbf{p}}$, only conveys the information of the effect of the true Hessian in a specific $\mathbf{p}$ direction. The reduced Hessian with respect to a given vector $\mathbf{p}$ satisfies

[TABLE]

Before introducing how we can compute the reduced Hessian, we present a method to estimate the true Hessian in the following lemma.

Lemma 3.5.

Suppose $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a function taking as input a vector $\mathbf{x}\in\mathbb{R}^{n}$. Given a set of small $n$ orthogonal vectors $\{\mathbf{d}_{1},\mathbf{d}_{2},\cdots,\mathbf{d}_{n}\}$, then the true Hessian at a point $\mathbf{x}$ is given by

[TABLE]

where

[TABLE]

The proof of this lemma can be found in the appendix (see Appendix).

Definition 3.6.

The reduced Hessian with respect to a given vector $\mathbf{d}$ is defined as:

[TABLE]

Theorem 3.7.

The reduced Hessian defined in (13) satisfies the following condition:

[TABLE]

Proof.

Suppose $\{\mathbf{d}_{1},\mathbf{d}_{2},\cdots,\mathbf{d}_{n-1}\}$ are a set of orthogonal vectors such that $\{\mathbf{d},\mathbf{d}_{1},\mathbf{d}_{2},\cdots,\mathbf{d}_{n-1}\}$ becomes a set of $n$ orthogonal vectors which spans the entire $\mathbb{R}^{n}$. Then using (12), the Hessian can be estimated as:

[TABLE]

Since all $\mathbf{d}_{i},i=1,\cdots,n-1$ are perpendicular to $\mathbf{d}$, we have

[TABLE]

On the other hand, from the definition in (13)

[TABLE]

∎

3.3. PSPO Algorithm

Algorithm 2 presents a step-by-step summary of the proposed approach in this paper. This algorithm is a conjugate gradient type algorithm which uses the PSP algorithm to estimate the gradient and Hessian at each iteration. The stopping criterion in this algorithm can be defined by the size of the gradient at the final point, the update size, or even the number of iterations.

Remark 3.1.

The PSPO algorithm gives the chance to set the error tolerance for each step. One is able to start a coarse search at first iterations and decrease the error tolerance as the algorithm gets closer to the optimal point (by increasing the number of the computation rounds).

Remark 3.2.

To estimate the Hessian in PSPO, at each iteration we compute a rank-one estimate of the Hessian with respect to the update direction in that iteration, which is called a reduced Hessian in this paper. In the gradient descent algorithms, the Hessian matrix scales the update vector based on the curvature of the objective function for each feature. The idea, here, is to estimate the curvature of the objective function for each feature at each iteration, but only in the direction of the update in that iteration. In fact, we only need to estimate the effect of the Hessian matrix in the update direction, and we are not interested on the effect of the Hessian matrix in other directions.

4. Simulation

4.1. Toy Example: Minimization of A Quadratic Function

To demonstrate the main contribution of this paper, we implemented the methodology on a variety of examples. We begin with a simple three-dimensional quadratic function to clearly illustrate the efficiency of the algorithm in the case we know the exact optimal solution. To empirically evaluate the proposed method, we investigated the efficiency of PSPO and SPSA algorithms on a stochastic optimization problem. We use the same initialization condition when comparing two optimization algorithms. We evaluate our proposed method on the following highly noisy objective function

[TABLE]

where, $w$ is a Gaussian noise $w\sim\mathcal{N}(0,3^{2})$. This convex objective function is suitable for comparison of different optimizers without worrying about local minimum issues. We rerun each optimization algorithm $200$ times for different initial conditions. According to Fig. 1a, we found that the PSPO algorithm converges faster than SPSA on average. Note that the maximum number of iteration set to $100$ for both algorithms. The efficiency of the two algorithms is also tested through running the algorithms for a fixed initial condition many times. Fig. 1b shows the result. In this example, the number of function evaluations per iteration for the PSPO algorithm is twice of the number of function evaluations per iteration for the SPSA algorithm. But, as it can be seen in Fig. 1b, in case we have access to supercomputers to run the stochastic simulations in parallel, then the PSPO algorithm converges faster.

4.2. Real Application: Epidemiological Model Calibration

To demonstrate the efficiency of the PSPO algorithm proposed in this paper, we implemented the methodology on a data set concerns a measles outbreak in a small town in Germany in 1861, which contains 188 infected individuals. This data set, called Hagelloch data set, is very popular in literature because of its completeness and depth of data. Fig. 2 gives the observed clinical data. The data are obtained from (Meyer et al., 2014). The objective, here, is finding a good model for the given epidemic data to investigate the properties of the disease spread. To estimate the model parameters, we use maximum likelihood estimate (MLE). Details on how we can formulate this problem as an optimization problem can be found in (Alaeddini and Klein, 2017).

We compare the efficiency of the PSPO algorithm with the conventional SPSA algorithm. To compare these two algorithms, we rerun each optimization algorithms 100 times. Fig. 3a shows the number of iterations required to converge. As it can be seen in this figure, PSPO algorithm on average needs fewer iterations compared with the conventional SPSA algorithm. The maximum number of iterations for both algorithms set to 30. Note that the number of iterations until convergence set to 30 if the algorithm does not converge after 30 iterations. In order to investigate the performance of the PSPO algorithm as the number of parallel rounds increases, we run the PSPO for this example with different number of parallel computation rounds. Fig. 3b shows the number of iterations for convergence as $M$ grows. In this example, again, the number of function evaluations per iteration for the PSPO is twice of the corresponding number for the SPSA algorithm. Having higher number of function evaluations per iteration in this example, we can still see in Fig. 3a that the number of total function evaluations for PSPO might be less than the total number of function evaluations for SPSA. However, in this paper we do not focus on reducing the number of objective function evaluations and only focus on the faster convergence in case of having access to parallel supercomputers.

5. Conclusions

In many stochastic optimization problems that noise presents, obtaining an analytical solution is hardly possible. In this paper, we described an algorithm for optimal parameter estimation in discretely observed stochastic model. We have introduced PSPO, which is a simple and computationally efficient algorithm for gradient-based optimization of stochastic objective functions. The comparison between SPSA and PSPO algorithm demonstrated the superiority of the PSPO algorithm in the number of iterations for convergence. Note that, comparing with SPSA, the PSPO algorithm with multiple parallel computations requires additional function evaluations per iteration. In case that we do not have access to powerful parallel computers, in order to prevent unnecessary long processing time per iteration, which results in long total processing time, the number of parallel computations need to be carefully computed using the level of the noise presented in the problem. We found the minimum number of parallel rounds of computation in order to have a bounded error in the estimated gradient. In some cases when we are dealing with a high SNR (signal to noise ratio) problem, a single round of computation might be enough for convergence, but in some other applications that the data are highly noisy the SPSA algorithm requires too many iterations to converge or even diverges. Our method is aimed towards highly uncertain objective functions.

Acknowledgements.

The authors thank Bill and Melinda Gates for their active support of this work and their sponsorship through the Global Good Fund. This work was performed at Institute for Disease Modeling with help from colleagues, especially Philip Eckhoff.

Appendix

Proof of Lemma 3.5

See 3.5

Proof.

We know that the Hessian matrix $H$ of $f$ is a square $n\times n$ matrix, defined as follows:

[TABLE]

where $\mathbf{g}_{i}$ is the value of the gradient in $\mathbf{e}_{i}$ direction, and can be approximated by

[TABLE]

The full gradient of $\mathbf{g}_{i}$ can also be approximated using the finite difference method as follows:

[TABLE]

Now, (15) can be written as

[TABLE]

Thus

[TABLE]

∎