0}L(r)r^{-\frac{1}{2}}(2-h^{2}r)^{-\frac{1}{2}}\mathbbm{1}_{[0,h^{-2})}(r)=\varpi(r)$
and $\int_{0}^{\infty}\varpi(r)dr<\infty$ by condition (D2), the Dominated
Convergence Theorem guarantees that

$\lim\limits_{h\to
0}\int_{0}^{\infty}L(r)r^{-\frac{1}{2}}(2-h^{2}r)^{-\frac{1}{2}}\mathbbm{1}_{[0,h^{-2})}(r)dr=\int_{0}^{\infty}\varpi(r)dr.$

For the second integral, as a consequence of condition (D2), $L$ must be
decrease faster than any power function in order for
$0<\int_{0}^{\infty}L^{k}(r)r^{\frac{q}{2}-1}dr<\infty$ for all $q\geq 1$ and
$k=1,2$. In particular, for some fixed $h_{0}>0$, $L(r)\leq r^{-1}$, $\forall
r\in[h^{-2},2h^{-2})$, $\forall h\in(0,h_{0})$. Using this, it results in:

$\lim_{h\to
0}\int_{h^{-2}}^{2h^{-2}}L(r)r^{-\frac{1}{2}}(2-h^{2}r)^{-\frac{1}{2}}dr\leq\lim_{h\to
0}\int_{h^{-2}}^{2h^{-2}}r^{-\frac{3}{2}}(2-h^{2}r)^{-\frac{1}{2}}dr=\lim_{h\to
0}h=0.$

This completes the proof.

The proofs of (b) and the first two integral results in (c) can be found in
García-Portugués et al. (2013) and thus omitted. We adopt some of the argument
of Lemma 3 in García-Portugués et al. (2013) to prove the last three integrals
in (c).
Recall that the $n$-dimensional spherical coordinates of
$\bm{x}=(x_{1},...,x_{n})^{T}$ with norm
$r:=\left|\left|\bm{x}\right|\right|_{2}$ are given by

$\begin{cases}x_{1}=r\cos\phi_{1},\\\
x_{j}=r\cos\phi_{j}\prod\limits_{k=1}^{j-1}\sin\phi_{k},\quad j=2,...,n-2,\\\
x_{n-1}=r\sin\theta\prod\limits_{k=1}^{n-2}\sin\phi_{k},\\\
x_{n}=r\cos\theta\prod\limits_{k=1}^{n-2}\sin\phi_{k},\end{cases}\quad
J=r^{n-1}\prod\limits_{k=1}^{n-2}\sin^{k}\phi_{n-1-k},$ (40)

where $0\leq\phi_{j}\leq\pi$, $j=1,...,n-2$, $0\leq\theta\leq 2\pi$, and
$0\leq r<\infty$. $J$ denotes the Jacobian of the transformation. Without loss
of generality, we assume, by the $q$-spherical coordinates (40), that
$x_{i}=\cos\phi_{1}$, $x_{j}=\cos\phi_{2}\sin\phi_{1}$, and
$x_{k}=\cos\phi_{3}\sin\phi_{2}\sin\phi_{1}$. Then,

$\displaystyle\int_{\Omega_{q}}x_{i}^{3}\,\omega_{q}(d\bm{x})$
$\displaystyle=\int_{0}^{2\pi}\int_{0}^{\pi}\times\stackrel{{\scriptstyle(q-1)}}{{\cdots}}\times\int_{0}^{\pi}\cos^{3}\phi_{1}\prod_{k=1}^{q-2}\sin^{k}\phi_{q-k}\sin^{q-1}\phi_{1}\prod_{j=q-1}^{1}d\phi_{j}d\theta$
$\displaystyle=\int_{0}^{2\pi}\int_{0}^{\pi}\times\stackrel{{\scriptstyle(q-2)}}{{\cdots}}\times\int_{0}^{\pi}\prod_{k=1}^{q-2}\sin^{k}\phi_{q-k}\prod_{j=q-1}^{2}d\phi_{j}d\theta\times\int_{0}^{\pi}\cos^{3}\phi_{1}\sin^{q-1}\phi_{1}d\phi_{1}$
$\displaystyle=\int_{0}^{2\pi}\int_{0}^{\pi}\times\stackrel{{\scriptstyle(q-2)}}{{\cdots}}\times\int_{0}^{\pi}\prod_{k=1}^{q-2}\sin^{k}\phi_{q-k}\prod_{j=q-1}^{2}d\phi_{j}d\theta\times\int_{0}^{\pi}(1-\sin^{2}\phi_{1})\sin^{q-1}\phi_{1}d(\sin\phi_{1})$
$\displaystyle=\bar{\omega}_{q-1}\times 0=0,$
$\displaystyle\int_{\Omega_{q}}x_{i}^{2}x_{j}\omega_{q}(d\bm{x})$
$\displaystyle=\int_{0}^{2\pi}\int_{0}^{\pi}\times\stackrel{{\scriptstyle(q-1)}}{{\cdots}}\times\int_{0}^{\pi}\cos^{2}\phi_{1}\cos\phi_{2}\sin\phi_{1}\prod_{k=1}^{q-3}\sin^{k}\phi_{q-k}\sin^{q-2}\phi_{2}\sin^{q-1}\phi_{1}\prod_{j=q-1}^{1}d\phi_{j}d\theta$
$\displaystyle=\int_{0}^{2\pi}\int_{0}^{\pi}\times\stackrel{{\scriptstyle(q-3)}}{{\cdots}}\times\int_{0}^{\pi}\prod_{k=1}^{q-3}\sin^{k}\phi_{q-k}\prod_{j=q-1}^{3}d\phi_{j}d\theta$
$\displaystyle\quad\times\int_{0}^{\pi}\cos^{2}\phi_{1}\sin^{q}\phi_{1}d\phi_{1}\int_{0}^{\pi}\cos\phi_{2}\sin^{q-2}\phi_{2}d\phi_{2}$
$\displaystyle=\bar{\omega}_{q-2}\times\int_{0}^{\pi}\cos^{2}\phi_{1}\sin^{q}\phi_{1}d\phi_{1}\times
0=0,$

and

$\displaystyle\int_{\Omega_{q}}x_{i}x_{j}x_{k}\omega_{q}(d\bm{x})$
$\displaystyle=\int_{0}^{2\pi}\int_{0}^{\pi}\times\stackrel{{\scriptstyle(q-1)}}{{\cdots}}\times\int_{0}^{\pi}\cos\phi_{1}\cos\phi_{2}\sin\phi_{1}\cos\phi_{3}\sin\phi_{2}\sin\phi_{1}$
$\displaystyle\quad\times\prod_{k=1}^{q-4}\sin^{k}\phi_{q-k}\sin^{q-3}\phi_{3}\sin^{q-2}\phi_{2}\sin^{q-1}\phi_{1}\prod_{j=q-1}^{1}d\phi_{j}d\theta$
$\displaystyle=\int_{0}^{2\pi}\int_{0}^{\pi}\times\stackrel{{\scriptstyle(q-4)}}{{\cdots}}\times\int_{0}^{\pi}\prod_{k=1}^{q-4}\sin^{k}\phi_{q-k}\prod_{j=q-1}^{4}d\phi_{j}d\theta$
$\displaystyle\quad\times\int_{0}^{\pi}\cos\phi_{1}\sin^{q}\phi_{1}d\phi_{1}\int_{0}^{\pi}\cos\phi_{2}\sin^{q-1}\phi_{2}d\phi_{2}\int_{0}^{\pi}\cos\phi_{3}\sin^{q-2}\phi_{3}d\phi_{3}$
$\displaystyle=\bar{\omega}_{q-3}\times 0\times 0\times 0=0.$

The preceding argument teaches us that

$\int_{\Omega_{q}}x_{i}x_{j}x_{k}x_{m}\,\omega_{q}(d\bm{x})=\int_{\Omega_{q}}x_{i}x_{j}x_{k}x_{m}x_{\ell}\,\omega_{q}(d\bm{x})=0$

as long as one of the unique factors in the integrand has an odd multiplicity.
(Indeed, any integration of a monomial with an odd degree on $\Omega_{q}$ will
yield 0.) Thus, the only nonzero integrals in
$\int_{\Omega_{q}}x_{i}x_{j}x_{k}x_{m}\,\omega_{q}(d\bm{x})$ and
$\int_{\Omega_{q}}x_{i}x_{j}x_{k}x_{m}x_{\ell}\,\omega_{q}(d\bm{x})$ are

$\int_{\Omega_{q}}x_{i}^{4}\,\omega_{q}(d\bm{x})\quad\text{ and
}\quad\int_{\Omega_{q}}x_{i}^{2}x_{j}^{2}\,\omega_{q}(d\bm{x})$

with $i\neq j$. To compute the first integral, we define a vector field as

$\bm{F}(\bm{x})=\left(F_{1}(\bm{x}),...,F_{q+1}(\bm{x})\right)=(x_{1}^{3},...,x_{q+1}^{3})$

with $\bm{x}=(x_{1},...,x_{q+1})\in\Omega_{q}$. By the divergence theorem
(Theorem 10.51 in Rudin 1976),

$\displaystyle\int_{\Omega_{q}}x_{i}^{4}\,\omega_{q}(d\bm{x})$
$\displaystyle=\frac{1}{q+1}\int_{\Omega_{q}}\left(\sum_{i=1}^{q+1}x_{i}^{4}\right)\omega_{q}(d\bm{x})$
$\displaystyle=\frac{1}{q+1}\int_{\Omega_{q}}\langle\bm{F},\bm{x}\rangle\,\omega_{q}(d\bm{x})$
$\displaystyle=\frac{1}{q+1}\int_{V_{q}}\mathtt{div}\,\bm{F}\,dV$
$\displaystyle=\frac{3}{q+1}\int_{0}^{1}r^{2}\cdot
r^{q}\bar{\omega}_{q}dr=\frac{3\bar{\omega}_{q}}{(q+1)(q+3)},$

where $\langle\cdot,\cdot\rangle$ is the usual inner product in
$\mathbb{R}^{q+1}$,
$\mathtt{div}\,\bm{F}=\sum\limits_{i=1}^{q+1}\frac{\partial F_{i}}{\partial
x_{i}}$, and $\int_{V_{q}}\cdots dV$ is integrating the solid $q$-dimensional
sphere $V_{q}$ in $\mathbb{R}^{q+1}$. The second integral can be evaluated
based on the preceding results as

$\displaystyle\int_{\Omega_{q}}x_{i}^{2}x_{j}^{2}\,\omega_{q}(d\bm{x})$
$\displaystyle=\frac{1}{q}\int_{\Omega_{q}}x_{i}^{2}\left(\sum_{j\neq
i}x_{j}^{2}\right)\,\omega_{q}(d\bm{x})$
$\displaystyle=\frac{1}{q}\int_{\Omega_{q}}(x_{i}^{2}-x_{i}^{4})\,\omega_{q}(d\bm{x})$
$\displaystyle=\frac{1}{q}\left[\frac{\bar{\omega}_{q}}{q+1}-\frac{3\bar{\omega}_{q}}{(q+1)(q+3)}\right]=\frac{\bar{\omega}_{q}}{(q+1)(q+3)}.$

As a specific application of our above results, we know that
$\int_{\Omega_{q-1}}\bm{B_{x}}\bm{\xi}\,\omega_{q-1}(d\bm{\xi})=0$.

###### Remark 22

García-Portugués et al. (2013) also provided a key remark about how to
generalize the arguments in (a) of Lemma 21. Under condition (D2’), one can
apply the same techniques in (a) to prove the result with the functions

$\begin{cases}\varpi_{h,i,j,k}(r)=L^{k}(r)r^{\frac{q}{2}+i}(2-h^{2}r)^{\frac{q}{2}-j}\mathbbm{1}_{[0,2h^{-2})}(r),\\\
\varpi_{i,j,k}(r)=\lim\limits_{h\to
0}\varpi_{h,i,j,k}(r)=L^{k}(r)r^{\frac{q}{2}+i}2^{\frac{q}{2}-j}\mathbbm{1}_{[0,\infty)}(r);\end{cases}$
$\begin{cases}\varpi_{h,i,j,k}^{\prime}(r)=[L^{\prime}(r)]^{k}r^{\frac{q}{2}+i}(2-h^{2}r)^{\frac{q}{2}-j}\mathbbm{1}_{[0,2h^{-2})}(r),\\\
\varpi_{i,j,k}^{\prime}(r)=\lim\limits_{h\to
0}\varpi_{h,i,j,k}^{\prime}(r)=[L^{\prime}(r)]^{k}r^{\frac{q}{2}+i}2^{\frac{q}{2}-j}\mathbbm{1}_{[0,\infty)}(r);\end{cases}$
$\begin{cases}\varpi_{h,i,j,k}^{\prime\prime}(r)=[L^{\prime\prime}(r)]^{k}r^{\frac{q}{2}+i}(2-h^{2}r)^{\frac{q}{2}-j}\mathbbm{1}_{[0,2h^{-2})}(r),\\\
\varpi_{i,j,k}^{\prime\prime}(r)=\lim\limits_{h\to
0}\varpi_{h,i,j,k}^{\prime\prime}(r)=[L^{\prime\prime}(r)]^{k}r^{\frac{q}{2}+i}2^{\frac{q}{2}-j}\mathbbm{1}_{[0,\infty)}(r)\end{cases}$

with $i\geq-1$, $j\leq 1$, and $k=1,2$. For the case where $\frac{q}{2}-j\geq
0$, use the Dominated Convergence Theorem. For the other cases, subdivide the
integral over $[0,2h^{-2})$ into the intervals $[0,h^{-2})$ and
$[h^{-2},2h^{-2})$. Then apply the Dominated Convergence Theorem in the former
and use a suitable power function to make the latter tend to 0 in the same way
as described in the proof of (a) in Lemma 21.

###### Theorem 2

Assume conditions (D1) and (D2’). For any fixed $\bm{x}\in\Omega_{q}$, we have

$\mathtt{grad}\,\widehat{f}_{h}(\bm{x})-\mathtt{grad}\,f(\bm{x})=O(h^{2})+O_{P}\left(\sqrt{\frac{1}{nh^{q+2}}}\right)$

as $h\to 0$ and $nh^{q+2}\to\infty$.
Under the same condition, for any fixed $\bm{x}\in\Omega_{q}$, we have

$\mathcal{H}\widehat{f}_{h}(\bm{x})-\mathcal{H}f(\bm{x})=O(h^{2})+O_{P}\left(\sqrt{\frac{1}{nh^{q+4}}}\right)$

as $h\to 0$ and $nh^{q+4}\to\infty$.

Proof  Part A: Pointwise convergence rate of the Riemannian gradient estimator
$\mathtt{grad}\,\widehat{f}_{h}(\bm{x})$. Recall from Section 4.1 that the
tangent/Riemannian gradient estimator of a directional KDE is uniquely defined
under a given kernel function $L$. Thus, we can establish the pointwise
convergence rate under any total gradient (or differential) estimator, that
is,
$\mathtt{grad}\,\widehat{f}_{h}(\bm{x})=\mathtt{grad}\,\widetilde{f}_{h}(\bm{x})\equiv\mathtt{Tang}\left(\nabla\widetilde{f}_{h}(\bm{x})\right)$.
Here, we stick on the differential form (21),
$\nabla\widetilde{f}_{h}(\bm{x})$. (One may also prove Theorem 2 with
$\nabla\widehat{f}_{h}(\bm{x})$. The proof of Lemma 10 provides a starting
point for this direction.)

$\bullet$ Result 1: The expectation of the Riemannian gradient estimator,
$\mathbb{E}\left[\mathtt{grad}\,\widehat{f}_{h}(\bm{x})\right]$, has the
following asymptotic behavior as $h\to 0$:

$\displaystyle\mathbb{E}\left[\mathtt{grad}\,\widetilde{f}_{h}(\bm{x})\right]$
$\displaystyle=\mathbb{E}\left[\mathtt{Tang}\left(\nabla\widetilde{f}_{h}(\bm{x})\right)\right]=(I_{q+1}-\bm{x}\bm{x}^{T})\mathbb{E}\left[\nabla\widetilde{f}_{h}(\bm{x})\right]$
$\displaystyle=\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})+\frac{h^{2}}{2}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})\cdot\frac{\int_{0}^{\infty}L(r)r^{\frac{q}{2}}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}$
$\displaystyle\quad+\frac{2h^{2}}{q}\sum_{i=1}^{q}\left(\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}\right)\bm{b}_{i}\cdot\frac{\int_{0}^{\infty}L^{\prime}(r)r^{\frac{q}{2}+1}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}+O(h^{2})+o(h^{2})$
$\displaystyle=\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})+O(h^{2}).$

Derivation of Result 1. With the definition of $\bm{B_{x}}$ from Lemma 21, the
expected value of $\nabla\widetilde{f}_{h}(\bm{x})$ is

$\displaystyle\begin{split}\mathbb{E}\left[\nabla\widetilde{f}_{h}(\bm{x})\right]&=\frac{c_{h,q}(L)}{h^{2}}\int_{\Omega_{q}}(\bm{x}-\bm{y})\cdot
L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{y}}{h^{2}}\right)f(\bm{y})\,\omega_{q}(d\bm{y})\\\
&=\frac{c_{h,q}(L)}{h^{2}}\int_{-1}^{1}\int_{\Omega_{q-1}}\left(\bm{x}-t\bm{x}-\sqrt{1-t^{2}}\bm{B_{x}}\bm{\xi}\right)L^{\prime}\left(\frac{1-t}{h^{2}}\right)\\\
&\quad\times
f\left(t\bm{x}+\sqrt{1-t^{2}}\bm{B_{x}}\bm{\xi}\right)(1-t^{2})^{\frac{q}{2}-1}\,\omega_{q-1}(d\bm{\xi})dt\\\
&=c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\left(rh^{2}\bm{x}-h\sqrt{r(2-h^{2}r)}\bm{B_{x}}\bm{\xi}\right)L^{\prime}(r)\\\
&\quad\times f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})\cdot
r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr\end{split}$
(41)

by (a) in Lemma 21 and a change of variable $r=\frac{1-t}{h^{2}}$, where
$\alpha_{\bm{x},\bm{\xi}}=-rh^{2}\bm{x}+h\sqrt{r(2-h^{2}r)}\bm{B_{x}}\bm{\xi}$.
By condition (D1), the Taylor’s expansion of $f$ at $\bm{x}$ is

$\displaystyle f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})$
$\displaystyle=f(\bm{x})+\alpha_{\bm{x},\bm{\xi}}^{T}\nabla
f(\bm{x})+\frac{1}{2}\alpha_{\bm{x},\bm{\xi}}^{T}\nabla\nabla
f(\bm{x})\alpha_{\bm{x},\bm{\xi}}+\frac{1}{6}\left(\sum_{i=1}^{q+1}(\alpha_{\bm{x},\bm{\xi}})_{i}\cdot\frac{\partial}{\partial
x_{i}}\right)^{3}f(\bm{x})+o\left(||\alpha_{\bm{x},\bm{\xi}}||_{2}^{3}\right)$
$\displaystyle\equiv\text{(I) + (II) + (III) + (IV)}+o(h^{3}),$

where
$\left|\left|\alpha_{\bm{x},\bm{\xi}}\right|\right|_{2}^{2}=r^{2}h^{4}+h^{2}r(2-h^{2}r)=2rh^{2}$
by the orthogonality of $\bm{x}$ and columns of $\bm{B_{x}}$, and
$(\alpha_{\bm{x},\bm{\xi}})_{i}$ stands for the $i^{th}$ entry of the vector
$\alpha_{\bm{x},\bm{\xi}}$. Now we plug (I), (II), (III), (IV), and $o(h^{3})$
back into (41) respectively to compute the dominating term of
$\mathbb{E}\left[\nabla\widetilde{f}_{h}(\bm{x})\right]$.

Plug in (I)
$\displaystyle=c_{h,q}(L)h^{q-2}f(\bm{x})\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\left(rh^{2}\bm{x}-h\sqrt{r(2-h^{2}r)}\bm{B_{x}}\bm{\xi}\right)$
$\displaystyle\quad\quad\times
L^{\prime}(r)\,r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\stackrel{{\scriptstyle\text{(i)}}}{{=}}\bar{\omega}_{q-1}\cdot\bm{x}f(\bm{x})\int_{0}^{2h^{-2}}c_{h,q}(L)h^{q}L^{\prime}(r)\cdot
r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-1}dr+0$
$\displaystyle\stackrel{{\scriptstyle\text{(ii)}}}{{=}}\bar{\omega}_{q-1}\cdot\bm{x}f(\bm{x})\cdot
c_{h,q}(L)h^{q}\bigg{\\{}L(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-1}\Big{|}_{0}^{2h^{-2}}$
$\displaystyle\hskip 85.35826pt-\int_{0}^{2h^{-2}}L(r)\left[\frac{q}{2}\cdot
r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}-h^{2}\left(\frac{q}{2}-1\right)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-2}\right]dr\bigg{\\}}$
$\displaystyle=-\frac{q}{2}\cdot\bar{\omega}_{q-1}\cdot\bm{x}f(\bm{x})\cdot
c_{h,q}(L)h^{q}\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}dr$
$\displaystyle\quad+\left(\frac{q-2}{2}\right)\bar{\omega}_{q-1}\cdot\bm{x}f(\bm{x})\cdot
c_{h,q}(L)h^{q+2}\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-2}dr$
$\displaystyle\stackrel{{\scriptstyle\text{(iii)}}}{{=}}-\frac{q}{2}\cdot\bm{x}f(\bm{x})+\left(\frac{q-2}{2}\right)\bm{x}f(\bm{x})h^{2}\cdot\frac{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-2}dr}{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}dr}$
$\displaystyle\stackrel{{\scriptstyle\text{(iv)}}}{{=}}-\frac{q}{2}\cdot\bm{x}f(\bm{x})+\left(\frac{q-2}{4}\right)\bm{x}f(\bm{x})h^{2}\cdot\frac{\int_{0}^{\infty}L(r)r^{\frac{q}{2}}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}+o(h^{2}),$

as $h\to 0$, where we use (c) of Lemma 21 with
$\bm{B_{x}}\bm{\xi}=\sum_{i=1}^{q}\xi_{i}\bm{b}_{i}$ in (i), conduct
integration by parts in (ii), plug in the expression (3) of $c_{h,q}(L)$ in
(iii), and take $h\to 0$ with our argument in (a) of Lemma 21 and Remark 22 to
obtain (iv). The $o(h^{2})$-term in (iv) takes into account those small error
terms as $h\to 0$. Likewise,

Plug in (II)
$\displaystyle=c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\left(rh^{2}\bm{x}-h\sqrt{r(2-h^{2}r)}\bm{B_{x}}\bm{\xi}\right)\alpha_{\bm{x},\bm{\xi}}^{T}f(\bm{x})$
$\displaystyle\quad\quad\times L^{\prime}(r)\cdot
r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle=-c_{h,q}(L)h^{q+2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\bm{x}\bm{x}^{T}\nabla
f(\bm{x})L^{\prime}(r)r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\quad+c_{h,q}(L)h^{q+1}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\bm{x}\bm{\xi}^{T}\bm{B_{x}}^{T}\nabla
f(\bm{x})L^{\prime}(r)r^{\frac{q+1}{2}}(2-h^{2}r)^{\frac{q-1}{2}}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\quad+c_{h,q}(L)h^{q+1}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\bm{B_{x}}\bm{\xi}\cdot\bm{x}^{T}\nabla
f(\bm{x})L^{\prime}(r)r^{\frac{q+1}{2}}(2-h^{2}r)^{\frac{q-1}{2}}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\quad-
c_{h,q}(L)h^{q}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\bm{B_{x}}\bm{\xi}\cdot\bm{\xi}^{T}\bm{B_{x}}^{T}\nabla
f(\bm{x})L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\stackrel{{\scriptstyle\text{(i)}}}{{=}}-c_{h,q}(L)h^{q+2}\cdot\bm{x}\bm{x}^{T}\nabla
f(\bm{x})\cdot\bar{\omega}_{q-1}\int_{0}^{2h^{-2}}L^{\prime}(r)r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}-1}dr+0+0$
$\displaystyle\quad-
c_{h,q}(L)h^{q}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\xi_{i}\bm{b}_{i}\right)\left(\sum_{i=1}^{q}\xi_{i}\bm{b}_{i}^{T}\nabla
f(\bm{x})\right)L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\stackrel{{\scriptstyle\text{(ii)}}}{{=}}-c_{h,q}(L)h^{q+2}\cdot\bm{x}\bm{x}^{T}\nabla
f(\bm{x})\cdot\bar{\omega}_{q-1}\Bigg{\\{}r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}-1}L(r)\Big{|}_{0}^{2h^{-2}}$
$\displaystyle\hskip
56.9055pt-\int_{0}^{2h^{-2}}L(r)\left[\left(\frac{q+2}{2}\right)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-1}-h^{2}\left(\frac{q-2}{2}\right)r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}-2}\right]dr\Bigg{\\}}$
$\displaystyle\quad-\frac{\bar{\omega}_{q-1}}{q}\left(\sum_{i=1}^{q}\bm{b}_{i}\bm{b}_{i}^{T}\right)\nabla
f(\bm{x})\cdot
c_{h,q}(L)h^{q}\int_{0}^{2h^{-2}}L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}dr$
$\displaystyle\stackrel{{\scriptstyle\text{(iii)}}}{{=}}c_{h,q}(L)h^{q+2}\cdot\bm{x}\bm{x}^{T}\nabla
f(\bm{x})\cdot\bar{\omega}_{q-1}\left(\frac{q+2}{2}\right)\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-1}dr$
$\displaystyle\quad-c_{h,q}(L)h^{q+4}\cdot\bm{x}\bm{x}^{T}\nabla
f(\bm{x})\cdot\bar{\omega}_{q-1}\left(\frac{q-2}{2}\right)\int_{0}^{2h^{-2}}r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}-2}dr$
$\displaystyle\quad-\frac{\bar{\omega}_{q-1}}{q}(I_{q+1}-\bm{x}\bm{x}^{T})\nabla
f(\bm{x})\cdot
c_{h,q}(L)h^{q}\Bigg{[}L(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}\Big{|}_{0}^{2h^{-2}}$
$\displaystyle\hskip
113.81102pt-\frac{q}{2}\int_{0}^{2h^{-2}}L(r)\left(r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}}-h^{2}r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-1}\right)dr\Bigg{]}$
$\displaystyle\stackrel{{\scriptstyle\text{(iv)}}}{{=}}\left(\frac{q+2}{2}\right)h^{2}\bm{x}\bm{x}^{T}\nabla
f(\bm{x})\cdot\frac{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-1}dr}{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}dr}$
$\displaystyle\quad-\left(\frac{q-2}{2}\right)h^{4}\bm{x}\bm{x}^{T}\nabla
f(\bm{x})\cdot\frac{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}-2}dr}{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}dr}$
$\displaystyle\quad+\frac{1}{2}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})\Bigg{[}\frac{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}}dr}{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}dr}$
$\displaystyle\hskip
142.26378pt-h^{2}\cdot\frac{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-1}dr}{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}dr}\Bigg{]}$
$\displaystyle\stackrel{{\scriptstyle\text{(v)}}}{{=}}\left(\frac{q+2}{2}\right)h^{2}\cdot\bm{x}\bm{x}^{T}\nabla
f(\bm{x})\cdot\frac{\int_{0}^{\infty}L(r)r^{\frac{q}{2}}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}-\left(\frac{q-2}{4}\right)h^{4}\cdot\bm{x}\bm{x}^{T}\nabla
f(\bm{x})\cdot\frac{\int_{0}^{\infty}L(r)r^{\frac{q}{2}+1}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}$
$\displaystyle\quad+\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})+O(h^{2})-\frac{h^{2}}{2}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})\cdot\frac{\int_{0}^{\infty}L(r)r^{\frac{q}{2}}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}+o(h^{2})$
$\displaystyle=(I_{q+1}-\bm{x}\bm{x}^{T})\nabla
f(\bm{x})+h^{2}\left(\frac{q+2}{2}\right)\bm{x}\bm{x}^{T}\nabla
f(\bm{x})\cdot\frac{\int_{0}^{\infty}L(r)r^{\frac{q}{2}}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}$
$\displaystyle\quad+\frac{h^{2}}{2}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})\cdot\frac{\int_{0}^{\infty}L(r)r^{\frac{q}{2}}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}+O(h^{2})+o(h^{2}),$

where we use (c) of Lemma 21 in (i) and (ii), leverage the fact that
$\sum_{i=1}^{q}\bm{b}_{i}\bm{b}_{i}^{T}=\bm{B_{x}}\bm{B_{x}}^{T}=I_{q+1}-\bm{x}\bm{x}^{T}$
in (iii), plug in the expression (3) of $c_{h,q}(L)$ in (iv), and take $h\to
0$ with arguments in Lemma 21 and Remark 22 to obtain (v). The $o(h^{2})$-term
incorporates higher-order error terms, while the $O(h^{2})$-term in (v) comes
from the following arguments:

$\frac{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}}dr}{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}dr}-2=-h^{2}\cdot\frac{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-1}dr}{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}dr}=O(h^{2}).$
(42)

We now move on to the calculation of (III), which is more complicated.

$\displaystyle\begin{split}\text{Plug in
(III)}&=c_{h,q}(L)h^{q}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\frac{1}{2}\alpha_{\bm{x},\bm{\xi}}^{T}\nabla\nabla
f(\bm{x})\alpha_{\bm{x},\bm{\xi}}\cdot\bm{x}L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr\\\
&-c_{h,q}(L)h^{q-1}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\frac{1}{2}\alpha_{\bm{x},\bm{\xi}}^{T}\nabla\nabla
f(\bm{x})\alpha_{\bm{x},\bm{\xi}}\cdot\bm{B_{x}}\bm{\xi}L^{\prime}(r)r^{\frac{q-1}{2}}(2-h^{2}r)^{\frac{q-1}{2}}\omega_{q-1}(d\bm{\xi})dr.\end{split}$
(43)

Notice that

$\displaystyle\begin{split}&\int_{\Omega_{q-1}}\alpha_{\bm{x},\bm{\xi}}^{T}\nabla\nabla
f(\bm{x})\alpha_{\bm{x},\bm{\xi}}\cdot\bm{x}\,\omega_{q-1}(\bm{\xi})\\\
&=r^{2}h^{4}\int_{\Omega_{q-1}}\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\cdot\bm{x}\,\omega_{q-1}(d\bm{\xi})\\\
&\quad-2rh^{3}\sqrt{r(2-h^{2}r)}\int_{\Omega_{q-1}}\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{B_{x}}\bm{\xi}\cdot\bm{x}\,\omega_{q-1}(d\bm{\xi})\\\
&\quad+h^{2}r(2-h^{2}r)\int_{\Omega_{q-1}}\bm{\xi}^{T}\bm{B_{x}}^{T}\nabla\nabla
f(\bm{x})\bm{B_{x}}\bm{\xi}\cdot\bm{x}\,\omega_{q-1}(d\bm{\xi})\\\
&=r^{2}h^{4}\bar{\omega}_{q-1}\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\cdot\bm{x}+h^{2}r(2-h^{2}r)\int_{\Omega_{q-1}}\left(\sum_{i,j=1}^{q}\bm{b}_{i}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{j}\xi_{i}\xi_{j}\right)\bm{x}\,\omega_{q-1}(d\bm{\xi})\\\
&=r^{2}h^{4}\bar{\omega}_{q-1}\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\cdot\bm{x}+h^{2}r(2-h^{2}r)\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\bm{b}_{i}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}\xi_{i}^{2}\right)\bm{x}\,\omega_{q-1}(d\bm{\xi})\\\
&=r^{2}h^{4}\bar{\omega}_{q-1}\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\cdot\bm{x}+h^{2}r(2-h^{2}r)\cdot\frac{\bar{\omega}_{q-1}}{q}\left[\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\right]\cdot\bm{x},\end{split}$ (44)

where we use (c) of Lemma 21 in the second, third, and fourth equations and
the fact that

$\displaystyle\sum_{i=1}^{q}\bm{b}_{i}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}=\text{tr}\left[\nabla\nabla
f(\bm{x})\sum_{i=1}^{q}\bm{b}_{i}\bm{b}_{i}^{T}\right]$
$\displaystyle=\text{tr}\left[\nabla\nabla
f(\bm{x})(I_{q+1}-\bm{x}\bm{x}^{T})\right]$ $\displaystyle=\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla f(\bm{x})\bm{x}.$

Here, $\Delta f(\bm{x})=\sum_{i=1}^{q+1}\frac{\partial^{2}}{\partial
x_{i}^{2}}f(\bm{x})$ is the Laplace of function $f$. At the same time,

$\displaystyle\int_{\Omega_{q-1}}\alpha_{\bm{x},\bm{\xi}}^{T}\nabla\nabla
f(\bm{x})\alpha_{\bm{x},\bm{\xi}}\cdot\bm{B_{x}}\bm{\xi}\,\omega_{q-1}(d\bm{\xi})$
$\displaystyle=r^{2}h^{4}\int_{\Omega_{q-1}}\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\cdot\bm{B_{x}}\bm{\xi}\,\omega_{q-1}(d\bm{\xi})$
$\displaystyle\quad-2rh^{3}\sqrt{r(2-h^{2}r)}\int_{\Omega_{q-1}}\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{B_{x}}\bm{\xi}\cdot\bm{B_{x}}\bm{\xi}\,\omega_{q-1}(d\bm{\xi})$
$\displaystyle\quad+h^{2}r(2-h^{2}r)\int_{\Omega_{q-1}}\bm{\xi}^{T}\bm{B_{x}}^{T}\nabla\nabla
f(\bm{x})\bm{B_{x}}\bm{\xi}\cdot\bm{B_{x}}\bm{\xi}\,\omega_{q-1}(d\bm{\xi})$
$\displaystyle=-2rh^{3}\sqrt{r(2-h^{2}r)}\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}\cdot\bm{b}_{i}\xi_{i}^{2}\right)\omega_{q-1}(d\bm{\xi})$
$\displaystyle\quad+h^{2}r(2-h^{2}r)\int_{\Omega_{q-1}}\left(\sum_{i,j=1}^{q}\bm{b}_{i}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{j}\xi_{i}\xi_{j}\right)\left(\sum_{k=1}^{q}\bm{b}_{k}\xi_{k}\right)\,\omega_{q-1}(d\bm{\xi})$
$\displaystyle=-2rh^{3}\sqrt{r(2-h^{2}r)}\cdot\frac{\bar{\omega}_{q-1}}{q}\sum_{i=1}^{q}\left(\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}\right)\bm{b}_{i},$

where we apply (c) of Lemma 21 in the last two equations. That is,

$\displaystyle\begin{split}\int_{\Omega_{q-1}}\alpha_{\bm{x},\bm{\xi}}^{T}\nabla\nabla
f(\bm{x})\alpha_{\bm{x},\bm{\xi}}\cdot\bm{B_{x}}\bm{\xi}\,\omega_{q-1}(d\bm{\xi})&=-2rh^{3}\sqrt{r(2-h^{2}r)}\cdot\frac{\bar{\omega}_{q-1}}{q}\sum_{i=1}^{q}\left(\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}\right)\bm{b}_{i}.\end{split}$ (45)

Plugging (44) and (45) back into (43), we proceed “Plug in (III)” as

Plug in (III)
$\displaystyle=\frac{c_{h,q}(L)h^{q+4}}{2}\cdot\bar{\omega}_{q-1}\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\cdot\bm{x}\int_{0}^{2h^{-2}}r^{\frac{q}{2}+2}(2-h^{2}r)^{\frac{q}{2}-1}L^{\prime}(r)dr$
$\displaystyle\quad+\frac{c_{h,q}(L)h^{q+2}}{2q}\cdot\bar{\omega}_{q-1}\left[\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\right]\bm{x}\int_{0}^{2h^{-2}}L^{\prime}(r)r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}}dr$
$\displaystyle\quad+\frac{c_{h,q}(L)h^{q+2}}{q}\cdot\bar{\omega}_{q-1}\sum_{i=1}^{q}\left(\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}\right)\bm{b}_{i}\int_{0}^{2h^{-2}}L^{\prime}(r)r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}}dr$
$\displaystyle\stackrel{{\scriptstyle\text{(i)}}}{{=}}\frac{h^{4}}{2}\cdot\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\cdot\bm{x}\cdot\frac{\int_{0}^{\infty}L^{\prime}(r)r^{\frac{q}{2}+2}(2-h^{2}r)^{\frac{q}{2}-1}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}dr}$
$\displaystyle\quad+\frac{h^{2}}{2q}\left[\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\right]\bm{x}\cdot\frac{\int_{0}^{\infty}L^{\prime}(r)r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}dr}$
$\displaystyle\quad+\frac{h^{2}}{q}\sum_{i=1}^{q}\left(\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}\right)\bm{b}_{i}\cdot\frac{\int_{0}^{\infty}L^{\prime}(r)r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}dr}$
$\displaystyle\stackrel{{\scriptstyle\text{(ii)}}}{{=}}\frac{h^{2}}{q}\left[\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\right]\bm{x}\cdot\frac{\int_{0}^{\infty}L^{\prime}(r)r^{\frac{q}{2}+1}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}$
$\displaystyle\quad+\frac{2h^{2}}{q}\sum_{i=1}^{q}\left(\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}\right)\bm{b}_{i}\cdot\frac{\int_{0}^{\infty}L^{\prime}(r)r^{\frac{q}{2}+1}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}+o(h^{2}),$

where we plug in the expression (3) of $c_{h,q}(L)$ in (i) and take $h\to 0$
with arguments in Lemma 21 and Remark 22 in (ii).

We argue that after plugging (IV)+$o(h^{3})$ back into (41), it yields a
$o(h^{2})$ term.

$\displaystyle\text{Plug in (IV)}+o(h^{3})$
$\displaystyle=c_{h,q}(L)h^{q}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\frac{1}{6}\left[\left(\sum_{i=1}^{q+1}(\alpha_{\bm{x},\bm{\xi}})_{i}\cdot\frac{\partial}{\partial
x_{i}}\right)^{3}f(\bm{x})\right]\bm{x}L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle+c_{h,q}(L)h^{q-1}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\frac{1}{6}\left[\left(\sum_{i=1}^{q+1}(\alpha_{\bm{x},\bm{\xi}})_{i}\cdot\frac{\partial}{\partial
x_{i}}\right)^{3}f(\bm{x})\right]\bm{B_{x}}\bm{\xi}$ $\displaystyle\hskip
170.71652pt\times
L^{\prime}(r)r^{\frac{q-1}{2}}(2-h^{2}r)^{\frac{q-1}{2}}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle+c_{h,q}(L)\cdot
o(h^{q+1})\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\left(rh^{2}\bm{x}-h\sqrt{r(2-h^{2}r)}\bm{B_{x}}\bm{\xi}\right)L^{\prime}(r)r^{\frac{q-1}{2}}(2-h^{2}r)^{\frac{q-1}{2}}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle=c_{h,q}(L)h^{q}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\frac{1}{6}\left[\left(\sum_{i=1}^{q+1}(\alpha_{\bm{x},\bm{\xi}})_{i}\cdot\frac{\partial}{\partial
x_{i}}\right)^{3}f(\bm{x})\right]\bm{x}L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle+c_{h,q}(L)h^{q-1}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\frac{1}{6}\left[\left(\sum_{i=1}^{q+1}(\alpha_{\bm{x},\bm{\xi}})_{i}\cdot\frac{\partial}{\partial
x_{i}}\right)^{3}f(\bm{x})\right]\bm{B_{x}}\bm{\xi}$ $\displaystyle\hskip
170.71652pt\times
L^{\prime}(r)r^{\frac{q-1}{2}}(2-h^{2}r)^{\frac{q-1}{2}}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle+c_{h,q}(L)\cdot
o(h^{q+3})\cdot\bar{\omega}_{q-1}\bm{x}\int_{0}^{2h^{-2}}L^{\prime}(r)r^{\frac{q+1}{2}}(2-h^{2}r)^{\frac{q-1}{2}}dr$

by (c) of Lemma 21 in the last equality. As
$c_{h,q}(L)=h^{q}\lambda_{h,q}(L)=O(h^{q})$ by (3) and (a) of Lemma 21, we
know that the third integral is of the order $o(h^{3})$. Since
$\bm{B_{x}}\bm{\xi}=\sum_{i=1}^{q}\xi_{i}\bm{b}_{i}$ and
$\alpha_{\bm{x},\bm{\xi}}=-rh^{2}\bm{x}+h\sqrt{r(2-h^{2}r)}\bm{B_{x}}\bm{\xi}$,
we derive that

$\displaystyle\left(\sum_{i=1}^{q+1}(\alpha_{\bm{x},\bm{\xi}})_{i}\cdot\frac{\partial}{\partial
x_{i}}\right)^{3}f(\bm{x})$ $\displaystyle=A_{f,1}\cdot
h^{3}r^{\frac{3}{2}}(2-h^{2}r)^{\frac{3}{2}}\sum_{i,j,k}\bm{b}_{i}\bm{b}_{j}\bm{b}_{k}\xi_{i}\xi_{j}\xi_{k}+A_{f,2}\cdot
h^{4}r^{2}(2-h^{2}r)\sum_{i,j}\bm{b}_{i}\bm{b}_{j}\xi_{i}\xi_{j}$
$\displaystyle\quad+A_{f,3}\cdot
h^{5}r^{\frac{5}{2}}(2-h^{2}r)^{\frac{1}{2}}\sum_{i}\bm{b}_{i}\xi_{i}+A_{f,4}h^{6}r^{3}$

and

$\displaystyle\left[\left(\sum_{i=1}^{q+1}(\alpha_{\bm{x},\bm{\xi}})_{i}\cdot\frac{\partial}{\partial
x_{i}}\right)^{3}f(\bm{x})\right]\bm{B_{x}}\bm{\xi}$
$\displaystyle=\widetilde{A}_{f,1}\cdot
h^{3}r^{\frac{3}{2}}(2-h^{2}r)^{\frac{3}{2}}\sum_{i,j,k,\ell}\bm{b}_{i}\bm{b}_{j}\bm{b}_{k}\bm{b}_{\ell}\xi_{i}\xi_{j}\xi_{k}\xi_{\ell}+\widetilde{A}_{f,2}\cdot
h^{4}r^{2}(2-h^{2}r)\sum_{i,j,k}\bm{b}_{i}\bm{b}_{j}\bm{b}_{k}\xi_{i}\xi_{j}\xi_{k}$
$\displaystyle\quad+\widetilde{A}_{f,3}\cdot
h^{5}r^{\frac{5}{2}}(2-h^{2}r)^{\frac{1}{2}}\sum_{i,j}\bm{b}_{i}\bm{b}_{j}\xi_{i}\xi_{j}+\widetilde{A}_{f,4}h^{6}r^{3}\sum_{i}\bm{b}_{i}\xi_{i},$

where $A_{f,i},i=1,...,4$ and $\widetilde{A}_{f,i},i=1,...,4$ are some
“constants” that depends on the partial derivatives of $f(\bm{x})$. Thus, by
(c) of Lemma 21 (that is, any integration of a monomial of $\bm{\xi}$ with an
odd degree on $\Omega_{q-1}$ will yield 0), we know that

$\int_{\Omega_{q-1}}\left[\left(\sum_{i=1}^{q+1}(\alpha_{\bm{x},\bm{\xi}})_{i}\cdot\frac{\partial}{\partial
x_{i}}\right)^{3}f(\bm{x})\right]\omega_{q-1}(d\bm{\xi})\asymp
h^{4}r^{2}(2-h^{2}r)+o(h^{4}),$
$\int_{\Omega_{q-1}}\left[\left(\sum_{i=1}^{q+1}(\alpha_{\bm{x},\bm{\xi}})_{i}\cdot\frac{\partial}{\partial
x_{i}}\right)^{3}f(\bm{x})\right]\bm{B_{x}}\bm{\xi}\,\omega_{q-1}(d\bm{\xi})\asymp
h^{3}r^{\frac{3}{2}}(2-h^{2}r)^{\frac{3}{2}}+h^{5}r\sqrt{2-h^{2}r}+o(h^{3}),$

where “$\asymp$” means an asymptotic equivalence. With condition (D2’) and our
arguments of (a) in Lemma 21 and Remark 22, we obtain that “Plug in
(IV)+$o(h^{3})$” yields a $o(h^{2})$ term. Therefore,

$\begin{split}\mathbb{E}\left[\nabla\widetilde{f}_{h}(\bm{x})\right]&=-\frac{q}{2}\cdot\bm{x}f(\bm{x})+\left(\frac{q-2}{4}\right)\bm{x}f(\bm{x})h^{2}\cdot\frac{\int_{0}^{\infty}L(r)r^{\frac{q}{2}}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}\\\
&\quad+(I_{q+1}-\bm{x}\bm{x}^{T})\nabla
f(\bm{x})+h^{2}\left(\frac{q+2}{2}\right)\bm{x}\bm{x}^{T}\nabla
f(\bm{x})\cdot\frac{\int_{0}^{\infty}L(r)r^{\frac{q}{2}}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}\\\
&\quad+\frac{h^{2}}{2}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})\cdot\frac{\int_{0}^{\infty}L(r)r^{\frac{q}{2}}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}\\\
&\quad+\frac{h^{2}}{q}\left[\Delta f(\bm{x})-\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\right]\bm{x}\cdot\frac{\int_{0}^{\infty}L^{\prime}(r)r^{\frac{q}{2}+1}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}\\\
&\quad+\frac{2h^{2}}{q}\sum_{i=1}^{q}\left(\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}\right)\bm{b}_{i}\cdot\frac{\int_{0}^{\infty}L^{\prime}(r)r^{\frac{q}{2}+1}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}+O(h^{2})+o(h^{2}),\end{split}$
(46)

which in turn shows that

$\displaystyle\mathbb{E}\left[\mathtt{grad}\,\widetilde{f}_{h}(\bm{x})\right]$
$\displaystyle=\mathbb{E}\left[\mathtt{Tang}\left(\nabla\widetilde{f}_{h}(\bm{x})\right)\right]=(I_{q+1}-\bm{x}\bm{x}^{T})\mathbb{E}\left[\nabla\widetilde{f}_{h}(\bm{x})\right]$
$\displaystyle=\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})+\frac{h^{2}}{2}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})\cdot\frac{\int_{0}^{\infty}L(r)r^{\frac{q}{2}}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}$
$\displaystyle\quad+\frac{2h^{2}}{q}\sum_{i=1}^{q}\left(\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}\right)\bm{b}_{i}\cdot\frac{\int_{0}^{\infty}L^{\prime}(r)r^{\frac{q}{2}+1}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}+O(h^{2})+o(h^{2})$
$\displaystyle=\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla f(\bm{x})+O(h^{2})$

as $h\to 0$. Result 1 thus follows and the Riemannian gradient estimator is
unbiased.

$\bullet$ Result 2 The covariance matrix of $\nabla\widetilde{f}_{h}(\bm{x})$
has the following asymptotic rate as $h\to 0$ and $nh^{q+2}\to\infty$:

$\text{Cov}\left[\nabla\widetilde{f}_{h}(\bm{x})\right]=\frac{1}{nh^{q+2}}\cdot
R(f,L)+o\left(\frac{1}{nh^{q+2}}\right),$

where
$R(f,L)=\frac{f(\bm{x})\int_{0}^{\infty}L^{\prime}(r)^{2}r^{\frac{q}{2}}dr}{2^{\frac{q}{2}-2}q\cdot\bar{\omega}_{q-1}\left(\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr\right)^{2}}\cdot(I_{q+1}-\bm{x}\bm{x}^{T})$.

_Derivation of Result 2_. By (46), the covariance matrix of
$\nabla\widetilde{f}_{h}(\bm{x})$ can be calculated as

$\displaystyle\text{Cov}\left[\nabla\widetilde{f}_{h}(\bm{x})\right]$
$\displaystyle=\frac{c_{h,q}(L)^{2}}{nh^{4}}\cdot\text{Cov}\left[(\bm{x}-\bm{X}_{1})L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{X}_{1}}{h^{2}}\right)\right]$
$\displaystyle=\frac{c_{h,q}(L)^{2}}{nh^{4}}\cdot\mathbb{E}\left[(\bm{x}-\bm{X}_{1})(\bm{x}-\bm{X}_{1})^{T}L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{X}_{1}}{h^{2}}\right)^{2}\right]-\frac{1}{n}\cdot\mathbb{E}\left[\nabla\widetilde{f}_{h}(\bm{x})\right]\cdot\mathbb{E}\left[\nabla\widetilde{f}_{h}(\bm{x})\right]^{T}$
$\displaystyle=\frac{c_{h,q}(L)^{2}}{nh^{4}}\int_{\Omega_{q}}(\bm{x}-\bm{y})(\bm{x}-\bm{y})^{T}L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{y}}{h^{2}}\right)^{2}f(\bm{y})\,\omega_{q}(d\bm{y})+O\left(\frac{1}{n}\right)$
$\displaystyle=\frac{c_{h,q}(L)^{2}}{n}h^{q-4}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\alpha_{\bm{x},\bm{\xi}}\alpha_{\bm{x},\bm{\xi}}^{T}L^{\prime}(r)^{2}f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\quad+O\left(\frac{1}{n}\right),$

where
$\alpha_{\bm{x},\bm{\xi}}=-rh^{2}\bm{x}+h\sqrt{r(2-h^{2}r)}\bm{B_{x}}\bm{\xi}$.
By condition (D1), the first-order Taylor’s expansion of $f$ at
$\bm{x}\in\Omega_{q}$ is

$f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})=f(\bm{x})+O(||\alpha_{\bm{x},\bm{\xi}}||_{2})=f(\bm{x})+O(h).$

Thus,

$\displaystyle\text{Cov}\left[\nabla\widetilde{f}_{h}(\bm{x})\right]$
$\displaystyle=\frac{c_{h,q}(L)^{2}}{n}h^{q-4}f(\bm{x})\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\Big{[}r^{2}h^{4}\bm{x}\bm{x}^{T}-rh^{3}\sqrt{r(2-h^{2}r)}\bm{x}(\bm{B_{x}}\bm{\xi})^{T}$
$\displaystyle\quad-
rh^{3}\sqrt{r(2-h^{2}r)}(\bm{B_{x}}\bm{\xi})\bm{x}^{T}+h^{2}r(2-h^{2}r)\bm{B_{x}}\bm{\xi}(\bm{B_{x}}\bm{\xi})^{T}\Big{]}L^{\prime}(r)^{2}r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\quad+o\left(\frac{1}{nh^{q+2}}\right)$
$\displaystyle\stackrel{{\scriptstyle\text{(i)}}}{{=}}\frac{c_{h,q}(L)^{2}}{n}h^{q}f(\bm{x})\bm{x}\bm{x}^{T}\bar{\omega}_{q-1}\int_{0}^{2h^{-2}}L^{\prime}(r)^{2}r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}-1}dr$
$\displaystyle\quad+\frac{c_{h,q}(L)^{2}}{n}h^{q-2}f(\bm{x})\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\xi_{i}\bm{b}_{i}\right)\left(\sum_{i=1}^{q}\xi_{i}\bm{b}_{i}^{T}\right)L^{\prime}(r)^{2}r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\quad+o\left(\frac{1}{nh^{q+2}}\right)$
$\displaystyle\stackrel{{\scriptstyle\text{(ii)}}}{{=}}\frac{c_{h,q}(L)^{2}}{n}h^{q}f(\bm{x})\bm{x}\bm{x}^{T}\bar{\omega}_{q-1}\int_{0}^{2h^{-2}}L^{\prime}(r)^{2}r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}-1}dr$
$\displaystyle\quad+\frac{c_{h,q}(L)^{2}\bar{\omega}_{q-1}}{nq}\cdot
h^{q-2}f(\bm{x})\int_{0}^{2h^{-2}}\left(\sum_{i=1}^{q}\bm{b}_{i}\bm{b}_{i}^{T}\right)L^{\prime}(r)^{2}r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}dr+o\left(\frac{1}{nh^{q+2}}\right)$
$\displaystyle\stackrel{{\scriptstyle\text{(iii)}}}{{=}}\frac{c_{h,q}(L)}{n}f(\bm{x})\bm{x}\bm{x}^{T}\cdot\frac{\int_{0}^{2h^{-2}}L^{\prime}(r)^{2}r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}-1}dr}{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}dr}$
$\displaystyle\quad+\frac{c_{h,q}(L)}{nq}\cdot
h^{q-2}f(\bm{x})\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\frac{\int_{0}^{2h^{-2}}L^{\prime}(r)^{2}r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}dr}{\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}dr}+o\left(\frac{1}{nh^{q+2}}\right)$
$\displaystyle\stackrel{{\scriptstyle\text{(iv)}}}{{=}}\frac{1}{nh^{q+2}}\cdot
R(f,L)+o\left(\frac{1}{nh^{q+2}}\right),$

where
$R(f,L)=\frac{f(\bm{x})\int_{0}^{\infty}L^{\prime}(r)^{2}r^{\frac{q}{2}}dr}{2^{\frac{q}{2}-2}q\cdot\bar{\omega}_{q-1}\left(\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr\right)^{2}}\cdot(I_{q+1}-\bm{x}\bm{x}^{T})$
is a matrix whose columns lie in the tangent space of $\Omega_{q}$ at
$\bm{x}$. During the derivation, we use (c) of Lemma 21 in (i) and (ii), plug
in the expression (3) of $c_{h,q}(L)$ in (iii), and take $h\to 0$ with
arguments in (a) of Lemma 21 and Remark 22 in (iv). Result 2 thus follows.

By the central limit theorem,

$\displaystyle\nabla\widetilde{f}_{h}(\bm{x})-\mathbb{E}\left[\nabla\widetilde{f}_{h}(\bm{x})\right]$
$\displaystyle=\text{Cov}\left[\nabla\widetilde{f}_{h}(\bm{x})\right]^{\frac{1}{2}}\cdot\text{Cov}\left[\nabla\widetilde{f}_{h}(\bm{x})\right]^{-\frac{1}{2}}\left\\{\nabla\widetilde{f}_{h}(\bm{x})-\mathbb{E}\left[\nabla\widetilde{f}_{h}(\bm{x})\right]\right\\}$
$\displaystyle=\left[\frac{1}{nh^{q+2}}\cdot
R(f,L)+o\left(\frac{1}{nh^{q+2}}\right)\right]^{\frac{1}{2}}\cdot\bm{Z}_{n}(\bm{x})$
$\displaystyle=O_{P}\left(\sqrt{\frac{1}{nh^{q+2}}}\right),$

where $\bm{Z}_{n}(\bm{x})\stackrel{{\scriptstyle
d}}{{\to}}N_{q+1}(\bm{0},I_{q+1})$.
The asymptotic rate for
$\mathtt{grad}\,\widetilde{f}_{h}(\bm{x})-\mathbb{E}\left[\mathtt{grad}\,\widetilde{f}_{h}(\bm{x})\right]=\mathtt{Tang}\left(\nabla\widetilde{f}_{h}(\bm{x})\right)-\mathbb{E}\left[\mathtt{Tang}\left(\nabla\widetilde{f}_{h}(\bm{x})\right)\right]$
remains unchanged, since the dominating constant $R(f,L)$ are within the
tangent space of $\Omega_{q}$ at $\bm{x}$.
In a nutshell, we conclude with bias (Result 1) and variance (Result 2)
estimation that

$\mathtt{grad}\,\widehat{f}_{h}(\bm{x})-\mathtt{grad}\,f(\bm{x})=\mathtt{Tang}\left(\nabla\widehat{f}_{h}(\bm{x})\right)-\mathtt{Tang}\left(\nabla
f(\bm{x})\right)=O(h^{2})+O_{P}\left(\sqrt{\frac{1}{nh^{q+2}}}\right)$

for any fixed $\bm{x}\in\Omega_{q}$, as $h\to 0$ and $nh^{q+2}\to\infty$.

Part B: Pointwise convergence rate of the Riemannian Hessian estimator
$\mathcal{H}\widehat{f}_{h}(\bm{x})$. As shown in Lemma 1,
$\mathcal{H}\widehat{f}_{h}(\bm{x})=\mathcal{H}\widetilde{f}_{h}(\bm{x})$ for
any $\bm{x}\in\Omega_{q}$ and we can establish the pointwise convergence rate
using either Riemannian Hessian estimator. Here, we stick to the Riemannian
Hessian estimator $\mathcal{H}\widehat{f}_{h}(\bm{x})$ in (38).

$\bullet$ Result 3. The expectation of the Riemannian Hessian estimator,
$\mathbb{E}\left[\mathcal{H}\widehat{f}_{h}(\bm{x})\right]$, has the following
asymptotic behavior as $h\to 0$:

$\displaystyle\mathbb{E}\left[\mathcal{H}\widehat{f}_{h}(\bm{x})\right]$
$\displaystyle=\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\mathbb{E}\left[\mathcal{A}\widehat{f}_{h}(\bm{x})\right]\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)$
$\displaystyle=\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\left[\nabla\nabla
f(\bm{x})-\bm{x}^{T}\nabla
f(\bm{x})\right]\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)+O(h^{2}),$

where
$\mathcal{A}\widehat{f}_{h}(\bm{x})=\frac{c_{h,q}(L)}{nh^{4}}\sum\limits_{i=1}^{n}\bm{X}_{i}\bm{X}_{i}^{T}L^{\prime\prime}\left(\frac{1-\bm{x}^{T}\bm{X}_{i}}{h^{2}}\right)+\frac{c_{h,q}(L)}{nh^{2}}\sum\limits_{i=1}^{n}\bm{x}^{T}\bm{X}_{i}I_{q+1}\cdot
L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{X}_{i}}{h^{2}}\right)$.

_Derivation of Result 3_. We first compute the expectation of
$\mathcal{A}\widehat{f}_{h}(\bm{x})$ and apply the left and right
multiplications of $\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)$ to simplify our
calculation. Notice that

$\displaystyle\mathbb{E}\left[\mathcal{A}\widehat{f}_{h}(\bm{x})\right]$
$\displaystyle=\frac{c_{h,q}(L)}{h^{4}}\int_{\Omega_{q}}\bm{y}\bm{y}^{T}L^{\prime\prime}\left(\frac{1-\bm{x}^{T}\bm{y}}{h^{2}}\right)f(\bm{y})\,\omega_{q}(d\bm{y})$
$\displaystyle\quad+\frac{c_{h,q}(L)}{h^{2}}\int_{\Omega_{q}}\bm{x}^{T}\bm{y}\cdot
I_{q+1}\cdot
L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{y}}{h^{2}}\right)f(\bm{y})\,\omega_{q}(d\bm{y})$
$\displaystyle=\frac{c_{h,q}(L)}{h^{4}}\int_{-1}^{1}\int_{\Omega_{q-1}}\left(t\bm{x}+\sqrt{1-t^{2}}\bm{B_{x}}\bm{\xi}\right)\left(t\bm{x}+\sqrt{1-t^{2}}\bm{B_{x}}\bm{\xi}\right)^{T}$
$\displaystyle\hskip 113.81102pt\times
L^{\prime\prime}\left(\frac{1-t}{h^{2}}\right)\cdot
f\left(t\bm{x}+\sqrt{1-t^{2}}\bm{B_{x}}\bm{\xi}\right)(1-t^{2})^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dt$
$\displaystyle\quad+\frac{c_{h,q}(L)}{h^{2}}\int_{-1}^{1}\int_{\Omega_{q-1}}tI_{q+1}\cdot
L^{\prime}\left(\frac{1-t}{h^{2}}\right)\cdot
f\left(t\bm{x}+\sqrt{1-t^{2}}\bm{B_{x}}\bm{\xi}\right)(1-t^{2})^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dt$
$\displaystyle=c_{h,q}(L)h^{q-4}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\left(\bm{x}+\alpha_{\bm{x},\bm{\xi}}\right)\left(\bm{x}+\alpha_{\bm{x},\bm{\xi}}\right)^{T}L^{\prime\prime}(r)$
$\displaystyle\hskip 170.71652pt\times
f\left(\bm{x}+\alpha_{\bm{x},\bm{\xi}}\right)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\quad+c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}(1-h^{2}r)I_{q+1}\cdot
L^{\prime}(r)f\left(\bm{x}+\alpha_{\bm{x},\bm{\xi}}\right)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr$

by $r=\frac{1-t}{h^{2}}$ and
$\alpha_{\bm{x},\bm{\xi}}=-rh^{2}\bm{x}+h\sqrt{r(2-h^{2}r)}\bm{B_{x}}\bm{\xi}$.
Since

$\displaystyle\left(\bm{x}+\alpha_{\bm{x},\bm{\xi}}\right)\left(\bm{x}+\alpha_{\bm{x},\bm{\xi}}\right)^{T}$
$\displaystyle=(1-rh^{2})^{2}\bm{x}\bm{x}^{T}+h(1-rh^{2})\sqrt{r(2-h^{2}r)}\left[\bm{x}(\bm{B_{x}}\bm{\xi})^{T}+(\bm{B_{x}}\bm{\xi})\bm{x}^{T}\right]$
$\displaystyle\quad+h^{2}r(2-h^{2}r)\cdot(\bm{B_{x}}\bm{\xi})(\bm{B_{x}}\bm{\xi})^{T},$

the preceding calculation proceeds as

$\displaystyle\begin{split}&\mathbb{E}\left[{\mathcal{A}}\widehat{f}_{h}(\bm{x})\right]\\\
&=c_{h,q}(L)h^{q-4}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\bm{x}\bm{x}^{T}L^{\prime\prime}(r)\cdot
f\left(\bm{x}+\alpha_{\bm{x},\bm{\xi}}\right)r^{\frac{q}{2}-1}(1-rh^{2})^{2}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr\\\
&\quad+c_{h,q}(L)h^{q-3}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\left[\bm{x}(\bm{B_{x}}\bm{\xi})^{T}+(\bm{B_{x}}\bm{\xi})\bm{x}^{T}\right]L^{\prime\prime}(r)\\\
&\hskip 142.26378pt\times
f\left(\bm{x}+\alpha_{\bm{x},\bm{\xi}}\right)r^{\frac{q-1}{2}}(1-rh^{2})(2-h^{2}r)^{\frac{q-1}{2}}\omega_{q-1}(d\bm{\xi})dr\\\
&\quad+c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}(\bm{B_{x}}\bm{\xi})(\bm{B_{x}}\bm{\xi})^{T}L^{\prime\prime}(r)\cdot
f\left(\bm{x}+\alpha_{\bm{x},\bm{\xi}}\right)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}\,\omega_{q-1}(d\bm{\xi})dr\\\
&\quad+c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}(1-h^{2}r)I_{q+1}\cdot
L^{\prime}(r)f\left(\bm{x}+\alpha_{\bm{x},\bm{\xi}}\right)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr\\\
&\equiv\text{(I) + (II) +(III) + (IV)}.\end{split}$ (47)

The above terms (I) and (II), after we apply the congruence operation

$\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\mathbb{E}\left[{\mathcal{A}}\widehat{f}_{h}(\bm{x})\right]\left(I_{q+1}-\bm{x}\bm{x}^{T}\right),$

yield zero. Hence, we will not continue to compute them. By condition (D1),
the Taylor’s expansion of $f$ at $\bm{x}$ is

$\displaystyle f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})$
$\displaystyle=f(\bm{x})+\alpha_{\bm{x},\bm{\xi}}^{T}\nabla
f(\bm{x})+\frac{1}{2}\alpha_{\bm{x},\bm{\xi}}^{T}\nabla\nabla
f(\bm{x})\alpha_{\bm{x},\bm{\xi}}+\frac{1}{6}\left(\sum_{i=1}^{q+1}(\alpha_{\bm{x},\bm{\xi}})_{i}\cdot\frac{\partial}{\partial
x_{i}}\right)^{3}f(\bm{x})+O\left(||\alpha_{\bm{x},\bm{\xi}}||_{2}^{4}\right)$
$\displaystyle=f(\bm{x})+\alpha_{\bm{x},\bm{\xi}}^{T}\nabla
f(\bm{x})+\frac{1}{2}\alpha_{\bm{x},\bm{\xi}}^{T}\nabla\nabla
f(\bm{x})\alpha_{\bm{x},\bm{\xi}}+\frac{1}{6}\left(\sum_{i=1}^{q+1}(\alpha_{\bm{x},\bm{\xi}})_{i}\cdot\frac{\partial}{\partial
x_{i}}\right)^{3}f(\bm{x})+O(h^{4}),$

where
$\left|\left|\alpha_{\bm{x},\bm{\xi}}\right|\right|_{2}^{2}=r^{2}h^{4}+h^{2}r(2-h^{2}r)=2rh^{2}$
by the orthogonality of $\bm{x}$ and columns of $\bm{B_{x}}$, and
$(\alpha_{\bm{x},\bm{\xi}})_{i}$ stands for the $i^{th}$ entry of the vector
$\alpha_{\bm{x},\bm{\xi}}$. Note that plugging $O(h^{4})$ into (III) or (IV)
in (47) both leads to a $O(h^{2})$ after integration, since with condition
(D2’) and our arguments in (a) of Lemma 21 and Remark 22,

$O(h^{4})\cdot
c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\phi(r,\bm{\xi})L^{\prime}(r)\,\omega_{q-1}(d\bm{\xi})dr\asymp
O(h^{2}),$

where $\phi(r,\bm{\xi})$ is a square integrable function of $(r,\bm{\xi})$ and
“$\asymp$” stands for the asymptotic equivalence. It shows that carrying out
the Taylor’s expansion of $f$ at $\bm{x}$ to the third order is sufficient in
our context.
More importantly, plugging the term
$\frac{1}{6}\left(\sum_{i=1}^{q+1}(\alpha_{\bm{x},\bm{\xi}})_{i}\cdot\frac{\partial}{\partial
x_{i}}\right)^{3}f(\bm{x})$ into (II) and (III) also gives rise to a
$O(h^{2})$ term. Because $\bm{B_{x}}\bm{\xi}=\sum_{i=1}^{q}\xi_{i}\bm{b}_{i}$
and
$\alpha_{\bm{x},\bm{\xi}}=-rh^{2}\bm{x}+h\sqrt{r(2-h^{2}r)}\bm{B_{x}}\bm{\xi}$,

$\displaystyle\begin{split}&\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)^{T}\left[\frac{1}{6}\left(\sum_{i=1}^{q+1}(\alpha_{\bm{x},\bm{\xi}})_{i}\cdot\frac{\partial}{\partial
x_{i}}\right)^{3}f(\bm{x})\right]\omega_{q-1}(d\bm{\xi})\\\
&=\int_{\Omega_{q-1}}\left[P_{r}(\bm{\xi},2)h^{6}+P_{r}(\bm{\xi},3)h^{5}+P_{r}(\bm{\xi},4)h^{4}+P_{r}(\bm{\xi},5)h^{3}\right]\omega_{q-1}(d\bm{\xi})\\\
&=O(h^{4})+\underbrace{\int_{\Omega_{q-1}}h^{3}P_{r}(\bm{\xi},5)\,\omega_{q-1}(d\bm{\xi})}_{=0},\end{split}$
(48)
$\displaystyle\begin{split}&\int_{\Omega_{q-1}}\left[\frac{1}{6}\left(\sum_{i=1}^{q+1}(\alpha_{\bm{x},\bm{\xi}})_{i}\cdot\frac{\partial}{\partial
x_{i}}\right)^{3}f(\bm{x})\right]\omega_{q-1}(d\bm{\xi})\\\
&=\int_{\Omega_{q-1}}\left[P_{r}(\bm{\xi},0)h^{6}+P_{r}(\bm{\xi},1)h^{5}+P_{r}(\bm{\xi},2)h^{4}+P_{r}(\bm{\xi},3)h^{3}\right]\omega_{q-1}(d\bm{\xi})\\\
&=O(h^{4})+\underbrace{\int_{\Omega_{q-1}}h^{3}P_{r}(\bm{\xi},3)\,\omega_{q-1}(d\bm{\xi})}_{=0},\end{split}$
(49)

where $P_{r}(\bm{\xi},n)$ is a polynomial of elements of
$\bm{\xi}=(\xi_{1},...,\xi_{q})$ with only degree $n$ terms, whose
coefficients may involve the variable $r$. The integral
$\int_{\Omega_{q-1}}h^{3}P_{r}(\bm{\xi},5)\,\omega_{q-1}(d\bm{\xi})=\int_{\Omega_{q-1}}h^{3}P_{r}(\bm{\xi},3)\,\omega_{q-1}(d\bm{\xi})=0$
is due to (c) of Lemma 21 and the fact that the integrand is a linear
combination of degree 5 or 3 monomials of elements of $\bm{\xi}$. With
condition (D2’) and our arguments in (a) of Lemma 21 and Remark 22, the final
$O(h^{4})$ terms in (48) and (49) both yield $O(h^{2})$ terms after being
plugged into (III) and (IV).
We now plug the Taylor’s expansion of $f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})$
back into (III) and obtain that

$\displaystyle\begin{split}&\text{Plug in (III) in \eqref{Hess_Expect1}}\\\
&=c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)^{T}L^{\prime\prime}(r)f(\bm{x})r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}\,\omega_{q-1}(d\bm{\xi})dr\\\
&\quad+c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)^{T}L^{\prime\prime}(r)\nabla
f(\bm{x})^{T}\alpha_{\bm{x},\bm{\xi}}\\\ &\hskip 227.62204pt\times
r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}\,\omega_{q-1}(d\bm{\xi})dr\\\
&\quad+c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)^{T}L^{\prime\prime}(r)\\\
&\hskip 113.81102pt\times\frac{1}{2}\alpha_{\bm{x},\bm{\xi}}^{T}\nabla\nabla
f(\bm{x})\alpha_{\bm{x},\bm{\xi}}\cdot
r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}\,\omega_{q-1}(d\bm{\xi})dr+O(h^{2})\\\
&=c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\frac{\bar{\omega}_{q-1}}{q}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)L^{\prime\prime}(r)f(\bm{x})r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}\,dr\\\
&\quad-
c_{h,q}(L)h^{q}\int_{0}^{2h^{-2}}\frac{\bar{\omega}_{q-1}}{q}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)L^{\prime\prime}(r)\nabla
f(\bm{x})^{T}\bm{x}r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}}\,dr\\\
&\quad+c_{h,q}(L)h^{q+2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)^{T}L^{\prime\prime}(r)\\\
&\hskip 113.81102pt\times\frac{1}{2}\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\cdot
r^{\frac{q}{2}+2}(2-h^{2}r)^{\frac{q}{2}}\,\omega_{q-1}(d\bm{\xi})dr\\\
&\quad-
c_{h,q}(L)h^{q+1}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)^{T}L^{\prime\prime}(r)\\\
&\hskip 113.81102pt\times\bm{x}^{T}\nabla\nabla
f(\bm{x})(\bm{B_{x}}\bm{\xi})\cdot
r^{\frac{q+3}{2}}(2-h^{2}r)^{\frac{q+1}{2}}\,\omega_{q-1}(d\bm{\xi})dr\\\
&\quad+c_{h,q}(L)h^{q}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)^{T}L^{\prime\prime}(r)\\\
&\hskip 113.81102pt\times\frac{1}{2}(\bm{B_{x}}\bm{\xi})^{T}\nabla\nabla
f(\bm{x})(\bm{B_{x}}\bm{\xi})\cdot
r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}+1}\,\omega_{q-1}(d\bm{\xi})dr+O(h^{2})\\\
&=c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\frac{\bar{\omega}_{q-1}}{q}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)L^{\prime\prime}(r)f(\bm{x})\cdot
r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}\,dr\\\ &\quad-
c_{h,q}(L)h^{q}\int_{0}^{2h^{-2}}\frac{\bar{\omega}_{q-1}}{q}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)L^{\prime\prime}(r)\nabla
f(\bm{x})^{T}\bm{x}\cdot
r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}}\,dr\quad-\quad 0\\\
&\quad+c_{h,q}(L)h^{q}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)^{T}L^{\prime\prime}(r)\\\
&\hskip 85.35826pt\times\frac{1}{2}(\bm{B_{x}}\bm{\xi})^{T}\nabla\nabla
f(\bm{x})(\bm{B_{x}}\bm{\xi})\cdot
r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}+1}\,\omega_{q-1}(d\bm{\xi})dr+O(h^{2}),\end{split}$
(50)

where we apply (c) of Lemma 21 and the fact that
$\sum_{i=1}^{q}\bm{b}_{i}\bm{b}_{i}^{T}=I_{q+1}-\bm{x}\bm{x}^{T}$. We also
absorb the third integral in the second equality into $O(h^{2})$ to obtain the
third equality, given condition (D2’) and our arguments in (a) of Lemma 21.
The “0” term in the third equality is due to the fact that

$\displaystyle\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)^{T}\bm{x}^{T}\nabla\nabla
f(\bm{x})(\bm{B_{x}}\bm{\xi})\,\omega_{q-1}(d\bm{\xi})=\int_{\Omega_{q-1}}P_{r}(\bm{\xi},3)\,\omega_{q-1}(d\bm{\xi})=0$

by (c) of Lemma 21, where the notation $P_{r}(\bm{\xi},3)$ is defined in (49).
Now, we consider the inner integral inside the last integration.

$\displaystyle\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)^{T}\cdot\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)^{T}\nabla\nabla
f(\bm{x})\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)\omega_{q-1}(d\bm{\xi})$
$\displaystyle=\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\bm{b}_{i}\bm{b}_{i}^{T}\xi_{i}^{2}+\sum_{i\neq
j}\bm{b}_{i}\bm{b}_{j}^{T}\xi_{i}\xi_{j}\right)\left(\sum_{i=1}^{q}\bm{b}_{i}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}\xi_{i}^{2}+\sum_{i\neq j}\bm{b}_{i}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{j}\xi_{i}\xi_{j}\right)\omega_{q-1}(d\bm{\xi})$
$\displaystyle\stackrel{{\scriptstyle(*)}}{{=}}\int_{\Omega_{q-1}}\Bigg{[}\left(\sum_{i=1}^{q}\bm{b}_{i}\bm{b}_{i}^{T}\xi_{i}^{2}\right)\left(\sum_{i=1}^{q}\bm{b}_{i}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}\xi_{i}^{2}\right)$ $\displaystyle\quad+\left(\sum_{i\neq
j}\bm{b}_{i}\bm{b}_{j}^{T}\xi_{i}\xi_{j}\right)\left(\sum_{i\neq
j}\bm{b}_{i}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{j}\xi_{i}\xi_{j}\right)\Bigg{]}\omega_{q-1}(d\bm{\xi})$
$\displaystyle=\int_{\Omega_{q-1}}\Bigg{[}\sum_{i=1}^{q}\bm{b}_{i}\bm{b}_{i}^{T}\cdot\bm{b}_{i}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}\xi_{i}^{4}+\sum_{i\neq
j}\bm{b}_{i}\bm{b}_{i}^{T}\cdot\bm{b}_{j}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{j}\xi_{i}^{2}\xi_{j}^{2}$ $\displaystyle\quad+2\sum_{i\neq
j}\bm{b}_{i}\bm{b}_{j}^{T}\cdot\bm{b}_{i}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{j}\xi_{i}^{2}\xi_{j}^{2}\Bigg{]}\omega_{q-1}(d\bm{\xi}),$

where the cross product terms vanish after integration by Lemma 21 and the
factor 2 in front of the last summation emerges because any fixed $(i,j)$ term
($i\neq j$) in the first factor of the second product in equality $(*)$ can be
matched up with both $(i,j)$ and $(j,i)$ terms in the second factor to yield a
summand. Using Lemma 21 and the facts that
$\sum_{i=1}^{q}\bm{b}_{i}\bm{b}_{i}^{T}=I_{q+1}-\bm{x}\bm{x}^{T}$ and

$\sum_{i=1}^{q}\bm{b}_{i}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}=\text{tr}\left[\nabla\nabla
f(\bm{x})\sum_{i=1}^{q}\bm{b}_{i}\bm{b}_{i}^{T}\right]=\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla f(\bm{x})\bm{x},$

the preceding display continues as

$\displaystyle\int_{\Omega_{q-1}}\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)^{T}\cdot\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)^{T}\nabla\nabla
f(\bm{x})\left(\sum_{i=1}^{q}\bm{b}_{i}\xi_{i}\right)\omega_{q-1}(d\bm{\xi})$
$\displaystyle=\frac{3\bar{\omega}_{q-1}}{q(q+2)}\left(\sum_{i=1}^{q}\bm{b}_{i}\bm{b}_{i}^{T}\cdot\bm{b}_{i}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}\right)+\frac{\bar{\omega}_{q-1}}{q(q+2)}\left(\sum_{i\neq
j}\bm{b}_{i}\bm{b}_{i}^{T}\cdot\bm{b}_{j}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{j}\right)$
$\displaystyle\quad+\frac{2\bar{\omega}_{q-1}}{q(q+2)}\left(\sum_{i\neq
j}\bm{b}_{i}\bm{b}_{j}^{T}\cdot\bm{b}_{i}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{j}\right)$
$\displaystyle=\frac{\bar{\omega}_{q-1}}{q(q+2)}\left[\sum_{i=1}^{q}\bm{b}_{i}\bm{b}_{i}^{T}\left(\sum_{j=1}^{q}\bm{b}_{j}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{j}\right)\right]$
$\displaystyle\quad+\frac{2\bar{\omega}_{q-1}}{q(q+2)}\left[\sum_{i=1}^{q}\bm{b}_{i}\left(\bm{b}_{i}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{i}\right)\bm{b}_{i}^{T}+\sum_{i\neq
j}\bm{b}_{i}\left(\bm{b}_{i}^{T}\nabla\nabla
f(\bm{x})\bm{b}_{j}\right)\bm{b}_{j}^{T}\right]$
$\displaystyle=\frac{\bar{\omega}_{q-1}}{q(q+2)}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\left[\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla f(\bm{x})\bm{x}\right]$
$\displaystyle\quad+\frac{2\bar{\omega}_{q-1}}{q(q+2)}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla\nabla
f(\bm{x})\left(I_{q+1}-\bm{x}\bm{x}^{T}\right).$

Plugging this result back into (50) and conduct some integration by parts, we
obtain that

Plug in (III) in (47)
$\displaystyle=c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\frac{\bar{\omega}_{q-1}}{q}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)f(\bm{x})L^{\prime\prime}(r)\cdot
r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}\,dr$ $\displaystyle\quad-
c_{h,q}(L)h^{q}\int_{0}^{2h^{-2}}\frac{\bar{\omega}_{q-1}}{q}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})^{T}\bm{x}L^{\prime\prime}(r)\cdot
r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}}\,dr$
$\displaystyle\quad+c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{2q(q+2)}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\left[\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\right]L^{\prime\prime}(r)r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}+1}dr$
$\displaystyle\quad+c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{q(q+2)}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla\nabla
f(\bm{x})\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)L^{\prime\prime}(r)r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}+1}dr$
$\displaystyle\quad+O(h^{2})$
$\displaystyle=c_{h,q}(L)h^{q-2}\cdot\frac{\bar{\omega}_{q-1}}{q}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)L^{\prime}(r)f(\bm{x})\cdot
r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}\Big{|}_{0}^{2h^{-2}}$
$\displaystyle\quad-
c_{h,q}(L)h^{q-2}\cdot\frac{\bar{\omega}_{q-1}}{2}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)f(\bm{x})L^{\prime}(r)\left[r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}}-h^{2}r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-1}\right]dr$
$\displaystyle\quad-
c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{q}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})^{T}\bm{x}L^{\prime}(r)\cdot
r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}}\Big{|}_{0}^{2h^{-2}}$
$\displaystyle\quad+c_{h,q}(L)h^{q}\int_{0}^{2h^{-2}}\frac{\bar{\omega}_{q-1}}{q}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})^{T}\bm{x}L^{\prime}(r)$ $\displaystyle\hskip
142.26378pt\times\left[\left(\frac{q+2}{2}\right)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}-\frac{qh^{2}}{2}r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}-1}\right]dr$
$\displaystyle\quad+c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{2q(q+2)}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\left[\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\right]L^{\prime}(r)r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}+1}\Big{|}_{0}^{2h^{-2}}$
$\displaystyle\quad-
c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{4q}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\left[\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla f(\bm{x})\bm{x}\right]L^{\prime}(r)$
$\displaystyle\hskip
142.26378pt\times\left[r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}+1}-h^{2}r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}}\right]dr$
$\displaystyle\quad+c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{q(q+2)}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla\nabla
f(\bm{x})\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)L^{\prime}(r)r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}+1}\Big{|}_{0}^{2h^{-2}}$
$\displaystyle\quad-
c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{2q}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla\nabla
f(\bm{x})\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)L^{\prime}(r)$
$\displaystyle\hskip
142.26378pt\times\left[r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}+1}-h^{2}r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}}\right]dr$
$\displaystyle\quad+O(h^{2})$
$\displaystyle=-c_{h,q}(L)h^{q-2}\bar{\omega}_{q-1}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)f(\bm{x})L^{\prime}(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}(1-h^{2}r)dr$
$\displaystyle\quad+c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{2}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})^{T}\bm{x}\cdot
L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}dr$
$\displaystyle\quad+c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{q}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})^{T}\bm{x}\cdot
L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}dr$ $\displaystyle\quad-
c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{4q}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\left[\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\right]L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}+1}dr$
$\displaystyle\quad-
c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{2q}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla\nabla
f(\bm{x})\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}+1}dr$
$\displaystyle\quad+O(h^{2}),$

where we use the fact that

$c_{h,q}(L)h^{q+2}\int_{0}^{2h^{-2}}L^{\prime}(r)r^{j}(2-h^{2}r)^{k}dr=O(h^{2})$
(51)

for any $k,j>0$ via Remark 22. With extra integration by parts on the third
and fifth term in the preceding display, we obtain that

$\displaystyle\begin{split}&\text{Plug in (III) in \eqref{Hess_Expect1}}\\\
&=-c_{h,q}(L)h^{q-2}\bar{\omega}_{q-1}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)f(\bm{x})L^{\prime}(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}(1-h^{2}r)dr\\\
&\quad+c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{2}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})^{T}\bm{x}L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}dr\\\
&\quad-
c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{2}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})^{T}\bm{x}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}}dr\\\ &\quad-
c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{4q}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\left[\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\right]L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}+1}dr\\\
&\quad+c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{4}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla\nabla
f(\bm{x})\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}+1}dr\\\
&\quad+O(h^{2}).\end{split}$ (52)

We now plug the Taylor’s expansion of $f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})$ at
$\bm{x}$ back into (IV) in (47) and deduce that

$\displaystyle\begin{split}&\text{Plug in (IV) in \eqref{Hess_Expect1}}\\\
&=c_{h,q}(L)h^{q-2}\bar{\omega}_{q-1}\int_{0}^{2h^{-2}}I_{q+1}f(\bm{x})L^{\prime}(r)\cdot
r^{\frac{q}{2}-1}(1-h^{2}r)(2-h^{2}r)^{\frac{q}{2}-1}dr\\\
&\quad+c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}I_{q+1}\nabla
f(\bm{x})^{T}\alpha_{\bm{x},\bm{\xi}}L^{\prime}(r)\cdot
r^{\frac{q}{2}-1}(1-h^{2}r)(2-h^{2}r)^{\frac{q}{2}-1}\,\omega_{q-1}(d\bm{\xi})dr\\\
&\quad+c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}I_{q+1}\cdot\frac{1}{2}\alpha_{\bm{x},\bm{\xi}}^{T}\nabla\nabla
f(\bm{x})\alpha_{\bm{x},\bm{\xi}}r^{\frac{q}{2}-1}(1-h^{2}r)(2-h^{2}r)^{\frac{q}{2}-1}\,\omega_{q-1}(d\bm{\xi})dr\\\
&\quad+O(h^{2})\\\
&=c_{h,q}(L)h^{q-2}\bar{\omega}_{q-1}\int_{0}^{2h^{-2}}I_{q+1}f(\bm{x})L^{\prime}(r)\cdot
r^{\frac{q}{2}-1}(1-h^{2}r)(2-h^{2}r)^{\frac{q}{2}-1}dr\\\ &\quad-
c_{h,q}(L)h^{q}\bar{\omega}_{q-1}\int_{0}^{2h^{-2}}I_{q+1}\nabla
f(\bm{x})^{T}\bm{x}L^{\prime}(r)r^{\frac{q}{2}}(1-h^{2}r)(2-h^{2}r)^{\frac{q}{2}-1}dr\\\
&\quad+c_{h,q}(L)h^{q-1}\int_{0}^{2h^{-2}}\underbrace{\int_{\Omega_{q-1}}I_{q+1}\cdot\nabla
f(\bm{x})^{T}(\bm{B_{x}}\bm{\xi})r^{\frac{q-1}{2}}(2-h^{2}r)^{\frac{q-1}{2}}(1-h^{2}r)\,\omega_{q-1}(d\bm{\xi})}_{=0}dr\\\
&\quad+c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}I_{q+1}\cdot\frac{1}{2}(\bm{B_{x}}\bm{\xi})^{T}\nabla\nabla
f(\bm{x})(\bm{B_{x}}\bm{\xi})L^{\prime}(r)r^{\frac{q}{2}}(1-h^{2}r)(2-h^{2}r)^{\frac{q}{2}}\\\
&\quad+O(h^{2})\\\
&=c_{h,q}(L)h^{q-2}\bar{\omega}_{q-1}\int_{0}^{2h^{-2}}I_{q+1}f(\bm{x})L^{\prime}(r)\cdot
r^{\frac{q}{2}-1}(1-h^{2}r)(2-h^{2}r)^{\frac{q}{2}-1}dr\\\ &\quad-
c_{h,q}(L)h^{q}\bar{\omega}_{q-1}\int_{0}^{2h^{-2}}I_{q+1}\nabla
f(\bm{x})^{T}\bm{x}L^{\prime}(r)r^{\frac{q}{2}}(1-h^{2}r)(2-h^{2}r)^{\frac{q}{2}-1}dr\\\
&\quad+c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{2q}\int_{0}^{2h^{-2}}I_{q+1}\left[\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\right]L^{\prime}(r)r^{\frac{q}{2}}(1-h^{2}r)(2-h^{2}r)^{\frac{q}{2}}dr\\\
&\quad+O(h^{2}),\end{split}$ (53)

where we expand
$\alpha_{\bm{x},\bm{\xi}}=-rh^{2}\bm{x}+h\sqrt{r(2-h^{2}r)}\bm{B_{x}}\bm{\xi}$,
absorb $O(h^{2})$ terms via (51), make use of (c) in Lemma 21, and leverage
our argument in (44). Combining (47), (52), and (53), we conclude that

$\displaystyle\mathbb{E}\left[\mathcal{H}\widehat{f}_{h}(\bm{x})\right]$
$\displaystyle=\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\mathbb{E}\left[{\mathcal{A}}\widehat{f}_{h}(\bm{x})\right]\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)$
$\displaystyle=\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\cdot\text{(III)}\cdot\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)+\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\cdot\text{(IV)}\cdot\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)$
$\displaystyle=-c_{h,q}(L)h^{q-2}\bar{\omega}_{q-1}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)f(\bm{x})L^{\prime}(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}(1-h^{2}r)dr$
$\displaystyle\quad+c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{2}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})^{T}\bm{x}L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}dr$
$\displaystyle\quad-
c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{2}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})^{T}\bm{x}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}}dr$
$\displaystyle\quad-
c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{4q}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\left[\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla
f(\bm{x})\bm{x}\right]L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}+1}dr$
$\displaystyle\quad+c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{4}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla\nabla
f(\bm{x})\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}+1}dr$
$\displaystyle\quad+c_{h,q}(L)h^{q-2}\cdot\bar{\omega}_{q-1}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)f(\bm{x})L^{\prime}(r)\cdot
r^{\frac{q}{2}-1}(1-h^{2}r)(2-h^{2}r)^{\frac{q}{2}-1}dr$ $\displaystyle\quad-
c_{h,q}(L)h^{q}\cdot\bar{\omega}_{q-1}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})^{T}\bm{x}L^{\prime}(r)r^{\frac{q}{2}}(1-h^{2}r)(2-h^{2}r)^{\frac{q}{2}-1}dr$
$\displaystyle\quad+c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{2q}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\left[\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla f(\bm{x})\bm{x}\right]$ $\displaystyle\hskip
170.71652pt\times
L^{\prime}(r)r^{\frac{q}{2}}(1-h^{2}r)(2-h^{2}r)^{\frac{q}{2}}dr$
$\displaystyle\quad+O(h^{2})$
$\displaystyle\stackrel{{\scriptstyle(**)}}{{=}}0+c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{2}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})^{T}\bm{x}$ $\displaystyle\hskip 142.26378pt\times
L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-1}(2-h^{2}r-2+2rh^{2})dr$
$\displaystyle\quad-\lambda_{h,q}(L)^{-1}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})^{T}\bm{x}\cdot\frac{\bar{\omega}_{q-1}}{2}\int_{0}^{2h^{-2}}L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}}dr$
$\displaystyle\quad+c_{h,q}(L)h^{q}\cdot\frac{\bar{\omega}_{q-1}}{4q}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\left[\Delta
f(\bm{x})-\bm{x}^{T}\nabla\nabla f(\bm{x})\bm{x}\right]$ $\displaystyle\hskip
142.26378pt\times\int_{0}^{2h^{-2}}L^{\prime}(r)r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}}(2-2h^{2}r-2+h^{2}r)dr$
$\displaystyle\quad+\lambda_{h,q}(L)^{-1}\cdot\frac{\bar{\omega}_{q-1}}{4}\int_{0}^{2h^{-2}}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla\nabla
f(\bm{x})\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)$ $\displaystyle\hskip
142.26378pt\times L(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}+1}dr$
$\displaystyle\quad+O(h^{2})$
$\displaystyle\stackrel{{\scriptstyle(***)}}{{=}}O(h^{2})-\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})^{T}\bm{x}\cdot\frac{\bar{\omega}_{q-1}}{2}\lambda_{q}(L)^{-1}\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}2^{\frac{q}{2}}dr+O(h^{2})$
$\displaystyle\quad+\lambda_{q}(L)^{-1}\cdot\frac{\bar{\omega}_{q-1}}{4}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla\nabla
f(\bm{x})\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}2^{\frac{q}{2}+1}dr+O(h^{2})$
$\displaystyle=-\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla
f(\bm{x})^{T}\bm{x}+\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla\nabla
f(\bm{x})\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)+O(h^{2}),$

where the first term matches up with the sixth term, the second term with the
seventh term, the fourth term with the eighth term, and (3) is applied to the
rest terms when $h\to 0$ in $(**)$. In addition, we leverage the asymptotic
rates (51) and (42) as well as recall that
$\lambda_{q}(L)=2^{\frac{q}{2}-1}\bar{\omega}_{q-1}\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr$
from (a) of Lemma 21 in $(***)$. Result 3 thus follows. It implies that the
bias
$\mathbb{E}\left[\mathcal{H}\widehat{f}_{h}(\bm{x})\right]-\mathcal{H}f(\bm{x})$
is of the rate $O(h^{2})$ and the Riemannian Hessian estimator is
asymptotically unbiased.

Now, we proceed to bound

$\mathcal{H}\widehat{f}_{h}(\bm{x})-\mathbb{E}\left[\mathcal{H}\widehat{f}_{h}(\bm{x})\right]=\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\left[{\mathcal{A}}\widehat{f}_{h}(\bm{x})-\mathbb{E}\left({\mathcal{A}}\widehat{f}_{h}(\bm{x})\right)\right]\left(I_{q+1}-\bm{x}\bm{x}^{T}\right).$

$\bullet$ Result 4. The covariance matrix
$\text{Cov}\left[\mathtt{vec}\left(\mathcal{H}\widehat{f}_{h}(\bm{x})\right)\right]$
has the following asymptotic rate as $h\to 0$ and $nh^{q+4}\to\infty$:

$\text{Cov}\left[\mathtt{vec}\left(\mathcal{H}\widehat{f}_{h}(\bm{x})\right)\right]=O\left(\frac{1}{nh^{q+4}}\right),$

where we define the matrix $\mathtt{vec}$ operator, which converts a matrix
into a vector by stacking the columns. That is, given a matrix
$A\in\mathbb{R}^{m\times n}$, $\mathtt{vec}(A)$ is a vector of length $mn$.

_Derivation of Result 4_. We first calculate the covariance matrix of
$\mathtt{vec}\left({\mathcal{A}}\widehat{f}_{h}(\bm{x})\right)$ as

$\displaystyle\text{Cov}\left[\mathtt{vec}\left({\mathcal{A}}\widehat{f}_{h}(\bm{x})\right)\right]$
$\displaystyle=\frac{c_{h,q}(L)^{2}}{nh^{8}}\cdot\mathbb{E}\left[\mathtt{vec}(\bm{X}_{1}\bm{X}_{1}^{T})\cdot\mathtt{vec}(\bm{X}_{1}\bm{X}_{1}^{T})^{T}L^{\prime\prime}\left(\frac{1-\bm{x}^{T}\bm{X}_{i}}{h^{2}}\right)^{2}\right]$
$\displaystyle\quad+\frac{2c_{h,q}(L)^{2}}{nh^{6}}\cdot\mathbb{E}\left[\mathtt{vec}(\bm{X}_{i}\bm{X}_{i}^{T})\cdot\mathtt{vec}(\bm{x}^{T}\bm{X}_{1}I_{q+1})^{T}L^{\prime\prime}\left(\frac{1-\bm{x}^{T}\bm{X}_{i}}{h^{2}}\right)L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{X}_{i}}{h^{2}}\right)\right]$
$\displaystyle\quad+\frac{c_{h,q}(L)^{2}}{nh^{4}}\cdot\mathbb{E}\left[\mathtt{vec}(\bm{x}^{T}\bm{X}_{1}I_{q+1})\cdot\mathtt{vec}(\bm{x}^{T}\bm{X}_{1}I_{q+1})^{T}L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{X}_{i}}{h^{2}}\right)^{2}\right]$
$\displaystyle\quad-\frac{1}{n}\cdot\mathbb{E}\left[\mathtt{vec}\left({\mathcal{A}}\widehat{f}_{h}(\bm{x})\right)\right]\mathbb{E}\left[\mathtt{vec}\left({\mathcal{A}}\widehat{f}_{h}(\bm{x})\right)\right]^{T}$
$\displaystyle=\frac{c_{h,q}(L)^{2}}{nh^{8}}\int_{\Omega_{q}}\mathtt{vec}(\bm{y}\bm{y}^{T})\cdot\mathtt{vec}(\bm{y}\bm{y}^{T})^{T}L^{\prime\prime}\left(\frac{1-\bm{x}^{T}\bm{y}}{h^{2}}\right)^{2}f(\bm{y})\,\omega_{q}(d\bm{y})$
$\displaystyle\quad+\frac{2c_{h,q}(L)^{2}}{nh^{6}}\int_{\Omega_{q}}\mathtt{vec}(\bm{y}\bm{y}^{T})\cdot\mathtt{vec}(\bm{x}^{T}\bm{y}I_{q+1})^{T}L^{\prime\prime}\left(\frac{1-\bm{x}^{T}\bm{y}}{h^{2}}\right)L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{y}}{h^{2}}\right)f(\bm{y})\,\omega_{q}(d\bm{y})$
$\displaystyle\quad+\frac{c_{h,q}(L)^{2}}{nh^{4}}\int_{\Omega_{q}}\mathtt{vec}(\bm{x}^{T}\bm{y}I_{q+1})\cdot\mathtt{vec}(\bm{x}^{T}\bm{y}I_{q+1})^{T}L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{y}}{h^{2}}\right)^{2}f(\bm{y})\,\omega_{q}(d\bm{y})$
$\displaystyle\quad-\frac{1}{n}\cdot\mathbb{E}\left[\mathtt{vec}\left({\mathcal{A}}\widehat{f}_{h}(\bm{x})\right)\right]\mathbb{E}\left[\mathtt{vec}\left({\mathcal{A}}\widehat{f}_{h}(\bm{x})\right)\right]^{T}$
$\displaystyle=\frac{c_{h,q}(L)^{2}h^{q}}{nh^{8}}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\mathtt{vec}\left[(\bm{x}+\alpha_{\bm{x},\bm{\xi}})(\bm{x}+\alpha_{\bm{x},\bm{\xi}})^{T}\right]\mathtt{vec}\left[(\bm{x}+\alpha_{\bm{x},\bm{\xi}})(\bm{x}+\alpha_{\bm{x},\bm{\xi}})^{T}\right]^{T}$
$\displaystyle\hskip 142.26378pt\times
L^{\prime\prime}(r)^{2}f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})\cdot
r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\,\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\quad+\frac{2c_{h,q}(L)^{2}h^{q}}{nh^{6}}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\mathtt{vec}\left[(\bm{x}+\alpha_{\bm{x},\bm{\xi}})(\bm{x}+\alpha_{\bm{x},\bm{\xi}})^{T}\right]\mathtt{vec}\left[I_{q+1}(1-rh^{2})\right]^{T}$
$\displaystyle\hskip 142.26378pt\times
L^{\prime\prime}(r)L^{\prime}(r)f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})\cdot
r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\,\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\quad+\frac{c_{h,q}(L)^{2}h^{q}}{nh^{4}}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\mathtt{vec}\left[I_{q+1}(1-rh^{2})\right]\mathtt{vec}\left[I_{q+1}(1-rh^{2})\right]^{T}$
$\displaystyle\hskip 142.26378pt\times
L^{\prime}(r)^{2}f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})\cdot
r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\,\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\quad-\frac{1}{n}\cdot\mathbb{E}\left[\mathtt{vec}\left({\mathcal{A}}\widehat{f}_{h}(\bm{x})\right)\right]\mathbb{E}\left[\mathtt{vec}\left({\mathcal{A}}\widehat{f}_{h}(\bm{x})\right)\right]^{T},$

where
$\alpha_{\bm{x},\bm{\xi}}=-rh^{2}\bm{x}+h\sqrt{r(2-h^{2}r)}\bm{B_{x}}\bm{\xi}$.
Note that by condition (D1), the first-order Taylor’s expansion of $f$ at
$\bm{x}\in\Omega_{q}$ is

$f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})=f(\bm{x})+O(||\alpha_{\bm{x},\bm{\xi}}||_{2})=f(\bm{x})+O(h).$

In addition,

$\displaystyle(\bm{x}+\alpha_{\bm{x},\bm{\xi}})(\bm{x}+\alpha_{\bm{x},\bm{\xi}})^{T}$
$\displaystyle=(1-rh^{2})^{2}\bm{x}\bm{x}^{T}+h(1-r^{2}h)\sqrt{r(2-h^{2}r)}\left[\bm{x}(\bm{B_{x}}\bm{\xi})^{T}+(\bm{B_{x}})\bm{x}^{T}\right]$
$\displaystyle\quad+h^{2}r(2-h^{2}r)(\bm{B_{x}}\bm{\xi})(\bm{B_{x}}\bm{\xi})^{T}.$

Moreover, when the congruence operation
$\left(I_{q+1}-\bm{x}\bm{x}^{T}\right){\mathcal{A}}\widehat{f}_{h}(\bm{x})\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)$
is introduced, it will be applied inside the $\mathtt{vec}$ operation. Thus,
after applying the congruence operation,

$\displaystyle\mathtt{vec}\left[\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)(\bm{x}+\alpha_{\bm{x},\bm{\xi}})(\bm{x}+\alpha_{\bm{x},\bm{\xi}})^{T}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\right]$
$\displaystyle\quad\times\mathtt{vec}\left[\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)(\bm{x}+\alpha_{\bm{x},\bm{\xi}})(\bm{x}+\alpha_{\bm{x},\bm{\xi}})^{T}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\right]^{T}=O(h^{4})$

and

$\displaystyle\mathtt{vec}\left[\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)(\bm{x}+\alpha_{\bm{x},\bm{\xi}})(\bm{x}+\alpha_{\bm{x},\bm{\xi}})^{T}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\right]$
$\displaystyle\quad\times\mathtt{vec}\left[\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)(1-rh^{2})\right]^{T}=O(h^{2}).$

Together with condition (D2’), (3), (51), and the bias bound
$\mathbb{E}\left[\mathtt{vec}\left(\mathcal{H}\widehat{f}_{h}(\bm{x})\right)\right]=\mathcal{H}f(\bm{x})+O(h^{2})$,
we conclude that

$\text{Cov}\left[\mathtt{vec}\left(\mathcal{H}\widehat{f}_{h}(\bm{x})\right)\right]=O\left(\frac{1}{nh^{q+4}}\right).$

Result 4 is thus proved. Finally, by the central limit theorem,

$\displaystyle\mathtt{vec}\left\\{\mathcal{H}\widehat{f}_{h}(\bm{x})-\mathbb{E}\left[\mathcal{H}\widehat{f}_{h}(\bm{x})\right]\right\\}$
$\displaystyle=\text{Cov}\left[\mathtt{vec}\left(\mathcal{H}\widehat{f}_{h}(\bm{x})\right)\right]^{\frac{1}{2}}\text{Cov}\left[\mathtt{vec}\left(\mathcal{H}\widehat{f}_{h}(\bm{x})\right)\right]^{-\frac{1}{2}}\mathtt{vec}\left\\{\mathcal{H}\widehat{f}_{h}(\bm{x})-\mathbb{E}\left[\mathcal{H}\widehat{f}_{h}(\bm{x})\right]\right\\}$
$\displaystyle=O\left(\sqrt{\frac{1}{nh^{q+4}}}\right)\cdot\widetilde{\bm{Z}}_{n}(\bm{x})$
$\displaystyle=O_{P}\left(\sqrt{\frac{1}{nh^{q+4}}}\right),$

where $\widetilde{\bm{Z}}_{n}(\bm{x})\stackrel{{\scriptstyle
d}}{{\to}}N_{(q+1)^{2}}\left(\bm{0},I_{(q+1)^{2}}\right)$. In total, we
conclude with our bias and stochastic variation bounds that

$\mathcal{H}\widehat{f}_{h}(\bm{x})-\mathcal{H}f(\bm{x})=O(h^{2})+O_{P}\left(\sqrt{\frac{1}{nh^{q+4}}}\right)$

for any fixed $\bm{x}\in\Omega_{q}$ as $h\to 0$ and $nh^{q+4}\to\infty$, where

$\mathcal{H}f(\bm{x})=\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)\nabla\nabla
f(\bm{x})\left(I_{q+1}-\bm{x}\bm{x}^{T}\right)-\nabla
f(\bm{x})^{T}\bm{x}\left(I_{q+1}-\bm{x}\bm{x}^{T}\right).$

###  D.3 Proof of Theorem 4

###### Theorem 4

Assume (D1), (D2’), and (K1). The uniform convergence rate of
$\widehat{f}_{h}$ is given by

$\sup_{\bm{x}\in\Omega_{q}}|\widehat{f}_{h}(\bm{x})-f(\bm{x})|=O(h^{2})+O_{P}\left(\sqrt{\frac{|\log
h|}{nh^{q}}}\right)$

as $h\to 0$ and $\frac{nh^{q}}{|\log h|}\to\infty$.
Furthermore, the uniform convergence rate of
$\mathtt{grad}\,\widehat{f}_{h}(\bm{x})$ on $\Omega_{q}$ is

$\sup_{\bm{x}\in\Omega_{q}}\left|\left|\mathtt{grad}\,\widehat{f}_{h}(\bm{x})-\mathtt{grad}\,f(\bm{x})\right|\right|_{\max}=O(h^{2})+O_{P}\left(\sqrt{\frac{|\log
h|}{nh^{q+2}}}\right),$

as $h\to 0$ and $\frac{nh^{q+2}}{|\log h|}\to\infty$. Finally, the uniform
convergence rate of $\mathcal{H}\widehat{f}_{h}(\bm{x})$ on $\Omega_{q}$ is

$\sup_{\bm{x}\in\Omega_{q}}\left|\left|\mathcal{H}\widehat{f}_{h}(\bm{x})-\mathcal{H}f(\bm{x})\right|\right|_{\max}=O(h^{2})+O_{P}\left(\sqrt{\frac{|\log
h|}{nh^{q+4}}}\right),$

as $h\to 0$ and $\frac{nh^{q+4}}{|\log h|}\to\infty$, where
$\left|\left|\cdot\right|\right|_{\max}$ is the elementwise maximum norm for a
vector in $\mathbb{R}^{q+1}$ or a matrix in $\mathbb{R}^{(q+1)\times(q+1)}$.

Proof  Note that with the directional KDE form (15), we have that

$D^{\tau_{j}}\widehat{f}_{h}(\bm{x})=\frac{c_{h,q}(L)}{nh}\sum_{i=1}^{n}\left(\frac{x_{\tau_{j}}-X_{\tau_{j}}}{h}\right)L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{X}_{i}}{h^{2}}\right),$
$D^{[\tau_{j},\tau_{k}]}\widehat{f}_{h}(\bm{x})=\begin{cases}\frac{c_{h,q}(L)}{nh^{2}}\sum_{i=1}^{n}\left(\frac{x_{\tau_{j}}-X_{\tau_{j}}}{h}\right)\left(\frac{x_{\tau_{k}}-X_{\tau_{k}}}{h}\right)L^{\prime\prime}\left(\frac{1-\bm{x}^{T}\bm{X}_{i}}{h^{2}}\right)&j\neq
k,\\\
\frac{c_{h,q}(L)}{nh^{2}}\sum_{i=1}^{n}\left(\frac{x_{\tau_{j}}-X_{\tau_{j}}}{h}\right)^{2}L^{\prime\prime}\left(\frac{1-\bm{x}^{T}\bm{X}_{i}}{h^{2}}\right)+\frac{c_{h,q}(L)}{nh^{2}}\sum_{i=1}^{n}L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{X}_{i}}{h^{2}}\right)&j=k,\end{cases}$

and

$\left|\left|D^{[\tau]}\widehat{f}_{h}-D^{[\tau]}f\right|\right|_{\infty}\leq\left|\left|\mathbb{E}\left[D^{[\tau]}\widehat{f}_{h}\right]-D^{[\tau]}f\right|\right|_{\infty}+\left|\left|D^{[\tau]}\widehat{f}_{h}-\mathbb{E}\left[D^{[\tau]}\widehat{f}_{h}\right]\right|\right|_{\infty}$

for $[[\tau]]=0,1,2$. The first term in the preceding display is of order
$O(h^{2})$ inside the tangent space by Theorem 2 and the differentiability of
$f$ under condition (D1). The proof of
$\left|\left|D^{[\tau]}\widehat{f}_{h}-\mathbb{E}\left[D^{[\tau]}\widehat{f}_{h}\right]\right|\right|_{\infty}=O_{P}\left(\sqrt{\frac{|\log
h|}{nh^{q+2[[\tau]]}}}\right)$ follows directly from the argument of Theorem
2.3 in Giné and Guillou (2002) and the following calculations:

$\displaystyle\mathbb{E}\left[L^{2}\left(\frac{1}{2}\left|\left|\frac{\bm{x}-\bm{X}}{h}\right|\right|_{2}^{2}\right)\right]$
$\displaystyle=\int_{\Omega_{q}}L^{2}\left(\frac{1-\bm{x}^{T}\bm{y}}{h^{2}}\right)\cdot
f(\bm{y})\,\omega_{q}(d\bm{y})$
$\displaystyle=\int_{-1}^{1}\int_{\Omega_{q-1}}L^{2}\left(\frac{1-t}{h^{2}}\right)\cdot
f\left(t\bm{x}+\sqrt{1-t^{2}}\bm{B_{x}}\bm{\xi}\right)(1-t^{2})^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dt$
$\displaystyle=h^{q}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})\cdot
L^{2}(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dt$
$\displaystyle\leq
h^{q}||f||_{\infty}\omega_{q-1}2^{\frac{q}{2}-1}\int_{0}^{\infty}L^{2}(r)r^{\frac{q}{2}-1}dr,$
$\displaystyle\mathbb{E}\left[\left(\frac{x_{i}-X_{i}}{h}\right)^{2}\left|L^{\prime}\left(\frac{1}{2}\left|\left|\frac{\bm{x}-\bm{X}}{h}\right|\right|_{2}^{2}\right)\right|^{2}\right]$
$\displaystyle\leq\mathbb{E}\left[\left|\left|\frac{\bm{x}-\bm{X}}{h}\right|\right|_{2}^{2}\cdot\left|L^{\prime}\left(\frac{1}{2}\left|\left|\frac{\bm{x}-\bm{X}}{h}\right|\right|_{2}^{2}\right)\right|^{2}\right]$
$\displaystyle=2\int_{\Omega_{q}}\left(\frac{1-\bm{x}^{T}\bm{y}}{h^{2}}\right)\left|L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{y}}{h^{2}}\right)\right|^{2}f(\bm{y})\,\omega_{q}(d\bm{y})$
$\displaystyle=2\int_{-1}^{1}\int_{\Omega_{q-1}}\left(\frac{1-t}{h^{2}}\right)\left|L^{\prime}\left(\frac{1-t}{h^{2}}\right)\right|^{2}f\left(t\bm{x}+\sqrt{1-t^{2}}\bm{B_{x}}\bm{\xi}\right)(1-t^{2})^{\frac{q}{2}-1}\,\omega_{q-1}(d\bm{\xi})dt$
$\displaystyle=2h^{q}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}f(x+\alpha_{\bm{x},\bm{\xi}})\cdot|L^{\prime}(r)|^{2}\cdot
r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\leq
2h^{q}||f||_{\infty}\omega_{q-1}2^{\frac{q}{2}-1}\int_{0}^{\infty}|L^{\prime}(r)|^{2}r^{\frac{q}{2}}dr,$

and

$\displaystyle\mathbb{E}\left[\max\left\\{\left(\frac{x_{i}-X_{i}}{h}\right)^{4}\left|L^{\prime\prime}\left(\frac{1}{2}\left|\left|\frac{\bm{x}-\bm{X}}{h}\right|\right|_{2}^{2}\right)\right|^{2},\left(\frac{x_{i}-X_{i}}{h}\right)^{2}\left(\frac{x_{j}-X_{j}}{h}\right)^{2}\left|L^{\prime\prime}\left(\frac{1}{2}\left|\left|\frac{\bm{x}-\bm{X}}{h}\right|\right|_{2}^{2}\right)\right|^{2}\right\\}\right]$
$\displaystyle\leq\mathbb{E}\left[\left|\left|\frac{\bm{x}-\bm{X}}{h}\right|\right|_{2}^{4}L^{\prime\prime}\left(\frac{1}{2}\left|\left|\frac{\bm{x}-\bm{X}}{h}\right|\right|_{2}^{2}\right)^{2}\right]$
$\displaystyle=4\int_{\Omega_{q}}\left(\frac{1-\bm{x}^{T}\bm{y}}{h^{2}}\right)^{2}\left|L^{\prime\prime}\left(\frac{1-\bm{x}^{T}\bm{y}}{h^{2}}\right)\right|^{2}f(\bm{y})\,\omega_{q}(d\bm{y})$
$\displaystyle=4\int_{-1}^{1}\int_{\Omega_{q-1}}\left(\frac{1-t}{h^{2}}\right)^{2}\left|L^{\prime\prime}\left(\frac{1-t}{h^{2}}\right)\right|^{2}f\left(t\bm{x}+\sqrt{1-t^{2}}\bm{B_{x}}\bm{\xi}\right)(1-t^{2})^{\frac{q}{2}-1}\,\omega_{q-1}(d\bm{\xi})dt$
$\displaystyle=4h^{q}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})\cdot
L^{\prime\prime}(r)^{2}r^{\frac{q}{2}+1}(2-h^{2}r)^{\frac{q}{2}-1}\,\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\leq
4h^{q}||f||_{\infty}\omega_{q-1}2^{\frac{q}{2}-1}\int_{0}^{\infty}r^{\frac{q}{2}+1}L^{\prime\prime}(r)^{2}dr$

for $i=1,...,q+1$, where we apply (a) in Lemma 21, the change of variable
$r=\frac{1-t}{h^{2}}$, and
$\alpha_{\bm{x},\bm{\xi}}=-rh^{2}\bm{x}+h\sqrt{r(2-h^{2}r)}\bm{B_{x}}\bm{\xi}$
in the preceding three displays.

###  D.4 Proof of Theorem 6

###### Theorem 6

Assume (D1), (D2’), (K1), and (M1-2). For any $\delta\in(0,1)$, when $h$ is
sufficiently small and $n$ is sufficiently large,

  1. (a)

there must be at least one estimated local mode $\widehat{\bm{m}}_{k}$ within
$S_{k}=\bm{m}_{k}\oplus\rho_{*}$ for every $\bm{m}_{k}\in\mathcal{M}$, and

  2. (b)

the collection of estimated modes satisfies
$\widehat{\mathcal{M}}_{n}\subset\mathcal{M}\oplus\rho_{*}$ and there is a
unique estimated local mode $\widehat{\bm{m}}_{k}$ within
$S_{k}=\bm{m}_{k}\oplus\rho_{*}$

with probability at least $1-\delta$. In total, when $h$ is sufficiently small
and $n$ is sufficiently large, there exist some constants $A_{3},B_{3}>0$ such
that

$\mathbb{P}\left(\widehat{K}_{n}\neq K\right)\leq B_{3}e^{-A_{3}nh^{q+4}}.$

  1. (c)

The Hausdorff distance between the collection of local modes and its estimator
satisfies

$\mathtt{Haus}\left(\mathcal{M},\widehat{\mathcal{M}}_{n}\right)=O(h^{2})+O_{P}\left(\sqrt{\frac{1}{nh^{q+2}}}\right),$

as $h\to 0$ and $nh^{q+2}\to\infty$.

Proof  The proof is partially adopted from the proof of Theorem 1 in Chen et
al. (2016).

Statement (a). Without loss of generality, we consider the local mode
$\bm{m}_{k}$ and the set
$S_{k}=\left\\{\bm{x}\in\Omega_{q}:\left|\left|\bm{x}-\bm{m}_{k}\right|\right|\leq\rho_{*}\right\\}$.
With condition (D1), we can apply the Taylor’s expansion on the exponential
map $\mathtt{Exp}_{\bm{m}_{k}}:D_{\epsilon}\subset
T_{\bm{m}_{k}}(\Omega_{q})\to\Omega_{q}$ with
$\mathtt{Exp}_{\bm{m}_{k}}(0)=\bm{m}_{k}$, where $D_{\epsilon}$ is a disk of
radius $\epsilon$ in $T_{\bm{m}_{k}}(\Omega_{q})$ with center in the origin
and $\epsilon>\arccos\left(1-\frac{\rho_{*}^{2}}{2}\right)$. Here,
$\arccos\left(1-\frac{\rho_{*}^{2}}{2}\right)$ is the geodesic distance from
the center $\bm{m}_{k}$ to $\partial S_{k}$ on $\Omega_{q}$, where $\partial
S_{k}=\\{\bm{x}\in\Omega_{q}:||\bm{x}-\bm{m}_{k}||_{2}=\rho_{*}\\}$ is the
boundary of $S_{k}$. With (M1) and the fact that the third order partial
derivatives of $f$ are upper bounded by $C_{3}$,

$\displaystyle\begin{split}\sup_{\bm{x}\in\partial
S_{k}}f(\bm{x})&\leq\sup_{\bm{x}\in\partial
S_{k}}\Bigg{[}f(\bm{m}_{k})+\left[\mathtt{grad}\,f(\bm{m}_{k})\right]^{T}\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{x})\\\
&\quad+\frac{1}{2}\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{x})^{T}\left(\mathcal{H}_{\bm{m}_{k}}f\right)\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{x})+\frac{C_{3}}{6}||\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{x})||_{2}^{3}\Bigg{]}\\\
&\leq
f(\bm{m}_{k})-\frac{\lambda_{*}}{2}\left(\frac{3\lambda_{*}}{2C_{3}}\right)^{2}+\frac{C_{3}}{6}\left(\frac{3\lambda_{*}}{2C_{3}}\right)^{3}=f(\bm{m}_{k})-\frac{9\lambda_{*}^{3}}{8C_{3}^{2}},\end{split}$
(54)

where recall that $\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{x})\in
T_{\bm{m}_{k}}(\Omega_{q})$ is in the direction from $\bm{m}_{k}$ to $\bm{x}$
with the length equal to the great-circle (or geodesic) distance on
$\Omega_{q}$. (We indeed apply the Cauchy-Schwarz inequality implicitly to
obtain the first inequality in (54).) Then, by the uniform consistency of
$\widehat{f}_{h}$ (Theorem 4), when $h$ is sufficiently small and
$\frac{nh^{q}}{|\log h|}$ is large enough,

$\left|\left|\widehat{f}_{h}-f\right|\right|_{\infty}<\frac{9\lambda_{*}}{16C_{3}^{2}}$
(55)

with probability at least $1-\delta$ for any $0<\delta<1$. We thus conclude
that there must be at least one estimated local mode $\widehat{\bm{m}}_{k}$
within $S_{k}$. (If, on the contrary, there exists no
$\widehat{\bm{m}}_{k}\in\widehat{\mathcal{M}}_{n}$ within $S_{k}$, then the
maximum of $\widehat{f}_{h}$ on $S_{k}$ is attained at the boundary $\partial
S_{k}$, that is, $\max_{\bm{x}\in
S_{k}}\widehat{f}_{h}(\bm{x})=\max_{\bm{x}\in\partial
S_{k}}\widehat{f}_{h}(\bm{x})$. However, $\max_{\bm{x}\in\partial
S_{k}}\widehat{f}_{h}(\bm{x})<\max_{\bm{x}\in\partial
S_{k}}f(\bm{x})+\frac{9\lambda_{*}}{16C_{3}^{2}}\leq
f(\bm{m}_{k})-\frac{9\lambda_{*}}{16C_{3}^{2}}<\widehat{f}_{h}(\widehat{\bm{m}}_{k})$
by (54), contradiction.) Note that this argument can be generalized to each
$k=1,...,K$.

Statement (b). With (M2), we know that whenever

$\displaystyle\begin{split}\sup_{\bm{x}\in\Omega_{q}}\left|\left|\mathtt{grad}\,\widehat{f}_{h}(\bm{x})-\mathtt{grad}\,f(\bm{x})\right|\right|_{\max}&=\sup_{\bm{x}\in\Omega_{q}}\left|\left|\mathtt{Tang}\left(\nabla\widehat{f}_{h}(\bm{x})\right)-\mathtt{Tang}\left(\nabla
f(\bm{x})\right)\right|\right|_{\max}\leq\Theta_{1},\\\
\sup_{\bm{x}\in\Omega_{q}}\left|\left|\mathcal{H}\widehat{f}_{h}(\bm{x})-\mathcal{H}f(\bm{x})\right|\right|_{\max}&\leq\Theta_{2}\end{split}$
(56)

for some $\Theta_{2}>0$, the followings hold simultaneously:

  1. (i)

$\left|\left|\mathtt{grad}\,f(\widehat{\bm{m}}_{k})\right|\right|_{\max}=\left|\left|\mathtt{grad}\,f(\widehat{\bm{m}}_{k})-\underbrace{\mathtt{grad}\,\widehat{f}_{h}(\widehat{\bm{m}}_{k})}_{=0}\right|\right|_{\max}\leq\Theta_{1}$,

  2. (ii)

$\sup_{\bm{x}\in
S_{k}}\lambda_{1}\left(\mathcal{H}\widehat{f}_{h}(\bm{x})\right)<0$ and
$\lambda_{1}\left(\mathcal{H}\widehat{f}_{h}(\widehat{\bm{m}}_{k})\right)\leq-\frac{\lambda_{*}}{2}-(q+1)\Theta_{2}$
by choosing $\Theta_{2}>0$ properly,

  3. (iii)

and

$\displaystyle\lambda_{1}\left(\mathcal{H}f(\widehat{\bm{m}}_{k})\right)$
$\displaystyle\leq\lambda_{1}\left(\mathcal{H}\widehat{f}_{h}(\widehat{\bm{m}}_{k})\right)+\lambda_{q-1}\left(\mathcal{H}f(\widehat{\bm{m}}_{k})-\mathcal{H}\widehat{f}_{h}(\widehat{\bm{m}}_{k})\right)$
$\displaystyle\leq-\frac{\lambda_{*}}{2}-(q+1)\Theta_{2}+(q+1)\Theta_{2}=-\frac{\lambda_{*}}{2}$

by Weyl’s theorem (Theorem 4.3.1 in Horn and Johnson 2012) and the fact that

$\displaystyle\left|\lambda_{q-1}(\mathcal{H}f(\widehat{\bm{m}}_{k})-\mathcal{H}\widehat{f}_{h}(\widehat{\bm{m}}_{k})\right|$
$\displaystyle\leq\sup_{\left|\left|v\right|\right|_{2}=1}\left|\left|\left[\mathcal{H}f(\widehat{\bm{m}}_{k})-\mathcal{H}\widehat{f}_{h}(\widehat{\bm{m}}_{k})\right]v\right|\right|_{2}$
$\displaystyle\leq\sqrt{(q+1)\times(q+1)}\left|\left|\mathcal{H}f(\widehat{\bm{m}}_{k})-\mathcal{H}\widehat{f}_{h}(\widehat{\bm{m}}_{k})\right|\right|_{\max}$
$\displaystyle\leq(q+1)\Theta_{2}.$

See Section 3.3 in Genovese et al. (2014) for detailed relations between
different types of matrix norms.

Notice that (ii) is true because $\lambda_{1}(\bm{m}_{k})\leq\lambda_{*}$ by
(M1) and the difference between $\mathcal{H}\widehat{f}_{h}$ and
$\mathcal{H}f$ will be minute given a small $\Theta_{2}$. By (i) and (iii), we
conclude that $\widehat{\mathcal{M}}_{n}\subset\mathcal{M}\oplus\rho_{*}$. By
(ii) and Lemma 3.2 in Banyaga and Hurtubise (2004), there is only one
estimated local mode $\widehat{\bm{m}}_{k}$ within $S_{k}$. They both hold
with probability at least $1-\delta$ for any $\delta\in(0,1)$.

In total, a sufficient condition for the number of true local modes and
estimated local modes being the same is a combination of the inequalities in
(55) and (56). That is,

$\displaystyle\begin{split}\left|\left|\widehat{f}_{h}-f\right|\right|_{\infty}&<\frac{9\lambda_{*}}{16C_{3}^{2}}\\\
\left|\left|\mathtt{Tang}\left(\nabla\widehat{f}_{h}\right)-\mathtt{Tang}\left(\nabla
f\right)\right|\right|_{\max,\infty}&\leq\Theta_{1}\\\
\left|\left|\widehat{\mathcal{H}}f-\mathcal{H}f\right|\right|_{\max,\infty}&\leq\Theta_{2}.\end{split}$
(57)

By bias bounds in Theorem 2 or Theorem 4, as $h$ is sufficiently small, we
have

$\displaystyle\left|\left|\mathbb{E}\left[\widehat{f}_{h}\right]-f\right|\right|_{\infty}<\frac{9\lambda_{*}}{32C_{3}^{2}},$
$\displaystyle\quad\sup_{\bm{x}\in\Omega_{q}}\left|\left|\mathbb{E}\left[\mathtt{grad}\,\widehat{f}_{h}(\bm{x})\right]-\mathtt{grad}\,f(\bm{x})\right|\right|_{\max}\leq\frac{\Theta_{1}}{2},$
$\displaystyle\text{ and
}\sup_{\bm{x}\in\Omega_{q}}\left|\left|\mathbb{E}\left[\mathcal{H}\widehat{f}_{h}(\bm{x})\right]-\mathcal{H}f(\bm{x})\right|\right|_{\max}\leq\frac{\Theta_{2}}{2}.$

Therefore, (57) holds whenever

$\displaystyle\begin{split}\left|\left|\widehat{f}_{h}-\mathbb{E}\left[\widehat{f}_{h}\right]\right|\right|_{\infty}&<\frac{9\lambda_{*}}{32C_{3}^{2}}\\\
\sup_{\bm{x}\in\Omega_{q}}\left|\left|\mathtt{grad}\,\widehat{f}_{h}(\bm{x})-\mathbb{E}\left[\mathtt{grad}\,\widehat{f}_{h}(\bm{x})\right]\right|\right|_{\max}&\leq\frac{\Theta_{1}}{2}\\\
\sup_{\bm{x}\in\Omega_{q}}\left|\left|\mathcal{H}\widehat{f}_{h}(\bm{x})-\mathbb{E}\left[\mathcal{H}\widehat{f}_{h}(\bm{x})\right]\right|\right|_{\max}&\leq\frac{\Theta_{2}}{2}\end{split}$
(58)

and $h$ is sufficiently small. Now applying Talagrand’s inequality Talagrand
(1996); Giné and Guillou (2002), there exist constants $A_{0},A_{1},A_{2}>0$
and $B_{0},B_{1},B_{2}>0$ such that when $n$ is large enough,

$\displaystyle\begin{split}\mathbb{P}\left(\left|\left|\widehat{f}_{h}-\mathbb{E}\left[\widehat{f}_{h}\right]\right|\right|_{\infty}\geq\epsilon\right)&\leq
B_{0}e^{-A_{0}\epsilon^{2}nh^{q}}\\\
\mathbb{P}\left(\sup_{\bm{x}\in\Omega_{q}}\left|\left|\mathtt{grad}\,\widehat{f}_{h}(\bm{x})-\mathbb{E}\left[\mathtt{grad}\,\widehat{f}_{h}(\bm{x})\right]\right|\right|_{\max}\geq\epsilon\right)&\leq
B_{1}e^{-A_{1}\epsilon^{2}nh^{q+2}}\\\
\mathbb{P}\left(\sup_{\bm{x}\in\Omega_{q}}\left|\left|\mathcal{H}\widehat{f}_{h}(\bm{x})-\mathbb{E}\left[\mathcal{H}\widehat{f}_{h}(\bm{x})\right]\right|\right|_{\max}\geq\epsilon\right)&\leq
B_{2}e^{-A_{2}\epsilon^{2}nh^{q+4}}.\end{split}$ (59)

Combining (58) and (59), we conclude that there exist some constants
$A_{3},B_{3}>0$ such that

$\mathbb{P}\Big{(}\text{(\ref{mode_consist_cond_com})
holds}\,\Big{)}>1-B_{3}e^{-A_{3}nh^{q+4}}$

when $h$ is sufficiently small. Since the condition (57) implies
$\widehat{K}_{n}=K$, we conclude that

$\mathbb{P}\left(\widehat{K}_{n}\neq K\right)\leq B_{3}e^{-A_{3}nh^{q+4}}$

for some constants $A_{3},B_{3}>0$ as $h$ is sufficiently small. This proves
the so-called modal consistency.

Statement (c). To establish the convergence rate of the Hausdorff distance
between $\widehat{\mathcal{M}}_{n}$ and $\mathcal{M}$, we assume that (57)
holds so that $K=\widehat{K}_{n}$ and each local mode is approximating by an
unique estimated local mode. Notice that
$\left|\left|\bm{m}_{k}-\widehat{\bm{m}}_{k}\right|\right|_{2}$ is upper
bounded by the great-circle distance between these two points. Then,

$\displaystyle\begin{split}&\mathtt{grad}\,f(\widehat{\bm{m}}_{k})\\\
&=\mathtt{Tang}\left(\nabla
f(\widehat{\bm{m}}_{k})\right)-\underbrace{\mathtt{Tang}\left(\nabla
f(\bm{m}_{k})\right)}_{=0}\\\ &=\nabla\mathtt{Tang}\left(\nabla
f(\bm{m}_{k})\right)\cdot\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\widehat{\bm{m}}_{k})+o\left(\left|\left|\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\widehat{\bm{m}}_{k})\right|\right|_{2}\right)\\\
&=\left[(I_{q+1}-\bm{m}_{k}\bm{m}_{k}^{T})\nabla\nabla
f(\bm{m}_{k})-\bm{m}_{k}^{T}\nabla f(\bm{m}_{k})I_{q+1}-\bm{m}_{k}\nabla
f(\bm{m}_{k})^{T}\right]\cdot\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\widehat{\bm{m}}_{k})\\\
&\quad+o\left(\left|\left|\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\widehat{\bm{m}}_{k})\right|\right|_{2}\right)\\\
&=\left[\mathcal{H}f(\bm{m}_{k})\right]\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\widehat{\bm{m}}_{k})+o\left(\left|\left|\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\widehat{\bm{m}}_{k})\right|\right|_{2}\right),\end{split}$
(60)

because $\nabla f(\bm{m}_{k})=\left|\left|\nabla
f(\bm{m}_{k})\right|\right|_{2}\cdot\bm{m}_{k}$ when $\bm{m}_{k}$ is a local
mode and $\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\widehat{\bm{m}}_{k})\in
T_{\bm{m}_{k}}$ is orthogonal to $\bm{m}_{k}$. Under (M1), the matrices
$\mathcal{H}f{\bm{m}_{k}}$ are nonsingular for all $\bm{m}_{k}\in\mathcal{M}$
inside the tangent space $T_{\bm{m}_{k}}$, respectively. As the chord distance
between two points on $\Omega_{q}$ is bounded by their great-circle distance,
we multiply $\left[\mathcal{H}f(\bm{m}_{k})\right]^{-1}$ on both sides of (60)
and obtain that

$\left|\left|\widehat{\bm{m}}_{k}-\bm{m}_{k}\right|\right|_{2}\leq\left|\left|\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\widehat{\bm{m}}_{k})\right|\right|_{2}=\left[\mathcal{H}f(\bm{m}_{k})\right]^{-1}\mathtt{grad}\,f(\widehat{\bm{m}}_{k})+o\left(\left|\left|\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\widehat{\bm{m}}_{k})\right|\right|_{2}\right),$

where the matrix inverse, strictly speaking, is taken with respect to the
local coordinate system near $\bm{m}_{k}$. Note that
$\left|\left|\left[\mathcal{H}f(\bm{m}_{k})\right]^{-1}\right|\right|_{\max}$
is bounded within $T_{\bm{m}_{k}}$ for all $\bm{m}_{k}\in\mathcal{M}$ under
the assumption (57). Moreover, by Theorem 2,

$\displaystyle\mathtt{grad}\,f(\widehat{\bm{m}}_{k})$
$\displaystyle=\mathtt{grad}\,f(\widehat{\bm{m}}_{k})-\underbrace{\mathtt{grad}\,\widehat{f}_{h}(\widehat{\bm{m}}_{k})}_{=0}$
$\displaystyle=O(h^{2})+O_{P}\left(\sqrt{\frac{1}{nh^{q+2}}}\right).$

Now applying this rate of convergence to each local mode and using the fact
that

$\mathtt{Haus}\left(\widehat{\mathcal{M}}_{n},\mathcal{M}\right)=\max_{k=1,...,K}||\widehat{\bm{m}}_{k}-\bm{m}_{k}||_{2},$

we obtain the final conclusion.

###  D.5 Proofs of Theorem 8, Lemma 10, and Theorem 11

######  Theorem 8 (Ascending Property)

If kernel $L:[0,\infty)\to[0,\infty)$ is monotonically decreasing,
differentiable, and convex with $L(0)<\infty$, then the sequence
$\left\\{\widehat{f}_{h}(\widehat{\bm{y}}_{s})\right\\}_{s=0}^{\infty}$ is
monotonically increasing and thus converges.

Proof  Obviously, the sequence
$\left\\{\widehat{f}_{h}(\widehat{\bm{y}}_{s})\right\\}_{s=0,1,...}$ is
bounded if the kernel function $L$ is monotonically decreasing with
$L(0)<\infty$. Hence, it suffices to show that it is monotonically increasing.
The convexity and differentiability of kernel $L$ imply that

$L(x_{2})-L(x_{1})\geq L^{\prime}(x_{1})\cdot(x_{2}-x_{1})$ (61)

for all $x_{1},x_{2}\in[0,\infty),x_{1}\neq x_{2}$. Using the fact that
(rearrangement from Algorithm 1)

$\sum\limits_{i=1}^{n}\bm{X}_{i}L^{\prime}\left(\frac{1-\widehat{\bm{y}}_{s}^{T}\bm{X}_{i}}{h^{2}}\right)=-\widehat{\bm{y}}_{s+1}\left|\left|\sum\limits_{i=1}^{n}\bm{X}_{i}L^{\prime}\left(\frac{1-\bm{y}_{s}^{T}\bm{X}_{i}}{h^{2}}\right)\right|\right|_{2}$

we have that

$\displaystyle\begin{split}\widehat{f}_{h}(\bm{y}_{s+1})-\widehat{f}_{h}(\bm{y}_{s})&=\frac{c_{h,q}(L)}{n}\sum_{i=1}^{n}\left[L\left(\frac{1-\bm{y}_{s+1}^{T}\bm{X}_{i}}{h^{2}}\right)-L\left(\frac{1-\bm{y}_{s}^{T}\bm{X}_{i}}{h^{2}}\right)\right]\\\
&\geq\frac{c_{h,q}(L)}{nh^{2}}\sum_{i=1}^{n}L^{\prime}\left(\frac{1-\bm{y}_{s}^{T}\bm{X}_{i}}{h^{2}}\right)\cdot(\bm{y}_{s}-\bm{y}_{s+1})^{T}\bm{X}_{i}\\\
&=\frac{c_{h,q}(L)}{nh^{2}}\cdot(\bm{y}_{s+1}-\bm{y}_{s})^{T}\bm{y}_{s+1}\cdot\left|\left|\sum\limits_{i=1}^{n}\bm{X}_{i}L^{\prime}\left(\frac{1-\bm{y}_{s}^{T}\bm{X}_{i}}{h^{2}}\right)\right|\right|_{2}\\\
&=\frac{c_{h,q}(L)}{2nh^{2}}\left|\left|\bm{y}_{s+1}-\bm{y}_{s}\right|\right|_{2}^{2}\cdot\left|\left|\sum\limits_{i=1}^{n}\bm{X}_{i}L^{\prime}\left(\frac{1-\bm{y}_{s}^{T}\bm{X}_{i}}{h^{2}}\right)\right|\right|_{2}\\\
&\geq 0,\end{split}$ (62)

where we use the fact that
$2(\bm{y}_{s+1}-\bm{y}_{s})^{T}\bm{y}_{s+1}=2-2\bm{y}_{s}^{T}\bm{y}_{s+1}=\left|\left|\bm{y}_{s+1}-\bm{y}_{s}\right|\right|_{2}^{2}$
between the third and fourth lines, given that
$\left|\left|\bm{y}_{s}\right|\right|_{2}=\left|\left|\bm{y}_{s+1}\right|\right|_{2}=1$.

###### Lemma 10

Assume conditions (D1) and (D2’). For any fixed $\bm{x}\in\Omega_{q}$, we have

$h^{2}\cdot\mathtt{Rad}\left(\nabla\widehat{f}_{h}(\bm{x})\right)\asymp
h^{2}\cdot\nabla\widehat{f}_{h}(\bm{x})=\bm{x}f(\bm{x})C_{L,q}+o\left(1\right)+O_{P}\left(\sqrt{\frac{1}{nh^{q}}}\right)$

as $nh^{q}\to\infty$ and $h\to 0$, where
$C_{L,q}=-\frac{\int_{0}^{\infty}L^{\prime}(r)r^{\frac{q}{2}-1}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}>0$
is a constant depending only on kernel $L$ and dimension $q$ and “$\asymp$”
stands for an asymptotic equivalence.

Proof  The proof follows the same logic as the one for Theorem 2. Note that

$\nabla\widehat{f}_{h}(\bm{x})=\mathbb{E}\left[\nabla\widehat{f}_{h}(\bm{x})\right]+\nabla\widehat{f}_{h}(\bm{x})-\mathbb{E}\left[\nabla\widehat{f}_{h}(\bm{x})\right].$
(63)

Recall that
$\nabla\widehat{f}_{h}(\bm{x})=-\frac{c_{h,q}(L)}{nh^{2}}\sum\limits_{i=1}^{n}\bm{X}_{i}L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{X}_{i}}{h^{2}}\right)$.
The expectation of $\nabla\widehat{f}_{h}(\bm{x})$ is

$\displaystyle\begin{split}\mathbb{E}\left[\nabla\widehat{f}_{h}(\bm{x})\right]&=\frac{c_{h,q}(L)}{h^{2}}\int_{\Omega_{q}}(-\bm{y})L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{y}}{h^{2}}\right)f(\bm{y})\,\omega_{q}(d\bm{y})\\\
&=\frac{c_{h,q}(L)}{h^{2}}\int_{-1}^{1}\int_{\Omega_{q-1}}\left(-t\bm{x}-\sqrt{1-t^{2}}\bm{B_{x}}\bm{\xi}\right)L^{\prime}\left(\frac{1-t}{h^{2}}\right)\\\
&\quad\times
f\left(-t\bm{x}-\sqrt{1-t^{2}}\bm{B_{x}}\bm{\xi}\right)(1-t^{2})^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dt\\\
&=c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}(-\bm{x}-\alpha_{\bm{x},\bm{\xi}})\cdot
L^{\prime}(r)\\\ &\quad\times f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})\cdot
r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr\end{split}$
(64)

by (a) in Lemma 21 and a change of variable $r=\frac{1-t}{h^{2}}$, where
$\alpha_{\bm{x},\bm{\xi}}=-rh^{2}\bm{x}+h\sqrt{r(2-h^{2}r)}\bm{B_{x}}\bm{\xi}$.
By condition (D1), the first-order Taylor’s expansion of $f$ at
$\bm{x}\in\Omega_{q}$ is

$f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})=f(\bm{x})+O\left(\left|\left|\alpha_{\bm{x},\bm{\xi}}\right|\right|_{2}\right),$

where $\left|\left|\alpha_{\bm{x},\bm{\xi}}\right|\right|_{2}^{2}=2rh^{2}$ by
the orthogonality of $\bm{x}$ and columns of $\bm{B_{x}}$. Now we plug it back
into (64) respectively to compute the dominating term of
$\mathbb{E}\left[\nabla\widehat{f}_{h}(\bm{x})\right]$.

$\displaystyle\mathbb{E}\left[\nabla\widehat{f}_{h}(\bm{x})\right]$
$\displaystyle=-c_{h,q}(L)h^{q-2}\bm{x}f(\bm{x})\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}L^{\prime}(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\quad-
c_{h,q}(L)h^{q-2}f(\bm{x})\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\alpha_{\bm{x},\bm{\xi}}L^{\prime}(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\quad+O(h)\cdot
c_{h,q}(L)h^{q-2}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}(-\bm{x}-\alpha_{\bm{x},\bm{\xi}})L^{\prime}(r)r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr$
$\displaystyle\stackrel{{\scriptstyle\text{(i)}}}{{=}}-c_{h,q}(L)h^{q-2}\bm{x}f(\bm{x})\cdot\omega_{q-1}\int_{0}^{2h^{-2}}L^{\prime}(r)\cdot
r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\,dr$
$\displaystyle\quad+c_{h,q}(L)h^{q}\bm{x}f(\bm{x})\cdot\omega_{q-1}\int_{0}^{2h^{-2}}L^{\prime}(r)\cdot
r^{\frac{q}{2}}(2-h^{2}r)^{\frac{q}{2}-1}\,dr$ $\displaystyle\quad-
c_{h,q}(L)h^{q-1}f(\bm{x})\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}\bm{B_{x}}\bm{\xi}\cdot
L^{\prime}(r)r^{\frac{q-1}{2}}(2-h^{2}r)^{\frac{q-1}{2}}\omega_{q-1}(d\bm{\xi})dr+O(h^{-1}),$
$\displaystyle\stackrel{{\scriptstyle\text{(ii)}}}{{=}}-c_{h,q}(L)h^{q-2}\bm{x}f(\bm{x})\cdot\omega_{q-1}\int_{0}^{2h^{-2}}L^{\prime}(r)\cdot
r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\,dr+O(1)+0+O(h^{-1})$
$\displaystyle\stackrel{{\scriptstyle\text{(iii)}}}{{=}}-\frac{\bm{x}f(\bm{x})}{h^{2}}\cdot\frac{\int_{0}^{\infty}L^{\prime}(r)r^{\frac{q}{2}-1}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}+o(h^{-2})$
$\displaystyle\equiv C_{L,q}\cdot\bm{x}f(\bm{x})h^{-2}+o(h^{-2}),$

where we use (b) of Lemma 21 and the fact that
$\bm{B_{x}}\bm{\xi}=\sum_{i=1}^{q}\xi_{i}\bm{b}_{i}$ in (ii), and apply
condition (D2’) and (3) to argue that

$c_{h,q}(L)h^{q}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}L^{\prime}(r)\cdot\phi(r,\bm{\xi})\,\,\omega_{q-1}(d\bm{\xi})dr\asymp
O(1)$ (65)

in both (i) and (ii). We also use the asymptotic relation (3) in (iii) and
denote
$C_{L,q}=-\frac{\int_{0}^{\infty}L^{\prime}(r)r^{\frac{q}{2}-1}dr}{\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr}$
in the last equality. Thus,

$\mathbb{E}\left[\mathtt{Rad}\left(\nabla\widehat{f}_{h}(\bm{x})\right)\right]=\bm{x}\bm{x}^{T}\mathbb{E}\left[\nabla\widehat{f}_{h}(\bm{x})\right]=C_{L,q}\cdot\bm{x}f(\bm{x})h^{-2}+o(h^{-2}).$

Based on the asymptotic rate of
$\mathbb{E}\left[\nabla\widehat{f}_{h}(\bm{x})\right]$, we calculate the
covariance matrix of $\nabla\widehat{f}_{h}(\bm{x})$ as

$\displaystyle\text{Cov}\left[\nabla\widehat{f}_{h}(\bm{x})\right]$
$\displaystyle=\frac{c_{h,q}(L)^{2}}{nh^{4}}\cdot\text{Cov}\left[\bm{X}_{1}\cdot
L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{X}_{1}}{h^{2}}\right)\right]$
$\displaystyle=\frac{c_{h,q}(L)^{2}}{nh^{4}}\cdot\mathbb{E}\left[\bm{X}_{1}\bm{X}_{1}^{T}\cdot
L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{X}_{1}}{h^{2}}\right)^{2}\right]-\frac{1}{n}\cdot\mathbb{E}\left[\nabla\widehat{f}_{h}(\bm{x})\right]\mathbb{E}\left[\nabla\widehat{f}_{h}(\bm{x})\right]^{T}$
$\displaystyle=\frac{c_{h,q}(L)^{2}}{nh^{4}}\int_{\Omega_{q}}\bm{y}\bm{y}^{T}L^{\prime}\left(\frac{1-\bm{x}^{T}\bm{y}}{h^{2}}\right)^{2}f(\bm{y})\,\omega_{q}(d\bm{y})+O\left(\frac{1}{nh^{4}}\right)$
$\displaystyle=\frac{c_{h,q}(L)^{2}}{n}h^{q-4}\int_{0}^{2h^{-2}}\int_{\Omega_{q-1}}(\bm{x}+\alpha_{\bm{x},\bm{\xi}})(\bm{x}+\alpha_{\bm{x},\bm{\xi}})^{T}L^{\prime}(r)^{2}$
$\displaystyle\quad\times
f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}\omega_{q-1}(d\bm{\xi})dr+O\left(\frac{1}{nh^{4}}\right).$

With condition (D1), we carry out the first-order Taylor’s expansion of $f$ at
$\bm{x}\in\Omega_{q}$ as

$f(\bm{x}+\alpha_{\bm{x},\bm{\xi}})=f(\bm{x})+O\left(\left|\left|\alpha_{\bm{x},\bm{\xi}}\right|\right|_{2}\right)=f(\bm{x})+O(h).$

Therefore,

$\displaystyle\text{Cov}\left[\nabla\widehat{f}_{h}(\bm{x})\right]$
$\displaystyle=\frac{c_{h,q}(L)^{2}}{n}h^{q-4}\bm{x}\bm{x}^{T}f(\bm{x})\omega_{q-1}\int_{0}^{2h^{-2}}L^{\prime}(r)^{2}r^{\frac{q}{2}-1}(2-h^{2}r)^{\frac{q}{2}-1}dr+o\left(\frac{1}{nh^{q+4}}\right)$
$\displaystyle=\frac{\bm{x}\bm{x}^{T}f(\bm{x})}{nh^{q+4}}\cdot\frac{\int_{0}^{\infty}L^{\prime}(r)^{2}r^{\frac{q}{2}-1}dr}{\omega_{q-1}2^{\frac{q}{2}-1}\left(\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr\right)^{2}}+o\left(\frac{1}{nh^{q+4}}\right),$

where we use (b) of Lemma 21, asymptotic rate (3), and (65) to absorb some
higher order terms into $o\left(\frac{1}{nh^{q+4}}\right)$. The dominating
term of $\text{Cov}\left[\nabla\widehat{f}_{h}(\bm{x})\right]$ is in the
radial direction, so by the central limit theorem,

$\displaystyle\mathtt{Rad}\left(\nabla\widehat{f}_{h}(\bm{x})\right)-\mathbb{E}\left[\mathtt{Rad}\left(\nabla\widehat{f}_{h}(\bm{x})\right)\right]$
$\displaystyle\asymp\nabla\widehat{f}_{h}(\bm{x})-\mathbb{E}\left[\nabla\widehat{f}_{h}(\bm{x})\right]$
$\displaystyle=\text{Cov}\left[\nabla\widehat{f}_{h}(\bm{x})\right]^{\frac{1}{2}}\cdot\text{Cov}\left[\nabla\widehat{f}_{h}(\bm{x})\right]^{-\frac{1}{2}}\left\\{\nabla\widehat{f}_{h}(\bm{x})-\mathbb{E}\left[\nabla\widehat{f}_{h}(\bm{x})\right]\right\\}$
$\displaystyle=\left[\frac{\bm{x}\bm{x}^{T}f(\bm{x})}{nh^{q+4}}\cdot\frac{\int_{0}^{\infty}L^{\prime}(r)^{2}r^{\frac{q}{2}-1}dr}{\omega_{q-1}2^{\frac{q}{2}-1}\left(\int_{0}^{\infty}L(r)r^{\frac{q}{2}-1}dr\right)^{2}}+o\left(\frac{1}{nh^{q+4}}\right)\right]^{\frac{1}{2}}\cdot\widehat{\bm{Z}}_{n}(\bm{x})$
$\displaystyle=O_{P}\left(\sqrt{\frac{1}{nh^{q+4}}}\right),$

where $\widehat{\bm{Z}}_{n}(\bm{x})\stackrel{{\scriptstyle
d}}{{\to}}N_{q+1}(\bm{0},I_{q+1})$. In total, we conclude with (63) that

$\mathtt{Rad}\left(h^{2}\nabla\widehat{f}_{h}(\bm{x})\right)\asymp
h^{2}\nabla\widehat{f}_{h}(\bm{x})=\bm{x}f(\bm{x})C_{L,q}+o\left(1\right)+O_{P}\left(\sqrt{\frac{1}{nh^{q}}}\right)$

for any fixed $\bm{x}\in\Omega_{q}$, as $h\to 0$ and $nh^{q}\to\infty$.

Before we prove Theorem 11, we first note the following useful result.

###### Proposition 23

Assume (C1) and the conditions on the kernel $L$ in Theorem 8. Then for any
mode $\widehat{\bm{m}}_{k}\in\widehat{\mathcal{M}}_{n}$ satisfying (C2), we
have that
$\widehat{\bm{m}}_{k}^{T}\nabla\widehat{f}_{h}(\widehat{\bm{m}}_{k})>0$.

Proof  Suppose, on the contrary, that
$\widehat{\bm{m}}_{k}^{T}\nabla\widehat{f}_{h}(\widehat{\bm{m}}_{k})<0$. By
the definition of a local mode $\widehat{\bm{m}}_{k}$ of $\widehat{f}_{h}$ on
$\Omega_{q}$, we know that
$\left|\left|\mathtt{grad}\,\widehat{f}_{h}(\widehat{\bm{m}}_{k})\right|\right|_{2}=\left|\left|\mathtt{Tang}\left(\nabla\widehat{f}_{h}(\widehat{\bm{m}}_{k})\right)\right|\right|_{2}=0$.
Then

$\frac{\nabla\widehat{f}_{h}(\widehat{\bm{m}}_{k})}{\left|\left|\nabla\widehat{f}_{h}(\widehat{\bm{m}}_{k})\right|\right|_{2}}=-\widehat{\bm{m}}_{k}$

and by (C1), there exist a $\widehat{r}_{k}\in(0,2]$ such that
$\mathtt{Tang}\left(\nabla\widehat{f}_{h}(\bm{y})\right)\neq 0$ and
$\widehat{f}_{h}(\bm{y})\leq\widehat{f}_{h}(\widehat{\bm{m}}_{k})$ for any
$\bm{y}\in\left\\{\bm{z}\in\Omega_{q}:\bm{z}^{T}\widehat{\bm{m}}_{k}\geq
1-\frac{\widehat{r}_{k}^{2}}{2}\right\\}\setminus\left\\{\widehat{\bm{m}}_{k}\right\\}=\left\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\widehat{\bm{m}}_{k}\right|\right|_{2}\leq\widehat{r}_{k}\right\\}\setminus\left\\{\widehat{\bm{m}}_{k}\right\\}$.
That is, $\widehat{\bm{m}}_{k}$ is the unique mode inside its neighborhood
$\left\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\widehat{\bm{m}}_{k}\right|\right|_{2}\leq\widehat{r}_{k}\right\\}$.
Since the sum of convex functions is convex, $\widehat{f}_{h}$ is indeed
convex and we deduce that when
$\bm{y}\in\left\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\widehat{\bm{m}}_{k}\right|\right|_{2}\leq\widehat{r}_{k}\right\\}\setminus\left\\{\widehat{\bm{m}}_{k}\right\\}$,

$\displaystyle\widehat{f}_{h}(\widehat{\bm{m}}_{k})-\widehat{f}_{h}(\bm{y})$
$\displaystyle\leq\frac{c_{h,q}(L)}{nh^{2}}\sum_{i=1}^{n}L^{\prime}\left(\frac{1-\widehat{\bm{m}}_{k}^{T}\bm{X}_{i}}{h^{2}}\right)\bm{X}_{i}^{T}(\bm{y}-\widehat{\bm{m}}_{k})$
$\displaystyle=\left|\left|\nabla\widehat{f}_{h}(\widehat{\bm{m}}_{k})\right|\right|_{2}\cdot(-\widehat{\bm{m}}_{k})^{T}(\widehat{\bm{m}}_{k}-\bm{y})$
$\displaystyle=\left|\left|\nabla\widehat{f}_{h}(\widehat{\bm{m}}_{k})\right|\right|_{2}\cdot(\widehat{\bm{m}}_{k}^{T}\bm{y}-1)$
$\displaystyle<0$

contradicting to the fact that $\widehat{\bm{m}}_{k}$ is the unique local mode
in
$\left\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\widehat{\bm{m}}_{k}\right|\right|_{2}\leq\widehat{r}_{k}\right\\}$.
The result follows.

###### Theorem 11

Assume (C1) and (C2) and the conditions on kernel $L$ in Theorem 8. We further
assume that $L$ is continuously differentiable. Then, for each local mode
$\widehat{\bm{m}}_{k}\in\widehat{\mathcal{M}}_{n}$, there exists a
$\widehat{r}_{k}>0$ such that the sequence
$\\{\widehat{\bm{y}}_{s}\\}_{s=0}^{\infty}$ converges to
$\widehat{\bm{m}}_{k}$ whenever the initial point
$\widehat{\bm{y}}_{0}\in\Omega_{q}$ satisfies
$\left|\left|\widehat{\bm{y}}_{0}-\widehat{\bm{m}}_{k}\right|\right|_{2}\leq\widehat{r}_{k}$.
Moreover, under conditions (D1) and (D2’), there exists a fixed constant
$r^{*}>0$ such that $\mathbb{P}(\widehat{r}_{k}\geq r^{*})\to 1$ as $h\to 0$
and $nh^{q}\to\infty$.

Proof  By the definition of a local mode $\widehat{\bm{m}}_{k}$ of
$\widehat{f}_{h}$ on $\Omega_{q}$,

$\left|\left|\mathtt{grad}\,\widehat{f}_{h}(\widehat{\bm{m}}_{k})\right|\right|_{2}=\left|\left|\mathtt{Tang}\left(\nabla\widehat{f}_{h}(\widehat{\bm{m}}_{k})\right)\right|\right|_{2}=0.$

Hence, with condition (C2) imposed on $\widehat{\bm{m}}_{k}$ and Proposition
23,

$\frac{\nabla\widehat{f}_{h}(\widehat{\bm{m}}_{k})}{\left|\left|\nabla\widehat{f}_{h}(\widehat{\bm{m}}_{k})\right|\right|_{2}}=\widehat{\bm{m}}_{k}\quad\text{
and
}\quad\widehat{\bm{m}}_{k}^{T}\cdot\frac{\nabla\widehat{f}_{h}(\widehat{\bm{m}}_{k})}{\left|\left|\nabla\widehat{f}_{h}(\widehat{\bm{m}}_{k})\right|\right|_{2}}=1.$

It indicates that our one-step fixed-point iteration of Algorithm 1 on the
local mode $\widehat{\bm{m}}_{k}$ will yield $\widehat{\bm{m}}_{k}$ itself.
(This is the so-called consistency of fixed-point iterations.) Moreover, there
exists a $\widehat{r}_{k}>0$ such that $\widehat{\bm{m}}_{k}$ is the only
point in
$\\{\bm{y}\in\Omega_{q}:||\bm{y}-\widehat{\bm{m}}_{k}||_{2}\leq\widehat{r}_{k}\\}$
satisfying
$\left|\left|\mathtt{Tang}\left(\nabla\widehat{f}_{h}(\bm{y})\right)\right|\right|_{2}=0$.
See Figure 4 for a graphical illustration.
In addition, given that $L$ is continuously differentiable, we may shrink
$\widehat{r}_{k}>0$ if necessary so that

$\widehat{\bm{m}}_{k}^{T}\cdot\frac{\nabla\widehat{f}_{h}(\bm{y})}{\left|\left|\nabla\widehat{f}_{h}(\bm{y})\right|\right|_{2}}\geq
1-\frac{\widehat{r}_{k}^{2}}{2}\quad\text{ and
}\quad\left|\left|\sum\limits_{i=1}^{n}\bm{X}_{i}L^{\prime}\left(\frac{1-\bm{y}^{T}\bm{X}_{i}}{h^{2}}\right)\right|\right|_{2}=\frac{nh^{2}}{c_{h,q}(L)}\left|\left|\nabla\widehat{f}_{h}(\bm{y})\right|\right|_{2}\geq\widehat{C}_{k}$
(66)

for all
$\bm{y}\in\left\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\widehat{\bm{m}}_{k}\right|\right|_{2}\leq\widehat{r}_{k}\right\\}$
and some constant $\widehat{C}_{k}>0$. The first inequality in (66) ensures
that the sequence $\\{\widehat{\bm{y}}_{s}\\}_{s=0}^{\infty}$ yielded by our
fixed-point iteration will not jump outside of the set
$\\{\bm{y}\in\Omega_{q}:\left|\left|\bm{y}-\widehat{\bm{m}}_{k}\right|\right|_{2}\leq\widehat{r}_{k}\\}$
as long as the initial point $\widehat{\bm{y}}_{0}$ is in the set. It also
guarantees the correctness of the second inequality in (66) for the iterative
sequence $\left\\{\widehat{\bm{y}}_{s}\right\\}_{s=0}^{\infty}$. By (62) in
the proof of Theorem 8, we know that

$\displaystyle\widehat{f}_{h}(\widehat{\bm{y}}_{s+1})-\widehat{f}_{h}(\widehat{\bm{y}}_{s})$
$\displaystyle\geq\frac{c_{h,q}(L)}{2nh^{2}}||\widehat{\bm{y}}_{s+1}-\widehat{\bm{y}}_{s}||_{2}^{2}\cdot\left|\left|\sum\limits_{i=1}^{n}\bm{X}_{i}L^{\prime}\left(\frac{1-\widehat{\bm{y}}_{s}^{T}\bm{X}_{i}}{h^{2}}\right)\right|\right|_{2}$
$\displaystyle\geq\frac{c_{h,q}(L)\cdot\widehat{C}_{k}}{2nh^{2}}\cdot||\widehat{\bm{y}}_{s+1}-\widehat{\bm{y}}_{s}||_{2}^{2},$

where we used (66) in the last strict inequality. Since
$\left\\{\widehat{f}_{h}(\widehat{\bm{y}}_{s})\right\\}_{s=0}^{\infty}$
converges by Theorem 8 as $s\to\infty$, we conclude that

$\lim_{s\to\infty}||\widehat{\bm{y}}_{s+1}-\widehat{\bm{y}}_{s}||_{2}^{2}=0\quad\text{
or equivalently,
}\quad\lim_{s\to\infty}\widehat{\bm{y}}_{s+1}^{T}\widehat{\bm{y}}_{s}=1.$ (67)

Now with the expression (25),

$\displaystyle\left|\left|\mathtt{Tang}\left(\nabla\widehat{f}_{h}(\widehat{\bm{y}}_{s})\right)\right|\right|_{2}^{2}$
$\displaystyle=\left|\left|\nabla\widehat{f}_{h}(\widehat{\bm{y}}_{s})-\widehat{\bm{y}}_{s}^{T}\nabla\widehat{f}_{h}(\widehat{\bm{y}}_{s})\cdot\widehat{\bm{y}}_{s}\right|\right|_{2}^{2}$
$\displaystyle=\left|\left|\nabla\widehat{f}_{h}(\widehat{\bm{y}}_{s})\right|\right|_{2}^{2}\cdot\left|\left|\widehat{\bm{y}}_{s+1}-\widehat{\bm{y}}_{s+1}^{T}\widehat{\bm{y}}_{s}\cdot\widehat{\bm{y}}_{s}\right|\right|_{2}^{2}$
$\displaystyle=\left|\left|\nabla\widehat{f}_{h}(\widehat{\bm{y}}_{s})\right|\right|_{2}^{2}\cdot\left[1-\left(\widehat{\bm{y}}_{s+1}^{T}\widehat{\bm{y}}_{s}\right)^{2}\right],$

where we plug in (23) in the second equality. As the function
$\bm{u}\mapsto\left|\left|\nabla\widehat{f}_{h}(\bm{u})\right|\right|_{2}^{2}$
is continuous on a compact set $\Omega_{q}$, it is upper bounded on
$\Omega_{q}$. As $s\to\infty$,
$\left|\left|\mathtt{Tang}\left(\nabla\widehat{f}_{h}(\widehat{\bm{y}}_{s})\right)\right|\right|_{2}\to
0$ by (67). Given that $\widehat{\bm{m}}_{k}$ is the unique point in
$\left\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\widehat{\bm{m}}_{k}\right|\right|_{2}\leq\widehat{r}_{k}\right\\}$
satisfying this
$\left|\left|\mathtt{Tang}\left(\nabla\widehat{f}_{h}(\bm{y})\right)\right|\right|_{2}=0$,
we conclude that $\widehat{\bm{y}}_{s}\to\widehat{\bm{m}}_{k}$ as $s\to\infty$
and
$\widehat{\bm{y}}_{0}\in\left\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\widehat{\bm{m}}_{k}\right|\right|_{2}\leq\widehat{r}_{k}\right\\}$.
Now with Lemma 10, we know that
$\widehat{\bm{m}}_{k}^{T}\nabla\widehat{f}_{h}(\widehat{\bm{m}}_{k})>0$ for
any $\widehat{\bm{m}}_{k}\in\widehat{\mathcal{M}}_{n}$ with probability
tending to 1 as $h\to 0$ and $nh^{q}\to\infty$. Therefore, as $h$ is small
enough and $n$ is sufficiently large, there exists a fixed constant $r^{*}>0$
such that $r^{*}\leq\min_{k}\widehat{r}_{k}$ with high probability. The
results follow.

###  D.6 Proof of Theorem 12

Before proving Theorem 12, we introduce the following two useful results. As
pointed out in Zhang and Sra (2016), a main hurdle in analyzing non-asymptotic
convergence of first-order methods on smooth manifolds is that the Euclidean
law of cosines does not hold. Fortunately, there is a trigonometric distance
bound stated below for Alexandrov space (Burago et al., 1992) with curvature
bounded below.

######  Lemma 24 (Lemma 5 in Zhang and Sra 2016; see also Bonnabel 2013)

If $a,b,c$ are the sides (that is, side lengths) of a geodesic triangle in an
Alexandrov space with sectional curvature (see Appendix B) lower bounded by
$\kappa$, and $A$ is the angle between sides $b$ and $c$, then

$a^{2}\leq\frac{\sqrt{|\kappa|}c}{\tanh(\sqrt{|\kappa|}c)}b^{2}+c^{2}-2bc\cos(A).$
(68)

The sketching proof of Lemma 24 can be founded in Lemma 5 of Zhang and Sra
(2016). Note that $\kappa=1$ on $\Omega_{q}$. We inherit the notation in Zhang
and Sra (2016) and denote $\frac{\sqrt{|\kappa|}c}{\tanh(\sqrt{|\kappa|}c)}$
by $\zeta(\kappa,c)$ for the curvature dependent quantity from inequality
(68). One can show by differentiating $\zeta(\kappa,c)$ with respect to $c$
that $\zeta(\kappa,c)$ is strictly increasing and greater than 1 for any $c>0$
and fixed $\kappa\neq 0$. With Lemma 24 in hand, we are able to state a
straightforward corollary indicating an important relation between two
consecutive updates of a gradient ascent algorithm on $\Omega_{q}$.

###### Corollary 25

For any point $\bm{x},\bm{y}_{s}$ in a convex set on $\Omega_{q}$, the update
in (30) satisfies

$2\eta\langle\mathtt{grad}\,f(\bm{y}_{s}),\mathtt{Exp}_{\bm{y}_{s}}^{-1}(\bm{x})\rangle\leq
d^{2}(\bm{y}_{s},\bm{x})-d^{2}(\bm{y}_{s+1},\bm{x})+\zeta(1,d(\bm{y}_{s},\bm{x}))\cdot\eta^{2}||\mathtt{grad}\,f(\bm{y}_{s})||_{2}^{2},$

recalling that
$d(\bm{x},\bm{y})=\sqrt{\langle\mathtt{Exp}_{\bm{x}}^{-1}(\bm{y}),\mathtt{Exp}_{\bm{x}}^{-1}(\bm{y})\rangle}=\left|\left|\mathtt{Exp}_{\bm{x}}^{-1}(\bm{y})\right|\right|_{2}$
on $\Omega_{q}$.

The proof is similar to Corollary 8 in Zhang and Sra (2016) and thus omitted.

###### Theorem 12

Assume (D1) and (M1).

  1. (a)

Linear convergence of gradient ascent with $f$: Given a convergence radius
$r_{0}$ with
$0<r_{0}\leq\sqrt{2-2\cos\left[\frac{3\lambda_{*}}{2(q+1)^{\frac{3}{2}}C_{3}}\right]}$,
the iterative sequence $\left\\{\bm{y}_{s}\right\\}_{s=0}^{\infty}$ defined by
the population-level gradient ascent algorithm (30) satisfies

$d(\bm{y}_{s},\bm{m}_{k})\leq\Upsilon^{s}\cdot
d(\bm{y}_{0},\bm{m}_{k})\quad\text{ with
}\quad\Upsilon=\sqrt{1-\frac{\eta\lambda_{*}}{2}},$

whenever
$\eta\leq\min\left\\{\frac{2}{\lambda_{*}},\frac{1}{(q+1)C_{3}\zeta(1,r_{0})}\right\\}$
and the initial point
$\bm{y}_{0}\in\left\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\bm{m}_{k}\right|\right|_{2}\leq
r_{0}\right\\}$ for some $\bm{m}_{k}\in\mathcal{M}$. We recall from Section
4.5 that $C_{3}$ is an upper bound for the derivatives of the directional
density $f$ up to the third order, $\lambda_{*}>0$ is defined in (M1), and
$\mathcal{M}$ is the set of local modes of the directional density $f$.

We further assume (D2’) and (K1) in the sequel.

  1. (b)

Linear convergence of gradient ascent with $\widehat{f}_{h}$: Let the sample-
based gradient ascent update on $\Omega_{q}$ be
$\widehat{\bm{y}}_{s+1}=\mathtt{Exp}_{\bm{y}_{s}}\left(\eta\cdot\mathtt{grad}\,\widehat{f}_{h}(\widehat{\bm{y}}_{s})\right)$.
With the same choice of the convergence radius $r_{0}>0$ and
$\Upsilon=\sqrt{1-\frac{\eta\lambda_{*}}{2}}$ as in (a), if $h\to 0$ and
$\frac{nh^{q+2}}{|\log h|}\to\infty$, then for any $\delta\in(0,1)$,

$d\left(\widehat{\bm{y}}_{s},\bm{m}_{k}\right)\leq\Upsilon^{s}\cdot
d\left(\widehat{\bm{y}}_{0},\bm{m}_{k}\right)+O(h^{2})+O_{P}\left(\sqrt{\frac{|\log
h|}{nh^{q+2}}}\right)$

with probability at least $1-\delta$, whenever
$\eta\leq\min\left\\{\frac{2}{\lambda_{*}},\frac{1}{(q+1)C_{3}\cdot\zeta(1,r_{0})}\right\\}$
and the initial point
$\widehat{\bm{y}}_{0}\in\left\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\bm{m}_{k}\right|\right|_{2}\leq
r_{0}\right\\}$ for some $\bm{m}_{k}\in\mathcal{M}$.

Proof  (a) Linear convergence of gradient ascent with $f$: The proof of the
linear convergence of the population-level gradient ascent algorithm (30) is
similar to some standard results in optimization theory, except that we are
under the manifold context now. Recall from (30) that the iterative formula
reads
$\bm{y}_{s+1}=\mathtt{Exp}_{\bm{y}_{s}}\left(\eta\cdot\mathtt{grad}\,f(\bm{y}_{s})\right)$
for $s=0,1,...$. We begin by deriving the following three facts.
$\bullet$ _Fact 1_ : Given (M1), $f$ is geodesically strongly concave around
some small neighborhoods of $\mathcal{M}$. In particular, when
$0<r_{0}\leq\sqrt{2-2\cos\left[\frac{3\lambda_{*}}{2(q+1)^{\frac{3}{2}}C_{3}}\right]}$,

$f(\bm{y})-f(\bm{m}_{k})-\langle\mathtt{grad}\,f(\bm{m}_{k}),\,\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})\rangle\leq-\frac{\lambda_{*}}{4}\left|\left|\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})\right|\right|_{2}^{2}$
(69)

for any
$\bm{y}\in\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\bm{m}_{k}\right|\right|_{2}\leq
r_{0}\\}$ and any $\bm{m}_{k}\in\mathcal{M}$.
$\bullet$ _Fact 2_. Given (D1) and (M1), we know that
$\left|\left|\mathtt{grad}\,f(\bm{x})\right|\right|_{2}\equiv\left|\left|\mathtt{Tang}\left(\nabla
f(\bm{x})\right)\right|\right|_{2}>0$ and

$f(\bm{m}_{k})-f\left(\mathtt{Exp}_{\bm{x}}\left(\frac{1}{(q+1)C_{3}}\mathtt{grad}\,f(\bm{x})\right)\right)\geq
0$

for all
$\bm{x}\in\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\bm{m}_{k}\right|\right|_{2}\leq
r_{0}\\}\setminus\\{\bm{m}_{k}\\}$ and any $\bm{m}_{k}\in\mathcal{M}$.
$\bullet$ _Fact 3_. Given (D1), the directional density $f$ is
$(q+1)C_{3}$-smooth.

As for _Fact 1_ , it follows from the differentiability of $f$ guaranteed by
(D1) and the eigenvalue condition (M1). By Taylor’s expansion on manifolds
(Pennec, 2006) and (M1),

$\displaystyle\begin{split}&f(\bm{y})-f(\bm{m}_{k})\\\
&=\langle\mathtt{grad}\,f(\bm{m}_{k}),\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})\rangle+\frac{1}{2}\cdot\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})^{T}\left[\mathcal{H}f(\bm{m}_{k})\right]\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})+o\left(\left|\left|\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})\right|\right|_{2}^{2}\right)\\\
&\leq\langle\mathtt{grad}\,f(\bm{m}_{k}),\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})\rangle-\frac{\lambda_{*}}{2}\left|\left|\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})\right|\right|_{2}^{2}+\frac{(q+1)^{\frac{3}{2}}C_{3}}{6}\cdot\left|\left|\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})\right|\right|_{2}^{3}\end{split}$
(70)

for any
$\bm{y}\in\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\bm{m}_{k}\right|\right|_{2}\leq
r_{0}\\}$ and $\bm{m}_{k}\in\mathcal{M}$. Since
$\left|\left|\bm{y}-\bm{m}_{k}\right|\right|_{2}\leq r_{0}$, the geodesic
distance between $\bm{y}$ and $\bm{m}_{k}$ satisfies
$d_{g}(\bm{y},\bm{m}_{k})=\left|\left|\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})\right|\right|_{2}=\arccos(\bm{y}^{T}\bm{m}_{k})\leq\frac{3\lambda_{*}}{2(q+1)^{\frac{3}{2}}C_{3}}$.
Plugging this result back into (70) yields that

$f(\bm{y})-f(\bm{m}_{k})\leq\langle\mathtt{grad}\,f(\bm{m}_{k}),\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})\rangle-\frac{\lambda_{*}}{4}\left|\left|\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})\right|\right|_{2}^{2}.$

For our purpose, it suffices to only prove (69) as above. One can shrink the
upper bound of the convergence radius $r_{0}>0$ so that the geodesically
strong concavity is valid for any pair of points within
$\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\bm{m}_{k}\right|\right|_{2}\leq
r_{0}\\}$. Indeed, the local strong concavity of $f$ is a natural consequence
of Morse Lemma (Lemma 3.11 in Banyaga and Hurtubise (2004)) given (M1).

_Fact 2_ is an obvious result under the eigenvalue condition (M1) and
differentiable condition (D1). This is because $\bm{m}_{k}$ is an unique local
mode of $f$ within the neighborhood
$\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\bm{m}_{k}\right|\right|_{2}\leq
r_{0}\\}$ and the geodesic distance between $\bm{x}$ and one-step gradient
ascent iteration from $\bm{x}$ with the step size $\frac{1}{(q+1)C_{3}}$
satisfies

$\displaystyle
d\left(\mathtt{Exp}_{\bm{x}}\left(\frac{1}{(q+1)C_{3}}\mathtt{grad}\,f(\bm{x})\right),\bm{x}\right)$
$\displaystyle=\frac{1}{(q+1)C_{3}}\left|\left|\mathtt{grad}\,f(\bm{x})\right|\right|_{2}$
$\displaystyle=\frac{1}{(q+1)C_{3}}\left|\left|\mathtt{grad}\,f(\bm{x})-\underbrace{\Gamma_{\bm{m}_{k}}^{\bm{x}}\left(\mathtt{grad}\,f(\bm{m}_{k})\right)}_{=0}\right|\right|_{2}$
$\displaystyle\leq\frac{1}{(q+1)C_{3}}\left|\left|\mathcal{H}f(\bm{x})\right|\right|_{2}\cdot\left|\left|\mathtt{Exp}_{\bm{x}}^{-1}(\bm{m}_{k})\right|\right|_{2}$
$\displaystyle\leq\left|\left|\mathtt{Exp}_{\bm{x}}^{-1}(\bm{m}_{k})\right|\right|_{2}=d_{g}(\bm{x},\bm{m}_{k}),$

where we use the fact that
$\left|\left|\mathcal{H}f(\bm{x})\right|\right|_{2}\leq(q+1)\left|\left|\mathcal{H}f(\bm{x})\right|\right|_{\max}\leq(q+1)C_{3}$
to deduce the last inequality. This shows that the one-step gradient ascent
iteration
$\mathtt{Exp}_{\bm{x}}\left(\frac{1}{(q+1)C_{3}}\cdot\mathtt{grad}\,f(\bm{x})\right)$
on $\Omega_{q}$ will stay within the neighborhood
$\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\bm{m}_{k}\right|\right|_{2}\leq
r_{0}\\}$ whenever
$\bm{x}\in\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\bm{m}_{k}\right|\right|_{2}\leq
r_{0}\\}\setminus\\{\bm{m}_{k}\\}$. Therefore,
$f(\bm{m}_{k})-f\left(\mathtt{Exp}_{\bm{x}}\left(\frac{1}{(q+1)C_{3}}\cdot\mathtt{grad}\,f(\bm{x})\right)\right)\geq
0$.

As for _Fact 3_ , note that
$\left|\left|\mathcal{H}f(\bm{x})\right|\right|_{\max}\leq C_{3}$ for all
$\bm{x}\in\Omega_{q}$. Thus,

$\displaystyle\left|\left|\mathtt{grad}\,f(\bm{x})-\Gamma_{\bm{y}}^{\bm{x}}\left(\mathtt{grad}\,f(\bm{y})\right)\right|\right|_{2}$
$\displaystyle=\left|\left|\left(\mathcal{H}f(\bm{x})\right)\mathtt{Exp}_{\bm{x}}^{-1}(\bm{x}^{*})\right|\right|_{2}$
$\displaystyle\leq\left|\left|\mathcal{H}f(\bm{x})\right|\right|_{2}\cdot\left|\left|\mathtt{Exp}_{\bm{x}}^{-1}(\bm{y})\right|\right|_{2}$
$\displaystyle\leq(q+1)C_{3}\left|\left|\mathtt{Exp}_{\bm{x}}^{-1}(\bm{y})\right|\right|_{2},$

where
$\bm{x}^{*}\in\left\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\bm{x}\right|\right|_{2}\leq\left|\left|\bm{y}-\bm{x}\right|\right|_{2}\right\\}$,
and we use the fact that $||A||_{2}\leq\sqrt{mn}||A||_{\max}$ for any
$A\in\mathbb{R}^{m\times n}$. See Section 3.3 in Genovese et al. (2014) or
Section 5.6 in Horn and Johnson (2012) for detailed relations between
different types of matrix norms.

With _Fact 1_ and _Fact 3_ , we have that

$\displaystyle\begin{split}&-\frac{(q+1)C_{3}}{2}\left|\left|\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})\right|\right|_{2}^{2}\leq
f(\bm{y})-f(\bm{m}_{k})-\langle\mathtt{grad}\,f(\bm{m}_{k}),\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})\rangle,\\\
&f(\bm{y})-f(\bm{m}_{k})-\langle\mathtt{grad}\,f(\bm{m}_{k}),\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})\rangle\leq-\frac{\lambda_{*}}{4}\left|\left|\mathtt{Exp}_{\bm{m}_{k}}^{-1}(\bm{y})\right|\right|_{2}^{2}<0\end{split}$
(71)

for any
$\bm{y}\in\left\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\bm{m}_{k}\right|\right|_{2}\leq
r_{0}\right\\}$ and $\bm{m}_{k}\in\mathcal{M}$.
Hence, given a point
$\bm{y}\in\left\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\bm{m}_{k}\right|\right|_{2}\leq
r_{0}\right\\}$ and using _Fact 2_ ,

$\displaystyle f(\bm{y})-f(\bm{m}_{k})$ $\displaystyle\leq
f(\bm{y})-f(\bm{m}_{k})+f(\bm{m}_{k})-f\left(\mathtt{Exp}_{\bm{y}}\left(\frac{1}{(q+1)C_{3}}\cdot\mathtt{grad}\,f(\bm{y})\right)\right)$
$\displaystyle=-\left[f\left(\mathtt{Exp}_{\bm{y}}\left(\frac{1}{(q+1)C_{3}}\cdot\mathtt{grad}\,f(\bm{y})\right)\right)-f(\bm{y})\right]$
$\displaystyle\leq-\left[\langle\mathtt{grad}\,f(\bm{y}),\frac{1}{(q+1)C_{3}}\mathtt{grad}\,f(\bm{y})\rangle-\frac{(q+1)C_{3}}{2}\left|\left|\mathtt{Exp}_{\bm{y}}^{-1}\left(\frac{1}{(q+1)C_{3}}\cdot\mathtt{grad}\,f(\bm{y})\right)\right|\right|_{2}^{2}\right]$
$\displaystyle=-\frac{1}{2(q+1)C_{3}}\left|\left|\mathtt{grad}\,f(\bm{y})\right|\right|_{2}^{2},$

where we apply the first inequality in (71) to obtain the fourth line. Thus,
for any
$\bm{x}\in\left\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\bm{m}_{k}\right|\right|_{2}\leq
r_{0}\right\\}$,

$\left|\left|\mathtt{grad}\,f(\bm{x})\right|\right|_{2}^{2}\leq
2(q+1)C_{3}\left[f(\bm{m}_{k})-f(\bm{x})\right].$ (72)

With
$\bm{y}_{0}\in\left\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\bm{m}_{k}\right|\right|_{2}\leq
r_{0}\right\\}$ and Corollary 25, we deduce that

$\displaystyle d^{2}(\bm{y}_{s+1},\bm{m}_{k})$ $\displaystyle\leq
d^{2}(\bm{y}_{s},\bm{m}_{k})-2\eta\langle\mathtt{grad}\,f(\bm{y}_{s}),\mathtt{Exp}_{\bm{y}_{s}}^{-1}(\bm{m}_{k})\rangle+\zeta(1,d(\bm{y}_{s},\bm{m}_{k}))\cdot\eta^{2}\left|\left|\mathtt{grad}\,f(\bm{y}_{s})\right|\right|_{2}^{2}$
$\displaystyle\stackrel{{\scriptstyle\text{(i)}}}{{\leq}}d^{2}(\bm{y}_{s},\bm{m}_{k})+2\eta\left[f(\bm{y}_{s})-f(\bm{m}_{k})-\frac{\lambda_{*}}{4}d^{2}(\bm{y}_{s},\bm{m}_{k})\right]$
$\displaystyle\quad+\zeta(1,r_{0})\cdot\eta^{2}\cdot
2(q+1)C_{3}\left[f(\bm{m}_{k})-f(\bm{y}_{s})\right]$
$\displaystyle=\left(1-\frac{\eta\lambda_{*}}{2}\right)d^{2}(\bm{y}_{s},\bm{m}_{k})-2\eta\left[1-\zeta(1,r_{0})(q+1)C_{3}\eta\right]\cdot\underbrace{\left[f(\bm{m}_{k})-f(\bm{y}_{s})\right]}_{\geq
0}$
$\displaystyle\leq\left(1-\frac{\eta\lambda_{*}}{2}\right)d^{2}(\bm{y}_{s},\bm{m}_{k})$

whenever
$\eta\leq\min\left\\{\frac{2}{\lambda_{*}},\frac{1}{(q+1)C_{3}\cdot\zeta(1,r_{0})}\right\\}$,
where we use the second inequality of (71), the monotonicity of $\zeta(1,c)$
with respect to $c$, and (72) to obtain (i). By telescoping, we conclude that
when
$\eta\leq\min\left\\{\frac{2}{\lambda_{*}},\frac{1}{(q+1)C_{3}\cdot\zeta(1,r_{0})}\right\\}$
and
$\bm{y}_{0}\in\left\\{\bm{z}\in\Omega_{q}:\left|\left|\bm{z}-\bm{m}_{k}\right|\right|_{2}\leq
r_{0}\right\\}$,

$d(\bm{y}_{s},\bm{m}_{k})=\left|\left|\mathtt{Exp}_{\bm{y}_{s}}^{-1}(\bm{m}_{k})\right|\right|_{2}\leq\left(1-\frac{\eta\lambda_{*}}{2}\right)^{\frac{s}{2}}\cdot