A second-order iterated smoothing algorithm

Abstract

Simulation-based inference for partially observed stochastic dynamic models is currently receiving much attention due to the fact that direct computation of the likelihood is not possible in many practical situations. Iterated filtering methodologies enable maximization of the likelihood function using simulation-based sequential Monte Carlo filters. Doucet et al. (2013) developed an approximation for the first and second derivatives of the log likelihood via simulation-based sequential Monte Carlo smoothing and proved that the approximation has some attractive theoretical properties. We investigated an iterated smoothing algorithm carrying out likelihood maximization using these derivative approximations. Further, we developed a new iterated smoothing algorithm, using a modification of these derivative estimates, for which we establish both theoretical results and effective practical performance. On benchmark computational challenges, this method beat the first-order iterated filtering algorithm. The method’s performance was comparable to a recently developed iterated filtering algorithm based on an iterated Bayes map. Our iterated smoothing algorithm and its theoretical justification provide new directions for future developments in simulation-based inference for latent variable models such as partially observed Markov process models.

Appendix: Proofs

1.1 Proof of Theorem 3

Let

$$\begin{aligned} R=\left[ \begin{array}{cccc} \tau _{0}I_{d\times d} &{}\quad 0_{d\times d} &{}\quad \cdots &{}\quad 0_{d\times d}\\ \tau _{0}I_{d\times d} &{}\quad \tau _{1}I_{d\times d} &{}\quad \ddots &{}\quad \vdots \\ \vdots &{}\quad \vdots &{}\quad \ddots &{}\quad \vdots \\ \tau _{0}I_{d\times d} &{}\quad \tau _{1}I_{d\times d} &{}\quad \cdots &{}\quad \tau _{N}I_{d\times d} \end{array}\right] , \end{aligned}$$

(8)

where $I_{d\times d}$ is identity matrix of dimension d and $0_{d\times d}$ is zero matrix of dimension d, then a random walk noise will be $R\tau Z_{0:N}$. From Assumption 5, $\breve{\ell }$ is four times continuously differentiable for $\theta ^{[N+1]}$. Since N is fixed, we can apply Theorem 1, with $\Sigma =\mathrm {Cov}(RZ_{0:N})=\breve{\Psi }_N$, to obtain the existence of an $\eta$ and a $C_8$ such that for every $\tau < \eta$ we have,

$$\begin{aligned}&\left| \breve{\mathbb {E}}\left( \breve{\Theta }_{0:N}-\theta ^{\left[ N+1\right] }\left| \breve{Y}_{1:N}=y^*_{1:N}\right. \right) -\tau ^{2}\breve{\Psi }_N\nabla {\breve{\ell }}\left( \theta ^{\left[ N+1\right] }\right) \right| \nonumber \\&\quad <C_8\tau ^{4}, \end{aligned}$$

(9)

where

$$\begin{aligned} \breve{\Psi }_N=\left[ \begin{array}{cccc} \tau _{0}^2\Psi &{}\quad \tau _{0}^2\Psi &{}\quad \cdots &{}\quad \tau _{0}^2\Psi \\ \tau _{0}^2\Psi &{}\quad \tau _{0}^2+\tau _{1}^2\Psi &{}\quad \ddots &{}\quad \tau _{0}^2+\tau _{1}^2\Psi \\ \vdots &{}\quad \vdots &{}\quad \ddots &{}\quad \vdots \\ \tau _{0}^2\Psi &{}\quad \tau _{0}^2+\tau _{1}^2\Psi &{}\quad \cdots &{}\quad \sum _{i=0}^N\tau _{i}^2\Psi \end{array}\right] . \end{aligned}$$

(10)

Note that Assumptions 1 and 2 are automatically satisfied for the multivariate normal distribution with mean zero and variance $\breve{\Psi }_N$, corresponding to the random variable $RZ_{0:N}$. As a result, for fixed $\tau _{0},\dots ,\tau _{N}$ and for a random walk noise, we have

$$\begin{aligned}&\left| \nabla {\breve{\ell }}\left( \theta ^{\left[ N+1\right] }\right) -\tau ^{-2}\breve{\Psi }_N^{-1}\breve{\mathbb {E}}\left( \breve{\Theta }_{0:N}-\theta ^{\left[ N+1\right] }\left| \breve{Y}_{1:N}=y^*_{1:N}\right. \right) \right| \\&\quad <C_9\tau ^{2}. \end{aligned}$$

An application of the Gaussian-Jordan inverse method gives

$$\begin{aligned}&\breve{\Psi }_N^{-1}=\\&\quad \left[ \begin{array}{cccc} (\tau _{0}^{-2}+\tau _{1}^{-2})\Psi ^{-1} &{}\quad -\tau _{1}^{-2}\Psi ^{-1} &{}\quad \cdots &{}\quad 0\\ -\tau _{1}^{-2}\Psi ^{-1} &{}\quad (\tau _{1}^{-2}+\tau _{2}^{-2})\Psi ^{-1} &{}\quad \cdots &{}\quad \vdots \\ 0 &{}\quad -\tau _{2}^{-2}\Psi ^{-1}&{}\quad \cdots &{}\quad \vdots \\ \vdots &{}\quad \vdots &{}\quad \vdots &{}\quad \vdots \\ 0 &{}\quad 0 &{}\quad (\tau _{N-1}^{-2}+\tau _{N}^{-2})\Psi ^{-1} &{}\quad -\tau _{N}^{-2}\Psi ^{-1}\\ 0 &{}\quad 0 &{}\quad -\tau _{N}^{-2}\Psi ^{-1} &{}\quad \tau _{N}^{-2}\Psi ^{-1} \end{array}\right] . \end{aligned}$$

We write $\nabla _{n} \breve{\ell }(\theta ^{[N+1]})$ for the d-dimensional vector of partial derivatives of $\breve{\ell }(\theta ^{[N+1]})$ with respect to each of the d components of $\theta _{n}$. An application of the chain rule gives the identity

$$\begin{aligned} \nabla {{\ell }}\left( \theta \right) =\sum _{n=0}^{N}\nabla _{n} \breve{\ell }\left( \theta ^{[N+1]}\right) , \end{aligned}$$

giving rise to an inequality,

$$\begin{aligned}&\left| \nabla {{\ell }}\left( \theta \right) -\tau ^{-2}\sum _{n=0}^{N}\left\{ \breve{\Psi }_N^{-1}\breve{\mathbb {E}}\left( \breve{\Theta }_{0:N}-\theta ^{[N+1]}|\breve{Y}_{1:N}=y^*_{1:N}\right) \right\} _{n}\right| \\&\quad <C_{10}\tau ^{2}, \end{aligned}$$

where $\left\{ s\right\} _{n}$ is the entries $\left\{ dn+1,...,d(n+1)\right\}$ of a vector $s\in R^{d(N+1)}$. Decomposing the matrix multiplication by $\breve{\Psi }_N^{-1}$ into $d\times d$ blocks, we have

$$\begin{aligned}&\tau ^{-2}\sum _{n=0}^{N}\left\{ \breve{\Psi }_N^{-1}\breve{\mathbb {E}}\left( \breve{\Theta }_{0:N}-\theta ^{[N+1]}|\breve{Y}_{1:N}=y^*_{1:N}\right) \right\} _{n} \nonumber \\&\quad =\tau ^{-2}\sum _{n=0}^{N}\mathrm {SumCol}_{n}\left( \breve{\Psi }_N^{-1}\right) \breve{\mathbb {E}}\left( \breve{\Theta }_{n}-\theta |\breve{Y}_{1:N}=y^*_{1:N}\right) ,\nonumber \\ \end{aligned}$$

(11)

where $\mathrm {SumCol}_{n}$ is the sum of the nth column in the $d\times d$ block construction of $\breve{\Psi }_N^{-1}$. Every column of $\breve{\Psi }_N^{-1}$ except the first sums to 0, and this special structure of $\tilde{\Psi }_N^{-1}$ gives a simple form,

$$\begin{aligned}&\left| \sum _{n=0}^{N}\nabla _{n}{\breve{\ell }}\left( \theta ^{[N+1]}\right) -\tau ^{-2}\Psi ^{-1}\tau _{0}^{-2}\breve{\mathbb {E}}\left( \breve{\Theta }_{0}-\theta |\breve{Y}_{1:N}=y^*_{1:N}\right) \right| \\&\quad <C_{11}\tau ^{2}. \end{aligned}$$

This can be written as

$$\begin{aligned} \left| \nabla {\ell }\left( \theta \right) -\tau ^{-2}\Psi ^{-1}\tau _{0}^{-2}\breve{\mathbb {E}}\left( \breve{\Theta }_{0}-\theta |\breve{Y}_{1:N}=y^*_{1:N}\right) \right| <C_{11}\tau ^{2}. \end{aligned}$$

1.2 Proof of Theorem 4

Using similar set up as above, let the random walk noise be $R\tau Z_{0:N}$ with R defined as in Eq. (8). By selecting $p_{\Theta _{0:N}}$ to follow a multivariate normal distribution, Assumption 4 is also satisfied. From Theorem 2, for fixed $\tau _{0}, \dots , \tau _{N}$, there exist $\eta$ and $C_{12}$ such that for $0<\tau <\eta$,

$$\begin{aligned}&\left| \nabla ^2{\breve{\ell }}\left( \theta ^{[N+1]}\right) -\tau ^{-4}\right. \nonumber \\&\quad \left. \left[ \breve{\Psi }_N^{-1}\left( {\breve{\mathrm{C}}\mathrm {ov}}_{\theta ^{[N+1]},\tau }\left( \breve{\Theta }_{0:N}|\breve{Y}_{1:N}=y^*_{1:N}\right) -\tau ^{2}\breve{\Psi }_N\right) \breve{\Psi }_N^{-1}\right] \right| \nonumber \\&\qquad <C_{12}\tau ^{2}. \end{aligned}$$

(12)

Define $\nabla ^2_{s,n}{\breve{\ell }}\left( \theta ^{[N+1]}\right)$ as

$$\begin{aligned} \nabla ^{2}_{s,n}\breve{\ell }\left( \theta ^{[N+1]}\right) =\frac{{{\partial }}^{2}\breve{\ell }\left( \theta ^{[N+1]}\right) }{{{\partial \theta }}_{s}{{\partial \theta }}_{n}}. \end{aligned}$$

Applying the chain rule, we have

$$\begin{aligned} \nabla ^2{\ell }\left( \theta \right) =\sum _{s=0}^{N}\sum _{n=0}^{N}\nabla _{s,n}^2{\breve{\ell }}\left( \theta ^{[N+1]}\right) . \end{aligned}$$

Adding up term in Eq. (12), we get

$$\begin{aligned}&\left| \sum _{s=0}^{N}\sum _{n=0}^{N}\nabla _{s,n}^{2}{\breve{\ell }}(\theta ^{[N+1]})-\tau ^{-4}\right. \\&\left. \sum _{s=0}^{N}\sum _{n=0}^{N}\left[ \breve{\Psi }_N^{-1}\left( {\breve{\mathrm{C}}\mathrm {ov}}_{\theta ^{[N+1]},\tau }\left( \breve{\Theta }_{0:N}|\breve{Y}_{1:N}=y^*_{1:N}\right) -\tau ^{2}\breve{\Psi }_N\right) \breve{\Psi }_N^{-1}\right] _{s,n}\right| \\&\quad <C_{14}\tau ^{2}, \end{aligned}$$

where $\left\{ A \right\} _{s,n}$ is the entries of rows $\left\{ ds+1,...,d(s+1)\right\}$ and of columns $\left\{ dn+1,...,d(n+1)\right\}$ of a matrix $A\in R^{d(N+1)\times d(N+1)}$. Therefore,

$$\begin{aligned}&\Bigg |\nabla ^2{\ell }\left( \theta \right) -\tau ^{-4}\Bigg .\\&\left. \sum _{s=0}^{N}\sum _{n=0}^{N}\left[ \breve{\Psi }_N^{-1}\left( {\breve{\mathrm{C}}\mathrm {ov}}_{\theta ^{[N+1]},\tau }\left( \breve{\Theta }_{0:N}|\breve{Y}_{1:N}=y^*_{1:N}\right) -\tau ^{2}\breve{\Psi }_N\right) \breve{\Psi }_N^{-1}\right] _{s,n}\right| \\&\quad <C_{14}\tau ^{2}. \end{aligned}$$

Defining $\mathrm {SumCol}_n$ as in Eq. (11), we have

$$\begin{aligned}&\sum _{s=0}^{N}\sum _{n=0}^{N}\left[ \breve{\Psi }_N^{-1}\left( {\breve{\mathrm{C}}\mathrm {ov}}_{\theta ^{[N+1]},\tau }\left( \breve{\Theta }_{0:N}|\breve{Y}_{1:N}=y^*_{1:N}\right) -\tau ^{2}\breve{\Psi }_N\right) \breve{\Psi }_N^{-1}\right] _{s,n}\\&\quad =\sum _{s=0}^{N}\sum _{n=0}^{N}\mathrm{SumCol}_{\mathrm {s}}(\breve{\Psi }_N^{-1})\mathrm {SumCol}_{n}(\breve{\Psi }_N^{-1})\\&\qquad \times \left( {\breve{\mathrm{C}}\mathrm {ov}}_{\theta ^{[N+1]},\tau }\left( \breve{\Theta }_{s},\breve{\Theta }_{n}|\breve{Y}_{1:N}=y^*_{1:N}\right) -\sum _{k=0}^{s\wedge n}\tau _k^{2}\tau ^2\Psi \right) \Psi ^{-1}\\&\quad =\left( {\breve{\mathrm{V}}\mathrm {ar}}_{\theta ^{[N+1]},\tau }\left( \breve{\Theta }_0|\breve{Y}_{1:N}=y^*_{1:N}\right) -\tau _0^{2}\tau ^2\Psi \right) . \end{aligned}$$

The last equality follows since $\breve{\Psi }_N^{-1}$ is symmetric matrix with block of $d\times d$ for which each column except the first sums to 0. Thus, we obtain

$$\begin{aligned}&\left| \nabla ^2{\ell }\left( \theta \right) -\tau ^{-4}\Psi ^{-1}\left( {\breve{\mathrm{V}}\mathrm {ar}}_{\theta ^{[N+1]},\tau }\left( \breve{\Theta }_{0}|\breve{Y}_{1:N}=y^*_{1:N}\right) -\tau _0^{2}\tau ^2\Psi \right) \Psi ^{-1}\right| \\&\quad <C_{15}\tau ^{2}. \end{aligned}$$

1.3 Proof of Theorem 5

In order to prove Theorem 5, we use the following corollary to Theorem 1.

Corollary 1

Suppose the perturbation kernel takes a value $\kappa \in \mathcal {K}$ satisfying Assumptions 6 and 7. Suppose also assumption 3. There exists an $\eta$ and a constant $C_{16}$ such that for every $0<\tau <\eta$ and every $\kappa \in \mathcal {K}$,

$$\begin{aligned} \left| \breve{\mathbb {E}}\left( \breve{\Theta }-\theta \left| \breve{Y}=y^*\right. \right) -\tau ^{2}\Sigma \nabla \ell \left( \theta \right) \right| <C_{16}\tau ^{4}. \end{aligned}$$

(13)

Corollary 1 follows directly from applying Theorem 1, noting that Assumptions 6 and 7 imply a uniform bound on $C_2$ in Theorem 1. Applying Corollary 1, we obtain the existence of an $\eta$ and $C_{18}$ such that for every $\tau < \eta$ we have

$$\begin{aligned}&\left| \tau ^{2}\breve{\Psi }_N\nabla {\breve{\ell }}\left( \theta ^{\left[ N+1\right] }\right) -\breve{\mathbb {E}}\left( \breve{\Theta }_{0:N}-\theta ^{\left[ N+1\right] }\left| \breve{Y}_{1:N}=y^*_{1:N}\right. \right) \right| \nonumber \\&\quad <C_{18}\tau ^{4}. \end{aligned}$$

(14)

For compactness of notation, we write $E_n=\breve{\mathbb {E}}\left( \breve{\Theta }_{n}-\theta \left| \breve{Y}_{1:N}=y^*_{1:N}\right. \right)$ and $D_n=\Psi \nabla _n{\breve{\ell }}\left( \theta ^{\left[ N+1\right] }\right)$. With $\breve{\Psi }_N$ as in Eq. (10), writing out terms of the vector equation in (14) gives

$$\begin{aligned}&E_0 = \tau ^2 \tau _0^2 \sum _{n=0}^N D_n +O(\tau ^4), \end{aligned}$$

(15)

$$\begin{aligned}&E_1 = \tau ^2 \tau _0^2 \sum _{n=0}^N D_n + \tau ^2 \tau _1^2 \sum _{n=1}^N D_n+O(\tau ^4), \end{aligned}$$

(16)

$$\begin{aligned}&\vdots \end{aligned}$$

(17)

$$\begin{aligned}&E_{N-1} = \tau ^2 \tau _0^2 \sum _{n=0}^N D_n +\cdots + \tau ^2 \tau _{N-1}^2 \sum _{n=N-1}^N D_n +O(\tau ^4), \nonumber \\ \end{aligned}$$

(18)

$$\begin{aligned}&E_N = \tau ^2 \tau _0^2 \sum _{n=0}^N D_n +\cdots + \tau ^2 \tau _{N}^2 D_N +O(\tau ^4). \end{aligned}$$

(19)

Using our assumption that for all $n=1\ldots N$, $\tau _n=O(\tau ^2)$, we get that $E_n=E_0+O(\tau ^{4})$, from which we can conclude that

$$\begin{aligned} \frac{1}{N+1}\sum _{n=0}^NE_n=E_0+O(\tau ^{4}). \end{aligned}$$

From Eq. (14) and the chain rule, for fixed $\tau _0$ we have

$$\begin{aligned} \left| \nabla {\ell }\left( \theta \right) -\tau ^{-2}\Psi ^{-1}\tau _{0}^{-2}\breve{\mathbb {E}}\left( \breve{\Theta }_{0}-\theta |\breve{Y}_{1:N}=y^*_{1:N}\right) \right| <C_{19}\tau ^{2}, \end{aligned}$$

which then completes the proof.

1.4 Proof of Theorem 6

As for Theorem 5, we need a corollary to Theorem 2 over the kernel set $\mathcal {K}$, making use of Assumptions 6 and 7.

Corollary 2

Suppose Assumptions 3, 6, and 7. Suppose also that every kernel in $\mathcal {K}$ satisfies the mesokurtic condition in Assumptions 3. There exists an $\eta$ and a constant $C_{20}$ such that for every $0<\tau <\eta$ and every $\kappa \in \mathcal {K}$,

$$\begin{aligned}&\left| \breve{\mathbb {E}}\left[ \left( \breve{\Theta }-\theta \right) \left( \breve{\Theta }-\theta \right) ^{\top }\left| \breve{Y}=y^*\right. \right] -\tau ^{2}\Sigma -\tau ^{4}\Sigma \left( \nabla ^{2}\ell (\theta )\right) \Sigma \right| \\&\quad < C_{20}\tau ^{6}. \end{aligned}$$

Applying Corollary 2, we have

$$\begin{aligned}&\left| \tau ^{4}\breve{\Psi }_N\nabla ^2{\breve{\ell }}\left( \theta ^{[N+1]}\right) \breve{\Psi }_N\right. \\&\quad \left. -\left( {\breve{\mathrm{C}}\mathrm {ov}}_{\theta ^{[N+1]},\tau }\left( \breve{\Theta }_{0:N}|\breve{Y}_{1:N}{=}y^*_{1:N}\right) {-}\tau ^{2}\breve{\Psi }_N\right) \right| <C_{21}\tau ^{6}. \end{aligned}$$

For compact notation, we write

$$\begin{aligned} \breve{C}ov_{s,n}=\breve{\mathrm{C}}{\mathrm {ov}}\left( \breve{\Theta }_{s},\breve{\Theta }_{n}|\breve{Y}_{1:N}=y^*_{1:N}\right) -\sum _{k=0}^{s\wedge {n}}\tau _k^{2}\tau ^2\Psi \end{aligned}$$

and

$$\begin{aligned} \nabla ^2_{s,n}\breve{\ell }= \nabla ^2_{s,n}\breve{\ell }\left( \theta ^{[N+1]}\right) . \end{aligned}$$

From the diagonal terms of the above matrix norm inequality, we derive $N+1$ equations,

$$\begin{aligned}&\breve{\mathrm{C}}{\mathrm {ov}}_{0,0}=\tau ^4\left[ \breve{\Psi }_N\nabla ^2\breve{\ell }\breve{\Psi }_N\right] _{0,0}+O(\tau ^6), \end{aligned}$$

(20)

$$\begin{aligned}&{\breve{\mathrm{C}}\mathrm {ov}}_{1,1}=\tau ^4\left[ \breve{\Psi }_N\nabla ^2\breve{\ell }\breve{\Psi }_N\right] _{1,1}+O(\tau ^6), \end{aligned}$$

(21)

$$\begin{aligned}&\vdots \end{aligned}$$

(22)

$$\begin{aligned}&{\breve{\mathrm{C}}\mathrm {ov}}_{N-1,N-1}=\tau ^4\left[ \breve{\Psi }_N\nabla ^2\breve{\ell }\breve{\Psi }_N\right] _{N-1,N-1}+O(\tau ^6), \end{aligned}$$

(23)

$$\begin{aligned}&{\breve{\mathrm{C}}\mathrm {ov}}_{N,N}=\tau ^4\left[ \breve{\Psi }_N\nabla ^2\breve{\ell }\breve{\Psi }_N\right] _{N,N}+O(\tau ^6). \end{aligned}$$

(24)

Using (20) through (24), and expanding out a matrix multiplication, we get

$$\begin{aligned}&\left[ \breve{\Psi }_N\nabla ^2\breve{\ell }\breve{\Psi }_N\right] _{n,n} = \Psi ^2\sum _{j=0}^n\\&\quad \left( \sum _{k=0}^j\tau _k^2\right) \left[ \sum _{i=0}^n\left( \sum _{k=0}^i\tau _k^2\right) \nabla ^2_{i,j}\breve{\ell }+\sum _{i=n+1}^N\left( \sum _{k=0}^n\tau _k^2\right) \nabla ^2_{i,j}\breve{\ell }\right] \\&\qquad + \Psi ^2\sum _{j=n+1}^N\left( \sum _{k=0}^n\tau _k^2\right) \\&\quad \left[ \sum _{i=0}^n\left( \sum _{k=0}^i\tau _k^2\right) \nabla ^2_{i,j}\breve{\ell }+\sum _{i=n+1}^N\left( \sum _{k=0}^n\tau _k^2\right) \nabla ^2_{i,j}\breve{\ell }\right] . \end{aligned}$$

Using our assumption that for all $n=1\ldots N$, $\tau _n=O(\tau ^2)$, we get that

$$\begin{aligned} \breve{\mathrm{C}}\mathrm {ov}_{n,n}=\breve{\mathrm{C}}\mathrm {ov}_{0,0}+O(\tau ^{6}), \end{aligned}$$

from which we can conclude that

$$\begin{aligned} \frac{1}{N+1}\sum _{n=0}^N\breve{\mathrm{C}}\mathrm {ov}_{n,n}=\breve{\mathrm{C}}\mathrm {ov}_{0,0}+O(\tau ^{6}). \end{aligned}$$

For $n=0$, we have

$$\begin{aligned} \left[ \breve{\Psi }_N\nabla ^2\breve{\ell }\breve{\Psi }_N\right] _{0,0}=\Psi ^2\tau _0^4\sum _{j=0}^N\sum _{i=0}^N\nabla ^2_{i,j}\breve{\ell }, \end{aligned}$$

(25)

from which, applying the chain rule completes the proof.