T^{-\lambda}}\left|\frac{1}{T}\sum_{i=1}^{T}\Delta_{i}(x_{k}+\eta_{T})-F_{\epsilon}(x_{k}+\eta_{T})+F_{\epsilon}(x_{k}+\eta_{T})-F_{\epsilon}(x_{k})+F_{\epsilon}(x_{k})-\frac{1}{T}\sum_{i=1}^{T}\Delta_{i}(x_{k})\right|\\\
&\leq\sup_{k\leq N_{T}}\sup_{|\eta_{T}|\leq
T^{-\lambda}}\left\\{\left|\frac{1}{T}\sum_{i=1}^{T}\Delta_{i}(x_{k}+\eta_{T})-F_{\epsilon}(x_{k}+\eta_{T})\right|+\left|F_{\epsilon}(x_{k}+\eta_{T})-F_{\epsilon}(x_{k})\right|+\left|F_{\epsilon}(x_{k})-\frac{1}{T}\sum_{i=1}^{T}\Delta_{i}(x_{k})\right|\right\\}.\end{split}$
(68)

Applying the Glivenko–Cantelli theorem again, we can find the first and third
term in the r.h.s. of Eq. 68 converges to 0 a.s.. For the middle term:

$\begin{split}&\sup_{k\leq N_{T}}\sup_{|\eta_{T}|\leq
T^{-\lambda}}\left|F_{\epsilon}(x_{k}+\eta_{T})-F_{\epsilon}(x_{k})\right|\\\
&=\sup_{k\leq N_{T}}\sup_{|\eta_{T}|\leq
T^{-\lambda}}\left|F_{\epsilon}(x_{k})+F_{\epsilon}^{\prime}(\tilde{x})\eta_{T}-F_{\epsilon}(x_{k})\right|\\\
&\leq\sup_{x}\sup_{|\eta_{T}|\leq
T^{-\lambda}}\left|F_{\epsilon}^{\prime}(\tilde{x})\eta_{T}\right|\overset{}{\to}0,~{}\text{Under
A11.}\end{split}$ (69)

We can do a similar analysis for the first term of the r.h.s of Eq. 66.
Combining all parts together, we prove Eq. 16. The proof of Eq. 15 with
heteroscedastic errors is similar. ∎

######  Proof of Theorem 4.1.

For simplifying the notations, we consider the NLAR model with order 1 and
homoscedastic errors. For NLAR models with higher order and heteroscedastic
errors, the proof can be written similarly. We show the proof of $h=2$ as an
example. The proof of higher steps or one-step prediction can be written
similarly.

First, by the tower property, we can show that $F_{X_{T+2}|X_{T}}(x)$ is
equivalent to:

$\begin{split}F_{X_{T+2}|X_{T}}(x)&=\mathbb{P}(X_{T+2}\leq x|X_{T})\\\
&=\mathbb{P}(\phi(X_{T+1},\theta_{1})+\epsilon_{T+2}\leq x|X_{T})\\\
&=\mathbb{P}\left(\epsilon_{T+2}\leq
x-\phi(\phi(X_{T},\theta_{1})+\epsilon_{T+1},\theta_{1})\bigg{|}X_{T}\right)\\\
&=\mathbb{E}\left[\mathbb{P}\left(\epsilon_{T+2}\leq
x-\phi(\phi(X_{T},\theta_{1})+\epsilon_{T+1},\theta_{1})\bigg{|}\epsilon_{T+1},X_{T}\right)\bigg{|}X_{T}\right]\\\
&=\mathbb{E}\left[F_{\epsilon}\left(x-\phi(\phi(X_{T},\theta_{1})+\epsilon_{T+1},\theta_{1})\right)\bigg{|}X_{T}\right]\\\
&=\mathbb{E}\left[F_{\epsilon}\left(\mathscr{L}(x,X_{T},\epsilon_{T+1})\right)\bigg{|}X_{T}\right];\end{split}$
(70)

we use $\mathscr{L}(x,X_{T},\epsilon_{T+1})$ to represent
$x-\phi(\phi(X_{T},\theta_{1})+\epsilon_{T+1},\theta_{1})$ to simplify
notations. Similarly, we can analyze $F_{X^{*}_{T+2}|X_{T},\ldots,X_{0}}(x)$,
it has below equivalent expressions:

$\begin{split}F_{X^{*}_{T+2}|X_{T},\ldots,X_{0}}(x)&=\mathbb{P}(X^{*}_{T+2}\leq
x|X_{T},\ldots,X_{0})\\\
&=\mathbb{E}\left[\mathbb{P}\left(\hat{\epsilon}^{*}_{T+2}\leq\widehat{\mathscr{L}}(x,X_{T},\hat{\epsilon}^{*}_{T+1})\bigg{|}\hat{\epsilon}^{*}_{T+1},X_{T},\ldots,X_{0}\right)\bigg{|}X_{T},\ldots,X_{0}\right]\\\
&=\mathbb{E}^{*}\left[\widehat{F}_{\epsilon}\left(\widehat{\mathscr{L}}(x,X_{T},\hat{\epsilon}^{*}_{T+1})\right)\right],\end{split}$
(71)

where $\widehat{\mathscr{L}}(x,X_{T},\hat{\epsilon}^{*}_{T+1})$ represents
$x-\phi(\phi(X_{T},\widehat{\theta}_{1})+\hat{\epsilon}^{*}_{T+1},\widehat{\theta}_{1})$
and $\mathbb{E}^{*}(\cdot)$ represents the expectation in the bootstrap world,
i.e., $\mathbb{E}(\cdot|X_{T},\ldots,X_{0})$. Thus, we hope to show:

$\sup_{|x|\leq
c_{T}}\bigg{|}\mathbb{E}^{*}\left[\widehat{F}_{\epsilon}(\widehat{\mathscr{L}}(x,X_{T},\hat{\epsilon}^{*}_{T+1}))\right]-\mathbb{E}\left[F_{\epsilon}\left(\mathscr{L}(x,X_{T},\epsilon_{T+1})\right)\bigg{|}X_{T}\right]\bigg{|}\overset{p}{\to}0.$
(72)

From here, we first bound the region of $X_{t}$ by Lemma 1 of Franke et al.,
(2004) under A1–A7:

$\mathbb{P}(|X_{T}|>\gamma_{T})\to 0,$ (73)

where $\\{\gamma_{T}\\}$ is a sequence of sets, such that
$\gamma_{1}\subseteq\cdots\subseteq\gamma_{T}\subseteq\gamma_{T+1}\subseteq\cdots$
with the form $\gamma_{T}=[-T^{c},T^{c}]$; $c$ is some appropriate constant.
For deriving such a result for a time series model with heteroscedastic
errors, we need the additional assumption of variance function in A1; then the
proof is referred to Lemma 1 of Franke et al., (2004). In addition, we have a
relationship:

$\begin{split}&\mathbb{P}\left(\sup_{|x|\leq
c_{T}}\bigg{|}\mathbb{E}^{*}\left[\widehat{F}_{\epsilon}(\widehat{\mathscr{L}}(x,X_{T},\hat{\epsilon}^{*}_{T+1}))\right]-\mathbb{E}\left[F_{\epsilon}\left(\mathscr{L}(x,X_{T},\epsilon_{T+1})\right)\bigg{|}X_{T}\right]\bigg{|}>\varepsilon\right)\\\
&\leq\mathbb{P}(|X_{T}|>\gamma_{T})+\mathbb{P}\left((|X_{T}|\leq\gamma_{T})\bigcap\left(\sup_{|x|\leq
c_{T}}\bigg{|}\mathbb{E}^{*}\left[\widehat{F}_{\epsilon}(\widehat{\mathscr{L}}(x,X_{T},\hat{\epsilon}^{*}_{T+1}))\right]-\mathbb{E}\left[F_{\epsilon}\left(\mathscr{L}(x,X_{T},\epsilon_{T+1})\right)\bigg{|}X_{T}\right]\bigg{|}>\varepsilon\right)\right).\end{split}$
(74)

Thus, to verify Eq. 72, we just need to show that the second term of the
r.h.s. of Eq. 74 converges to 0. We can take the sequence $c_{T}$ and
$\gamma_{T}$ to be the same sequence which converges to infinity slowly
enough. Then, it is enough for us to analyze the asymptotic probability of the
below expression:

$\sup_{|x|,|y|\leq
c_{T}}\bigg{|}\mathbb{E}^{*}\left[\widehat{F}_{\epsilon}(\widehat{\mathscr{L}}(x,y,\hat{\epsilon}^{*}_{T+1}))\right]-\mathbb{E}\left[F_{\epsilon}\left(\mathscr{L}(x,y,\epsilon_{T+1})\right)\right]\bigg{|}>\varepsilon.$
(75)

Decompose the l.h.s. of Eq. 75 as:

$\begin{split}&\sup_{|x|,|y|\leq
c_{T}}\bigg{|}\mathbb{E}^{*}\left[\widehat{F}_{\epsilon}(\widehat{\mathscr{L}}(x,y,\hat{\epsilon}^{*}_{T+1}))\right]-\mathbb{E}\left[F_{\epsilon}\left(\mathscr{L}(x,y,\epsilon_{T+1})\right)\right]\bigg{|}\\\
&\leq\sup_{|x|,|y|\leq
c_{T}}\bigg{|}\mathbb{E}^{*}\left[\widehat{F}_{\epsilon}(\widehat{\mathscr{L}}(x,y,\hat{\epsilon}^{*}_{T+1}))\right]-\mathbb{E}^{*}\left[F_{\epsilon}(\widehat{\mathscr{L}}(x,y,\hat{\epsilon}^{*}_{T+1}))\right]\bigg{|}\\\
&+\sup_{|x|,|y|\leq
c_{T}}\bigg{|}\mathbb{E}^{*}\left[F_{\epsilon}(\widehat{\mathscr{L}}(x,y,\hat{\epsilon}^{*}_{T+1}))\right]-\mathbb{E}\left[F_{\epsilon}\left(\mathscr{L}(x,y,\epsilon_{T+1})\right)\right]\bigg{|}.\end{split}$
(76)

Then, we analyze two terms on the r.h.s. of Eq. 76 separately. For the first
term, we have:

$\begin{split}&\sup_{|x|,|y|\leq
c_{T}}\bigg{|}\mathbb{E}^{*}\left[\widehat{F}_{\epsilon}(\widehat{\mathscr{L}}(x,y,\hat{\epsilon}^{*}_{T+1}))\right]-\mathbb{E}^{*}\left[F_{\epsilon}(\widehat{\mathscr{L}}(x,y,\hat{\epsilon}^{*}_{T+1}))\right]\bigg{|}\\\
&\leq\sup_{|x|,|y|\leq
c_{T}}\mathbb{E}^{*}\bigg{|}\widehat{F}_{\epsilon}(\widehat{\mathscr{L}}(x,y,\hat{\epsilon}^{*}_{T+1}))-F_{\epsilon}(\widehat{\mathscr{L}}(x,y,\hat{\epsilon}^{*}_{T+1}))\bigg{|}\\\
&\leq\sup_{|x|,|y|\leq
c_{T},z}\bigg{|}\widehat{F}_{\epsilon}(\widehat{\mathscr{L}}(x,y,z))-F_{\epsilon}(\widehat{\mathscr{L}}(x,y,z))\bigg{|}\overset{p}{\to}0,~{}\text{under
<EMAIL_ADDRESS>(77)

For the second term on the r.h.s. of Eq. 76, we have:

$\begin{split}&\sup_{|x|,|y|\leq
c_{T}}\bigg{|}\mathbb{E}^{*}\left[F_{\epsilon}(\widehat{\mathscr{L}}(x,y,\hat{\epsilon}^{*}_{T+1}))\right]-\mathbb{E}\left[F_{\epsilon}\left(\mathscr{L}(x,y,\epsilon_{T+1})\right)\right]\bigg{|}\\\
&\leq\sup_{|x|,|y|\leq
c_{T}}\bigg{|}\frac{1}{T}\sum_{i=1}^{T}F_{\epsilon}(\widehat{\mathscr{L}}(x,y,\hat{\epsilon}_{i}))-\frac{1}{T}\sum_{i=1}^{T}F_{\epsilon}(\mathscr{L}(x,y,\epsilon_{i}))\bigg{|}\\\
&+\sup_{|x|,|y|\leq
c_{T}}\bigg{|}\frac{1}{T}\sum_{i=1}^{T}F_{\epsilon}(\mathscr{L}(x,y,\epsilon_{i}))-\mathbb{E}\left[F_{\epsilon}\left(\mathscr{L}(x,y,\epsilon_{T+1})\right)\right]\bigg{|},\end{split}$
(78)

where $\\{\epsilon_{j}\\}_{j=1}^{T}$ are taken as
$X_{i}-\phi(X_{i-1},\theta_{1})$ for $i=1,\ldots,T$. We can show:

$\begin{split}&\mathbb{P}\left(\max_{i=1,\ldots,T}\bigg{|}\epsilon_{i}-\hat{\epsilon}_{i}\bigg{|}>\varepsilon\right)\\\
&=\mathbb{P}\left(\max_{i=1,\ldots,T}\bigg{|}X_{i}-\phi(X_{i-1},\theta_{1})-X_{i}+\phi(X_{i-1},\widehat{\theta}_{1})\bigg{|}>\varepsilon\right)\\\
&\leq\mathbb{P}\left((\max_{i=1,\ldots,T}|X_{i-1}|>c_{T})\right)\\\
&+\mathbb{P}\left(\left(\max_{i=1,\ldots,T}|X_{i-1}|<c_{T}\right)\bigcap\left(\max_{i=1,\ldots,T}\bigg{|}\phi(X_{i-1},\widehat{\theta}_{1})-\phi(X_{i-1},\theta_{1})\bigg{|}>\varepsilon\right)\right)\\\
&\leq o(1)+\mathbb{P}\left(\sup_{|x|\leq
c_{T}}\bigg{|}\phi(x,\widehat{\theta}_{1})-\phi(x,\theta_{1})\bigg{|}>\varepsilon\right)\\\
&\to 0,~{}\text{under A9 and A10.}\end{split}$ (79)

We further consider two terms on the r.h.s. of Eq. 78 separately. For the
first term, by Taylor expansion, we have:

$\begin{split}&\sup_{|x|,|y|\leq
c_{T}}\bigg{|}\frac{1}{T}\sum_{i=1}^{T}F_{\epsilon}(\widehat{\mathscr{L}}(x,y,\hat{\epsilon}_{i}))-\frac{1}{T}\sum_{i=1}^{T}F_{\epsilon}(\mathscr{L}(x,y,\epsilon_{i}))\bigg{|}\\\
&=\sup_{|x|,|y|\leq
c_{T}}\bigg{|}\frac{1}{T}\sum_{i=1}^{T}\left(F_{\epsilon}(\mathscr{L}(x,y,\epsilon_{i}))+f_{\epsilon}(o_{i})(\widehat{\mathscr{L}}(x,y,\hat{\epsilon}_{i})-\mathscr{L}(x,y,\epsilon_{i}))\right)-\frac{1}{T}\sum_{i=1}^{T}F_{\epsilon}(\mathscr{L}(x,y,\epsilon_{i}))\bigg{|}\\\
&=\sup_{|x|,|y|\leq
c_{T}}\bigg{|}\frac{1}{T}\sum_{i=1}^{T}f_{\epsilon}(o_{i})(\widehat{\mathscr{L}}(x,y,\hat{\epsilon}_{i})-\mathscr{L}(x,y,\epsilon_{i}))\bigg{|}\\\
&\leq\sup_{|x|,|y|\leq
c_{T}}\frac{1}{T}\sum_{i=1}^{T}\bigg{|}f_{\epsilon}(o_{i})(\widehat{\mathscr{L}}(x,y,\hat{\epsilon}_{i})-\mathscr{L}(x,y,\epsilon_{i}))\bigg{|}\\\
&\leq\sup_{|x|,|y|\leq
c_{T}}\sup_{z}|f_{\epsilon}(z)|\cdot\frac{1}{T}\sum_{i=1}^{T}\bigg{|}\widehat{\mathscr{L}}(x,y,\hat{\epsilon}_{i})-\mathscr{L}(x,y,\epsilon_{i})\bigg{|}\\\
&\leq\sup_{|x|,|y|\leq
c_{T}}C\cdot\frac{1}{T}\sum_{i=1}^{T}\bigg{|}\widehat{\mathscr{L}}(x,y,\hat{\epsilon}_{i})-\mathscr{L}(x,y,\epsilon_{i})\bigg{|}~{}\text{(under
A11)}\\\ &\leq\sup_{|x|,|y|\leq
c_{T},j\in\\{1,\ldots,T\\}}C\cdot\bigg{|}\widehat{\mathscr{L}}(x,y,\hat{\epsilon}_{j})-\mathscr{L}(x,y,\epsilon_{j})\bigg{|}.\end{split}$
(80)

From Eq. 79 and A9–A12, we have Eq. 80 converges to 0 in probability. For the
second term on the r.h.s. of Eq. 78, by the uniform law of large numbers, we
have:

$\sup_{|x|,|y|\leq
c_{T}}\bigg{|}\frac{1}{T}\sum_{i=1}^{T}F_{\epsilon}(\mathscr{L}(x,y,\epsilon_{i}))-\mathbb{E}\left[F_{\epsilon}\left(\mathscr{L}(x,y,\epsilon_{T+1})\right)\right]\bigg{|}\overset{p}{\to}0.$
(81)

Combine all pieces, Eq. 75 converges to 0 in probability, which implies 4.1. ∎

######  Proof of Lemma 4.2.

Under A4-A7, the time series $\\{X_{t}\\}$ is ergodic. Besides, under A2 and
A10, we can show $L(\vartheta)$ is uniformly finite for
$\vartheta\in\Theta_{1}$. With A13, by the uniform law of larger numbers for
the ergodic series, see Theorem 6 of Kirch and Kamgaing, (2012) for a
reference, we have:

$\sup_{\vartheta\in\Theta_{1}}|L_{T}(\vartheta)-L(\vartheta)|\overset{p}{\to}0.$
(82)

Under A13 and A14, we can easily show that:

$\inf_{|\vartheta-\theta_{1}|>\varepsilon}L(\vartheta)>L(\theta_{1}),\text{for}~{}\forall\epsilon,$
(83)

Eq. 83 implies that given $\forall\varepsilon>0$, $\exists C>0$ such that
$|\vartheta-\theta_{1}|>\epsilon\Rightarrow L(\vartheta)-L(\theta_{1})\geq
C>0$, thus we have:

$\begin{split}\mathbb{P}(|\widehat{\theta}_{1}-\theta_{1}|>\epsilon)&\leq\mathbb{P}(L(\widehat{\theta}_{1})-L(\theta_{1})\geq
C)\\\
&=\mathbb{P}(L(\widehat{\theta}_{1})-L_{T}(\widehat{\theta}_{1})+L_{T}(\widehat{\theta}_{1})-L(\theta_{1})\geq
C)\\\
&\leq\mathbb{P}(L(\widehat{\theta}_{1})-L_{T}(\widehat{\theta}_{1})+L_{T}(\theta_{1})-L(\theta_{1})\geq
C)\\\
&\leq\mathbb{P}(2\sup_{\vartheta\in\Theta_{1}}|L_{T}(\vartheta)-L(\vartheta)|>C)\overset{}{\to}0.\end{split}$
(84)

The last inequality of Eq. 84 is guaranteed by Eq. 82. ∎

######  Proof of Lemma 4.3.

The proof of the consistency of $\widehat{\theta}_{1}$ to $\theta_{1}$ is the
same as the proof of 4.2. Similar with the proof of 4.2, by the ergodic
property of the series $\\{X_{t}\\}$, we know:

$\sup_{\vartheta\in\Theta_{2}}\left|\frac{1}{T}\sum_{t=1}^{T}\left(\frac{X_{t}-\phi(X_{t-1},\theta_{1})}{h(X_{t-1},\vartheta)}\right)^{2}-\mathbb{E}\left(\frac{X_{t}-\phi(X_{t-1},\theta_{1})}{h(X_{t-1},\vartheta)}\right)^{2}\right|\overset{p}{\to}0.$
(85)

Actually, after the first step estimation, we have got the
$\widehat{\theta}_{1}$ which is consistent to the true parameter $\theta_{1}$.
Thus, we need to find the below convergence relationship:

$\sup_{\vartheta\in\Theta_{2}}\left|\frac{1}{T}\sum_{t=1}^{T}\left(\frac{X_{t}-\phi(X_{t-1},\widehat{\theta}_{1})}{h(X_{t-1},\vartheta)}\right)^{2}-\mathbb{E}\left(\frac{X_{t}-\phi(X_{t-1},\theta_{1})}{h(X_{t-1},\vartheta)}\right)^{2}\right|\overset{p}{\to}0.$
(86)

Since $\widehat{\theta}_{1}$ converges to $\theta_{1}$ in probability, it is
easily to find:

$\sup_{\vartheta\in\Theta_{2}}\left|\frac{1}{T}\sum_{t=1}^{T}\left(\frac{X_{t}-\phi(X_{t-1},\widehat{\theta}_{1})}{h(X_{t-1},\vartheta)}\right)^{2}-\frac{1}{T}\sum_{t=1}^{T}\left(\frac{X_{t}-\phi(X_{t-1},\theta_{1})}{h(X_{t-1},\vartheta)}\right)^{2}\right|\overset{p}{\to}0,$
(87)

which implies Eq. 86 for a compact set $\Theta_{2}$ in conjunction with Eq.
85. Then, since $|~{}\cdot~{}-a|$ is a uniformly continuous function for a
constant $a$, by applying the uniform continuous mapping theorem which is the
Theorem 1 of Kasy, (2019) on Eq. 86, we can get:

$\sup_{\vartheta\in\Theta_{2}}|K_{T}(\vartheta,\widehat{\theta}_{1})-K(\vartheta,\theta_{1})|\overset{p}{\to}0,$
(88)

where $K(\vartheta,\theta_{1})$ and $K_{T}(\vartheta,\widehat{\theta}_{1})$
are:

$K(\vartheta,\theta_{1})=\left|\mathbb{E}\left(\frac{X_{t}-\phi(X_{t-1},\theta_{1})}{h(X_{t-1},\vartheta)}\right)^{2}-1\right|~{};~{}K_{T}(\vartheta,\widehat{\theta}_{1})=\left|\frac{1}{T}\sum_{t=1}^{T}\left(\frac{X_{t}-\phi(X_{t-1},\widehat{\theta}_{1})}{h(X_{t-1},\vartheta)}\right)^{2}-1\right|.$
(89)

Then, we can repeat the procedure in the proof of 4.2 to show:

$\mathbb{P}(|\widehat{\theta}_{2}-\theta_{2}|>\varepsilon)\overset{}{\to}0,~{}\forall\varepsilon>0.$
(90)

∎

######  Proof of Theorem 4.2.

Throughout this proof, we focus on a sequence of sets
$\Omega_{T+p}\subseteq\mathbb{R}^{T+p}$, such that
$\mathbb{P}\left((X_{-p+1},\ldots,X_{T})\in\Omega_{T+p}\right)=o(1)$.We first
explain the truth that the bootstrap series is ergodic for
$(X_{-p+1},\ldots,X_{T})\in\Omega_{T+p}$. Particularly, we want to check
analogous A4-A7 in the bootstrap world. With the consistency property of
parameter estimators, i.e., A9 and the continuous density function of
residuals after the convolutional manipulation, we can easily find A4-A7 are
satisfied in the bootstrap world; see Theorem 2 of Franke et al., (2004) for
more discussions.

To show the closeness of the stationary distribution of bootstrap and real
series. Based on the proof of Theorem 3 of Franke et al., 2002b  and Theorem 3
of Franke et al., (2004), we can get:

$\sup_{B}|\Pi(B)-\Pi^{*}(B)|=o(1),$ (91)

which holds for all measurable sets $B$, where $\Pi(B)$ and $\Pi^{*}(B)$
represent stationary distribution for real series and bootstrap series,
respectively. This implies, in conjunction with Eq. 33 and the condition that
$\mathbb{P}(X_{t}\notin\gamma_{T})=o(1)$, we can get:

$\begin{split}&\int_{\mathbb{R}^{2}}(x_{1}-\phi(x_{0},\Xi_{\theta}))^{2}\hat{f}_{\epsilon}(x_{1}-\phi(x_{0},\hat{\theta}))\pi^{*}(x_{0})dx_{1}dx_{0}\\\
&=\int_{\gamma_{T}^{2}}(x_{1}-\phi(x_{0},\Xi_{\theta}))^{2}\hat{f}_{\epsilon}(x_{1}-\phi(x_{0},\hat{\theta}))\pi^{*}(x_{0})dx_{1}dx_{0}+o(1)\\\
&\overset{}{\to}\mathbb{E}[X_{1}-\phi(X_{0},\Xi_{\theta})]^{2},\end{split}$
(92)

which implies Eq. 34 directly.

∎

######  Proof of Theorem 4.3.

To simplify notation, we consider the case $p=q=1$. Higher order cases can be
handled similarly. We first build the estimation inference for
$\widehat{\theta}_{1}$. Starting from A19, we have:

$op(T^{-1/2})=\nabla L_{T}(\widehat{\theta}_{1})=\nabla
L_{T}(\theta_{1})+\nabla^{2}L_{T}(\widetilde{\theta}_{1})(\widehat{\theta}_{1}-\theta_{1}),$
(93)

where $\widetilde{\theta}_{1}$ is between $\widehat{\theta}_{1}$ and
$\theta_{1}$; hence $\widetilde{\theta}_{1}$ also converges to $\theta_{1}$ in
probability. First, we consider $\nabla^{2}L_{T}(\widetilde{\theta}_{1})$,
which has a form as below:

$\nabla^{2}L_{T}(\widetilde{\theta}_{1})=\frac{1}{T}\nabla^{2}\sum_{t=1}^{T}\left(X_{t}-\phi(X_{t-1},\widetilde{\theta}_{1})\right)^{2}=\frac{1}{T}\sum_{t=1}^{T}\nabla^{2}q_{t}(\widetilde{\theta}_{1}).$
(94)

Under A18, since it is easy to check
$\mathbb{E}\sup_{\vartheta\in\Theta_{1}^{0}}$ $||\nabla^{2}q_{1}(\vartheta)$
$||<\infty$, by combining the uniform law of larger numbers for the ergodic
series, dominated convergence theorem, the consistency between
$\widetilde{\theta}_{1}$ and $\theta_{1}$ and the continuity of $L_{T}(\cdot)$
w.r.t. $\theta_{1}$, we can get:

$\nabla^{2}L_{T}(\widetilde{\theta}_{1})\overset{p}{\to}B_{1},$ (95)

where
$B_{1}=2\cdot\mathbb{E}\left(\nabla\phi(X_{0},\theta_{1})(\nabla\phi(X_{0},\theta_{1}))^{\top}\right)$.
Thus, we can multiply both side of Eq. 93 by $\sqrt{T}$ to get:

$op(1)=\sqrt{T}\nabla
L_{T}(\theta_{1})+(B_{1}+op(1))\sqrt{T}(\widehat{\theta}_{1}-\theta_{1}),$
(96)

which further implies that:

$\sqrt{T}(\widehat{\theta}_{1}-\theta_{1})=-B_{1}^{-1}\sqrt{T}\nabla
L_{T}(\theta_{1})+op(1),$ (97)

where $\sqrt{T}\nabla L_{T}(\theta_{1})$ has a concrete form as below:

$-\sqrt{T}\nabla
L_{T}(\theta_{1})=\frac{2}{\sqrt{T}}\sum_{t=1}^{T}(X_{t}-\phi(X_{t-1},\theta_{1}))\nabla\phi(X_{t-1},\theta_{1}).$
(98)

We need the CLT for strongly mixing processes to show the asymptotic normality
of Eq. 98. Based on Theorem 1.7 of Bosq, (2012) for univariate case, we can
show the normality with the Cramér–Wold theorem. Also, Kirch and Kamgaing,
(2012) applied the strong invariance principle from Kuelbs and Philipp, (1980)
to show the asymptotic distribution of multivariate cases directly. Since we
need to analyze the asymptotic distribution of parameters in the bootstrap
world, we take the first approach which is more clear. We consider:

$-\sqrt{T}\bm{a}^{\top}\nabla
L_{T}(\theta_{1})=\frac{2}{\sqrt{T}}\sum_{t=1}^{T}(X_{t}-\phi(X_{t-1},\theta_{1}))\bm{a}^{\top}\nabla\phi(X_{t-1},\theta_{1}),$
(99)

where $\bm{a}^{\top}$ is any real vector that has the same dimension as
$\theta_{1}$. Since we assume we can correctly specify the model, we have:

$-\sqrt{T}\bm{a}^{\top}\nabla
L_{T}(\theta_{1})\overset{d}{\to}N(0,\tau_{1}^{2}),$ (100)

where $\tau_{1}^{2}$ has the form as follows:

$\begin{split}\tau_{1}^{2}&=\sum_{i=-\infty}^{\infty}\text{Cov}\left(2\cdot(X_{1}-\phi(X_{0},\theta_{1}))\bm{a}^{\top}\nabla\phi(X_{0},\theta_{1}),2\cdot(X_{i+1}-\phi(X_{i},\theta_{1}))\bm{a}^{\top}\nabla\phi(X_{i},\theta_{1})\right)\\\
&=4\cdot\mathbb{E}\left((\bm{a}^{\top}\nabla\phi(X_{0},\theta_{1}))(X_{1}-\phi(X_{0},\theta_{1}))^{2}(\nabla\phi(X_{0},\theta_{1})^{\top}\bm{a})\right)\\\
&+2\sum_{i=1}^{\infty}\text{Cov}\left(2\cdot(X_{1}-\phi(X_{0},\theta_{1}))\bm{a}^{\top}\nabla\phi(X_{0},\theta_{1}),2\cdot(X_{i+1}-\phi(X_{i},\theta_{1}))\bm{a}^{\top}\nabla\phi(X_{i},\theta_{1})\right)\\\
&=4\cdot\mathbb{E}(\bm{a}^{\top}\sigma(X_{0},\theta_{2})\nabla\phi(X_{0},\theta_{1})\nabla\phi(X_{0},\theta_{1})^{\top}\sigma(X_{0},\theta_{2})\bm{a})\\\
&+2\sum_{i=1}^{\infty}\text{Cov}\left(2\cdot\epsilon_{1}\bm{a}^{\top}\nabla\phi(X_{0},\theta_{1}),2\cdot\epsilon_{i+1}\bm{a}^{\top}\nabla\phi(X_{i},\theta_{1})\right)\\\
&=4\cdot\mathbb{E}(\bm{a}^{\top}\sigma(X_{0},\theta_{2})R_{1}\sigma(X_{0},\theta_{2})\bm{a}),\end{split}$
(101)

where $R_{1}=\nabla\phi(X_{0},\theta_{1})\nabla\phi(X_{0},\theta_{1})^{\top}$.
Thus, applying Cramér–Wold theorem, we have:

$-\sqrt{T}\nabla L_{T}(\theta_{1})\overset{d}{\to}N(0,\Omega_{1}),$ (102)

where $\Omega_{1}$ is
$4\cdot\mathbb{E}(\sigma(X_{0},\theta_{2})R_{1}\sigma(X_{0},\theta_{2}))$
Thus, we can conclude from Eq. 97 that:

$\sqrt{T}(\widehat{\theta}_{1}-\theta_{1})\overset{d}{\to}N(0,B_{1}^{-1}\Omega_{1}B_{1}^{-1}).$
(103)

Compared to the result in Kirch and Kamgaing, (2012), they got another form of
$\Omega_{1}$:

$\Omega_{1}=\lim_{T\to\infty}\frac{1}{T}\mathbb{E}\left[\left(\nabla\sum_{t=1}^{T}\left(X_{t}-\phi(X_{t-1},\theta_{1})\right)^{2}\right)\left(\nabla\sum_{t=1}^{T}\left(X_{t}-\phi(X_{t-1},\theta_{1})\right)^{2}\right)^{T}\right].$
(104)

With a correctly specified model, it is equivalent to our form.

We can also analyze the distribution of parameter estimation
$\widehat{\theta}_{2}$. By A19, we first have:

$op(T^{-1/2})=\nabla
K_{T}(\theta_{2},\widehat{\theta}_{1})+\nabla^{2}K_{T}(\widetilde{\theta}_{2},\widehat{\theta}_{1})(\widehat{\theta}_{2}-\theta_{2}).$
(105)

For simplifying the notation, we write:

$\begin{split}K_{T}(\widetilde{\theta}_{2},\widehat{\theta}_{1})&=\left(\frac{1}{T}\sum_{t=1}^{T}\left(\frac{X_{t}-\phi(X_{t-1},\widehat{\theta}_{1})}{h(X_{t-1},\widetilde{\theta}_{2})}\right)^{2}-1\right)^{2}\\\
&=\left(\frac{1}{T}\sum_{t=1}^{T}g(X_{t},X_{t-1},\widetilde{\theta}_{2},\widehat{\theta}_{1})-1\right)^{2}.\\\
\end{split}$ (106)

We can find $\nabla K_{T}(\theta_{2},\theta_{1})$ and
$\nabla^{2}K_{T}(\theta_{2},\theta_{1})$ w.r.t $\theta_{2}$ have a following
forms respectively:

$\begin{split}&\nabla
K_{T}(\theta_{2},\theta_{1})=2\cdot\left(\frac{1}{T}\sum_{t=1}^{T}g(X_{t},X_{t-1},\theta_{2},\theta_{1})-1\right)\cdot\left(\frac{1}{T}\sum_{t=1}^{T}\nabla
g(X_{t},X_{t-1},\theta_{2},\theta_{1})\right)\\\
&\nabla^{2}K_{T}(\theta_{2},\theta_{1})=2\cdot\left(\frac{1}{T}\sum_{t=1}^{T}\nabla
g(X_{t},X_{t-1},\theta_{2},\theta_{1})\right)\cdot\left(\frac{1}{T}\sum_{t=1}^{T}\nabla
g(X_{t},X_{t-1},\theta_{2},\theta_{1})\right)^{\top}\\\
&+2\cdot\left(\frac{1}{T}\sum_{t=1}^{T}g(X_{t},X_{t-1},\theta_{2},\theta_{1})-1\right)\cdot\left(\frac{1}{T}\sum_{t=1}^{T}\nabla^{2}g(X_{t},X_{t-1},\theta_{2},\theta_{1})\right).\end{split}$
(107)

Similarly to analyze $\nabla^{2}L_{T}(\widetilde{\theta}_{1})$, under the
consistence relationship between $\widetilde{\theta}_{2}$ and $\theta_{2}$,
$\widehat{\theta}_{1}$ and $\theta_{1}$, we can get:

$\nabla^{2}K_{T}(\widetilde{\theta}_{1},\widehat{\theta}_{1})\overset{p}{\to}B_{2}.$
(108)

where $B_{2}=2\cdot\mathbb{E}(\nabla
g(X_{1},X_{0},\theta_{2},\theta_{1}))\cdot\mathbb{E}(\nabla
g(X_{1},X_{0},\theta_{2},\theta_{1})^{\top})$. Looking back to Eq. 105, it is
left to analyze $\sqrt{T}\nabla K_{T}(\theta_{2},\widehat{\theta}_{1})$.
Through Slutsky’s theorem, it is asymptotically equivalent to analyze
$\sqrt{T}\left(\frac{2}{T}\sum_{t=1}^{T}[g(X_{t},X_{t-1},\theta_{2},\theta_{1})-1]\right)\cdot
B_{3}$; $B_{3}=\mathbb{E}(\nabla g(X_{1},X_{0},\theta_{2},\theta_{1}))$. Apply
the same technique as we analyzed the distribution of $\widehat{\theta}_{1}$,
we can get:

$\sqrt{T}\bm{a}^{\top}\nabla
K_{T}(\theta_{2},\widehat{\theta}_{1})\overset{d}{\to}N(0,\tau_{2}^{2}),$
(109)

where
$\tau_{2}^{2}=4\cdot\mathbb{E}(\bm{a}^{\top}B_{3}R_{2}B_{3}^{\top}\bm{a})$;
$R_{2}=(g(X_{1},X_{0},\theta_{2},\theta_{1})-1)^{2}$. Thus, we have:

$\sqrt{T}\nabla
K_{T}(\theta_{2},\widehat{\theta}_{1})\overset{d}{\to}N(0,\Omega_{2}),$ (110)

where $\Omega_{2}=4\cdot\mathbb{E}(B_{3}R_{2}B_{3}^{\top})$. This further
implies that:

$\sqrt{T}(\widehat{\theta}_{2}-\theta_{2})\overset{d}{\to}N(0,B_{2}^{-1}\Omega_{2}B_{2}^{-1}).$
(111)

∎

## References

  * An and Huang, (1996) An, H. and Huang, F. (1996).  The geometrical ergodicity of nonlinear autoregressive models.  Statistica Sinica, 6(4):943–956.
  * Berg et al., (2012) Berg, A., McMurry, T., and Politis, D. N. (2012).  Testing time series linearity: traditional and bootstrap methods.  Handbook of Statistics, 30:27–42.
  * Boldin, (1983) Boldin, M. (1983).  Estimation of the distribution of noise in an autoregression scheme.  Theory of Probability & Its Applications, 27(4):866–871.
  * Bosq, (2012) Bosq, D. (2012).  Nonparametric statistics for stochastic processes: estimation and prediction, volume 110.  Springer Science & Business Media.
  * Box et al., (2015) Box, G. E., Jenkins, G. M., Reinsel, G. C., and Ljung, G. M. (2015).  Time Series Analysis: Forecasting and Control.  John Wiley & Sons.
  * Bradley, (2005) Bradley, R. C. (2005).  Basic properties of strong mixing conditions. a survey and some open questions.  Probability surveys, 2:107–144.
  * Chen and Politis, (2019) Chen, J. and Politis, D. N. (2019).  Optimal multi-step-ahead prediction of arch/garch models and novas transformation.  Econometrics, 7(3):1–34.
  * Chen et al., (2004) Chen, R., Yang, L., and Hafner, C. (2004).  Nonparametric multistep-ahead prediction in time series analysis.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 66(3):669–686.
  * De Gooijer and Kumar, (1992) De Gooijer, J. G. and Kumar, K. (1992).  Some recent developments in non-linear time series modelling, testing, and forecasting.  International Journal of Forecasting, 8(2):135–156.
  * Doob, (1953) Doob, J. L. (1953).  Stochastic processes.  Wiley New York.
  * Efron, (1979) Efron, B. (1979).  Bootstrap methods: Another look at the jackknife.  The Annals of Statistics, 7(1):1–26.
  * Eğrioğlu and Fildes, (2020) Eğrioğlu, E. and Fildes, R. (2020).  A new bootstrapped hybrid artificial neural network approach for time series forecasting.  Computational Economics, pages 1–29.
  * (13) Franke, J., Kreiss, J.-P., and Mammen, E. (2002a).  Bootstrap of kernel smoothing in nonlinear time series.  Bernoulli, 8(1):1–37.
  * (14) Franke, J., Kreiss, J.-P., Mammen, E., and Neumann, M. H. (2002b).  Properties of the nonparametric autoregressive bootstrap.  Journal of Time Series Analysis, 23(5):555–585.
  * Franke and Neumann, (2000) Franke, J. and Neumann, M. H. (2000).  Bootstrapping neural networks.  Neural Computation, 12(8):1929–1949.
  * Franke et al., (2004) Franke, J., Neumann, M. H., and Stockis, J.-P. (2004).  Bootstrapping nonparametric estimators of the volatility function.  Journal of Econometrics, 118(1-2):189–218.
  * Guo et al., (1999) Guo, M., Bai, Z., and An, H. Z. (1999).  Multi-step prediction for nonlinear autoregressive models based on empirical distributions.  Statistica Sinica, 9(2):559–570.
  * Jones and Cox, (1978) Jones, D. A. and Cox, D. R. (1978).  Nonlinear autoregressive processes.  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 360(1700):71–95.
  * Kasy, (2019) Kasy, M. (2019).  Uniformity and the delta method.  Journal of Econometric Methods, 8(1):1–19.
  * Kirch and Kamgaing, (2012) Kirch, C. and Kamgaing, J. T. (2012).  Testing for parameter stability in nonlinear autoregressive models.  Journal of Time Series Analysis, 33(3):365–385.
  * Kreiss and Paparoditis, (2023) Kreiss, J.-P. and Paparoditis, E. (2023).  Bootstrap for Time Series: Theory and Methods.  Springer-Verlag, New York.
  * Kuelbs and Philipp, (1980) Kuelbs, J. and Philipp, W. (1980).  Almost sure invariance principles for partial sums of mixing b-valued random variables.  The Annals of Probability, pages 1003–1036.
  * Lee and Billings, (2003) Lee, K. and Billings, S. (2003).  A new direct approach of computing multi-step ahead predictions for non-linear models.  International Journal of Control, 76(8):810–822.
  * Majerek et al., (2005) Majerek, D., Nowak, W., and Zieba, W. (2005).  Conditional strong law of large number.  Int. J. Pure Appl. Math, 20(2):143–156.
  * Manzan and Zerom, (2008) Manzan, S. and Zerom, D. (2008).  A bootstrap-based non-parametric forecast density.  International Journal of Forecasting, 24(3):535–550.
  * Min and Hongzhi, (1999) Min, C. and Hongzhi, A. (1999).  The probabilistic properties of the nonlinear autoregressive model with conditional heteroskedasticity.  Acta Mathematicae Applicatae Sinica, 15(1):9–17.
  * Pan and Politis, (2016) Pan, L. and Politis, D. N. (2016).  Bootstrap prediction intervals for linear, nonlinear and nonparametric autoregressions.  Journal of Statistical Planning and Inference, 177:1–27.
  * Paparoditis and Politis, (1998) Paparoditis, E. and Politis, D. N. (1998).  The backward local bootstrap for conditional predictive inference in nonlinear time series.  In Proceedings of the 4th Hellenic-European conference on computer mathematics and its applications (HERCMA’98). Lea Publishing, Athens, pages 467–470.
  * Paparoditis and Politis, (2002) Paparoditis, E. and Politis, D. N. (2002).  The local bootstrap for markov processes.  Journal of Statistical Planning and Inference, 108(1-2):301–328.
  * Pemberton, (1987) Pemberton, J. (1987).  Exact least squares multi-step prediction from nonlinear autoregressive models.  Journal of Time Series Analysis, 8(4):443–448.
  * Politis, (2009) Politis, D. N. (2009).  Financial time series.  Wiley Interdisciplinary Reviews: Computational Statistics, 1(2):157–166.
  * Politis, (2013) Politis, D. N. (2013).  Model-free model-fitting and predictive distributions.  Test, 22(2):183–221.
  * Politis, (2015) Politis, D. N. (2015).  Model-free prediction in regression.  In Model-Free Prediction and Regression, pages 57–80. Springer.
  * Politis et al., (1997) Politis, D. N., Romano, J. P., and Wolf, M. (1997).  Subsampling for heteroskedastic time series.  Journal of Econometrics, 81(2):281–317.
  * Tjøstheim, (1990) Tjøstheim, D. (1990).  Non-linear time series and markov chains.  Advances in Applied Probability, 22(3):587–611.
  * Tjøstheim, (1994) Tjøstheim, D. (1994).  Non-linear time series: a selective review.  Scandinavian Journal of Statistics, 21(2):97–130.
  * Tong, (1990) Tong, H. (1990).  Non-linear Time Series: A Dynamical System Approach.  Oxford university press.
  * Trucíos et al., (2017) Trucíos, C., Hotta, L. K., and Ruiz, E. (2017).  Robust bootstrap forecast densities for garch returns and volatilities.  Journal of Statistical Computation and Simulation, 87(16):3152–3174.
  * Tweedie, (1975) Tweedie, R. L. (1975).  Sufficient conditions for ergodicity and recurrence of markov chains on a general state space.  Stochastic Processes and Their Applications, 3(4):385–403.
  * Wang and Politis, (2021) Wang, Y. and Politis, D. N. (2021).  Model-free bootstrap and conformal prediction in regression: Conditionality, conjecture testing, and pertinent prediction intervals.  arXiv preprint arXiv:2109.12156.
  * Zhang et al., (1998) Zhang, G., Patuwo, B. E., and Hu, M. Y. (1998).  Forecasting with artificial neural networks:: The state of the art.  International Journal of Forecasting, 14(1):35–62.


# The Heat is On: Thermal Facial Landmark Tracking

James Baker
Department of Computer Science
University of Houston

(25/4/15)

###### Abstract

Facial landmark tracking for thermal images requires tracking certain
important regions of subjects’ faces, using images from thermal images, which
omit lighting and shading, but show the temperatures of their subjects. The
fluctuations of heat in particular places reflect physiological changes like
bloodflow and perspiration, which can be used to remotely gauge things like
anxiety and excitement. Past work in this domain has been limited to only a
very limited set of architectures and techniques. This work goes further by
trying a comprehensive suit of various models with different components, such
as residual connections, channel and feature-wise attention, as well as the
practice of ensembling components of the network to work in parallel. The best
model integrated convolutional and residual layers followed by a channel-wise
self-attention layer, requiring less than 100K parameters.

##  1 Introduction

Detecting physiological changes in human faces, like bloodflow or
perspiration, can measure stress [42, 43], empathy [12], and even deceit [46].
Thermal images are ideal for conveying physiological changes. They are non-
invasive, unlike sensors, and are sensitive to the heat changes caused by
physiological changes, unlike normal images. However, using thermal images
requires consistently measuring the same region of interest on a human
subject, which is complicated by the fact that humans move, tilt and rotate
their bodies. It is imperative to use a robust, reliable framework for
landmark tracking for these tasks.

Previous work on thermal image facial landmark tracking did not implement
techniques such as non-locality, residual layers and ensembling or ”wisdom-of-
crowds”. This paper rectifies that by using hyperparameter optimization to
construct models with residual and attentional elements. Multiple, parallel
components were added to the top models, in order to test if ensembling would
improve performance.

###  1.1 Contributions of This Work

After experimentation, the key insight from this work was that Convolutional
Neural Networks that used attention and residual components were not only
accurate, but also very lightweight in their number of parameters, which will
assist in their use in smaller, mobile devices. In fact, increasing their
complexity by adding parallel components made the best models overfit.

##  2 Related Work

In affective computing, researchers need to measure physiology of their
subjects. Unfortunately, this often requires attaching sensors to people,
which is burdensome to the subjects and does not scale economically. Using
cameras solves this, but non-thermal images do not provide any information on
temperature. Dusty, foggy or dark environments may make rob normal images of
any information, and normal cameras will barely pick up anything, whereas
thermal imaging is still robust under bad conditions.

The specific task we are interested in is landmark tracking. Given a face, we
want to reliably find specific objects or regions (landmarks), like the eyes,
nose or lips. There have been many datasets made for this task for normal
images [44, 40]. What makes thermal facial images unique is skin temperature
changes in response to things like bloodflow or breathing, which means that
regions of interest may be prone to suddenly blend in or stand out or change
color over time. Thermal images are also difficult as they tend to lack the
textural information in normal images [37]. While they offer more information
in some ways, they are challenging in others, thus using thermal images is a
”higher risk, higher reward” version of working with normal images.

###  2.1 Statistical Methods

An adjacent task to landmark detection is region-of-interest tracking. Using
videos of faces, [60] used a particle filter tracker in order to follow
regions of interest. Similarly, [30] used a particle filter and an object
detector that follows the position of the eyes to build features for
classification of regions using a random fern algorithm. Other landmark
detection has been in the domain of normal images. Many methods use a
regression framework, which can be formulated like so. Given a set of p
coordinates of landmarks at time t $S_{t}$:

$S_{t}=\\{(x_{0},y_{0}),(x_{1},y_{1})...(x_{p},y_{p})$ (1)

Given regressor $r_{t}$ and image $I$, we can model:

$S_{t+1}=S_{t}+r_{t}(I,S_{t})$ (2)

Given predicted $\hat{S_{t}}$ and actual $S_{t}^{*}$, we choose $r_{t}$ as the
solution to the optimization problem

$\operatorname*{arg\,min}_{r_{t}}=\sum^{N}||S_{t}^{*}-\hat{S_{t}}-r_{t}||_{2}$
(3)

In [41], the authors used a continuous regression method, where the input
space was approximated using Taylor Expansions, which proved to be robust to
translation, tilts and other shifts, as well as computationally efficient to
solve for. [19] also used a cascaded regression to find landmarks. Then, using
the location of landmarks, they found similar faces in the dataset and used
those to fit a regression for each face. [27] used an approach where at each
step, the sampling region for each landmark was updated, and a Markov Random
Field ensured that the new sampling regions were consistently in a face shape.
[54] used a random regression forest, followed by a cascade of sieves to
filter out votes that are too inconsistent or distant from the hypothesis for
the predicted landmark.

###  2.2 Deep Learning Methods

[21] used the U-Net, introduced by [39], which is a deep convolutional neural
network (CNN) that first downsamples and then upsamples the input images, with
a fully connected head. The CNN was trained, and then the fully connected
layer. They used one U-Net to find a square enclosure of the face, and passed
the enclosed region to another U-Net that predicted the coordinates of the
landmarks. [10] expanded on the sequential U-Net model by replacing the final
fully connected layer with two parallel layers, one for landmark detection and
one for emotion classification. [37] compared three different CNN models for
landmark detection: a Multitask Cascaded Convolutional Network [59], a Deep
Alignment Network [22], and a Multiclass Patch Based Classifier, a CNN which
for every 60 x 60 patch of the input image classifies it as one of six
regions.

##  3 Proposed Approach

###  3.1 Singular Models

All models consisted of 4 basic components. The first was a “root”, which was
the same in all models, that consisted of 2 convolutional layers for
downsampling the input while also increasing the number of feature channels,
translating the images $\mathbb{R}^{480\times 640\times
3}\rightarrow\mathbb{R}^{120\times 160\times 64}$. This was also done because
models that were initially convolutional but then implemented attentional
layers later on had shown to perform better [38] and are used in other
successful experiments [51]. Following the root was a ”stem”, then an optional
”branch”. Following the stem, or branch if the model had one, was a ”head”,
which flattened the outputs of the preceding layer,applied a fully connected
layer with 2N nodes and dropout, where N is the amount of points to be
predicted, and then reshaped the outputs from
$\mathbb{R}^{2N}\rightarrow\mathbb{R}^{2\times N}$.

A stem was made by repeating one of five different layers. More detail
describing :

  * •

Convolutional layers (referrred to as Conv Stem)

  * •

ResNeXt layers (referred to as ResNeXt stem)

  * •

Alternating Convolutional layers and Bahdanau feature-wise attention layers
(referred to as Alternating Conv-Bahdanau Stem)

  * •

Alternating Convolutional layers and Luong feature-wise attention layers
(referred to as Alternating Conv-Luong Stem)

  * •

Alternating Convolutional layers and ResNeXt layers (referred to as
Alternating Conv-ResNeXt stem)

A branch was made by repeating one of four different layers

  * •

Nothing; the stem was directly connected to the head (referred to as No
Branch)

  * •

Luong feature-wise attention layers (referred to as Luong branch)

  * •

Bahdanau feature-wise attention layers (referred to as Bahdanau branch)

  * •

A single Patch Encoder layer (not repeated), and then Transformer spatial-wise
attention layers (referred to as Vision Transformer branch)

Figure 1: Singular Model (left), and Ensemble Model (right), in testing mode,
with Image 8282 from subject 64

####  3.1.1 Optimization

Due to the massive search space of hyperparameters (depth, kernel size, etc.)
for these models, we used the Optuna library [1] to efficiently and
automatically run trials to find the optimal hyperparameters. For searching
for which hyperparameters to test for each trial, we used the Tree-Structured
Parzen Estimator (TPE), given it had performed well experimentally [32]. The
TPE algorithm is a greedy algorithm that samples hyperparameters for new
trials from a distribution estimated using the hyperparameters from past
previous trials [7]. For pruning unpromising trials, we used the Asynchronous
Successive Halving (ASHA), which also performed well in experiments [24]. ASHA
works by asynchronously promoting the best trials based on some performance
threshold, and then canceling the rest.

For each stem, for each branch, we ran 40 trials that ran for 20 epochs each,
and chose the model configurations that returned the lowest loss to be the
optimal, leading to 20 different models. The loss function that was to be
minimized was wing loss [13], which penalizes small and medium weight errors,
and had been shown to perform better for facial landmark localisation tasks.
The size (in quantity of parameters) and wing loss of the final best models is
given in Tab. 1. The singular models have been indexed with arbitrary letters
of the alphabet, and in later tables, models will be referred to using their
alphabetical letter.

Model | Params | Wing Loss
---|---|---
C-1: Conv Stem + No Branch | 383100 | 2.0686
C-2: Conv Stem + Luong Branch | 139404 | 4.1657
C-3: Conv Stem + Bahdanau Branch | 533948 | 1.4898
C-4: Conv Stem + Vision Transformer Branch | 215868 | 1.7841
R-1: ResNeXt Stem + No Branch | 287792 | 3.3708
R-2: ResNeXt Stem + Luong Branch | 294844 | 3.5428
R-3: ResNext Stem + Bahdanau Branch | 2275528 | 1.6349
R-4: ResNeXt Stem + Vision Transformer Branch | 159420 | 2.4904
L-1: Alternating Conv-Luong Stem + No Branch | 618844 | 1.7184
L-2: Alternating Conv-Luong + Luong Branch | 303836 | 4.6675
L-3: Alternating Conv-Luong + Bahdanau Branch | 283772 | 1.5763
L-4: Alternating Conv-Luong + Vision Transformer Branch | 911292 | 2.306
B-1: Alternating Conv-Bahdanau Stem + No Branch | 592748 | 2.7187
B-2: Alternating Conv-Bahdanau + Luong Branch | 262028 | 5.092
B-3: Alternating Conv-Bahdanau + Bahdanau Branch | 599276 | 2.5331
B-4: Alternating Conv-Bahdanau + Vision Transformer Branch | 747948 | 2.5257
A-1: Alternating Conv-ResNeXt Stem + No Branch | 762052 | 1.3416
A-2: Alternating Conv-ResNeXt + Luong Branch | 91280 | 2.5564
A-3: Alternating Conv-ResNeXt + Bahdanau Branch | 541380 | 1.1292
A-4: Alternating Conv-ResNeXt + Vision Transformer Branch | 2070140 | 1.3802
Table 1: Optimal Models from Optuna

###  3.2 Ensemble Models

The success of ensemble learning, using multiple models together for the same
task, has proven successful in fields such as robotics [58], finance [53] and
medicine [3], as ensembling reduces variance [8] and evades local optima [14].
Specifically, this paper used stacked generalizations [49], where a final
algorithm is trained to use the outputs of the component models ”stacked”
together. In this case, this was implemented by concatenating the flattened
outputs of all the component models and applying a fully-connected head. Each
of the component models also shared a convolutional root, as using the same
root would add parameters, and each root would likely be supplying redundant
information.

##  4 Experiments

###  4.1 Implementation Details

The data consisted of roughly 2000 thermal images collected from videos with
73 different subjects, each of which had been annotated with the locations of
16 landmarks. Ten of the videos were collected as part of an IRB approved
study. The subjects were recruited through email and personal
communications—the subjects work/study at a public research university. They
were filmed in an office setting doing work like grant proposals and writing
papers. They were filmed visually using a Tau 640 long-wave infrared (LWIR)
camera (FLIR Systems, Wilsonville, OR. The camera features 50◦ mK thermal
resolution, 640 × 480 pixels spatial resolution, and an auto-focus mechanism.
The camera was located under the participant’s desktop screen, attached to a
Bescor MP- 101 Motorized Pan; Tilt Head (Bescor, Farmingdale, NY) to
facilitate face tracking. Thermal facial data were collected at a frame rate
of approximately 30 fps using a 35 mm lens. The other 63 subjects were
knowledge workers and were recorded using a Tau 640 long-wave infrared (LWIR)
camera (FLIR Systems, Wilsonville, OR), featuring a small size (44 × 44 × 30
mm) and adequate thermal (50 mK) and spatial resolution (640×512 pixels),
though the images were cropped to 640 × 480 pixels for this experiment. A LWIR
35mm lens f/1.2, controlled by a custom auto-focus mechanism was fitted on the
camera. The thermal camera is located under the participant’s computer screen,
attached to Bescor MP-101 Motorized Pan and Tilt Head (Bescor, Farmingdale,
NY) to facilitate face tracking. The 63 images were previously published [56].
As the upper lip region is the most useful for stress studies [42, 43], we
only trained the Model to find the 6 points surrounding that region. In order
to make the data more robust, each image and the relevant landmarks were
rotated by a random angle between 20 and 30 degrees around the center of the
drawing, to the left and the right, as shown in Fig. 2.

Figure 2: Original and Augmented Images, Annotated with locations of landmark
points

###  4.2 Results

####  4.2.1 Singular Results

A total of 20 different models singular models were tested. The exact
hyperparameters were found using the methods discussed in Sec. 3.1. Each of
these models was trained for 100 epochs with a batch size of 32. Tab. 2
reports the Accuracy, Wing Loss, Mean Absolute Error (MAE) and Mean Squared
Error (MSE) of each model when tested.

As Tab. 2 shows, the best models were those that alternated ResNeXt and
Convolutional layers.

| metrics
---|---
Model | Accuracy | Wing Loss | MAE | MSE
Model C-1 | 0.9389 | 1.7195 | 0.0292 | 0.0018
Model C-2 | 0.9182 | 2.2061 | 0.0376 | 0.0026
Model C-3 | 0.9482 | 1.4323 | 0.0244 | 0.0012
Model C-4 | 0.9517 | 1.3282 | 0.0226 | 0.0011
Model R-1 | 0.6848 | 21.9716 | 0.4085 | 0.1873
Model R-2 | 0.5723 | 21.9296 | 0.4054 | 0.1846
Model R-3 | 0.6188 | 22.1513 | 0.4096 | 0.1885
Model R-4 | 0.679 | 22.1508 | 0.4103 | 0.19
Model L-1 | 0.9541 | 1.2664 | 0.0209 | 0.001
Model L-2 | 0.8572 | 3.7362 | 0.0649 | 0.0075
Model L-3 | 0.9532 | 1.277 | 0.0214 | 0.001
Model L-4 | 0.9446 | 1.5254 | 0.0259 | 0.0014
Model B-1 | 0.9193 | 1.9792 | 0.0339 | 0.0024
Model B-2 | 0.9078 | 2.4484 | 0.0418 | 0.0033
Model B-3 | 0.9146 | 2.2086 | 0.0371 | 0.0028
Model B-4 | 0.9262 | 1.9015 | 0.0325 | 0.0022
Model A-1 | 0.966 | 0.9175 | 0.0153 | 0.0005
Model A-2 | 0.96 | 1.0997 | 0.0188 | 0.0007
Model A-3 | 0.9715 | 0.7628 | 0.0128 | 0.0003
Model A-4 | 0.9538 | 1.0896 | 0.0186 | 0.0007
Table 2: Singular Model Results

Model A-3 stands out as the best one. It achieves the lowest wing loss of
0.7628, while using only 91,280 parameters, while most state of the art models
use millions.

####  4.2.2 Ensemble Results

Given that the Alternating Conv-ResNeXt models (Models Q,R,S,T) performed
best, we then tried to improve them by using an ensemble of 3 parallel stem-
branch components connecting the root to the head. Results are shown in Tab.
3. Performance consistently decreased with the use of ensembling.

| metrics
---|---
Model | Accuracy | Wing Loss | MAE | MSE
Model A-1 | 0.9643 | 0.9921 | 0.0167 | 0.0008
Model A-2 | 0.9592 | 1.153 | 0.0197 | 0.0009
Model A-3 | 0.6812 | 21.9525 | 0.407 | 0.1864
Model A-4 | 0.5726 | 22.2677 | 0.4137 | 0.1914
Table 3: Ensemble Model Results

###  4.3 Ablation Study

####  4.3.1 Data Without Rotations

In order to test if rotating the data was actually useful, we then repeated
the experiments with the alternating residual and convolutional singular
models. Given that the data set was one third of the size without the rotated
images, we trained them for 300 epochs instead of the usual 100. Results are
shown in Tab. 5. There were significant decreases in performance, as
predicted.

| metrics
---|---
Model | Accuracy | Wing Loss | MAE | MSE
Model A-1 | 0.9474 | 1.3288 | 0.0225 | 0.0012’
Model A-2 | 0.9389 | 1.4316 | 0.0243 | 0.0014’
Model A-3 | 0.9599 | 1.0501 | 0.0176 | 0.0008
Model A-4 | 0.9523 | 1.5673 | 0.0263 | 0.0016
Table 4: Models Trained Without Rotations

####  4.3.2 Removing Layers

To test if the models were too complex, we removed the last convolutional and
ResNeXt layer. As Tab. 5 shows, wing loss increased in all cases but the Model
A-4 (alternating convolutional and ResNeXt stem and a Vision Transformer
branch), in which case performance stayed nearly the same.

| metrics
---|---
Model | Accuracy | Wing Loss | MAE | MSE
Model A-1 | 0.7341 | 18.0931 | 0.3301 | 0.1339
Model A-2 | 0.9553 | 1.1807 | 0.0192 | 0.0008
Model A-3 | 0.9654 | 0.9057 | 0.0152 | 0.0005
Model A-4 | 0.9632 | 1.0893 | 0.0185 | 0.0007
Table 5: Models with Layers Removed

###  4.4 Activation Maximization

Activation Maximization consists of searching for an input that produces the
maximum possible value for some particular output of the model [29]. This can
mean maximizing the output of a particular unit or even a whole layer in a
network. Activation Maximization is useful for visualizing what features
hidden layers are responsive to.

####  4.4.1 Gradient Ascent

In this case, we used gradient ascent [26] on a random noise image to maximize
the output of particular channels of output layers. The outputs were extremely
abstract and surreal. In Appendix A, Fig. 3 and Fig. 4 show samples of images
that were made to maximize different output channels in the second to last
layer, and last layer, respectively, of the alternating convolutional and
residual singular model.

##  5 Conclusions and Recommendations

In this paper, we demonstrated that models that had alternating convolutional
and residual components in early layers, with attentional components in later
models, were very effective in facial landmark tracking for thermal images.
Specifically, the best model (Model A-3) used alternating residual and
convolutional layers, and then a luong channel-wise attention layer, while
using less than 100,000 parameters. This is not surprising given that one of
the selling points of the ResNeXt layer [50] was its high performance relative
to its small number of parameters.

The failure of ensembling is likely due to the ability of models that are too
complicated learning to overfit the training data, resulting in high
generalization error [6]. Even standard forms of regularization employed here
like dropout and data augmentation may not sufficiently correct this [57].

This particular task of locating this particular region of the face is useful
for monitoring stress, which can affect whether we operate machinery [33],
drive automobiles [34, 15] safely, and even perform surgery [35]. Ideally,
this work can be used for further research on how humans perform under stress
and monitoring stress as we perform tasks to gauge when we are prone to make
potentially dangerous mistakes. Sleep studies also require unobtrusive facial
measurements, and this work may also be useful for those as well.

###  5.1 Limitations

The experiments in this paper were done on only one dataset. Our model may not
robust to facial attributes not present in this dataset (some people only have
one eye, but none of the subjects did). The images were also very high-
quality; the assumption is that this would be used for scenarios where humans
were very close to the cameras. For example, when operating an automobile, or
in a lab, it would be feasible to keep a heat camera a few feet from the
humans face. The subjects were also not filmed doing physically strenuous
tasks like exercising or manual labor. The increased blood flow and
perspiration may create very different heat patterns in subjects faces.

## References

  * [1] Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama.  Optuna: A next-generation hyperparameter optimization framework.  CoRR, abs/1907.10902, 2019.
  * [2] Jay Alammar.  The illustrated transformer.  2020\.
  * [3] Farman Ali, Shaker El-Sappagh, S.M. Riazul Islam, Daehan Kwak, Amjad Ali, Muhammad Imran, and Kyung-Sup Kwak.  A smart healthcare monitoring system for heart disease prediction based on ensemble deep learning and feature fusion.  Information Fusion, 63:208–222, 2020.
  * [4] Jimmy Ba, Volodymyr Mnih, and Koray Kavukcuoglu.  Multiple object recognition with visual attention, 2015.
  * [5] Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio.  Neural machine translation by jointly learning to align and translate, 2016.
  * [6] Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal.  Reconciling modern machine learning practice and the bias-variance trade-off, 2019.
  * [7] James Bergstra, R. Bardenet, Balázs Kégl, and Y. Bengio.  Algorithms for hyper-parameter optimization.  12 2011.
  * [8] Leo Breiman.  Bias, variance, and arcing classifiers.  Technical report, 1996.
  * [9] Long Chen, Hanwang Zhang, Jun Xiao, Liqiang Nie, Jian Shao, Wei Liu, and Tat-Seng Chua.  Sca-cnn: Spatial and channel-wise attention in convolutional networks for image captioning.  In 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 6298–6306, 2017.
  * [10] Wei-Ta Chu and Yu-Hui Liu.  Thermal facial landmark detection by deep multi-task learning.  In 2019 IEEE 21st International Workshop on Multimedia Signal Processing (MMSP), pages 1–6, 2019.
  * [11] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale.  CoRR, abs/2010.11929, 2020.
  * [12] Sjoerd J. Ebisch, Tiziana Aureli, Daniela Bafunno, Daniela Cardone, Gian Luca Romani, and Arcangelo Merla.  Mother and child in synchrony: Thermal facial imprints of autonomic contagion.  Biological Psychology, 89(1):123–129, 2012.
  * [13] Zhen-Hua Feng, Josef Kittler, Muhammad Awais, Patrik Huber, and Xiaojun Wu.  Wing loss for robust facial landmark localisation with convolutional neural networks.  CoRR, abs/1711.06753, 2017.
  * [14] Mudasir A. Ganaie, Minghui Hu, Mohammad Tanveer, and Ponnuthurai N. Suganthan.  Ensemble deep learning: A review.  CoRR, abs/2104.02395, 2021.
  * [15] Juan Pablo Gomez, Derya Akleman, Ergun Akleman, and Ioannis Pavlidis.  Causality effects of interventions and stressors on driving behaviors under typical conditions.  Mathematics, 6(8), 2018.
  * [16] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep Learning.  MIT Press, 2016.  http://www.deeplearningbook.org.
  * [17] Alex Graves, Greg Wayne, and Ivo Danihelka.  Neural turing machines.  CoRR, abs/1410.5401, 2014.
  * [18] Karol Gregor, Ivo Danihelka, Alex Graves, and Daan Wierstra.  DRAW: A recurrent neural network for image generation.  CoRR, abs/1502.04623, 2015.
  * [19] Yangyang Hao, Hengliang Zhu, Kai Wu, Xiao Lin, and Lizhuang Ma.  Salient-points-guided face alignment.  Multimedia Systems, pages 1–11, 2017.
  * [20] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition, 2015.
  * [21] Jin Keong, Xingbo Dong, Zhe Jin, Khawla Mallat, and J. Dugelay.  Multi-spectral facial landmark detection.  2020 IEEE International Workshop on Information Forensics and Security (WIFS), pages 1–6, 2020.
  * [22] Marek Kowalski, Jacek Naruniec, and Tomasz Trzcinski.  Deep alignment network: A convolutional neural network for robust face alignment.  CoRR, abs/1706.01789, 2017.
  * [23] A. Krizhevsky, Ilya Sutskever, and Geoffrey E. Hinton.  Imagenet classification with deep convolutional neural networks.  Communications of the ACM, 60:84 – 90, 2012.
  * [24] Liam Li, Kevin Jamieson, Afshin Rostamizadeh, Ekaterina Gonina, Moritz Hardt, Benjamin Recht, and Ameet Talwalkar.  A system for massively parallel hyperparameter tuning, 2020.
  * [25] Minh-Thang Luong, Hieu Pham, and Christopher D. Manning.  Effective approaches to attention-based neural machine translation.  CoRR, abs/1508.04025, 2015.
  * [26] Aravindh Mahendran and Andrea Vedaldi.  Understanding deep image representations by inverting them, 2014.
  * [27] Brais Martinez, Michel F. Valstar, Xavier Binefa, and Maja Pantic.  Local evidence aggregation for regression-based facial point detection.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(5):1149–1163, 2013.
  * [28] Volodymyr Mnih, Nicolas Heess, Alex Graves, and Koray Kavukcuoglu.  Recurrent models of visual attention.  CoRR, abs/1406.6247, 2014.
  * [29] Grégoire Montavon, Wojciech Samek, and Klaus-Robert Müller.  Methods for interpreting and understanding deep neural networks.  Digital Signal Processing, 73:1–15, 2018.
  * [30] Eslam Mostafa, Aly Farag, Ahmed Shalaby, Asem Ali, Travis Gault, and Ali Mahmoud.  Long term facial parts tracking in thermal imaging for uncooperative emotion recognition.  In 2013 IEEE Sixth International Conference on Biometrics: Theory, Applications and Systems (BTAS), pages 1–6, 2013.
  * [31] Chigozie Nwankpa, Winifred Ijomah, Anthony Gachagan, and Stephen Marshall.  Activation functions: Comparison of trends in practice and research for deep learning.  CoRR, abs/1811.03378, 2018.
  * [32] Skogby Steinholtz Olof.  A comparative study of black-box optimization algorithms for tuning of hyper-parameters in deep neural networks.  2018\.
  * [33] George Panagopoulos and Ioannis Pavlidis.  Forecasting markers of habitual driving behaviors associated with crash risk.  IEEE Transactions on Intelligent Transportation Systems, 21(2):841–851, 2020.
  * [34] Ioannis Pavlidis, Ashik Khatri, Pradeep Buddharaju, Michael Manser, Robert Wunderlich, Ergun Akleman, and Panagiotis Tsiamyrtzis.  Biofeedback arrests sympathetic and behavioral effects in distracted driving.  IEEE Transactions on Affective Computing, 12(2):453–465, 2021.
  * [35] I. Pavlidis, Dmitry Zavlin, Ashik Khatri, Amanveer Wesley, G. Panagopoulos, and A. Echo.  Absence of stressful conditions accelerates dexterous skill acquisition in surgery.  Scientific Reports, 9, 2019.
  * [36] George Philipp, D. Song, and J. Carbonell.  Gradients explode - deep networks are shallow - resnet explained.  ArXiv, abs/1712.05577, 2018.
  * [37] Domenick Poster, Shuowen Hu, Nasser Nasrabadi, and Benjamin Riggan.  An examination of deep-learning based landmark detection methods on thermal face imagery.  In 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), pages 980–987, 2019.
  * [38] Prajit Ramachandran, Niki Parmar, Ashish Vaswani, Irwan Bello, Anselm Levskaya, and Jonathon Shlens.  Stand-alone self-attention in vision models.  CoRR, abs/1906.05909, 2019.
  * [39] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  CoRR, abs/1505.04597, 2015.
  * [40] Christos Sagonas, Epameinondas Antonakos, Georgios Tzimiropoulos, Stefanos Zafeiriou, and Maja Pantic.  300 faces in-the-wild challenge.  Image Vision Comput., 47(C):3–18, Mar. 2016.
  * [41] Enrique Sanchez-Lozano, Georgios Tzimiropoulos, Brais Martinez, Fernando De la Torre, and Michel Valstar.  A functional regression approach to facial landmark tracking.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 40(9):2037–2050, Sep 2018.
  * [42] Dvijesh Shastri, Arcangelo Merla, Panagiotis Tsiamyrtzis, and Ioannis Pavlidis.  Imaging facial signs of neurophysiological responses.  IEEE Transactions on Biomedical Engineering, 56(2):477–484, 2009\.
  * [43] Dvijesh Shastri, Manos Papadakis, Panagiotis Tsiamyrtzis, Barbara Bass, and Ioannis Pavlidis.  Perinasal imaging of physiological stress and its affective potential.  IEEE Transactions on Affective Computing, 3(3):366–378, 2012.
  * [44] Jie Shen, S. Zafeiriou, Grigorios G. Chrysos, Jean Kossaifi, Georgios Tzimiropoulos, and M. Pantic.  The first facial landmark tracking in-the-wild challenge: Benchmark and results.  2015 IEEE International Conference on Computer Vision Workshop (ICCVW), pages 1003–1011, 2015.
  * [45] K. Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  CoRR, abs/1409.1556, 2015.
  * [46] Panagiotis Tsiamyrtzis, J. Dowdall, Dvijesh Shastri, Ioannis Pavlidis, Mark Frank, and P. Ekman.  Imaging facial physiology for the detection of deceit.  International Journal of Computer Vision, 71:197–214, 02 2007.
  * [47] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  CoRR, abs/1706.03762, 2017.
  * [48] Eric Weisstein.  Convolution.
  * [49] David H. Wolpert.  Stacked generalization.  Neural Networks, 5(2):241–259, 1992.
  * [50] Saining Xie, Ross B. Girshick, Piotr Dollár, Zhuowen Tu, and Kaiming He.  Aggregated residual transformations for deep neural networks.  2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 5987–5995, 2017.
  * [51] Kelvin Xu, Jimmy Ba, Ryan Kiros, Kyunghyun Cho, Aaron Courville, Ruslan Salakhutdinov, Richard Zemel, and Yoshua Bengio.  Show, attend and tell: Neural image caption generation with visual attention, 2016.
  * [52] Kelvin Xu, Jimmy Ba, Ryan Kiros, Kyunghyun Cho, Aaron C. Courville, Ruslan Salakhutdinov, Richard S. Zemel, and Yoshua Bengio.  Show, attend and tell: Neural image caption generation with visual attention.  CoRR, abs/1502.03044, 2015.
  * [53] Hongyang Yang, Xiao-Yang Liu, Shanli Zhong, and A. Walid.  Deep reinforcement learning for automated stock trading: An ensemble strategy.  2020\.
  * [54] Heng Yang and I. Patras.  Sieving regression forest votes for facial feature detection in the wild.  2013 IEEE International Conference on Computer Vision, pages 1936–1943, 2013.
  * [55] Minghao Yin, Zhuliang Yao, Yue Cao, Xiu Li, Zheng Zhang, Stephen Lin, and Han Hu.  Disentangled non-local neural networks, 2020.
  * [56] Shaila Zaman, Amanveer Wesley, Dennis Rodrigo Da Cunha Silva, Pradeep Buddharaju, Fatema Akbar, Ge Gao, Gloria Mark, Ricardo Gutierrez-Osuna, and Ioannis Pavlidis.  Stress and productivity patterns of interrupted, synergistic, and antagonistic office activities.  Scientific Data, 6(1):264, Nov 2019.
  * [57] Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals.  Understanding deep learning requires rethinking generalization, 2017.
  * [58] Junhao Zhang, Wei Zhang, Ran Song, Ling guo Ma, and Yibin Li.  Grasp for stacking via deep reinforcement learning.  2020 IEEE International Conference on Robotics and Automation (ICRA), pages 2543–2549, 2020.
  * [59] Kaipeng Zhang, Zhanpeng Zhang, Zhifeng Li, and Yu Qiao.  Joint face detection and alignment using multitask cascaded convolutional networks.  IEEE Signal Processing Letters, 23(10):1499–1503, 2016.
  * [60] Y. Zhou, P. Tsiamyrtzis, P. Lindner, I. Timofeyev, and Ioannis T. Pavlidis.  Spatiotemporal smoothing as a basis for facial tissue tracking in thermal imaging.  IEEE Transactions on Biomedical Engineering, 60:1280–1289, 2013\.

##  Appendix A Components

###  A.1 Convolution

Generally, a convolution of two functions is their product over a range [48].
In the continuous case, this can be expressed as an integral, like so:

$[f*g](t)=\int f(t)g(t-\tau)d\tau$ (4)

In the discrete case, this can be expressed as a sum, like so:

$[f*g](t)=\sum_{\tau}f(t)g(t-\tau)$ (5)

The sums can be across multiple axes

$[f*g](t,s,r)=\sum_{\tau}\sum_{\sigma}\sum_{\rho}f(t,s,r)g(t-\tau,s-\sigma,r-\rho)$
(6)

A 2D Convolution of an image $I$, that can be represented as an $H\times
W\times C$ tensor (usually for black and white images, C=1, for color images
C=3), and a kernel $K$, that can be represented as an $H_{k}\times W_{k}\times
C$ tensor, (where typically $H_{k}<<<H$ and $W_{k}<<<W$) is usually of the
form:

$\mathbf{Conv2D}(x,y)=\sum_{c}^{C}\sum_{h}^{H_{k}}\sum_{w}^{W_{k}}I_{x,y,c}K_{x-h,y-w,c}$
(7)

Often, convolutions do not sample every possible $x,y$ value. The ’stride’
denotes the amount of pixels between each pixel to apply the convolution
function to. For example, with a $stride=(2,2$, we would only perform:
$\mathbf{Conv2D}(x,y),\mathbf{Conv2D}(x+2,y),\mathbf{Conv2D}(x,y+2),\mathbf{Conv2D}(x+2,y+2)$,
skipping intermediate values of $x,y$. Thus, each 2D Convolution produces a
new $H^{\prime\prime}\times W^{\prime\prime}$ output map. If there are
$C^{\prime\prime}$ different kernels (all with the same shape), then we can
concatenate all of the output maps to create a $H^{\prime\prime}\times
W^{\prime\prime}\times C^{\prime\prime}$ tensor. The weights of each kernel
are tunable parameters that are optimized during training [16].

###  A.2 Attention

When humans process input from a sequence or an image, we contextualize each
part of the input using other parts of the input. However, we do not pay equal
attention to every other part of the input [28]. For a textual example, in the
sentence ”Alan said he was hungry”, the word ”alan” is more useful in
determining the meaning of ”he” than the word ”hungry”. For a visual example,
in order to guess the location of the right eye on someones face, the most
important information would be the location of the left eye, not the length of
their beard. Attention, also known as Non-Locality was first applied to text
sequences, for neural machine translation [25, 5, 17]. Later, networks with
attention were used for image generation [18], video captioning [52] and
object recognition [4]. Attention [55] is formulated as

$Attention(Q,K,V)=\mathrm{softmax}(\mathrm{score}(Q,K))V$ (8)

$Q,K,V$ are query, key and value matrices. When $Q=K=V$, this is called self-
attention. Softmax [31] is defined as:

$\mathrm{softmax}{(x_{i})}=\frac{\exp{x_{i}}}{\sum_{j}\exp{x_{j}}}$ (9)

In [25], they represent the hidden states of each sequence of the input,
derived from the output of the LSTM layers that precede the attentional layers
in the model. In other contexts, given an input vector $x_{i}$,
$Q_{i}=W^{Q}x_{i},K_{i}=W^{K}x_{i},V_{i}=W^{V}x_{i}$, where
$W^{Q},W^{K},W^{K}$ are trainable weight matrices [2]. Two score functions
were used, one from Bahdanau [5]:

$\mathrm{BahdanauScore}=V^{T}\tanh{(K+Q)}$ (10)

and one from Luong [25]

$\mathrm{LuongScore}=K^{T}Q$ (11)

Following the work of [9], we performed self-attention across channels. Given
an image $I\in\mathbb{R}^{H\times L\times D}$, where $H,L,D$ are the height,
length and depth of the input images after a few layers of convolutions, for
each $(h,w)\in\mathbb{R}^{H\times W}$, we performed feature-wise, or
equivalently channel-wise attention, by performing the self-attention on the
corresponding channel vector $\in\mathbb{R}^{D}$, for both Luong and Bahdanau
scoring.

To perform self-attention across spatial dimensions, we used Transformers
[47]. The Transformer consisted of alternating fully connected and multi-head
attention (MHA) layers, the latter being multiple attention layers in
parallel, whose outputs were concatenated to produce the output of the MHA
layer. Given weight matrix $W^{O}$:

$\mathrm{MultiHead}(Q,K,V)=\mathrm{concat}(h_{0},h_{1}...h_{i})W^{O}$ (12)

Where each $h_{j}$ is the output of an Attention layer with its own set of
weights:

$h_{j}=Attention(QW^{Q}_{j},KW^{K}_{j},VW^{V}_{j})$ (13)

In this paper, the Attention function used Luong scoring, as is standard in
the literature. [11] introduced the Vision Transformer (ViT), which performs
patch embedding by reshaping $x\in\mathbb{R}^{H\times L\times D}\rightarrow
x_{p}\in\mathbb{R}^{\frac{HL}{P^{2}}\times P^{2}D}$, where $P\times P$ is the
size of each patch in pixels. This patch embedded input is then passed to a
transformer.

###  A.3 Residual Networks

Increasingly, deeper and larger convolutional networks have been used for
vision tasks [45, 23]. However, this led to the degradation problem, where
accuracy for classification tasks would saturate after a certain level of
depth. To solve this, [20] proposed the residual connection, where the input
of past layers would be added to the output of past layers. [36] proved that
ResNets were also capable of circumventing the problem of exploding gradients.
[50] further improved on this by introducing the ResNeXt, which offered
superior accuracy to ResNet models with the same number of parameters. The
output of a ResNeXt layer $y$ with input $x$ is

$y=x+\sum_{i}^{C}T_{i}(x)$ (14)

Where each $T_{i}$ is a transformation, $C$, known as the cardinality, is the
amount of transformations to be applied. This is equivalently implemented as
the concatenation of the outputs of convolutional layers, an aggregation of
residual layers, or a grouped convolution [23].

##  Appendix B Visualization

Figure 3: Images generated for activation maximization of second to last layer
of Model A-1 using gradient ascent

Figure 4: Images generated for activation maximization of last layer using
gradient ascent