# Accelerate iterated filtering

Dao<EMAIL_ADDRESS>[ Some Universitym1 and Another
Universitym2  Departments of Mathematics,
University of Mississippi, Oxford, Mississippi, USA

###### Abstract

In simulation-based inferences for partially observed Markov process models
(POMP), the by-product of the Monte Carlo filtering is an approximation of the
log likelihood function. Recently, iterated filtering [14, 13] has originally
been introduced and it has been shown that the gradient of the log likelihood
can also be approximated. Consequently, different stochastic optimization
algorithm can be applied to estimate the parameters of the underlying models.
As accelerated gradient is an efficient approach in the optimization
literature, we show that we can accelerate iterated filtering in the same
manner and inherit that high convergence rate while relaxing the restricted
conditions of unbiased gradient approximation. We show that this novel
algorithm can be applied to both convex and nonconvex log likelihood
functions. In addition, this approach has substantially outperformed most of
other previous approaches in a toy example and in a challenging scientific
problem of modeling infectious diseases.

accelerate iterated filtering,

accelerate inexact gradient method,

sequential Monte Carlo,

state space model,

parameter estimation,

###### keywords:

arXiv:0000.0000

##  1 Introduction

The last decade has seen a great increase in the use of simulation-based
inference where numerical approximations are based on either Markov chain
Monte Carlo or sequential Monte Carlo sampling. These approaches have become
popularized, in part, because of the increasing computational power and the
emergence of efficient stochastic optimization algorithms. On the Bayesian
paradigm, particle Markov chain Monte Carlo has been introduced and
popularized by Doucet and collaborators [1, 2, 31]. Similar ideas have been
developed previously [21, 7, 10, 16] but in different contexts than
simulation-based inferences. On the frequentist paradigm, [14, 13] have
introduced an original approach to perform simulation-based parameter
inference in POMP models by combining stochastic gradient approximation and
particle filtering. In this paper, we will focus on improving one of the most
popular algorithm of this class, namely, iterated filtering (IF). Iterated
filtering uses an approximation estimate of the gradient of the log likelihood
computed from particle filters while proposing an artificial perturbation
moves to update the parameters. This class of algorithm is attractive because
it enables routine simulation-based parameter inferences in general POMP
model, even in the cases of intractable likelihoods. Due to some interesting
theoretical properties [13, 15, 29], its applications range in various fields
such as biology, ecology, economics and engineering [24, 23, 20, 3, 5, 4].

Iterated filtering was later theoretically developed by Ionides et al., [13].
Recently, Lindström et al., [24] extended it to improve on numerical
performance while Doucet et al., [8] expanded it to include
filtering/smoothing with quite attractive theoretical properties. Ionides et
al., [15] generalized Lindström et al., [24]’s approach and combined the idea
with data cloning [22], developing a Bayes map iterated filtering with an
entirely different theoretical approach. Nguyen and Ionides, [30] revisited
the approach of Doucet et al., [8], using a different perturbation noise and
computed both the gradient and the Hessian. Similar to intractable likelihood
in the context of iterated filtering, Poyiadjis et al., [32], Nemeth et al.,
[26], Doucet et al., [8] showed that the gradient and Hessian information can
also be computed from particle filter. In the same line, manifold Langevin
Monte Carlo (mMALA) [12] exploits the Hessian information to simplify the
tedious tuning method while improving on convergence rate. However, this
relies on rather strong assumptions that the gradient, and Hessian information
of transition density and observation density can be sampled from. This is
quite unrealistic in many real world applications. We, therefore, followed the
formal approaches, based solely on very weak assumptions of being able to
sample from transition density and evaluate from observation density.
Motivated from the fact that the gradient and Hessian information can be
approximated using the first and the second moments [13, 8], we propose to use
such approximations in the context of accelerate iterated filtering. Ionides
uses score vector merely while Doucet includes the Hessian information for the
independent white noise, which is not quite useful in the context of iterated
filtering with natural random walk noise. Nguyen and Ionides, [30] proposed to
approximate the gradient and Hessian using random walk noise to efficiently
explore the mode of the likelihood. Other than exploiting approximations of
the Hessian under weak assumption, we chose an alternative approach. That is,
we apply the accelerate gradient approach to the approximation of the gradient
of the log likelihood for an effective estimation approach.

The key contributions of this paper are three folds. Firstly, we developed and
showed that accelerate iterated filtering algorithm converges using a general
non-increasing step size with bias approximation of the gradient. It is
simple, elegant, and generalizable to faster algorithms. Secondly, we proved
that it has a higher convergence rate in general convex and non-convex
conditions of the objective log likelihood. Finally, we showed substantial
improvements of the method on a toy problem and on a real world challenge
problem of vivax malaria model compared to previous simulation-based inference
approaches.

The paper is organized as follows. In the next section we introduce some
notations and we develop the framework of accelerate iterated filtering. In
Sections 3, we state the convergence of this approximation method to the true
maximum likelihood estimation by iterating and accelerating noisy gradient of
the log likelihood. We validate the proposed methodology by a toy example and
a challenging inference problem of fitting a malaria transmission model to
time series data in Section 4, showing substantial gains for our methods over
current alternatives. We conclude in Section 5 with the suggesting of the
future works to be extended. The proofs are postponed to the Appendix.

##  2 Background of simulation-based inferences

We are interested in a general latent variable model since this is an
ubiquitous model for applied sciences. Let $\mathcal{X}$ be a latent state
space with a density $q_{\theta}(x)$ parameterized by
$\theta\in\Theta=\mathbb{R}^{d}$, and let $\mathcal{Y}$ be an observation
space equipped with a conditional density $f_{\theta}(y|x)$. The observation
$y\in\mathcal{Y}$ are considered as fixed and we write the log-likelihood
function of the data $\ell(\theta)\overset{\triangle}{=}\log\int
q_{\theta}(x)f_{\theta}(y|x)dx$. We work with the maximum likelihood
estimator, $\hat{\theta}=\arg\max\ell(\theta)$ where $\ell(\theta)$ is
intractable but $f_{\theta}(y|x)$ can be evaluated, by using samples where
$f_{\theta}(y|x)$ is also intractable. This process often uses the first order
stochastic approximation [19], which involves a Monte Carlo approximation to a
difference equation,
$\theta_{m}=\theta_{m-1}+\gamma_{m}\nabla\ell(\theta_{m-1}),$ where
$\theta_{0}\in\Theta$ is an arbitrary initial estimate and
$\\{\gamma_{m}\\}_{m\geq 1}$ is a sequence of step sizes with ${\sum_{m\geq
1}\gamma_{m}=\infty}$ and ${\sum_{m\geq 1}\gamma_{m}^{2}<\infty}$. The
algorithm converges to a local maximum of $\ell(\theta)$ under regularity
conditions. The term $\nabla\ell(\theta)$, also called the score function, is
shorthand for the $\mathbb{R}^{d}$-valued vector of partial derivatives,
$\nabla\ell(\theta)=\frac{\partial\ell(\theta)}{\partial\theta}$.

Sequential Monte Carlo (SMC) approaches have previously been developed to
estimate the score function [32, 25, 6]. However, under the simulation-based
setting, which does not require the ability to evaluate transition densities
and their derivatives, these approaches are not applicable. As a result, [13],
[8] used an artificial dynamics approach to estimate the derivatives.
Specifically, [30] considers a parametric model consisting of a density
$p_{Y}(y;\theta)$ with the log-likelihood of the data $y^{*}\in\mathcal{Y}$
given by $\ell(\theta)=\log p_{Y}(y^{*};\theta)$. A stochastically perturbed
model corresponding to a pair of random variables $(\breve{\Theta},\breve{Y})$
having a joint probability density on $\mathbb{R}^{d}\times\mathcal{Y}$ can be
defined as $p_{\breve{\Theta},\breve{Y}}(\breve{\vartheta},\ y;\theta,\
\tau)=\tau^{-d}\kappa\left\\{\tau^{-1}(\breve{\vartheta}-\theta)\right\\}p_{Y}(y;\breve{\vartheta}).$
Suppose the following regularity conditions, identical to the assumptions of
[8]:

######  Assumption 1.

There exists $C<\infty$ such that for any integer $k\geq 1,1\leq i_{1},\
\ldots,\ i_{k}\leq d$ and $\beta_{1},\ \ldots,\ \beta_{k}\geq 1$,
$\int\left|u_{i_{1}}^{\beta_{1}}u_{i_{2}}^{\beta_{2}}\cdots
u_{i_{k}}^{\beta_{k}}\right|\kappa(u)\ du\leq C,$ where $\kappa$ is a
symmetric probability density on $\mathbb{R}^{d}$ with respect to Lebesgue
measure and $\Sigma=(\sigma_{i,j})_{i,j=1}^{d}$ is the non-singular covariance
matrix associated to $\kappa$.

######  Assumption 2.

There exist $\gamma,\ \delta,\ M>0,$ such that for all $u\in\mathbb{R}^{d}$,
$|u|>M\Rightarrow\kappa(u)<e^{-\gamma|u|^{\delta}}.$

######  Assumption 3.

$\ell$ is four times continuously differentiable and $\delta$ defined as in
Assumption 2. For all $\theta\in\mathbb{R}^{d}$, there exists $0<\eta<\delta,\
\epsilon,\ D>0,$ such that for all $u\in\mathbb{R}^{d}$,
$\mathcal{L}(\theta+u)\leq De^{\epsilon|u|^{\eta}},$ where $\mathcal{L}$ :
$\mathbb{R}^{d}\rightarrow\mathbb{R}$ is the associated likelihood function
$\mathcal{L}=\exp\ell$.

Under these regularity assumptions, [8] show that

$\left|\tau^{-2}\Sigma^{-1}{\mathbb{E}}\left(\breve{\Theta}-\theta\left|\breve{Y}=y^{*}\right.\right)-\nabla\ell\left(\theta\right)\right|<C\tau^{2}.$
(2.1)

These approximations are useful for latent variable models, where the log-
likelihood of the model consists of marginalizing over a latent variable, $X$,

$\ell(\theta)=\log\int{p_{{X},{Y}}(x,y^{*};\theta)\,dx}.$

In this case, the expectations in equation 2.1 can be approximated by Monte
Carlo importance sampling, as proposed by [13] and [8]. In [30], the POMP
model is a specific latent variable model with ${X}=X_{0:N}$ and
${Y}=Y_{1:N}$. A perturbed POMP model is defined to have a similar
construction to our perturbed latent variable model with
${\breve{X}}=\breve{X}_{0:N}$, $\breve{Y}=\breve{Y}_{1:N}$ and
$\breve{\Theta}=\breve{\Theta}_{0:N}$. [13] perturbed the parameters by
setting $\breve{\Theta}_{0:N}$ to be a random walk starting at $\theta$,
whereas [8] took $\breve{\Theta}_{0:N}$ to be independent additive white noise
perturbations of $\theta$. We take advantage of the asymptotic developments of
[8] while maintaining some practical advantages of random walk perturbations
for finite computations, so we use the construct $\breve{\Theta}_{0:N}$ as in
[30] as follows.

Let $Z_{0},\ \ldots,\ Z_{N}$ be $N+1$ independent draws from a density $\psi$.
[30] introduces $N+2$ perturbation parameters, $\tau$ and
$\tau_{0},\ldots,\tau_{N}$, and construct a process $\breve{\Theta}_{0:N}$ by
setting ${\breve{\Theta}_{n}=\theta+\tau\sum_{i=0}^{n}\tau_{i}Z_{i}}$ for
$0\leq n\leq N$. We later consider a limit where $\tau_{0:N}$ as fixed and the
scale factor $\tau$ decreases toward zero, and subsequently another limit
where $\tau_{0}$ is fixed but $\tau_{1:N}$ decrease toward zero together with
$\tau$. Let $p_{\breve{\Theta}_{0:N}}(\breve{\vartheta}_{0:N};\theta,\ \tau,\
\tau_{0:N})$ be the probability density of $\breve{\Theta}_{0:N}$. We define
the artificial random variables $\breve{\Theta}_{0:N}$ via their density,

$p_{\breve{\Theta}_{0:N}}(\breve{\vartheta}_{0:N};\theta,\ \tau,\
\tau_{0:N})=\\\
(\tau\tau_{0})^{-d}\psi\left\\{(\tau\tau_{0})^{-1}(\breve{\vartheta_{0}}-\theta)\right\\}\times\prod_{n=1}^{N}(\tau\tau_{n})^{-d}\psi\left\\{(\tau\tau_{n})^{-1}(\breve{\vartheta}_{t}-\breve{\vartheta}_{t-1})\right\\}.$

We define the stochastically perturbed model with a Markov process
$\\{(\breve{X}_{n},\breve{\Theta}_{n}),\ 0\leq n\leq N\\}$, observation
process $\breve{Y}_{1:N}$ and parameter $(\theta,\ \tau,\ \tau_{0:N})$ by the
factorization of their joint probability density

$p_{\breve{X}_{0:N},\breve{Y}_{1:N},\breve{\Theta}_{0:N}}(x_{0:N},y_{1:N},\breve{\vartheta}_{0:N};\theta,\
\tau,\ \tau_{0:N})\\\
=p_{\breve{\Theta}_{0:N}}(\breve{\vartheta}_{0:N};\theta,\ \tau,\
\tau_{0:N})p_{\breve{X}_{0:N},\breve{Y}_{1:N}|\breve{\Theta}_{0:N}}(x_{0:N},\
y_{1:N}|\breve{\vartheta}_{0:N}),$

where

$p_{\breve{X}_{0:N},\breve{Y}_{1:N}|\breve{\Theta}_{0:N}}(x_{0:N},y_{1:N}|\breve{\vartheta}_{0:N};\theta,\tau,\
\tau_{0:N})=\\\
\mu(x_{0};\breve{\vartheta}_{0})\prod_{n=1}^{N}f_{n}(x_{n}|x_{n-1};\breve{\vartheta}_{n})\prod_{n=1}^{N}g_{n}(y_{n}|x_{n};\breve{\vartheta}_{n}).$

This extended model can be used to define a perturbed parameter log-likelihood
function, defined as

$\breve{\ell}(\breve{\vartheta}_{0:N})=\log
p_{\breve{Y}_{1:N}|\breve{\Theta}_{0:N}}(y_{1:N}^{*}|\breve{\vartheta}_{0:N};\theta,\tau,\tau_{0:N}).$
(2.2)

Here, the right hand side does not depend on $\theta$, $\tau$ or $\tau_{0:N}$.
We have designed (2.2) so that, setting
$\breve{\vartheta}^{[N+1]}=(\theta,\theta,\dots,\theta)\in\mathbb{R}^{d(N+1)},$
the log-likelihood of the unperturbed model can be written as
$\ell(\theta)=\breve{\ell}(\breve{\vartheta}^{[N+1]}).$ For the perturbed
likelihood, we need an additional assumption of the extended version.

######  Assumption 4.

$\breve{\ell}$ is four times continuously differentiable. For all
$\theta\in\mathbb{R}^{d}$, there exist $\epsilon>0$, $D>0$ and $\delta$
defined as in Assumption 2, such that for all $0<\eta<\delta$ and
$u_{0:N}\in\mathbb{R}^{d(N+1)}$,$\breve{\mathcal{L}}(\breve{\vartheta}^{[N+1]}+u_{0:N})\leq
De^{\epsilon\sum_{n=1}^{N}|u_{n}|^{\eta}},$ where
$\breve{\mathcal{L}}(\breve{\vartheta}_{0:N})=\exp\\{\breve{\ell}(\breve{\vartheta}_{0:N})\\}$
is the perturbed likelihood.

Let $\breve{\mathbb{E}}_{\theta,\tau,\tau_{0:N}}$,
${\mathrm{\breve{C}ov}}_{\theta,\tau,\tau_{0:N}}$,
${\mathrm{\breve{V}ar}}_{\theta,\tau,\tau_{0:N}}$ denote the expectation,
covariance and variance with respect to the associated posterior,
$p_{\breve{\Theta}_{0:N}|\breve{Y}_{1:N}}(\breve{\vartheta}_{0:N}|y_{1:N}^{*};\theta,\
\tau,\tau_{0:N}).$ By using $\breve{\mathbb{E}}$, ${\mathrm{\breve{C}ov}}$,
${\mathrm{\breve{V}ar}}$ instead of
$\breve{\mathbb{E}}_{\theta,\tau,\tau_{0:N}}$,
${\mathrm{\breve{C}ov}}_{\theta,\tau,\tau_{0:N}}$,
${\mathrm{\breve{V}ar}}_{\theta,\tau,\tau_{0:N}}$ respectively, a theorem
similar to theorem 4 of [8] but for random walk noise instead of independent
white noise is derived.

######  Theorem 1.

[Theorem 2 of [30]] Suppose Assumptions 1, 2 and 4, there exists a constant
$C$ independent of $\tau,\tau_{1},...\tau_{N}$ such that,

$\left|\nabla\ell\left(\theta\right)-\tau^{-2}\Psi^{-1}\left\\{\tau_{0}^{-2}\breve{\mathbb{E}}\left(\breve{\Theta}_{0}-\theta|\breve{Y}_{1:N}=y_{1:N}^{*}\right)\right\\}\right|<C\tau^{2},$

where $\Psi$ is the non-singular covariance matrix associated to $\psi$.

Theorem 1 formally allows an approximation of
$\nabla{\ell}\left(\theta\right)$. [30] also presents an alternative
variations on these results which lead to more stable Monte Carlo estimation.

######  Theorem 2.

[Theorem 3 of [30]] Suppose Assumption 1, 2 and 4 hold. In addition, assume
that $\tau_{n}=O(\tau^{2})$ for all $n=1\ldots N$, the following holds true,

$\left|\nabla\ell\left(\theta\right)-\frac{1}{N+1}\tau^{-2}\tau_{0}^{-2}\Psi^{-1}\sum_{n=0}^{N}\left\\{\breve{\mathbb{E}}\left(\breve{\Theta}_{n}-\theta|\breve{Y}_{1:N}=y_{1:N}^{*}\right)\right\\}\right|=O(\tau^{2}).$
(2.3)

These theorems are useful for our approaches because we can approximate the
gradient of the log-likelihood of the extended model to the second order of
$\tau$ which we will later show that it fits well with our accelerate
simulation based setup.

##  3 Proposed accelerate iterated filtering

Our motivation comes from the accelerated gradient method for smooth non-
linear stochastic programming literature. By using an approximation of the
score function, it is possible to use an accelerated gradient method as in
Nesterov acceleration scheme in optimization literature. One issue with the
accelerated gradient approach is that it is not clear how the technique can be
used in situations where both the likelihood and the gradient are intractable.
These sorts of examples are common in scientific applications of state space
models where the state process is a diffusion process or an ordinary
differential equation (ODE) with stochastic coefficients. However, in these
family of iterated filtering approaches, the score function can be
approximated with noise under control without affecting the convergence rate.
Specifically, applying an accelerated inexact gradient algorithm in the
iterated filtering approach can obtain an optimal rate of convergence.

In this paper, $\epsilon_{k}$ denotes the error in the approximation of the
gradient. Using the same notation as [11], denote the sequences of magnitudes
of the errors in the gradient approximations
$\left\\{\left\|\epsilon_{k}\right\|\right\\}$. Suppose the following
assumptions:

######  Assumption 5.

The function $\ell$ : $\Theta\rightarrow\mathbb{R}$ is differentiable, bounded
from above and has a L-Lipschitz-continuous gradient, i.e. for all $\theta,\
\vartheta\in\Theta$,
$\left\|\nabla\ell(\theta)-\nabla\ell(\vartheta)\right\|\leq
L\left\|\theta-\vartheta\right\|,$ where $\nabla\ell$ denotes the gradient of
$\ell$. The function $\ell$ attains its maximum at a certain
$\theta^{*}\in\Theta$.

In the sequel, $\Theta$ denotes a finite-dimensional Euclidean space with norm
$\left\|\cdot\right\|$ and inner product
$\left\langle\cdot,\cdot\right\rangle$. It can be shown that (e.g. in [27])
Assumption 5 is equivalent to

$\left|\ell(\vartheta)-\ell(\theta)-\left\langle\nabla\ell(\theta),\mathrm{\vartheta}-\theta\right\rangle\right|\leq\frac{L}{2}\left\|\vartheta-\theta\right\|^{2},\;{\displaystyle\forall\theta,\
\vartheta\in\Theta}$ (3.1)

It is well-known that the gradient descent method converges for a general non-
convex optimization problem but it does not achieve the optimal rate of
convergence, in terms of the functional optimality gap, when $\ell(\cdot)$ is
convex [11]. In contrast, the accelerated gradient method in [28] is optimal
for solving convex optimization problems, but does not necessarily converge
for solving nonconvex optimization problems. [11] proposed a modified
accelerated gradient method which can converge in both convex and non-convex
optimization problem. However, they assumed unbiased estimation of the
gradient which is not satisfied for most simulation-based inferences. Below,
we extend the approach of Ghadimi to an accelerated inexact gradient (AIG)
method in the context of accelerate iterated filtering. That is, we allow bias
in gradient approximation by properly specifying the stepsize policy. We prove
that it not only achieves the same optimal rate of convergence for both convex
and non-convex optimizations, but also exhibits the best-known rate of
convergence for simulation-based inference problems.

Algorithm 1 Accelerate Inexact Gradient (AIG)

1:

2:

3:$\theta_{0}\in\Theta.$

4:$\left\\{\beta_{\mathrm{k}}>0\right\\}$, $\left\\{\lambda_{k}>0\right\\}$
for any $k\geq 2$.

5:$\left\\{\alpha_{k}\right\\}\in\left(0,1\right)$ for $k>1$ and
$\alpha_{1}=1$.

6:$\theta_{0}^{ag}=\theta_{0}$. $\triangleright$ Initialize

7:for $k$ in $1...N$ do

8: $\theta_{k}^{md}=(1-\alpha_{k})\theta_{k-1}^{ag}+\alpha_{k}\theta_{k-1}$
(3.2)

9:
$\theta_{k}=\theta_{k-1}-\lambda_{k}\left(\widehat{\nabla\ell(\theta_{k}^{md})}\right)$
(3.3)

10:
$\theta_{k}^{ag}=\theta_{k-1}^{md}-\beta_{k}\left(\widehat{\nabla\ell(\theta_{k}^{md})}\right)$
(3.4) $\triangleright$ where $\widehat{\nabla\ell(\theta_{k}^{md})}$ is an
estimation of $\nabla\ell(\theta_{k}^{md})$ with error $\epsilon_{k}$.

11:end for

In addition to Assumption 5, we assume a noise control condition for Algorithm
1.

######  Assumption 6.

$\Theta$ is bounded. There exists an $A<\infty$ such that
$\sum_{k=1}^{N}\lambda_{k}\left\|\epsilon_{k}\right\|<A.$

Given some mild conditions often satisfied by controlling the artificial
noises, we have the following result.

######  Theorem 3.

(Extension of Theorem 1 of [11]).
Suppose Assumptions 5 and 6 hold. In addition, let $\\{\theta_{k},\
\theta_{k}^{ag}\\}$ $k\geq 1$ be computed by Algorithm 1.
a) If sequences $\left\\{\alpha_{k}\right\\},\left\\{\beta_{k}\right\\}$,
$\left\\{\lambda_{k}\right\\}$ and $\left\\{\Gamma_{k}\right\\}$ satisfy

$\Gamma_{k}:=\begin{cases}1&k=1\\\ (1-\alpha_{k})\Gamma_{k-1}&k\geq
2\end{cases},\ $ (3.5)
$C_{k}:=1-L{\displaystyle\lambda_{k}-\frac{L(\lambda_{k}-\beta_{k})^{2}}{2\lambda_{k}\alpha_{k}\Gamma_{k}}\left(\sum_{\tau=k}^{N}\frac{1}{\Gamma_{\tau}}\right)>0,\mbox{
for}\ 1\leq k\leq N},$ (3.6)

then for any $N\geq 1$, we have for some $B<\infty$,

${\displaystyle\min_{k=1,...,N}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}\leq\frac{\ell(\theta_{0})-\ell^{*}+B}{\sum_{k=1}^{N}\lambda_{k}C_{k}}}.$
(3.7)

b) Suppose that $\ell(\cdot)$ is convex. If sequences
$\left\\{\alpha_{k}\right\\},\left\\{\beta_{k}\right\\}$,$\left\\{\lambda_{k}\right\\}$
and $\left\\{\Gamma_{k}\right\\}$ satisfy

$\alpha_{k}\lambda_{k}\leq\beta_{k}<\frac{1}{L},$ (3.8)
$\frac{\alpha_{1}}{\lambda_{1}\Gamma_{1}}\geq\frac{\alpha_{2}}{\lambda_{2}\Gamma_{2}}\geq\ldots,$
(3.9)

then for any $N\geq 1$, we have

$\min_{k=1,...,N}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}\\\
\leq
2\frac{\frac{\|\theta^{*}-\theta_{0}\|^{2}}{2\lambda_{1}}+\sum_{k=1}^{N}\Gamma_{k}^{-1}\left[\beta_{k}\left\|\epsilon_{k}\right\|\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|+\alpha_{k}\left\|\epsilon_{k}\right\|\left\|\theta_{k-1}-\theta_{0}\right\|\right]}{\sum_{k=1}^{N}\Gamma_{k}^{-1}\beta_{k}(1-L\beta_{k})},$
(3.10) $\ell(\theta_{N}^{ag})-\ell(\theta^{*})\\\
\leq\Gamma_{N}\left[\frac{\left\|\theta_{0}-\theta^{*}\right\|^{2}}{\lambda_{1}}+\sum_{k=1}^{N}\Gamma_{k}^{-1}\left[\beta_{k}\left\|\epsilon_{k}\right\|\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|+\alpha_{k}\left\|\epsilon_{k}\right\|\left\|\theta_{k-1}-\theta_{0}\right\|\right]\right].$
(3.11)

There are various options for selecting
$\left\\{\alpha_{k}\right\\},\left\\{\beta_{k}\right\\}$,$\left\\{\lambda_{k}\right\\}$,$\left\\{\Gamma_{k}\right\\}$.
By controlling error $\epsilon_{k}$, we can provide some of these selections
below which guarantee the optimal convergence rate of the AIG algorithm for
both convex and nonconvex problems.

######  Theorem 4.

Suppose Assumptions 5 and 6 hold. In addition, suppose that
$\left\\{\beta_{k}\right\\}$ in the accelerated gradient method are set to
$\beta_{k}=\frac{1}{2L}$.

a) If sequences $\left\\{\alpha_{k}\right\\}$ and
$\left\\{\lambda_{k}\right\\}$ satisfy

$\lambda_{k}\in\left[\beta_{k},(1+\frac{1}{k})\beta_{k}\right],\mbox{ for
}{\displaystyle\;\forall k\geq 1},$ (3.12)

then for any $N\geq 1$, we have

${\displaystyle\min_{k=1,\ldots,N}\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\|^{2}\leq
O\left(\frac{1}{N}\right)}.$ (3.13)

Suppose that $\epsilon_{k}=O\left(\tau^{2}\right)\leq O(\frac{1}{k})$, then
the AIG method can find a solution $\bar{\theta}$ such that
$\left\|\nabla\ell(\bar{\theta})\right\|^{2}\leq\epsilon$ in at most
$O(1/\epsilon^{2})$ iterations.

b) Suppose that $\ell(\cdot)$ is convex and
$\epsilon_{k}=O\left(\tau^{2}\right)\leq O(\frac{1}{k^{2+\delta+\delta_{1}}})$
for some $\delta_{1}>0$. If $\left\\{\lambda_{k}\right\\}$ satisfies

${\displaystyle\lambda_{k}=\left(k^{1+\delta}-\left(k-1\right)^{1+\delta}\right)\>\forall
k\geq 1},$ (3.14)

then for any $N\geq 1$, we have

$\min_{k=1,...,N}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}{\displaystyle\leq}O\left(\frac{1}{N^{2+\delta}}\right),$
(3.15) ${\displaystyle\ell(\theta_{N}^{ag})-\ell(\theta^{*})\leq
O\left(\frac{1}{N^{1+\delta}}\right)},$ (3.16)

then the AIG method can find a solution $\bar{\theta}$ such that
$\left\|\nabla\ell(\bar{\theta})\right\|^{2}\leq\epsilon$ in
$O\left(1/\epsilon^{\frac{1}{2+\delta}}\right)$ at most.

Algorithm 2 Accelerate Iterated Filtering (AIF)

1:

2:

3:Starting parameter, $\theta_{0}=\theta^{ag}_{0}$ , sequences,
$\alpha_{n},\beta_{n},\lambda_{n},\Gamma_{n}$

4:simulator for $f_{X_{0}}(x_{0}|\theta)$,
$f_{X_{n}|X_{n-1}}(x_{n}|x_{n-1}|\theta)$, evaluator for
$f_{Y_{n}|X_{n}}(y_{n}|x_{n}|\theta)$

5:data, $y^{*}_{1:N}$, labels designating IVPs, $I\subset\\{1,\dots,p\\}$,
initial scale multiplier, $C>0$

6:number of particles, $J$, number of iterations, $M$, cooling rate, $0<a<1$,
perturbation scales, $\sigma_{1:p}$

7:

8:Maximum likelihood estimate $\theta_{MLE}$

9:$\theta^{md}_{0}=\theta_{0}$ $\triangleright$ Initialize

10:$[\Theta^{F}_{0,j}]_{i}\sim
N\left([\theta^{md}_{0}]_{i},(Ca^{m-1}\sigma_{i})^{2}\right)$ for $i$ in
$1..p$, $j$ in $1...J$. $\triangleright$ Initialize filter mean for parameters

11:simulate $X_{0,j}^{F}\sim f_{X_{0}}\big{(}\cdot;{\Theta^{F}_{0,j}}\big{)}$
for $j$ in $1..J$. $\triangleright$ Initialize states

12:for $m$ in $1...M$ do

13:
$\theta^{md}_{m}=(1-\alpha_{m})\theta^{ag}_{m-1}+\alpha_{m_{1}}\theta_{m-1}$.

14: for $n$ in $1...N$ do

15:
$\big{[}\Theta_{n,j}^{P}\big{]}_{i}\sim\mathcal{N}\big{(}\big{[}\Theta^{F}_{n-1,j}\big{]}_{i},(c^{m-1}\sigma_{i})^{2}\big{)}$
for $i\notin I$, $j$ in $1:J$. $\triangleright$ Perturb

16:
${X}_{n,j}^{P}\sim{f}_{n}\big{(}{x}_{n}|{X}_{n-1,j}^{F};{\Theta_{n,j}^{P}}\big{)}$
for $j$ in $1:J$. $\triangleright$ Simulate prediction particles

17: $w(n,j)=g_{n}(y_{n}^{*}|X_{n,j}^{P};\Theta_{n,j}^{P})$ for $j$ in $1:J$.
$\triangleright$ Evaluate weights

18: $\breve{w}(n,j)=w(n,j)/\sum_{u=1}^{J}w(n,u)$. $\triangleright$ Normalize
weights

19: $k_{1:J}$ with $P\left\\{k_{u}=j\right\\}=\breve{w}\left(n,j\right)$.
$\triangleright$ Apply systematic resampling to select indices

20: $X_{n,j}^{F}=X_{n,k_{j}}^{P}$ and $\Theta_{n,j}^{F}=\Theta_{n,k_{j}}^{P}$
for $j$ in $1:J$. $\triangleright$ Resample particles

21: end for

22:
$S_{m}=c^{-2(m-1)}\Psi^{-1}\sum_{n=1}^{N}\big{[}\left(\bar{\theta}_{n}-\theta^{md}_{m-1}\right)\big{]}$
$\triangleright$ Update Parameters

23:
$\big{[}\theta_{m}\big{]}_{i}=\theta_{m-1}-\lambda_{m-1}\big{[}S_{m}\big{]}_{i}$
for $i\notin I$.

24:
$\big{[}\theta^{ag}_{m}\big{]}_{i}=\theta^{md}_{m-1}-\beta_{m-1}\big{[}S_{m}\big{]}_{i}$
for $i\notin I$.

25:
$\big{[}\theta_{m}\big{]}_{i}=\frac{1}{J}\sum_{j=1}^{J}\big{[}\Theta^{F}_{L,j}\big{]}_{i}$
for $i\in I$.

26:end for

We now add a few remarks about the extension results obtained in Theorem 4.
First, if the problem is convex, by choosing more aggressive stepsizes
$\\{\lambda_{k}\\}$ in (3.14), the AIG method exhibits the optimal rate of
convergence in (3.16). It is also worth noting that with such a selection of
$\\{\lambda_{k}\\}$, the AIG method can find a solution $\bar{\theta}$ such
that $\left\|\nabla\ell(\bar{\theta})\right\|^{2}\leq\epsilon$ in at most
$O(1/\epsilon^{1/2+\delta})$ iterations. The latter result has been shown by
[27], [11] but only for the accelerate unbiased gradient method. Second,
observe that $\left\\{\lambda_{k}\right\\}$ in (3.12) for general nonconvex
problems is in the order of $O(1/L)$, while the one in (3.14) for convex
problems are more aggressive (in $O(k/L))$. The value $\delta$ is optimal at
$1$ for convergence rate. However, it may not be optimal for computation of
controlling the noises. Finally, we show that we can apply the stepsize policy
in (3.12) for solving general inexact gradient problems for both convex and
nonconvex optimization. The sequential Monte Carlo filter can be arbitrarily
approximated to the exact filter by choosing sufficiently large number of
particles [13]. It can be seen that we can choose the perturbation sequence so
that the gradient noise satisfies condition in Theorem 4. For completeness, we
present the pseudo code of the proposed algorithm as in Algorithm 2.

##  4 Numerical examples

To measure the performance of the new inference algorithm, we evaluate our
accelerate iterated filtering on some benchmark examples and compare it to the
existing simulation-based approaches. We make use of well tested and
maintained code of R [33] packages such as pomp [17]. Specifically, models are
coded using C snipet declarations [17]. New algorithm is written in R package
is2, which provides user friendly interfaces in R and efficient matrix
operations in the highly optimized Rcpp [9]. All the simulation-based
approaches mentioned above use sequential Monte Carlo algorithm (SMC),
implemented using bootstrap filter. Experiments were carried out on a cluster
of $32$ cores Intel Xeon E5-2680 $2.7$ Ghz with $256$ GB memory. For a fair
comparison, we try to use the same setup and assessment for every inference
method. A public Github repository containing scripts for reproducing our
results may be found at https://github.com/nxdao2000/AIFcomparisons.

###  4.1 Toy example: A linear, Gaussian model

In this subsection, we compare our accelerate iterated filtering algorithm to
the original iterated filtering algorithm IF1 [14], Bayes map iterated
filtering (IF2) [15] and the second-order iterated smoothing (IS2) [30]. It
has been shown in [30] and [15] that the second-order iterated smoothing with
white noise (IS1) [8] and particle Markov chain Monte Carlo (PMCMC) [1] do not
perform as well as Bayes map iterated filtering so we leave them out. For a
computationally convenient setting, simple models provide an opportunity to
test the basic features of inference algorithms. Therefore, we first consider
a bivariate discrete time Gaussian autoregressive process, a relatively simple
mechanistic model. This model is chosen so that the Monte Carlo calculations
can be verified using a Kalman filter. For this example, there are some
alternatives to iterated filtering class. For example, EM and MCMC algorithms
would be practical in this case although they do not scale well to large
dynamic models, so we do not include them here. The model is given by the
state space forms: $X_{n}|X_{n-1}=x_{n-1}\sim\mathcal{N}(\alpha
x_{n-1},\sigma^{\top}\sigma),$ $Y_{n}|X_{n}=x_{n}\sim\mathcal{N}(x_{n},I_{2})$
where $\alpha$, $\sigma$ are $2\times 2$ matrices and $I_{2}$ is $2\times 2$
identity matrix. The data are simulated from the following parameters:

${\displaystyle\alpha=\left[\begin{array}[]{cc}\alpha_{1}&\alpha_{2}\\\
\alpha_{3}&\alpha_{4}\end{array}\right]=\left[\begin{array}[]{cc}0.8&-0.5\\\
0.3&0.9\end{array}\right],\ \sigma=\left[\begin{array}[]{cc}3&0\\\
-0.5&2\end{array}\right].}$

The number of time points $N$ is set to $100$ and initial starting point
$X_{0}=\left(-3,4\right)$. For each method mentioned above, we estimate
parameters $\alpha_{2}$ and $\alpha_{3}$ for this model using $J=1000$
particles and run our estimation for $M=25$ iterations. We start the initial
search uniformly on a large rectangular region $[-1,1]\times[-1,1]$. As can be
seen from Fig. 1, all of the distributions of estimated maximized log
likelihoods touch the true MLE (computed from Kalman filter) at the vertical
broken line, implying that they all successfully converged. The results show
that AIF is the most efficient method of all because using AIF the results
have higher mean and smaller variance compared to other approaches, indicating
a higher empirical convergence rate. Algorithmically, AIF has similar
computational costs with the first order approaches IF1, IF2, and is cheaper
than the second order approach IS2. In deed, average computational time of
twenty independent runs of each approach is given in Table 1. Additional
overheads for estimating score make the computation time of AIF a bit larger
compared to computational time of IF2. However, with complex models and large
enough number of particles, these overheads become negligible and
computational time of AIF will be similar to other first order approaches. The
fact that it has the convergence rate of second order with computation
complexity of first-order shows that it is a very promising algorithm. In
addition, the results also imply that AIF is robust to initial starting
guesses.

Table 1: Computation times, in seconds, for the toy example. |  $J=100$ |  $J=1000$ |  $J=10000$
---|---|---|---
IF1 | 1.656 | 5.251 | 62.632
IF2 | 1.591 | 5.156 | 61.072
IS2 | 2.530 | 10.198 | 135.248
AIF | 2.729 | 10.278 | 132.016
Figure 1: Comparison of estimators for the linear, Gaussian toy example,
showing the densities of the MLEs estimated by the IF1, IF2, AIF and IS2
methods. The parameters $\alpha_{2}$ and $\alpha_{3}$ were estimated, started
from 200 randomly uniform initial values over a large rectangular region
$[-1,1]\times[-1,1]$.  Figure 2: Comparison of different estimators. The
likelihood surface for the linear model, with the location of the MLE is
marked with a green cross. The crosses show final points from 40 Monte Carlo
replications of the estimators: (A) Original iterated filtering method; (B)
Bayes map iterated filtering method; (C) Accelerate iterated filtering method;
(D) Second-order iterated filtering method; Each method, was started uniformly
over the rectangle shown, with $M=25$ iterations, $N=1000$ particles, and a
random walk standard deviation decreasing from $0.02$ geometrically to $0.011$
for both $\alpha_{2}$ and $\alpha_{3}$.

To see how the final MLEs clustered around the true MLE, we only show $40$
Monte Carlo replications for this toy example. As can be observed from Fig. 2,
most of the replications clustered near the true MLE for AIF approach, while
none of them stays in a lower likelihood region. It can be interpreted as a
statistical summary of Fig. 2, with $200$ Monte Carlo replications. These
results indicate that AIF is clearly the best of the investigated methods for
this test compared to others. Given additional computational resources, we
also checked how the results of each method compared. Specifically, we set
$M=100$ iterations and $J=10000$ particles, with the random walk standard
deviation decreasing geometrically from $0.02$ down to $0.0018$ for each
method. In this situation, we confirm that AIF is the best among other IF1,
IF2 and IS2. All methods have comparable computational demands for given $M$
and $J$.

Figure 3: Comparison of estimators for the linear, Gaussian toy example,
showing the densities of the MLEs estimated by the IF1, IF2, AIF and IS2
methods using $M=100$ iterations and $J=10000$ particles. The parameters
$\alpha_{2}$ and $\alpha_{3}$ were estimated, started from 200 randomly
uniform initial values over a large rectangular region $[-1,1]\times[-1,1]$.

###  4.2 Malaria benchmark

Many real world dynamic systems are highly nonlinear, partially observed and
even weakly identifiable. To demonstrate the capabilities of accelerate
iterated filtering for such situations, we apply it to evaluate the likelihood
in a stochastic differential equation for vivax malaria model of Roy et al.,
[34]. The reason to choose this challenging model is that it provides a
rigorous performance benchmark for our verification. The model $SEIH^{3}QS$ we
consider splits up the study population of size $P(t)$ into seven classes:
susceptible individuals, $S(t)$, exposure $E(t)$, infected individuals,
$I(t)$, dormant classes $H_{1}(t)$, $H_{2}(t)$, $H_{3}(t)$ and recovered
individuals, $Q(t)$. This strain of malaria characterized by relapse following
initial recovery from symptoms [30]. Therefore the the last $S$ in the model
name indicates the possibility that a recovered person can return to the class
of susceptible individuals. The data, denoted by $y_{1:N}^{*}$, are in the
form of monthly time series over a 20-year period, counting the malaria
morbidity. $\delta$ denotes the mortality rate, $\kappa(t)$ a delay stage,
$\mu_{SE}(t)$ the current force of infection, and $\tau_{D}$ the mean latency
time. The state process is

$X(t)=\big{(}S(t),E(t),I(t),Q(t),H_{1}(t),H_{2}(t),H_{3}(t),\kappa(t),\mu_{SE}(t)\big{)},$

where transition rates from stage $H_{1}$ to $H_{2}$, $H_{2}$ to $H_{3}$ and
$H_{3}$ to $Q$ are specified to be $3\mu_{HI}$ while infected population to
dormancy transition rate is $\mu_{IH}$. The model satisfies the following
stochastic differential equation system

$\displaystyle dS/dt$ $\displaystyle=$ $\displaystyle\delta
P+\mathrm{d}P/dt+\mu_{IS}I+\mu_{QS}Q$ $\displaystyle\hskip
14.22636pt+a\mu_{IH}I+b\mu_{EI}E-\mu_{SE}(t)S-\delta S,$ $\displaystyle dE/dt$
$\displaystyle=$ $\displaystyle\mu_{SE}(t)S-\mu_{EI}E-\delta E,$
$\displaystyle dI/dt$ $\displaystyle=$
$\displaystyle(1-b)\mu_{EI}E+3\mu_{HI}H_{n}-(\mu_{IH}+\mu_{IS}+\mu_{IQ})I-\delta
I,$ $\displaystyle dH_{1}/dt$ $\displaystyle=$
$\displaystyle(1-a)\mu_{IH}I-n\mu_{HI}H_{1}-\delta H_{1},$ $\displaystyle
dH_{i}/dt$ $\displaystyle=$ $\displaystyle
3\mu_{HI}H_{i-1}-3\mu_{HI}H_{i}-\delta H_{i}\hskip 14.22636pt\mbox{for
$i\in\\{2,3\\}$},$ $\displaystyle dQ/dt$ $\displaystyle=$
$\displaystyle\mu_{IQ}I-\mu_{QS}Q-\delta Q.$

In addition, the malaria pathogen reproduction within the mosquito vector is
given by

$\displaystyle\mathrm{d}\kappa/\mathrm{d}t$ $\displaystyle=$
$\displaystyle[\lambda(t)-\kappa(t)]/\tau_{D},$
$\displaystyle\mathrm{d}\mu_{SE}/\mathrm{d}t$ $\displaystyle=$
$\displaystyle[\kappa(t)-\mu_{SE}(t)]/\tau_{D},$

where $\lambda(t)$ is the latent force of infection and $\lambda(t)$,
$\kappa(t)$ and $\mu_{SE}(t)$ satisfies

$\mu_{SE}(t)=\int_{-\infty}^{t}\gamma(t-s)\lambda(s)\mathrm{d}s,$ (4.1)

with $\gamma(s)=\frac{(2/\tau_{D})^{2}s^{2-1}}{(2-1)!}\exp(-2s/\tau_{D})$, a
gamma distribution with shape parameter $2$. Since the latent force of
infection is constrained by rainfall covariate $R(t)$ and some Gamma white
noise, from Roy et al., [34] we have:

${\lambda(t)=\left(\frac{I+qQ}{P}\right)\times\exp\left\\{\sum_{i=1}^{N_{s}}b_{i}s_{i}(t)+b_{r}R(t)\right\\}{\times\left[\frac{\mathrm{d}\Gamma(t)}{\mathrm{d}t}\right]}}.$

In this equation, $q$ denotes a reduced infection risk from humans in the
$\mathrm{Q}$ class and $\\{s_{i}(t),i=1,\dots,N_{s}\\}$ is a periodic cubic
B-spline basis, with $N_{s}=6$. The observation model for $Y_{n}$ is a
negative binomial distribution with mean $M_{n}$ and variance
$M_{n}+M_{n}^{2}\sigma_{\mathrm{obs}}^{2}$ where
$M_{n}=\rho\int_{t_{n-1}}^{t_{n}}[\mu_{EI}E(s)+3\mu_{HI}H_{3}(s)]ds$ is the
number of new cases observed from time $t_{n-1}$to time $t_{n}$ and $\rho$ it
the mean age. The coupled system of stochastic differential equations is
solved using an Euler-Maruyama scheme [18] with a time step of $1/20$ month in
our case.

Figure 4: The density of the maximized log likelihood approximations estimated
by IF1, IF2, IS2, RIS1 and AIF for the malaria model when using J = 1000 and M
= 50. The log likelihood at a previously computed MLE is shown as a dashed
vertical line.

Given the data obtained from National Institutes of Malaria Research [34], we
carried out simulation-based inference via the original iterated filtering
(IF1), the perturbed Bayes map iterated filtering (IF2), the second order
iterated smoothing (IS2), and the new accelerate iterated filtering (AIF). The
inference goal used to assess all of these methods is to find high likelihood
parameter values starting from randomly drawn values in a large
hyperrectangle. In the presence of possible multi-modality, weak
identifiability, and considerable Monte Carlo error of this model, we start
$200$ random searches. The random walk standard deviation is initially set to
$0.1$ for estimated parameters while the cooling rate $c$ is set to
$0.1^{0.02}\approx 0.95$. These corresponding quantities for initial value
parameters are $2$ and $0.1^{0.02}$, respectively, but they are applied only
at time zero. We run our experiment on a cluster computers with $M=50$
iterations and with $J=1000$ particles. The reason to choose these values for
this model is that increasing the iterations to $100$ and the number of
particles to $10000$ does not improve the results much but it takes
significant longer time. Figure 4 shows the distribution of the MLEs estimated
by IF1, IF2, IS2 and AIF. All distributions touch the global maximum as
expected and the higher mean and smaller variance of IF2, AIF estimation
clearly demonstrate that they are considerably more effective than IF1. Note
that the computational times for IF1, IF2, IS2 and AIF are 44.86, 43.92, 53.10
and 52.25 minutes respectively, confirming that accelerate iterated filtering
has essentially the same computational cost as first order methods IF1, IF2
and is cheaper a bit than IS2, for a given Monte Carlo sample size and number
of iterations. In this hard problem, while IF1 reveal their limitations, we
have shown that IF2 and AIF can still offer a substantial improvement. A
natural heuristic idea to further improve the method is hybridizing IF2 and
AIF but we leave it for the future work.

##  5 Conclusion

In this paper, we have proposed a novel class of iterated filtering theory
using an accelerated inexact gradient approach. We have shown that choosing
perturbation sequence and number of particles carefully results in an
algorithm which has led to many advances including the statistical and
computational efficiency. This is also very fruitful as it is extendable to a
more generalized class of algorithm, based on proximal theory. Previous proof
of iterated filtering class require some difficult conditions, which is not
easily verifiable. However, in this article, we use only general standard
gradient conditions. We are going further down the road of a more systematic
approach which could be easily generalized to the state of the art algorithm
in the optimization literatures. The convergence rate is also explicitly
stated and it is better than standard theory. From a theoretical point of
view, it could be an interesting perspective and insight.

In addition, from practical point of view, we have provided an efficient
framework, applicable to a general class of nonlinear, non-Gaussian non-
standard POMP models, especially suitable in the control feedback system.
There are a lot of such systems, which are not well-treated by current
available modeling framework. We simultaneously present the performance of our
open source software package is2 to facilitate the needs of the community. The
performance of this new approaches surpass the other frameworks by a large
margin of magnitude.

It may be surprising that this simple accelerated inexact gradient approach
has the needed convergence properties, and can easily be generalized, at least
in some asymptotic sense. It is not hard to show that the accelerated inexact
proximal gradient iterated filtering theory can be adapted to apply with
iterated smoothing and with either independent white noise or random walk
perturbations while our empirical results still show strong evidences of the
improvements. In principle, different simulation-based inference methods can
readily be hybridized to build on the strongest features of multiple
algorithms. Our results could also be applied to develop other simulation-
based methodologies which can take advantage of proximal map. For example, it
may be possible to use our approach to help design efficient proposal
distributions for particle Markov chain Monte Carlo algorithms. The
theoretical and algorithmic innovations of this paper will help to build a new
direction for future developments on this frontier. Applying this approach to
methodologies like Approximate Bayesian Computation (ABC), Liu-West Particle
Filter (LW-PF), Particle Markov chain Monte Carlo (PMCMC), with different
samplers scheme, e.g. forward backward particle filter, forward smoothing or
forward backward smoothing are foreseeable extensions.

##  Appendix A Proofs

We first need a simple technical result (see Lemma 1 of [11]). The proof is
the same as that of Lemma 1 of [11] but we provide it here for completeness.

######  Lemma 1.

(Lemma1 of [11].

Assume sequences $\\{\alpha_{k}\\}\in\left(0,1\right)$ for $k>1$ and
$\alpha_{1}=1$ and sequences $\\{a_{k}\\},$ $\\{\eta_{k}\\}$ satisfy

$a_{k}\leq(1\ -\alpha_{k})a_{k-1}+\eta_{k},\ k=1,2,\ldots$ (A.1)

If we define a positive sequence $\\{\Gamma_{k}\\}$ as in 3.5 then for any
$k\geq 1$, we have

${\displaystyle a_{k}\leq\Gamma_{k}\sum_{i=1}^{k}(\eta_{i}/\Gamma_{i})}.$

###### Proof.

Since $\alpha_{1}=1$ and $\Gamma_{1}=1$, from 3.5 we have

$a_{1}\leq\eta_{1}$

or

$\frac{a_{1}}{\Gamma_{1}}\leq\frac{\eta_{1}}{\Gamma_{1}}.$

Since $\Gamma_{k}>0$ for every $k>1$, dividing both sides of A.1 by
$\Gamma_{k}$,

$\frac{a_{k}}{\Gamma_{k}}\leq\frac{(1-\alpha_{k})a_{k-1}+\eta_{k}}{\Gamma_{k}}=\frac{a_{k-1}}{\Gamma_{k-1}}+\frac{\eta_{k}}{\Gamma_{k}},\
\forall k\geq 2.$

Summing up the above inequalities and rearranging the terms, the conclusion
follows. ∎

######  Lemma 2.

$\sum_{\tau=1}^{k}\frac{\alpha_{\tau}}{\Gamma_{\tau}}=\frac{1}{\Gamma_{k}}.$
(A.2)

###### Proof.

We have

$\sum_{\tau=1}^{k}\frac{\alpha_{\tau}}{\Gamma_{\tau}}=\frac{\alpha_{1}}{\Gamma_{1}}+\sum_{\tau=2}^{k}\frac{1}{\Gamma_{\tau}}\left(1-(1-\alpha_{\tau})\right)$
$={\displaystyle\frac{1}{\Gamma_{1}}+\sum_{\tau=2}^{k}(\frac{1}{\Gamma_{\tau}}-\frac{1}{\Gamma_{\tau-1}})=\frac{1}{\Gamma_{k}}}.$

∎

###  A.1 Proof of Theorem 3

###### Proof.

The proof follows closely to the proof of theorem 1 of [11] except we consider
bias estimate of the gradient. We first prove part a.

By 3.1 and 3.3, we have

$\ell(\theta_{k})\leq\ell(\theta_{k-1})+\left\langle\nabla\ell(\theta_{k-1}),\
\theta_{k}-\theta_{k-1}\right\rangle{\displaystyle+\frac{L}{2}\|\theta_{k}-\theta_{k-1}\|^{2}}$
$=\ell(\theta_{k-1})+\left\langle\left(\nabla\ell(\theta_{k-1})-\nabla\ell(\theta_{k}^{md})-\epsilon_{k}\right)+\left(\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right),\
-\lambda_{k}\left(\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right)\right\rangle$
$+\frac{L\lambda_{k}^{2}}{2}\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\|^{2}$
$=\ell(\theta_{k-1})-\lambda_{k}(1-\frac{L\lambda_{k}}{2})\|\nabla\ell\
(\theta_{k}^{md})+\epsilon_{k}\|^{2}-\lambda_{k}\left\langle\left(\nabla\ell(\theta_{k-1})-\nabla\ell(\theta_{k}^{md})-\epsilon_{k}\right),\
\left(\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right)\right\rangle$
$\leq\ell(\theta_{k-1})-\lambda_{k}(1-\frac{L\lambda_{k}}{2})\|\nabla\ell\
(\theta_{k}^{md})+\epsilon_{k}\|^{2}+\lambda_{k}\left(\|\nabla\ell(\theta_{k-1})-\nabla\ell(\theta_{k}^{md})\|+\left\|\epsilon_{k}\right\|\right)\cdot\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\|,$
$\leq\ell(\theta_{k-1})-\lambda_{k}(1-\frac{L\lambda_{k}}{2})\|\nabla\ell\
(\theta_{k}^{md})+\epsilon_{k}\|^{2}+\lambda_{k}\left(L\|\theta_{k-1}-\theta_{k}^{md}\|+\left\|\epsilon_{k}\right\|\right)\cdot\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\|,$
$=\ell(\theta_{k-1})-\lambda_{k}(1-\frac{L\lambda_{k}}{2})\|\nabla\ell\
(\theta_{k}^{md})+\epsilon_{k}\|^{2}+\lambda_{k}\left(L(1-\alpha_{k})\|\theta_{k-1}^{ag}-\theta_{k-1}\|+\left\|\epsilon_{k}\right\|\right)\cdot\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\|,$
$=\ell(\theta_{k-1})-\lambda_{k}\left(1-\frac{L\lambda_{k}}{2}\right)\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|{}^{2}$
$+L(1-\alpha_{k})\lambda_{k}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|\cdot\left\|\theta_{k-1}^{ag}-\theta_{k-1}\right\|+\lambda_{k}\left\|\epsilon_{k}\right\|\cdot\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|$
$\leq\ell\
(\theta_{k-1})-\lambda_{k}\left(1-\frac{L\lambda_{k}}{2}\right)\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|{}^{2}$
$+\frac{L\lambda_{k}^{2}}{2}\left\|\nabla\ell\
(\theta_{k}^{md})+\epsilon_{k}\right\|{}^{2}+\frac{L(1-\alpha_{k})^{2}}{2}\left\|\theta_{k-1}^{ag}-\theta_{k-1}\right\|{}^{2}+\lambda_{k}\left\|\epsilon_{k}\right\|\cdot\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|$
$=\ell(\theta_{k-1})-\lambda_{k}(1-L\lambda_{k})\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|{}^{2}$
$+{\displaystyle\frac{L(1-\alpha_{k})^{2}}{2}\left\|\theta_{k-1}^{ag}-\theta_{k-1}\right\|{}^{2}+\lambda_{k}\left\|\epsilon_{k}\right\|\cdot\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|}$
(A.3)

The second inequality is from triangular inequality and the Cauchy-Schwarz
inequality while the second inequality is due to the Lipschitz of gradient
assumption (1.2) and last equality comes from 3.2. We have the last inequality
follows from $ab\leq(a^{2}+b^{2})/2$. From 3.2, 3.3, and 3.4, it follows that

$\theta_{k}^{ag}-\theta_{k}=(1-\alpha_{k})\theta_{k-1}^{ag}+\alpha_{k}\theta_{k-1}-\beta_{k}\left(\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right)-\left(\theta_{k-1}-\lambda_{k}\left(\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right)\right)$
$=(1-\alpha_{k})(\theta_{k-1}^{ag}-\theta_{k-1})+(\lambda_{k}-\beta_{k})\left(\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right).$

Applying Lemma 1 where $\theta_{k}^{ag}-\theta_{k}:=a_{k}$ and
$\eta_{k}:=(\lambda_{k}-\beta_{k})\left(\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right)$,
we obtain

$\theta_{k}^{ag}-\theta_{k}=\Gamma_{k}\sum_{\tau=1}^{k}(\frac{\lambda_{\tau}-\beta_{\tau}}{\Gamma_{\tau}})\left(\nabla\ell(\theta_{\tau}^{md})+\epsilon_{\tau}\right).$

Since$\left\|\cdot\right\|{}^{2}$ is convex, using Jensen’s inequality and
Lemma 2 we have

$\left\|\theta_{k}^{ag}-\theta_{k}\right\|{}^{2}=\left\|\Gamma_{k}\sum_{\tau=1}^{k}(\frac{\lambda_{\tau}-\beta_{\tau}}{\Gamma_{\tau}})\left(\nabla\ell\
(\theta_{\tau}^{md})+\epsilon_{k}\right)\right\|^{2}$
$=\left\|\Gamma_{k}\sum_{\tau=1}^{k}\frac{\alpha_{\tau}}{\Gamma_{\tau}}\left[\left(\frac{\lambda_{\tau}-\beta_{\tau}}{\alpha_{\tau}}\right)\left(\nabla\ell(\theta_{\tau}^{md})+\epsilon_{k}\right)\right]\right\|^{2}$
$\leq\Gamma_{k}\sum_{\tau=1}^{k}\frac{\alpha_{\tau}}{\Gamma_{\tau}}\left\|\left(\frac{\lambda_{\tau}-\beta_{\tau}}{\alpha_{\tau}}\right)\left(\nabla\ell(\theta_{\tau}^{md})+\epsilon_{k}\right)\right\|^{2}$
$={\displaystyle\Gamma_{k}\sum_{\tau=1}^{k}\frac{(\lambda_{\tau}-\beta_{\tau})^{2}}{\Gamma_{\tau}\alpha_{\tau}}\left\|\nabla\ell(\theta_{\tau}^{md})+\epsilon_{\tau}\right\|^{2}}.$
(A.4)

Replacing the above bound in A.3, and the fact that
$\Gamma_{k}=\Gamma_{k-1}(1-\alpha_{k})$ as in 3.5 and that $\alpha_{k}\in(0,\
1]$ for all $k\geq 1$we obtain

$\ell(\theta_{k})\leq\ell(\theta_{k-1})-\lambda_{k}(1-L\lambda_{k})\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}$
$+\frac{L\Gamma_{k-1}(1-\alpha_{k})^{2}}{2}\sum_{\tau=1}^{k-1}\frac{(\lambda_{\tau}-\beta_{\tau})^{2}}{\Gamma_{\tau}\alpha_{\tau}}\left\|\nabla\ell(\theta_{\tau}^{md})+\epsilon_{\tau}\right\|^{2}+\lambda_{k}\left\|\epsilon_{k}\right\|\cdot\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|$
$\leq\ell(\theta_{k-1})-\lambda_{k}(1-L\lambda_{k})\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}$
$+{\displaystyle\frac{L\Gamma_{k}}{2}\sum_{\tau=1}^{k}\frac{(\lambda_{\tau}-\beta_{\tau})^{2}}{\Gamma_{\tau}\alpha_{\tau}}\left\|\nabla\ell(\theta_{\tau}^{md})+\epsilon_{\tau}\right\|^{2}+\lambda_{k}\left\|\epsilon_{k}\right\|\cdot\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|}$
(A.5)

for every $k\geq 1$. Using the definition of $C_{k}$ in 3.6 and summing up the
above inequalities, we have

$\ell(\theta_{N})\ \leq\ell(\theta_{0})-\sum_{k=1}^{N}\lambda_{k}(1\
-L\lambda_{k})\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}$
$+\frac{L}{2}\sum_{k=1}^{N}\Gamma_{k}\sum_{\tau=1}^{k}\frac{(\lambda_{\tau}-\beta_{\tau})^{2}}{\Gamma_{\tau}\alpha_{\tau}}\left\|\nabla\ell(\theta_{\tau}^{m})+\epsilon_{k}\right\|^{2}+\sum_{k=1}^{N}\lambda_{k}\epsilon_{k}\cdot\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|$
$=\ell(\theta_{0})-\sum_{k=1}^{N}\lambda_{k}(1\
-L\lambda_{k})\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}$
$+\frac{L}{2}\sum_{k=1}^{N}\frac{(\lambda_{k}-\beta_{k})^{2}}{\Gamma_{k}\alpha_{k}}(\sum_{\tau=k}^{N}\Gamma_{\tau})\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\|^{2}+\sum_{k=1}^{N}\lambda_{k}\left\|\epsilon_{k}\right\|\cdot\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|$
$=\ell(\theta_{0})-{\displaystyle\sum\lambda_{k}C_{k}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}}+\sum_{k=1}^{N}\lambda_{k}\left\|\epsilon_{k}\right\|\cdot\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|$
(A.6)

Rearranging the terms in the above inequality

${\displaystyle\sum_{k=1}^{N}\lambda_{k}C_{k}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}\leq\ell(\theta_{0})-\ell(\theta^{*})}+\sum_{k=1}^{N}\lambda_{k}\left\|\epsilon_{k}\right\|\cdot\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|$

By assumption 4 that $\left\|\nabla\ell(\cdot)\right\|$ and
$\sum_{k=1}^{N}\lambda_{k}\left\|\epsilon_{k}\right\|$ are bounded. Since
$\ell(\theta_{N})\geq\ell(\theta^{*})$ and in view of the assumption that
$C_{k}>0$, we obtain for some constant $B$,

${\displaystyle\mathrm{min_{k=1,.N}}\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\|^{2}\leq\frac{\ell(\theta_{0})-\ell(\theta^{*})+B}{\sum_{k=1}^{N}\lambda_{k}C_{k}}}$

which clearly implies 3.7.

We now prove part b).

First, from L-Lipschitz-continuous gradient property 3.4, we have

$\ell(\theta_{k}^{ag})\leq\ell(\theta_{k}^{md})+\left\langle\nabla\ell(\theta_{k}^{md}),\
\theta_{k}^{ag}-\theta_{k}^{md}\right\rangle+\frac{L}{2}\left\|\theta_{k}^{ag}-\theta_{k}^{md}\right\|^{2}$
$\leq\ell(\theta_{k}^{md})-{\displaystyle\beta_{k}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}+\left\|\epsilon_{k}\right\|\beta_{k}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|+\frac{L\beta_{k}^{2}}{2}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}.}$
(A.7)

By the assumption that $\ell(\cdot)$ is convex and 3.2,

$\ell(\theta_{k}^{md})-\left[(1-\alpha_{k})\ell(\theta_{k-1}^{ag})+\alpha_{k}\ell(\theta)\right]$
$=\alpha_{k}\left[\ell(\theta_{k}^{md})-\ell(\theta)\right]+(1-\alpha_{k})\left[\ell(\theta_{k}^{md})-\ell(\theta_{k-1}^{ag})\right]$
$\leq\alpha_{k}\left\langle\nabla\ell(\theta_{k}^{md}),\
\theta_{k}^{md}-\theta\right\rangle+(1-\alpha_{k})\left\langle\nabla\ell(\theta_{k}^{md}),\
\theta_{k}^{md}-\theta_{k-1}^{ag}\right\rangle$
$=\langle\nabla\ell(\theta_{k}^{md}),\ \alpha_{k}(\theta_{k}^{md}-\theta)\
+(1-\alpha_{k})(\theta_{k}^{md}-\theta_{k-1}^{ag})\rangle$
$=\alpha_{k}\left\langle\nabla\ell(\theta_{k}^{md}),\
\theta_{k-1}-\theta\right\rangle.$ (A.8)

From 3.3, we have

$\left\|\theta_{k}-\theta\right\|^{2}=\left\|\theta_{k-1}-\lambda_{k}\widehat{\nabla\ell(\theta_{k}^{md})}-\theta\right\|^{2}$
$=\left\|\theta_{k-1}-\theta\right\|^{2}-2\lambda_{k}\langle\widehat{\nabla\ell(\theta_{k}^{md})},\
\theta_{k-1}-\theta\rangle+\lambda_{k}^{2}\left\|\widehat{\nabla\ell(\theta_{k}^{md})}\right\|^{2},$
$=\left\|\theta_{k-1}-\theta\right\|^{2}-2\lambda_{k}\langle\nabla\ell(\theta_{k}^{md})+\epsilon_{k},\
\theta_{k-1}-\theta\rangle+\lambda_{k}^{2}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2},$

which implies

$\alpha_{k}\left\langle\nabla\ell(\theta_{k}^{md})+\epsilon_{k},\
\theta_{k-1}-\theta\right\rangle\
=\frac{\alpha_{k}}{2\lambda_{k}}\left[\left\|\theta_{k-1}-\theta\right\|^{2}-\left\|\theta_{k}-\theta\right\|^{2}\right]$
$+{\displaystyle\frac{\alpha_{k}\lambda_{k}}{2}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}}.$

Hence we obtain

$\alpha_{k}\left\langle\nabla\ell(\theta_{k}^{md}),\
\theta_{k-1}-\theta\right\rangle\
\leq\frac{\alpha_{k}}{2\lambda_{k}}\left[\left\|\theta_{k-1}-\theta\right\|^{2}-\left\|\theta_{k}-\theta\right\|^{2}\right]$
$+{\displaystyle\frac{\alpha_{k}\lambda_{k}}{2}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}}+\alpha_{k}\left\|\epsilon_{k}\right\|\left\|\theta_{k-1}-\theta\right\|$
(A.9)

Using the results of A.7, A.8, and A.9, we get

${\displaystyle\ell(\theta_{k}^{ag})\leq(1-\alpha_{k})\ell(\theta_{k-1}^{ag})+\alpha_{k}\ell(\theta)+\frac{\alpha_{k}}{2\lambda_{k}}\left[\left\|\theta_{k-1}-\theta\right\|^{2}-\left\|\theta_{k}-\theta\right\|^{2}\right]+\alpha_{k}\left\|\epsilon_{k}\right\|\left\|\theta_{k-1}-\theta\right\|}$
$-\beta_{k}(1-\frac{L\beta_{k}}{2}-\frac{\alpha_{k}\lambda_{k}}{2\beta_{k}})\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}+\left\|\epsilon_{k}\right\|\beta_{k}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|$
$\leq\ (1\
-\alpha_{k})\ell(\theta_{k-1}^{ag})+\alpha_{k}\ell(\theta)+\frac{\alpha_{k}}{2\lambda_{k}}\left[\left\|\theta_{k-1}-\theta\right\|^{2}-\left\|\theta_{k}-\theta\right\|^{2}\right]$
$-{\displaystyle\frac{\beta_{k}}{2}(1\
-L\beta_{k})\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}}+\left\|\epsilon_{k}\right\|\beta_{k}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|+\alpha_{k}\left\|\epsilon_{k}\right\|\left\|\theta_{k-1}-\theta\right\|,$
(A.10)

where the last inequality follows from the assumption in 3.8. Subtracting
$\ell(\theta)$ from both sides of the above inequality and using Lemma 1, we
conclude that

${\displaystyle\ell(\theta_{N}^{ag})-\ell(\theta)\leq\Gamma_{N}\left[\sum_{k=1}^{N}\frac{\alpha_{k}}{2\lambda_{k}\Gamma_{k}}\left[\left\|\theta_{k-1}-\theta\right\|^{2}-\left\|\theta_{k}-\theta\right\|^{2}\right]\right.}$
$\left.-\sum_{k=1}^{N}\frac{\beta_{k}}{2\Gamma_{k}}(1-L\beta_{k})\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}+\sum_{k=1}^{N}\frac{1}{\Gamma_{k}}\left[\left\|\epsilon_{k}\right\|\beta_{k}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|+\alpha_{k}\left\|\epsilon_{k}\right\|\left\|\theta_{k-1}-\theta\right\|\right]\right]$
$\leq\Gamma_{N}\frac{\left\|\theta_{0}-\theta\right\|^{2}}{2\lambda_{1}}-\Gamma_{N}\sum_{k=1}^{N}\frac{\beta_{k}}{2\Gamma_{k}}\
(1\ -L\beta_{k})\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}$
$+\Gamma_{N}\sum_{k=1}^{N}\frac{1}{\Gamma_{k}}\left[\left\|\epsilon_{k}\right\|\beta_{k}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|+\alpha_{k}\left\|\epsilon_{k}\right\|\left\|\theta_{k-1}-\theta\right\|\right]$
(A.11)

for every $\theta\in\mathbb{R}^{n}.$ By our contruction 3.9 that sequence
$\left\\{\frac{\alpha_{k}}{\lambda_{k}\Gamma_{k}}\right\\}$is decreasing and
the fact that $\alpha_{1}=\Gamma_{1}=1$, we have

${\displaystyle\sum_{k=1}^{N}\frac{\alpha_{k}}{\lambda_{k}\Gamma_{k}}\left[\left\|\theta_{k-1}-\theta\right\|^{2}-\left\|\theta_{k}-\theta\right\|^{2}\right]\leq\frac{\alpha_{1}\left\|\theta_{0}-\theta\right\|^{2}}{\lambda_{1}\Gamma_{1}}=\frac{\left\|\theta_{0}-\theta\right\|^{2}}{\lambda_{1}}}$
(A.12)

which immediately implies the last inequality of A.11.

Hence, we can conclude 3.11 from the above inequality and the assumption in
3.8:

${\displaystyle\ell(\theta_{N}^{ag})-\ell(\theta^{*})\leq\Gamma_{N}\left[\frac{\left\|\theta_{0}-\theta^{*}\right\|^{2}}{\lambda_{1}}+\sum_{k=1}^{N}\Gamma_{k}^{-1}\left[\left\|\epsilon_{k}\right\|\beta_{k}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|+\alpha_{k}\left\|\epsilon_{k}\right\|\left\|\theta_{k-1}-\theta\right\|\right]\right]}$

Finally, noting the fact that $\ell(\theta_{N}^{ag})\geq\ell(\theta^{*})$,
substitute $\theta:=\theta^{*}$, re-arranging the terms in A.11 we obtain

${\displaystyle\sum_{k=1}^{N}\frac{\beta_{k}}{2\Gamma_{k}}(1\
-L\beta_{k})\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}\mbox{\,\
}k=1,\ldots,N}$
$\leq\frac{\|\theta^{*}-\theta_{0}\|^{2}}{2\lambda_{1}}+\sum_{k=1}^{N}\frac{1}{\Gamma_{k}}\left[\left\|\epsilon_{k}\right\|\beta_{k}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|+\alpha_{k}\left\|\epsilon_{k}\right\|\left\|\theta_{k-1}-\theta\right\|\right],$

or

$\mathrm{min_{k=1,.N}}\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}{\displaystyle\|^{2}\leq
2\frac{\frac{\|\theta^{*}-\theta_{0}\|^{2}}{2\lambda_{1}}+\sum_{k=1}^{N}\frac{1}{\Gamma_{k}}\left[\left\|\epsilon_{k}\right\|\beta_{k}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|+\alpha_{k}\left\|\epsilon_{k}\right\|\left\|\theta_{k-1}-\theta\right\|\right]}{\sum_{k=1}^{N}\Gamma_{k}^{-1}\beta_{k}(1-L\beta_{k})}}$

which together with 3.8, clearly imply 3.10. ∎

###  A.2 Proof of Theorem 4

###### Proof.

We first prove part a). Note that by choosing

$\beta_{k}=\frac{1}{2L}$
${\displaystyle{\displaystyle\Gamma_{k}=\frac{1}{k^{1+\delta}}}},$ (A.13)

which implies that for sufficient large $k$

$\sum_{\tau=k}^{N}\Gamma_{\tau}=\sum_{\tau=k}^{N}\frac{1}{\tau^{1+\delta}}=O(\frac{1}{k^{\delta}})$

We also have

$1-\alpha_{k}=\frac{(k-1)^{1+\delta}}{k^{1+\delta}}$ (A.14)

for every $k>1,$ or
$\alpha_{k}=\frac{\left(k^{1+\delta}-\left(k-1\right)^{1+\delta}\right)}{k^{1+\delta}}=O(\frac{\left(1+\delta\right)k^{\delta}}{k^{1+\delta}})=O(\frac{1}{k})$.
If we choose $\lambda_{k}$ such that $\lambda_{k}-\beta_{k}=o(k^{-1})$ then

$\frac{(\lambda_{k}-\beta_{k})^{2}}{2\alpha_{k}\Gamma_{k}\lambda_{k}}(\sum_{\tau=k}^{N}\Gamma_{\tau})=\frac{o(k^{-2})}{k^{-1}k^{-(1+\delta)}}\frac{1}{k^{\delta}}=o(1)$

so for sufficiently large $k$ we have

$C_{k}=1-L[\lambda_{k}+\frac{(\lambda_{k}-\beta_{k})^{2}}{2\alpha_{k}\Gamma_{k}\lambda_{k}}(\sum_{\tau=k}^{N}\Gamma_{\tau})]>\frac{1}{4}$

Hence, it can also be seen from 3.7 that for some positive bounded constant
$B_{2}$,

${\displaystyle\mathrm{min_{k=1,.N}}\left\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}\right\|^{2}\leq\frac{\ell(\theta_{0})-\ell^{*}+B}{NB_{2}}=O(\frac{1}{N})},$

which concludes the first part of the proof. Since
$\left\|\epsilon_{k}\right\|=O\left(\tau^{2}\right)\leq O(\frac{1}{k})$, we
have $\nabla\ell(\theta_{k}^{md})$ converge to $0$ at the rate of

$\min\left\\{O(\frac{1}{\sqrt{N}}),O\left(\left\|\epsilon_{k}\right\|\right)\right\\}=O(\frac{1}{\sqrt{N}}),$

which gives us the desired result.

We now show part b). Let
$\lambda_{k}=\left(k^{1+\delta}-\left(k-1\right)^{1+\delta}\right)c$ for some
constant $c$ then

$\frac{\alpha_{1}}{\lambda_{1}\Gamma_{1}}=\frac{\alpha_{2}}{\lambda_{2}\Gamma_{2}}=\cdots=\frac{\alpha_{k}}{\lambda_{k}\Gamma_{k}}.$

Observe that

$\alpha_{k}\lambda_{k}=\frac{c^{2}\left(k^{1+\delta}-\left(k-1\right)^{1+\delta}\right)^{2}}{k^{1+\delta}}=\frac{c^{2}\left(1+\delta\right)^{2}O(k^{2\delta})}{k^{1+\delta}}\rightarrow
0$

for $\delta<1$ so
$\frac{\left(1+\delta\right)^{2}k^{2\delta}}{k^{1+\delta}}<\beta_{k}=\frac{1}{2L}$
for sufficient large $k$, which implies that conditions 3.8 and 3.9 hold.
Moreover, it can also be easily seen from A.3 that

$\mathrm{min_{k=1,.N}}\|\nabla\ell(\theta_{k}^{md})+\epsilon_{k}{\displaystyle\|^{2}\leq\frac{\frac{\|\theta^{*}-\theta_{0}\|^{2}}{2\lambda_{1}}+C\sum_{k=1}^{N}\Gamma_{k}^{-1}\left[\left\|\epsilon_{k}\right\|+O(\frac{1}{k})\left\|\epsilon_{k}\right\|\right]}{\sum_{k=1}^{N}\Gamma_{k}^{-1}}}=O(N^{-2-\delta}).$

The last equality is due to the fact that
$\sum_{k=1}^{N}\Gamma_{k}^{-1}=\sum_{k=1}^{N}k^{(1+\delta)}=O(N^{2+\delta})$.
Combining the above relation with 3.7, and since
$\left\|\epsilon_{k}\right\|=O\left(\tau^{2}\right)\leq
O(\frac{1}{k^{2+\delta+\delta_{1}}})$ for some $\delta_{1}>0$, we have
$\nabla\ell(\theta_{k}^{md})$ converge to $0$ at the rate of
$O\left(\sqrt{\frac{1}{N^{2+\delta}}}\right).$

Since $\alpha_{k}\lambda_{k}<\beta_{k}=\frac{1}{2L}$, we have $\delta\leq 1$
which implies that the best convergence rate is $O(N^{-3})$. ∎

## Acknowledgements

This research was funded in part by National Science Foundation grant
DMS-1308919.

## References

  * Andrieu et al., [2010] Andrieu, C., Doucet, A., and Holenstein, R. (2010).  Particle Markov chain Monte Carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(3):269–342.
  * Andrieu and Vihola, [2015] Andrieu, C. and Vihola, M. (2015).  Convergence properties of pseudo-marginal Markov chain Monte Carlo algorithms.  The Annals of Applied Probability, 25(2):1030–1077.
  * Bhadra, [2010] Bhadra, A. (2010).  Discussion of ‘particle Markov chain Monte Carlo methods’ by C. Andrieu, A. Doucet and R. Holenstein.  Journal of the Royal Statistical Society, Series B (Statistical Methodology), 72:314–315.
  * Bretó, [2014] Bretó, C. (2014).  On idiosyncratic stochasticity of financial leverage effects.  Statistics & Probability Letters, 91:20–26.
  * Bretó and Ionides, [2011] Bretó, C. and Ionides, E. L. (2011).  Compound Markov counting processes and their applications to modeling infinitesimally over-dispersed systems.  Stoch. Process. Their Appl., 121(11):2571–2591.
  * Dahlin et al., [2015] Dahlin, J., Lindsten, F., and Schön, T. B. (2015).  Particle Metropolis-Hastings using gradient and Hessian information.  Statistics and Computing, 25(1):81–92.
  * Doucet et al., [2002] Doucet, A., Godsill, S. J., and Robert, C. P. (2002).  Marginal maximum a posteriori estimation using Markov chain Monte Carlo.  Statistics and Computing, 12:77–84.
  * Doucet et al., [2013] Doucet, A., Jacob, P. E., and Rubenthaler, S. (2013).  Derivative-free estimation of the score vector and observed information matrix with application to state-space models.  ArXiv:1304.5768.
  * Eddelbuettel et al., [2011] Eddelbuettel, D., François, R., Allaire, J., Ushey, K., Kou, Q., Russel, N., Chambers, J., and Bates, D. (2011).  Rcpp: Seamless r and c++ integration.  Journal of Statistical Software, 40(8):1–18.
  * Gaetan and Yao, [2003] Gaetan, C. and Yao, J.-F. (2003).  A multiple-imputation Metropolis version of the EM algorithm.  Biometrika, 90(3):643–654.
  * Ghadimi and Lan, [2016] Ghadimi, S. and Lan, G. (2016).  Accelerated gradient methods for nonconvex nonlinear and stochastic programming.  Mathematical Programming, 156(1-2):59–99.
  * Girolami and Calderhead, [2011] Girolami, M. and Calderhead, B. (2011).  Riemann manifold Langevin and Hamiltonian Monte Carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(2):123–214.
  * Ionides et al., [2011] Ionides, E. L., Bhadra, A., Atchadé, Y., and King, A. (2011).  Iterated filtering.  Annals of Statistics, 39:1776–1802.
  * Ionides et al., [2006] Ionides, E. L., Bretó, C., and King, A. A. (2006).  Inference for nonlinear dynamical systems.  Proceedings of the National Academy of Sciences of the USA, 103:18438–18443.
  * Ionides et al., [2015] Ionides, E. L., Nguyen, D., Atchadé, Y., Stoev, S., and King, A. A. (2015).  Inference for dynamic and latent variable models via iterated, perturbed Bayes maps.  Proceedings of the National Academy of Sciences of the USA, 112(3):719–724.
  * Jacquier et al., [2007] Jacquier, E., Johannes, M., and Polson, N. (2007).  MCMC maximum likelihood for latent state models.  Journal of Econometrics, 137(2):615–640.
  * King et al., [2016] King, A. A., Nguyen, D., and Ionides, E. L. (2016).  Statistical inference for partially observed Markov processes via the R package pomp.  Journal of Statistical Software, 69(12).
  * Kloeden and Platen, [1999] Kloeden, P. E. and Platen, E. (1999).  Numerical Soluion of Stochastic Differential Equations.  Springer, New York, 3rd edition.
  * Kushner and Clark, [1978] Kushner, H. J. and Clark, D. S. (1978).  Stochastic Approximation Methods for Constrained and Unconstrained Systems.  Springer-Verlag, New York.
  * Laneri et al., [2010] Laneri, K., Bhadra, A., Ionides, E. L., Bouma, M., Dhiman, R. C., Yadav, R. S., and Pascual, M. (2010).  Forcing versus feedback: Epidemic malaria and monsoon rains in Northwest India.  PLoS Computational Biology, 6(9):e1000898.
  * Leggetter and Woodland, [1995] Leggetter, C. J. and Woodland, P. C. (1995).  Maximum likelihood linear regression for speaker adaptation of continuous density hidden Markov models.  Computer Speech & Language, 9(2):171–185.
  * Lele et al., [2007] Lele, S. R., Dennis, B., and Lutscher, F. (2007).  Data cloning: Easy maximum likelihood estimation for complex ecological models using Bayesian Markov chain Monte Carlo methods.  Ecology Letters, 10(7):551–563.
  * Lindström, [2013] Lindström, E. (2013).  Tuned iterated filtering.  Statistics & Probability Letters, 83(9):2077–2080.
  * Lindström et al., [2012] Lindström, E., Ionides, E. L., Frydendall, J., and Madsen, H. (2012).  Efficient iterated filtering.  In 16th IFAC Symposium on System Identification.
  * Nemeth et al., [2013] Nemeth, C., Fearnhead, P., and Mihaylova, L. (2013).  Particle approximations of the score and observed information matrix for parameter estimation in state space models with linear computational cost.  ArXiv:1306.0735.
  * Nemeth et al., [2014] Nemeth, C., Fearnhead, P., and Mihaylova, L. (2014).  Sequential monte carlo methods for state and parameter estimation in abruptly changing environments.  IEEE Transactions on Signal Processing, 62(5):1245–1255.
  * Nesterov, [2005] Nesterov, Y. (2005).  Smooth minimization of non-smooth functions.  Mathematical Programming, 103(1):127–152.
  * Nesterov, [2013] Nesterov, Y. (2013).  Introductory Lectures on Convex Optimization: A Basic Course, volume 87.  Springer Science & Business Media.
  * Nguyen, [2016] Nguyen, D. (2016).  Another look at Bayes map iterated filtering.  Statistics & Probability Letters.
  * Nguyen and Ionides, [2017] Nguyen, D. and Ionides, E. L. (2017).  A second-order iterated smoothing algorithm.  Statistical Computing.
  * Pitt et al., [2012] Pitt, M. K., dos Santos Silva, R., Giordani, P., and Kohn, R. (2012).  On some properties of Markov chain Monte Carlo simulation methods based on the particle filter.  Journal of Econometrics, 171(2):134–151.
  * Poyiadjis et al., [2011] Poyiadjis, G., Doucet, A., and Singh, S. S. (2011).  Particle approximations of the score and observed information matrix in state space models with application to parameter estimation.  Biometrika, 98(1):65–80.
  * R Core Team, [2013] R Core Team (2013).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Roy et al., [2013] Roy, M., Bouma, M. J., Ionides, E. L., Dhiman, R. C., and Pascual, M. (2013).  The potential elimination of plasmodium vivax malaria by relapse treatment: Insights from a transmission model and surveillance data from NW India.  PLoS Neglected Tropical Diseases, 7(1):e1979.