$|\lambda_{i}-\gamma_{i}|\leq|\lambda_{i}-\mu_{i}|+|\mu_{i}-\gamma_{i}|.$

The first term of the right-hand side of the above equation can be controlled
by Theorem 2. For the second term, since $\mathbf{K}_{n}^{*}$ and
$\widetilde{\mathcal{K}}$ have the same bandwidth chosen as (20) in the main
paper, we can use a sample splitting argument on $\mathsf{h}$ and apply Lemma
5 to obtain the final upper bound. The proof then follows from (79). ∎

###### Eigenfunction convergence..

By Proposition 1, one can apply the same argument leading to (B.1) above, to
replace the bandwidth $h_{n}$ in $\mathbf{K}_{n}$ by $\mathsf{h}$ defined in
(20) of the main paper, with an negligible error compared with the final
result. In the following, we focus on the setting where the bandwidth is
$\mathsf{h}$. As before, without loss of generality, we will consider the
reduced data $\\{\mathbf{y}_{i}^{0}\\}_{1\leq i\leq n}$ and their related
quantities $\\{\mathbf{x}_{i}^{0}\\}_{1\leq i\leq n}$ and
$\\{\mathbf{z}_{i}^{0}\\}_{1\leq i\leq n}$. For simplicity, we omit their
superscripts and denote them as $\\{\mathbf{y}_{i}\\}$, $\\{\mathbf{x}_{i}\\}$
and $\\{\mathbf{z}_{i}\\}$.

We denote $\mathbf{y}\in\mathbb{R}^{p}$ as

$\mathbf{y}=(\bm{y},\bm{y}^{\perp}),\ \bm{y}\in\mathbb{R}^{r}\ \text{and}\
\bm{y}^{\perp}\in\mathbb{R}^{p-r},$

and set $\mathbf{z}_{i}=(\bm{z}_{i},\bm{z}_{i}^{\perp})^{\top}$ in a similar
fashion. Consequently, by the fact that $\mathbf{x}_{i}=(\bm{x}_{i},0,...,0)$,
we have that

$\|\mathbf{y}-\mathbf{y}_{i}\|_{2}^{2}=\|\bm{y}-\bm{x}_{i}\|_{2}^{2}+\|\mathbf{z}_{i}\|_{2}^{2}+\|\bm{y}^{\perp}\|_{2}^{2}-2(\bm{y}^{\perp})^{\top}\bm{z}_{i}^{\perp}-2(\bm{y}-\bm{x}_{i})^{\top}\bm{z}_{i}.$
(84)

This leads to

$\displaystyle\exp\left(-\frac{\|\mathbf{y}-\mathbf{y}_{i}\|_{2}^{2}}{\mathsf{h}}\right)=\exp\left(-\frac{\|\bm{y}-\bm{x}_{i}\|_{2}^{2}}{\mathsf{h}}\right)\cdot\ell,$

where $\ell$ is denoted as

$\ell:=\exp\left(-\frac{\|\mathbf{z}_{i}\|_{2}^{2}+\|\bm{y}^{\perp}\|_{2}^{2}-2(\bm{y}^{\perp})^{\top}\bm{z}_{i}^{\perp}-2(\bm{y}-\bm{x}_{i})^{\top}\bm{z}_{i}}{\mathsf{h}}\right).$

As we will examine the RKHS norm in $\mathcal{H}_{K}$, containing functions
defined only on the support of $\tilde{\mathsf{P}}$, without loss of
generality, we can ignore $\bm{y}^{\perp}$ (that is, set $\bm{y}^{\perp}=0$)
and only focus on $\bm{y}.$ Since $\bm{y}$ is sub-Gaussian, using the
elementary estimate $\exp(-x)\sim 1-x$ as $x\to 0,$ we obtain

$\ell=1+{\mathrm{O}_{\prec}\bigg{(}\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}+\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}\bigg{)},}$

uniformly in $\mathbf{y}$ on the support of $\widetilde{\mathsf{P}}$. This
immediately shows that

$\displaystyle\exp\left(-\frac{\|\mathbf{y}-\mathbf{y}_{i}\|_{2}^{2}}{\mathsf{h}}\right)$
$\displaystyle=\exp\left(-\frac{\|\bm{y}-\bm{x}_{i}\|_{2}^{2}}{\mathsf{h}}\right)$
$\displaystyle+{\mathrm{O}_{\prec}\bigg{(}\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}+\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}\bigg{)}.}$
(85)

With the above preparation, we start analyzing the eigenfunctions. Recall the
definition in (48). In the first step, we study the following term

$\left(\widehat{\phi}_{i}^{(n)}(\mathbf{y})-\phi_{i}^{(n)}(\bm{y})\right)^{2}=\left(\mathrm{E}_{1}+\mathrm{E}_{2}+\mathrm{E}_{3}\right)^{2},$
(86)

where $\mathrm{E}_{k},k=1,2,3,$ are denoted as follows

$\displaystyle\mathrm{E}_{1}:=\left(\frac{1}{\lambda_{i}}-\frac{1}{\mu_{i}}\right)\frac{1}{\sqrt{n}}\sum_{j=1}^{n}\exp\left(-\frac{\|\mathbf{y}-\mathbf{y}_{j}\|_{2}^{2}}{\mathsf{h}}\right)v_{ij},$
$\displaystyle\mathrm{E}_{2}:=\frac{1}{\sqrt{n}\mu_{i}}\sum_{j=1}^{n}\left[\exp\left(-\frac{\|\mathbf{y}-\mathbf{y}_{j}\|_{2}^{2}}{\mathsf{h}}\right)-\exp\left(-\frac{\|\bm{y}-\bm{x}_{j}\|_{2}^{2}}{\mathsf{h}}\right)\right]v_{ij},$
$\displaystyle\mathrm{E}_{3}:=\frac{1}{\sqrt{n}\mu_{i}}\sum_{j=1}^{n}\exp\left(-\frac{\|\bm{y}-\bm{x}_{j}\|_{2}^{2}}{\mathsf{h}}\right)(v_{ij}-u_{ij}).$

First, by Theorem 2, when $1\leq i\leq\mathsf{K},$ we have that
$\lambda_{i}\geq\gamma_{i}-|\lambda_{i}-\gamma_{i}|\geq c$ for sufficiently
large $n$ with high probability, and thus

$\left|\frac{1}{\lambda_{i}}-\frac{1}{\mu_{i}}\right|\prec\frac{|\lambda_{i}-\mu_{i}|}{\lambda_{i}\mu_{i}}\prec{\frac{1}{\mu_{i}}\bigg{(}\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}+\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}\bigg{)}.}$

Consequently, using the boundedness of the kernels and the fact that
$\|\mathbf{v}\|_{1}\leq\sqrt{n}\|\mathbf{v}\|_{2}=\sqrt{n}$ for any unit
vector $\mathbf{v}\in\mathbb{R}^{n}$, we conclude that for $i\leq\mathsf{K}$

$\mathrm{E}_{1}\prec{\frac{1}{\mu_{i}}\bigg{(}\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}+\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}\bigg{)}.}$

Second, by (B.2), we readily have that for $i\leq\mathsf{K}$

$\mathrm{E}_{2}\prec{\frac{1}{\mu_{i}}\bigg{(}\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}+\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}\bigg{)}.}$

Third, we control $\mathrm{E}_{3}.$ On one hand, for $i\leq\mathsf{K},$ by
Lemma 5, we have that $\mathrm{E}_{3}=\mathrm{O}(1).$ One the other hand, by
setting

$\mathbf{a}=(a_{1},\cdots,a_{n})^{\top},\
\mathbf{b}=(b_{1},\cdots,b_{n})^{\top},$

where
$a_{j}=\frac{1}{\sqrt{n}\mu_{i}}\exp\left(-\frac{\|\bm{y}-\bm{x}_{j}\|_{2}^{2}}{\mathsf{h}}\right)$
and $b_{j}=v_{ij}-u_{ij},$ we can rewrite

$\mathrm{E}_{3}^{2}=\mathbf{b}^{\top}\mathbf{a}\mathbf{a}^{\top}\mathbf{b}.$

Note that $\|\mathbf{b}\|_{2}^{2}=\|\mathbf{v}_{i}-\mathbf{u}_{i}\|_{2}^{2}.$
Then we have that

$\mathrm{E}_{3}^{2}\leq\|\mathbf{a}\|_{2}^{2}\|\mathbf{b}\|_{2}^{2}\prec{\frac{1}{\mu_{i}^{2}}\cdot\frac{1}{\mathsf{r}_{i}^{2}}\left(\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{\alpha_{i}-\beta})^{1/2}}\right),}$
(87)

where in the second step we used Part 2 of Theorem 2. Combining all the above
bounds, using (86), when $i\leq\mathsf{K},$ we arrive at

$\displaystyle\mu_{i}\left(\widehat{\phi}_{i}^{(n)}(\mathbf{y})-\phi_{i}^{(n)}(\bm{y})\right)^{2}$
$\displaystyle\prec{\bigg{[}\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}+\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}\bigg{]}^{2}+\frac{1}{\mathsf{r}_{i}^{2}}\left(\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{\alpha_{i}-\beta})^{1/2}}\right),}$
(88)

for any $\mathbf{y}$ on the support of $\widetilde{\mathsf{P}}$. Finally,
since by Part 1 of Theorem 2, we have for any $\mathbf{y}\in\mathcal{S},$

$\displaystyle\sqrt{\lambda_{i}}\phi_{i}^{(n)}(\mathbf{y})-\sqrt{\mu_{i}}\phi_{i}^{(n)}(\mathbf{y})\leq(\sqrt{\lambda_{i}}-\sqrt{\mu_{i}})\phi_{i}^{(n)}(\mathbf{y})\prec{\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}+\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}},}$
(89)

and similarly

$\sqrt{\lambda_{i}}\widehat{\phi}_{i}^{(n)}(\mathbf{y})-\sqrt{\mu_{i}}\widehat{\phi}_{i}^{(n)}(\mathbf{y})\prec{\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}+\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}.}$
(90)

Combining the above bounds, we have that

$\displaystyle\sqrt{\lambda_{i}}\widehat{\phi}_{i}^{(n)}(\mathbf{y})-\sqrt{\gamma_{i}}{\widetilde{\phi}_{i}}(\mathbf{y})$
$\displaystyle\leq\sqrt{\lambda_{i}}\widehat{\phi}_{i}^{(n)}(\mathbf{y})-\sqrt{\lambda_{i}}\phi_{i}^{(n)}(\mathbf{y})+\sqrt{\lambda_{i}}\phi_{i}^{(n)}(\mathbf{y})-\sqrt{\mu_{i}}\phi_{i}^{(n)}(\mathbf{y})+\sqrt{\mu_{i}}\phi_{i}^{(n)}-\sqrt{\gamma}_{i}\widetilde{\phi}_{i}(\mathbf{y})$
$\displaystyle\leq|\sqrt{\mu_{i}}\widehat{\phi}_{i}^{(n)}(\mathbf{y})-\sqrt{\mu_{i}}\phi_{i}^{(n)}(\mathbf{y})|+|\sqrt{\lambda_{i}}\widehat{\phi}_{i}^{(n)}(\mathbf{y})-\sqrt{\mu_{i}}\widehat{\phi}_{i}^{(n)}(\mathbf{y})|+2|\sqrt{\lambda_{i}}\phi_{i}^{(n)}-\sqrt{\mu_{i}}\phi_{i}^{(n)}(\mathbf{y})|$
$\displaystyle\quad+|\sqrt{\mu_{i}}\phi_{i}^{(n)}(\mathbf{y})-\sqrt{\gamma}_{i}\widetilde{\phi}_{i}(\mathbf{y})|$
$\displaystyle\prec{|\sqrt{\mu_{i}}\widehat{\phi}_{i}^{(n)}(\mathbf{y})-\sqrt{\mu_{i}}\phi_{i}^{(n)}(\mathbf{y})|+\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}+\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{\sqrt{n}\mathsf{r}_{i}},}$
(91)

where in the last inequality we used

$|f(\mathbf{y})|=|\langle
f,K(\mathbf{y},\cdot)\rangle_{K}|\leq\|f\|_{K}\sqrt{K(\mathbf{y},\mathbf{y})}\leq
C\|f\|_{K},$ (92)

and Lemma 5. Finally, plugging in (B.2), we obtain the final result as
$i\leq\mathsf{K}$. ∎

###  B.3 Guarantee of Bandwith Selection: Proof of Proposition 1

We first obtain upper and lower bounds for $\mathsf{h}$ defined in (20).
Recall from (63) that

$\max_{1\leq i,j\leq
n}\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}^{2}\prec\sum_{i=1}^{r}\theta_{i}.$
(93)

This implies the upper bound that

$\mathsf{h}\prec\sum_{i=1}^{r}\theta_{i}.$ (94)

For the lower bound, we show that, with high probability, for any sufficiently
small constant $\epsilon>0$, there exist at most
$\frac{n^{2}}{\log^{\epsilon}n}$ pairs of $(\mathbf{x}_{i},\mathbf{x}_{j})$ so
that

$\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}^{2}\leq\frac{C}{\log^{4\epsilon}n}\sum_{i=1}^{r}\theta_{i},$
(95)

where $C>0$ is some universal constant. To see this, we define the events

$E(i,j)=\bigg{\\{}\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}^{2}\leq\frac{C}{\log^{4\epsilon}n}\sum_{i=1}^{r}\theta_{i}\bigg{\\}},\
\ 1\leq i\neq j\leq n.$

We need the following lemma, whose proof will be given at the end of this
section.

######  Lemma 15.

For any given $i\in\\{1,2,...,n\\}$, with high probability, at most
$\frac{n}{\log^{\epsilon}n}$ events in $\\{E(i,j):1\leq j\leq n,j\neq i\\}$
are true.

By Lemma 15, we can apply the union bound over all $1\leq i\leq n$, and
conclude that, with high probability, at most $\frac{n^{2}}{\log^{\epsilon}n}$
pairs of $(\mathbf{x}_{i},\mathbf{x}_{j})$ satisfy (95). Now for any
percentile parameter $\omega\in(0,1)$ that is independent of $n$, when $n$ is
sufficiently large, we always have
$n(n-1)\omega>\frac{n^{2}}{\log^{\epsilon}n}$, which means the signal
bandwidth $\mathsf{h}$ selected according to (20) satisfies

$\mathsf{h}\geq\frac{C}{\log^{4\epsilon}n}\sum_{i=1}^{r}\theta_{i},$

with high probability. This leads to the lower bound

$\mathsf{h}\succ\sum_{i=1}^{r}\theta_{i}.$ (96)

Next we control $|\mathsf{h}/h_{n}-1|$. Note that

$\|\mathbf{y}_{i}-\mathbf{y}_{j}\|_{2}^{2}=\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}^{2}+\|\mathbf{z}_{i}-\mathbf{z}_{j}\|_{2}^{2}-2(\mathbf{x}_{i}-\mathbf{x}_{j})^{\top}(\mathbf{z}_{i}-\mathbf{z}_{j}).$
(97)

By Lemma 10 and Assumption 2, we have that

$\max_{1\leq i,j\leq
n}\|\mathbf{z}_{i}-\mathbf{z}_{j}\|_{2}^{2}\prec{\sigma^{2}p\asymp
n^{\beta+\eta}.}$

By an argument similar to that leading to (64) and a union bound, we have

$\max_{1\leq i,j\leq
n}|2(\mathbf{x}_{i}-\mathbf{x}_{j})^{\top}(\mathbf{z}_{i}-\mathbf{z}_{j})|\prec\sigma\bigg{(}\sum_{i=1}^{r}\theta_{i}\bigg{)}^{1/2}\asymp\bigg{(}\sum_{i=1}^{r}n^{\alpha_{i}+\beta}\bigg{)}^{1/2}.$
(98)

Combining the above arguments, we find that for any $(i,j),i\neq j$

$\displaystyle\big{|}\|\mathbf{y}_{i}-\mathbf{y}_{j}\|_{2}^{2}-\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}^{2}\big{|}$
$\displaystyle=\big{|}\|\mathbf{z}_{i}-\mathbf{z}_{j}\|_{2}^{2}-2(\mathbf{x}_{i}-\mathbf{x}_{j})^{\top}(\mathbf{z}_{i}-\mathbf{z}_{j})\big{|}$
$\displaystyle\prec{n^{\beta+\eta}+n^{\beta/2}\big{(}\sum_{i=1}^{r}\theta_{i}\big{)}^{1/2}.}$
(99)

Moreover, for any $k$ such that $k=\omega n$ for some $\omega\in(0,1).$
Without loss of generality, we assume that $k$ is an integer. Let $X_{(k)}$ be
the $k$-th largest element in
$\\{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|^{2}_{2}:i\neq j\\}$, then by Lemma 15,
we have that

$X_{(k)}\geq\frac{C}{\log^{4\epsilon}n}\sum_{i=1}^{r}\theta_{i}.$ (100)

Together with (99) and (97), we find that there are at least $k$ elements in
$\\{\|\mathbf{y}_{i}-\mathbf{y}_{j}\|^{2}_{2}:i\neq j\\}$ satisfying

$\|\mathbf{y}_{i}-\mathbf{y}_{j}\|_{2}^{2}\geq
X_{(k)}-{\mathrm{O}_{\prec}\bigg{(}n^{\beta+\eta}+n^{\beta/2}\big{(}\sum_{i=1}^{r}\theta_{i}\big{)}^{1/2}\bigg{)}.}$

which follows from the fact that the first $k$ elements $\\{X_{(i)}:1\leq
i\leq k\\}$ are no smaller than $X_{(k)}$ and hence they do not satisfy (95).
As a result, the $k$-th largest element $Y_{(k)}$ in
$\\{\|\mathbf{y}_{i}-\mathbf{y}_{j}\|^{2}_{2}:i\neq j\\}$ satisfies

$Y_{(k)}\geq
X_{(k)}-{\mathrm{O}_{\prec}\bigg{(}n^{\beta+\eta}+n^{\beta/2}\big{(}\sum_{i=1}^{r}\theta_{i}\big{)}^{1/2}\bigg{)}.}$
(101)

On the other hand, we claim that there are at least $n-k+1$ elements in
$\\{\|\mathbf{y}_{i}-\mathbf{y}_{j}\|_{2}:i\neq j\\}$ satisfying

$\|\mathbf{y}_{i}-\mathbf{y}_{j}\|_{2}^{2}\leq
X_{(k)}+{\mathrm{O}_{\prec}\bigg{(}n^{\beta+\eta}+n^{\beta/2}\big{(}\sum_{i=1}^{r}\theta_{i}\big{)}^{1/2}\bigg{)}.}$

This is because Lemma 15 and (99) imply that the last $n-k+1$ elements
$\\{X_{(i)}:k\leq i\leq n\\}$ are no greater than $X_{(k)}$. As a result, we
also have

$Y_{(k)}\leq
X_{(k)}+{\mathrm{O}_{\prec}\bigg{(}n^{\beta+\eta}+n^{\beta/2}\big{(}\sum_{i=1}^{r}\theta_{i}\big{)}^{1/2}\bigg{)}.}$
(102)

Hence, using the definitions of $\mathsf{h}$ and $h_{n},$ by (101), (102) and
(100), we have that

$\frac{\mathsf{h}}{h_{n}}=1+{\mathrm{O}_{\prec}\bigg{(}\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{\alpha_{i}-\beta})^{1/2}}\bigg{)}\equiv
1+\mathrm{O}_{\prec}(\psi_{n}).}$ (103)

Finally, we prove Lemma 15.

###### Proof of Lemma 15..

Define the indicator function $1_{E(i,j)}$ so that $1_{E(i,j)}=1$ if $E(i,j)$
is true and $1_{E(i,j)}=0$ otherwise. We shall prove that, for any fixed
$i\in\\{1,2,...,n\\}$, we have

$\mathbb{P}\bigg{(}\sum_{1\leq j\leq n,j\neq
i}1_{E(i,j)}\leq\frac{n}{\log^{\epsilon}n}\bigg{)}\geq 1-n^{-C},$ (104)

for any large constant $C>0$ and sufficiently large $n$. Note that

$\displaystyle\mathbb{P}\bigg{(}\sum_{1\leq j\leq n,j\neq
i}1_{E(i,j)}\leq\frac{n}{\log^{\epsilon}n}\bigg{)}$
$\displaystyle=\int\mathbb{P}\bigg{(}\sum_{1\leq j\leq n,j\neq
i}1_{E(i,j)}\leq\frac{n}{\log^{\epsilon}n}\bigg{|}\mathbf{x}_{i}\bigg{)}\mathrm{d}\mathbb{P}_{\mathbf{x}_{i}}.$

It suffices to show that, for any
$\mathbf{x}_{i}\in\text{supp}(\widetilde{\mathsf{P}})$, we have

$\mathbb{P}\bigg{(}\sum_{1\leq j\leq n,j\neq
i}1_{E(i,j)}\leq\frac{n}{\log^{\epsilon}n}\bigg{|}\mathbf{x}_{i}\bigg{)}\geq
1-n^{-C}.$ (105)

Now since conditional on $\mathbf{x}_{i}$, the events $\\{E(i,j):1\leq j\leq
n,j\neq i\\}$ are mutually independent, so that $\\{1_{E(i,j)}:1\leq j\leq
n,j\neq i\\}$ are $i.i.d.$ Bernoulli random variables, with mean

$\mu_{\mathbf{x}_{i}}=\mathbb{P}(E(i,j)|\mathbf{x}_{i})=\mathbb{E}[1_{E(i,j)}|\mathbf{x}_{i}].$

Note that when

$\frac{C\log n}{n}\leq\mu_{\mathbf{x}_{i}}\leq\frac{1}{\log^{\epsilon}n},$
(106)

using Lemma 8 with $t=en\mu_{\mathbf{x}_{i}}$ leads to

$\mathbb{P}\bigg{(}\sum_{1\leq j\leq n,j\neq i}1_{E(i,j)}\geq
en\mu_{\mathbf{x}_{i}}\bigg{|}\mathbf{x}_{i}\bigg{)}\leq\exp(-n\mu_{\mathbf{x}_{i}}),$

which implies (105) and therefore (104). Moreover, if

$0\leq\mu_{\mathbf{x}_{i}}\leq\frac{C\log n}{n},$

using Lemma 8 with $t=\frac{n}{\log^{\epsilon}n}$ leads to

$\mathbb{P}\bigg{(}\sum_{1\leq j\leq n,j\neq
i}1_{E(i,j)}\geq\frac{n}{\log^{\epsilon}n}\bigg{|}\mathbf{x}_{i}\bigg{)}\leq\exp(-n\mu_{\mathbf{x}_{i}})(e\mu_{\mathbf{x}_{i}}\log^{\epsilon}n)^{n/\log^{\epsilon}n}\leq
n^{-C},$

for any large constant $C>0$ and sufficiently large $n$. As a result, we only
need to show that

$\mu_{\mathbf{x}_{i}}\leq\frac{1}{\log^{\epsilon}n}.$ (107)

The proof of (107) is separated into two cases as follows.

Case I: $r=\mathrm{o}(\log^{\epsilon}n).$ For each $j=\\{1,2,...,n\\}$, we
have

$\displaystyle\mathbb{P}\bigg{(}\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}^{2}\geq\frac{C}{\log^{4\epsilon}n}\sum_{k=1}^{r}\theta_{i}\bigg{|}\mathbf{x}_{i}\bigg{)}$
$\displaystyle\geq\mathbb{P}\bigg{(}\bigcap_{k=1}^{r}\bigg{\\{}(x_{ik}-x_{jk})^{2}\geq\frac{C\theta_{k}}{\log^{4\epsilon}n}\bigg{|}\mathbf{x}_{i}\bigg{\\}}\bigg{)}$
$\displaystyle\geq
1-\mathbb{P}\bigg{(}\bigcup_{k=1}^{r}\bigg{\\{}(x_{ik}-x_{jk})^{2}<\frac{C\theta_{k}}{\log^{4\epsilon}n}\bigg{\\}}\bigg{|}\mathbf{x}_{i}\bigg{)}$
$\displaystyle\geq
1-\sum_{i=1}^{r}\mathbb{P}\bigg{(}(x_{ik}-x_{jk})^{2}<\frac{C\theta_{k}}{\log^{4\epsilon}n}\bigg{|}\mathbf{x}_{i}\bigg{)}.$

Now under the assumption that, for $1\leq k\leq r$, the marginal probability
density of $x_{jk}/\sqrt{\theta_{k}}$ is continuous, bounded, and strictly
bounded away from zero on
$\text{supp}(\widetilde{\mathsf{P}}_{k})/\sqrt{\theta_{k}}$, where
$\widetilde{\mathsf{P}}_{k}$ is the marginal distribution of
$x_{jk},j=1,2,...,n$, by the mean value theorem, we have

$\mathbb{P}\bigg{(}|x_{ik}-x_{jk}|\leq\frac{\sqrt{\theta_{k}}}{\log^{2\epsilon}n}\bigg{|}\mathbf{x}_{i}\bigg{)}\asymp\frac{1}{\log^{2\epsilon}n},\qquad
k=1,2,...,r.$ (108)

Thus, we have

$\mathbb{P}\bigg{(}\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}^{2}\geq\frac{C}{\log^{4\epsilon}n}\sum_{i=1}^{r}\theta_{i}\bigg{|}\mathbf{x}_{i}\bigg{)}\geq
1-\frac{Cr}{\log^{2\epsilon}n}.$

So that (107) holds whenever $r=\mathrm{o}(\log^{\epsilon}n)$.

Case II: $r\gtrsim\log^{\epsilon}n$. The argument depends on the magnitude of
the ratio $\sum_{i=1}^{r}\theta_{i}/\theta_{1}$. On the one hand, when
$\sum_{i=1}^{r}\theta_{i}/\theta_{1}\geq c_{1}\log\log^{\epsilon}n$ for some
large constant $c_{1}>0$, since $(x_{ik}-x_{jk})^{2},1\leq k\leq n$ are i.i.d.
sub-exponential random variables, we can apply Lemma 9 to obtain

$\displaystyle\mathbb{P}\bigg{(}\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}^{2}\leq\frac{C}{\log^{4\epsilon}n}\sum_{k=1}^{r}\theta_{i}\bigg{|}\mathbf{x}_{i}\bigg{)}$
$\displaystyle\leq
2\exp\bigg{(}-c\min\bigg{\\{}\frac{(\sum_{i=1}^{r}\theta_{i})^{2}+\|\mathbf{x}_{i}\|_{2}^{4}}{\sum_{i=1}^{r}\theta_{i}^{2}},\frac{\sum_{i=1}^{r}\theta_{i}+\|\mathbf{x}_{i}\|_{2}^{2}}{\theta_{1}}\bigg{\\}}\bigg{)}$
$\displaystyle\leq
2\exp\bigg{(}-c\min\bigg{\\{}\frac{(\sum_{i=1}^{r}\theta_{i})^{2}}{\sum_{i=1}^{r}\theta_{i}^{2}},\frac{\sum_{i=1}^{r}\theta_{i}}{\theta_{1}}\bigg{\\}}\bigg{)}$
$\displaystyle\leq
2\exp\bigg{(}-c\frac{\sum_{i=1}^{r}\theta_{i}}{\theta_{1}}\bigg{)}\leq\log^{-\epsilon
c_{2}}n,$

for some $c_{2}>1$ and sufficiently large $r$, where the third step follows
from the Hölder’s inequality
$\sum_{i=1}^{r}\theta_{i}^{2}\leq\theta_{1}\cdot\sum_{i=1}^{r}\theta_{i}$.
This implies (107) for sufficiently large $n$. On the other hand when
$\sum_{i=1}^{r}\theta_{i}/\theta_{1}<c_{1}\log\log^{\epsilon}n$, we have

$\displaystyle\mathbb{P}\bigg{(}\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}^{2}$
$\displaystyle\leq\frac{C}{\log^{4\epsilon}n}\sum_{k=1}^{r}\theta_{i}\bigg{|}\mathbf{x}_{i}\bigg{)}$
$\displaystyle\leq\mathbb{P}\bigg{(}\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}^{2}\leq\frac{C_{2}\log\log^{\epsilon}n}{\log^{4\epsilon}n}\theta_{1}\bigg{|}\mathbf{x}_{i}\bigg{)}$
$\displaystyle\leq\mathbb{P}\bigg{(}(x_{i1}-x_{j1})^{2}\leq\frac{C_{2}\theta_{1}}{\log^{3\epsilon}n}\bigg{|}\mathbf{x}_{i}\bigg{)}\lesssim\frac{1}{\log^{3\epsilon/2}n},$

where the last inequality follows from (108). This again implies (107) for
sufficiently large $n$. This proves (107) and therefore completes the proof of
Lemma 15. ∎

###  B.4 Extension to General Kernel Functions: Proof of Theorem 4

Similar as in the proof of Theorems 2 and 3, without loss of generality, we
consider the reduced random vectors $\\{\mathbf{y}_{i}^{0}\\}_{1\leq i\leq
n}$, $\\{\mathbf{x}_{i}^{0}\\}_{1\leq i\leq n}$ and
$\\{\mathbf{z}_{i}^{0}\\}_{1\leq i\leq n}$, and denote them for simplicity as
$\\{\mathbf{y}_{i}\\}$, $\\{\mathbf{x}_{i}\\}$ and $\\{\mathbf{z}_{i}\\}$.
Without loss of generality, we assume they are centered. Again, we denote the
kernel matrices for ${\mathbf{y}_{i}},{\mathbf{x}_{i}},$ and
${\mathbf{z}_{i}},$ with bandwidth $\mathsf{h}$ as
$\mathbf{K}_{y},\mathbf{K}_{x}$ and $\mathbf{K}_{z},$ respectively.

####  B.4.1 Convergence to Noiseless Kernel Matrix

###### Eigenvalue convergence..

Firstly, we study $n^{-1}\|\mathbf{K}_{y}-\mathbf{K}_{x}\|$. Our discussion
can be separated into three cases as follows.

Case I: $\inf_{0<x<1}\nu(x)\geq\tau_{0}$. By (iii) of Assumption 3, for any
$1\leq i,j\leq n$ and $i\neq j$, if we denote
$x=\frac{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}}{\mathsf{h}^{1/2}}$ and
$y=\frac{\|\mathbf{y}_{i}-\mathbf{y}_{j}\|_{2}}{\mathsf{h}^{1/2}}$, we have

$\displaystyle\mathbf{K}_{y}(i,j)$ $\displaystyle=f(y)\leq
f(x)+C\\{x^{\nu(x)}+y^{\nu(y)}\\}\cdot|x-y|^{\tau_{0}}$ (109)
$\displaystyle\leq
f(x)+C\\{x^{\nu(x)}+y^{\nu(y)}\\}\cdot\bigg{|}\frac{x^{2}-y^{2}}{x+y}\bigg{|}^{\tau_{0}}$
$\displaystyle\leq
f(x)+C\\{x^{\nu(x)-\tau_{0}}+y^{\nu(y)-\tau_{0}}\\}\cdot|{x^{2}-y^{2}}|^{\tau_{0}}$
$\displaystyle\leq
f(x)+\mathrm{O}_{\prec}\bigg{(}\bigg{|}\frac{\|\mathbf{y}_{i}-\mathbf{y}_{j}\|_{2}^{2}-\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}^{2}}{\mathsf{h}}\bigg{|}^{\tau_{0}}\bigg{)},$

where in the last step we used the assumption $\nu(x)\geq\tau_{0}$ for all
$0<x<1$, and $\nu(x)<0$ for all $x>1$, and the fact

$\|\mathbf{y}_{i}-\mathbf{y}_{j}\|_{2}\prec\mathsf{h}^{1/2},\quad\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}\prec\mathsf{h}^{1/2}.$

This shows that

$|\mathbf{K}_{y}(i,j)-\mathbf{K}_{x}(i,j)|\prec\bigg{(}\frac{\|\mathbf{z}_{i}-\mathbf{z}_{j}\|_{2}^{2}+2|(\mathbf{x}_{i}-\mathbf{x}_{j})^{\top}(\mathbf{z}_{i}-\mathbf{z}_{j})|}{\mathsf{h}}\bigg{)}^{\tau}.$

By (57), we have

$\frac{\|\mathbf{z}_{i}-\mathbf{z}_{j}\|_{2}^{2}}{\mathsf{h}}\prec{\frac{\sigma^{2}p}{\mathsf{h}}\prec\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}.}$

This along with (67) yields

$\max_{1\leq i,j\leq
n}|\mathbf{K}_{y}(i,j)-\mathbf{K}_{x}(i,j)|\prec{\bigg{(}\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}\bigg{)}^{\tau_{0}}.}$

Using a discussion similar to (B.1) with Lemma 13, we find that

$\displaystyle\|n^{-1}\mathbf{K}_{y}-n^{-1}\mathbf{K}_{x}\|$
$\displaystyle\leq\max_{1\leq i,j\leq
n}|\mathbf{K}_{y}(i,j)-\mathbf{K}_{x}(i,j)|$
$\displaystyle\prec{\bigg{(}\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}\bigg{)}^{\tau_{0}}}.$
(110)

Case II: $0\leq\inf_{0<x<1}\nu(x)<\tau_{0}$. In this case, we denote
$\nu_{0}=\inf_{0<x<1}\nu(x)$. Similar to the discussion of (109), we have

$\displaystyle\mathbf{K}_{y}(i,j)$ $\displaystyle=f(y)\leq
f(x)+C\\{x^{\nu(x)}+y^{\nu(y)}\\}\cdot\frac{|x^{2}-y^{2}|^{\nu_{0}}}{|x+y|^{\nu_{0}}}|x-y|^{\tau_{0}-{\nu_{0}}}$
$\displaystyle\leq
f(x)+C\\{x^{\nu(x)-\nu_{0}}+y^{\nu(y)-\nu_{0}}\\}\cdot{|x^{2}-y^{2}|^{\nu_{0}}}|x-y|^{\tau_{0}-{\nu_{0}}}$
$\displaystyle\leq
f(x)+C\\{x^{\nu(x)-\nu_{0}}+y^{\nu(y)-\nu_{0}}\\}\cdot{|x^{2}-y^{2}|^{\frac{\tau_{0}+\nu_{0}}{2}}}$
$\displaystyle\leq
f\bigg{(}\frac{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}}{\mathsf{h}^{1/2}}\bigg{)}+\mathrm{O}_{\prec}\bigg{(}\bigg{|}\frac{\|\mathbf{y}_{i}-\mathbf{y}_{j}\|_{2}^{2}-\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}^{2}}{\mathsf{h}}\bigg{|}^{\frac{\tau_{0}+\nu_{0}}{2}}\bigg{)},$

where in the third step we used the inequality

$x-y\leq\sqrt{x^{2}-y^{2}},\ \text{for all}\ x\geq y\geq 0.$

and in the last step we used the assumption that $\nu(x)>\nu_{0}$ for $0<x<1$
and $\nu(x)<0$ for $x\geq 1$. This leads to

$\displaystyle\|n^{-1}\mathbf{K}_{y}-n^{-1}\mathbf{K}_{x}\|$
$\displaystyle\leq\max_{1\leq i,j\leq
n}|\mathbf{K}_{y}(i,j)-\mathbf{K}_{x}(i,j)|$
$\displaystyle\prec{\bigg{(}\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}\bigg{)}^{\frac{\tau_{0}+\nu_{0}}{2}}}.$
(111)

Next, we study $\|\mathbf{K}_{y}-\mathbf{K}_{n}\|$. Note that for some $h^{*}$
between $h_{n}$ and $\mathsf{h}$, by (iii) of Assumption 3, it holds that

$\displaystyle|\mathbf{K}_{y}(i,j)-\mathbf{K}_{n}(i,j)|$ $\displaystyle\leq
L\bigg{|}\frac{\|\mathbf{y}_{i}-\mathbf{y}_{j}\|_{2}}{\mathsf{h}^{1/2}}-\frac{\|\mathbf{y}_{i}-\mathbf{y}_{j}\|_{2}}{h_{n}^{1/2}}\bigg{|}^{\tau_{0}}$
$\displaystyle\leq\bigg{|}\sqrt{\frac{\mathsf{h}}{h_{n}}}-1\bigg{|}\cdot\bigg{(}\frac{\|\mathbf{y}_{i}-\mathbf{y}_{j}\|_{2}^{2}}{\mathsf{h}}\bigg{)}^{\tau_{0}/2}$
$\displaystyle\prec{\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{\alpha_{i}-\beta})^{1/2}},}$

where the last step follows from (74) and Proposition 1. Similar to the
discussion of (B.1), we have

$\displaystyle\frac{1}{n}\|\mathbf{K}_{y}-\mathbf{K}_{n}\|$
$\displaystyle\leq\max_{i\neq j}|\mathbf{K}_{y}(i,j)-\mathbf{K}_{n}(i,j)|$
$\displaystyle\prec{\left(\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{\alpha_{i}-\beta})^{1/2}}\right).}$
(112)

Combining the above results, we have proven (39). ∎

###### Eigenvector convergence..

Based on the integral representation (77), for each $i$ satisfying (41), we
have

$\displaystyle\langle\mathbf{u}_{i},\mathbf{v}_{i}\rangle^{2}$
$\displaystyle=\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{i}}\mathbf{v}_{i}^{\top}(z{\bf
I}-n^{-1}\mathbf{K}_{x})^{-1}\mathbf{v}_{i}\mathrm{d}z$
$\displaystyle\quad+\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{i}}\mathbf{v}_{i}^{\top}\left[(z{\bf
I}-n^{-1}\mathbf{K}_{y})^{-1}-(z{\bf
I}-n^{-1}\mathbf{K}_{x})^{-1}\right]\mathbf{v}_{i}\mathrm{d}z$
$\displaystyle\quad+\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{i}}\mathbf{v}_{i}^{\top}\left[(z{\bf
I}-n^{-1}\mathbf{K}_{n})^{-1}-(z{\bf
I}-n^{-1}\mathbf{K}_{y})^{-1}\right]\mathbf{v}_{i}\mathrm{d}z$
$\displaystyle:=\mathsf{L}_{1}+\mathsf{L}_{2}+\mathsf{L}_{3},$ (113)

for some contour
$\Gamma_{i}:=\mathbb{B}(\gamma_{i},\frac{\mathsf{r}_{i}}{C}),$ for some large
constant $C>0,$ where $\mathbb{B}(\gamma_{i},\frac{\mathsf{r}_{i}}{C})$ is the
disk centered at $\gamma_{i}$ with radius $\mathsf{r}_{i}/C$ and we recall
$\mathsf{r}_{i}:=\min\\{\gamma_{i-1}-\gamma_{i},\gamma_{i}-\gamma_{i+1}\\}$.

On the one hand, by

$|\mu_{i}-\gamma_{i}|\prec\frac{1}{\sqrt{n}},$ (114)

and the assumption of (41) that
${\mathsf{r}_{i}\gg\big{[}\big{(}\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}\big{)}^{\xi}+\frac{1}{\sqrt{n}}\big{]}}$,
when $n$ is sufficiently large, $\lambda_{i}(\mathbf{K}_{x})$ is the only
simple pole of the resolvent $(z{\bf I}-n^{-1}{\mathbf{K}}_{x})^{-1}$ inside
$\Gamma_{i}$. Using the spectral decomposition of $(z{\bf
I}-n^{-1}\mathbf{K}_{x})^{-1}$ and residual theorem, we have that

$\mathsf{L}_{1}=1.$ (115)

On the other hand, by the definition of $\Gamma_{i},$ we have

$\displaystyle|\lambda_{i}(\mathbf{K}_{y})-\gamma_{i}|\leq|\lambda_{i}(\mathbf{K}_{y})-\mu_{i}|+|\gamma_{i}-\mu_{i}|$
$\displaystyle\prec{\bigg{(}\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}\bigg{)}^{\xi}+\frac{1}{\sqrt{n}},}$

so that under the same event, by (28), we have

$\inf_{z\in\Gamma_{i}}\min\\{|\lambda_{i}(\mathbf{K}_{y})-z|,|\lambda_{i-1}(\mathbf{K}_{y})-z|,|\lambda_{i+1}(\mathbf{K}_{y})-z|\\}\asymp\mathsf{r}_{i}.$

Moreover, by the definition of the resolvent, we have that

$\displaystyle\sup_{z\in\Gamma_{i}}\|(z{\bf I}-n^{-1}\mathbf{K}_{y})^{-1}\|$
$\displaystyle\leq\sup_{z\in\Gamma_{i}}\left(\frac{1}{|\lambda_{i}(\mathbf{K}_{y})-z|}+\frac{1}{|\lambda_{i-1}(\mathbf{K}_{y})-z|}+\frac{1}{|\lambda_{i+1}(\mathbf{K}_{y})-z|}\right)$
$\displaystyle\prec\frac{1}{\mathsf{r}_{i}}.$

Similarly, by (114), we also have $|\mu_{i}-\gamma_{i}|=o(\mathsf{r}_{i})$
with high probability, which implies

$\displaystyle\sup_{z\in\Gamma_{i}}\|(z{\bf
I}-n^{-1}\mathbf{K}_{x})^{-1}\|\prec\frac{1}{\mathsf{r}_{i}}.$

Now recall the resolvent identity

$(z{\bf I}-n^{-1}\mathbf{K}_{y})^{-1}-(z{\bf
I}-n^{-1}\mathbf{K}_{x})^{-1}=(z{\bf
I}-n^{-1}\mathbf{K}_{y})^{-1}\left[n^{-1}\mathbf{K}_{y}-n^{-1}\mathbf{K}_{x}\right](z{\bf
I}-n^{-1}\mathbf{K}_{x})^{-1}.$

Applying (110), we have

$\mathsf{L}_{2}\prec{\frac{1}{\mathsf{r}_{i}^{2}}\left[\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}\right]^{\xi}.}$

Similarly, by (112), we also have

$\displaystyle\sup_{z\in\Gamma_{i}}\|(z{\bf
I}-n^{-1}\mathbf{K}_{n})^{-1}\|\prec\frac{1}{\mathsf{r}_{i}},$

and therefore

$\mathsf{L}_{3}\prec{\frac{1}{\mathsf{r}_{i}^{2}}\left(\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{\alpha_{i}-\beta})^{1/2}}\right).}$

Insert the above bounds and (115) back into (B.4.1), we immediately obtain
that

$\displaystyle\left|\langle\mathbf{u}_{i},\mathbf{v}_{i}\rangle^{2}-1\right|\prec{\frac{1}{\mathsf{r}_{i}^{2}}\left[\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-1}}+\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}\right]^{\xi},}$
(116)

as $\xi\leq 1$. This completes our proof of (42).

∎

####  B.4.2 Convergence to Population Integral Operator

###### Eigenvalue convergence..

The proof of (40) is the same as Theorem 3. ∎

###### Eigenfunction convergence..

Again, in light of Proposition 1 and by the same argument as in the beginning
of Section B.2, it suffices to prove the result under the setting where the
bandwidth is chosen as (52). Recall (84). By an argument similar to Section
B.4.1, we have

$\displaystyle
f\left(\frac{\|\mathbf{y}-\mathbf{y}_{i}\|_{2}}{\mathsf{h}^{1/2}}\right)\leq
f\left(\frac{\|\bm{y}-\bm{x}_{i}\|_{2}}{\mathsf{h}^{1/2}}\right)+{\mathrm{O}_{\prec}\bigg{(}\left[\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}\right]^{\xi}\bigg{)},}$

uniformly in $\mathbf{y}$ on the support of $\widetilde{\mathsf{P}}$. Similar
to (86), we decompose

$\left(\widehat{\phi}_{i}^{(n)}(\mathbf{y})-\phi_{i}^{(n)}(\bm{y})\right)^{2}=\left(\mathrm{E}_{1}+\mathrm{E}_{2}+\mathrm{E}_{3}\right)^{2},$
(117)

where

$\displaystyle\mathrm{E}_{1}:=\left(\frac{1}{\lambda_{i}}-\frac{1}{\mu_{i}}\right)\frac{1}{\sqrt{n}}\sum_{j=1}^{n}f\left(\frac{\|\mathbf{y}-\mathbf{y}_{j}\|_{2}}{\mathsf{h}^{1/2}}\right)v_{ij},$
$\displaystyle\mathrm{E}_{2}:=\frac{1}{\sqrt{n}\mu_{i}}\sum_{j=1}^{n}\left[f\left(\frac{\|\mathbf{y}-\mathbf{y}_{j}\|_{2}}{\mathsf{h}^{1/2}}\right)-f\left(\frac{\|\bm{y}-\bm{x}_{j}\|_{2}}{\mathsf{h}^{1/2}}\right)\right]v_{ij},$
$\displaystyle\mathrm{E}_{3}:=\frac{1}{\sqrt{n}\mu_{i}}\sum_{j=1}^{n}f\left(\frac{\|\bm{y}-\bm{x}_{j}\|_{2}}{\mathsf{h}^{1/2}}\right)(v_{ij}-u_{ij}).$

For $\mathrm{E}_{1}$ and $\mathrm{E}_{2}$, following the same argument as in
proof of Theorem 3, along with (ii) of Assumption 3, we have

$\mathrm{E}_{1}\prec{\frac{1}{\mu_{i}}\bigg{(}\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}\bigg{)}^{\xi},}$
$\mathrm{E}_{2}\prec{\frac{1}{\mu_{i}}\bigg{(}\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}\bigg{)}^{\xi}.}$

For $\mathrm{E}_{3}$, we denote $\mathbf{a}=(a_{1},\cdots,a_{n})^{\top}$ and
$\mathbf{b}=(b_{1},\cdots,b_{n})^{\top}$, where
$a_{j}=\frac{1}{\sqrt{n}\mu_{i}}f\left(\frac{\|\bm{y}-\bm{x}_{j}\|_{2}}{\mathsf{h}^{1/2}}\right)$
and $b_{j}=v_{ij}-u_{ij}.$ Similar to (87), we have that

$\mathrm{E}_{3}^{2}\leq\|\mathbf{a}\|_{2}^{2}\|\mathbf{b}\|_{2}^{2}\prec{\frac{1}{\mu_{i}^{2}}\cdot\frac{1}{\mathsf{r}_{i}^{2}}\left[\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}\right]^{\xi},}$

where in the second step we used (116). Combining all the above bounds, using
(117), when $i\leq\mathsf{K},$ we arrive at

$\displaystyle\mu_{i}\left(\widehat{\phi}_{i}^{(n)}(\mathbf{y})-\phi_{i}^{(n)}(\bm{y})\right)^{2}$
$\displaystyle\prec{\bigg{[}\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}\bigg{]}^{2\xi}}$
$\displaystyle\quad+{\frac{1}{\mathsf{r}_{i}^{2}}\left[\frac{1}{\sum_{i=1}^{r}n^{\alpha_{i}-\beta-\eta}}+\frac{1}{(\sum_{i=1}^{r}n^{{\alpha_{i}-\beta}})^{1/2}}\right]^{\xi},}$
(118)

for any $\mathbf{y}$ on the support of $\widetilde{\mathsf{P}}$.Finally, by
Part 1 of Theorem 4, we have

$\displaystyle|\sqrt{\lambda_{i}}\phi_{i}^{(n)}(\mathbf{y})-\sqrt{\mu_{i}}\phi_{i}^{(n)}(\mathbf{y})|$
$\displaystyle\leq(\sqrt{\lambda_{i}}-\sqrt{\mu_{i}})|\phi_{i}^{(n)}(\mathbf{y})|\prec{\psi_{n}^{\xi},}$

and similarly

$|\sqrt{\lambda_{i}}\widehat{\phi}_{i}^{(n)}(\mathbf{y})-\sqrt{\mu_{i}}\widehat{\phi}_{i}^{(n)}(\mathbf{y})|\leq(\sqrt{\lambda_{i}}-\sqrt{\mu_{i}})|\widehat{\phi}_{i}^{(n)}(\mathbf{y})|\prec{\psi_{n}^{\xi}.}$
(119)

Then by a discussion similar to (91), we have

$\displaystyle|\sqrt{\lambda_{i}}\widehat{\phi}_{i}^{(n)}(\mathbf{y})-\sqrt{\gamma_{i}}{\widetilde{\phi}_{i}}(\mathbf{y})|\prec|\sqrt{\mu_{i}}\widehat{\phi}_{i}^{(n)}(\mathbf{y})-\sqrt{\mu_{i}}\phi_{i}^{(n)}(\mathbf{y})|+{\psi_{n}^{\xi}}+\frac{1}{\sqrt{n}\mathsf{r}_{i}}.$

Plugging the above bounds into (B.4.2), we obtain the final result. ∎

##  Appendix C Some Further Discussions

###  C.1 Extensions and Related Problems

We propose a kernel-spectral embedding algorithm for learning low-dimensional
nonlinear structure from noisy and high-dimensional datasets under a manifold
setup. A key component is to develop an adaptive bandwidth selection procedure
that is free from prior knowledge of the manifold. Our proposed method is
theoretically justified and our embeddings correspond to an integral operator
from the RKHS associated to the underlying manifold. While our results pave
the way towards statistical foundations for nonlinear dimension reduction,
manifold learning, among others, there are various problems to be investigated
further. We now highlight a few of them.

First, although this paper focuses on a general class of kernel functions for
the construction of the low-dimensional embeddings, there are still many other
interesting kernel functions, such as $f(x)=\frac{\sin x}{x}$ and
$f(x)=\frac{1}{e^{x}+e^{-x}}$, that are excluded from our discussion.
Moreover, the current paper only concerns the distance-type kernel matrices as
in (3), leaving the inner-product kernel matrices with
$K(i,j)=f(\mathbf{y}_{i}^{\top}\mathbf{y}_{j}/h_{n})$ unaddressed.
Nevertheless, we believe this is only an initial step towards interpretable
manifold learning under high-dimensional noisy data, and our results stand as
a vantage point for exploring and understanding other kernel-spectral
embedding algorithms.

Second, in Part 2 of Theorem 3, the convergence of eigenfunctions is obtained
in a pointwise manner. In some applications, a more precise local
characterization of the embeddings is needed. For example, one may need to
study the eigenfunction convergence under the $L_{\infty}$ norm, defined by
$\|f\|_{\infty}=\sup_{\mathbf{x}}|f(\mathbf{x})|$. However, to our best
knowledge, even in the noiseless setting, the $L_{\infty}$ convergence has
only been studied for the normalized integral operators [93, 92] and the
linear differential Laplace–Beltrami operator [21, 37, 98]. Less is known for
the unnormalized integral operator considered in our paper. Moreover, to
understand the $L_{\infty}$ convergence under noisy datasets, in contrast to
(29), we also need to study the $\ell_{\infty}$ convergence of the
eigenvectors. The analysis requires a more sophisticated argument from random
matrix theory as in [6, 41]. We will pursue these directions in future works.

Additionally, the current work mainly concerns the super-critical non-null
regime $n^{\beta+\eta}=\mathrm{o}(\sum_{i=1}^{r}n^{\alpha_{i}})$. When
$n^{\beta+\eta}\asymp\sum_{i=1}^{r}n^{\alpha_{i}}$, according to (2) of
Theorem 3.1 of [32], since the signal part has a comparable strength as the
noise part, neither of them dominates the other in terms of the behavior of
$n^{-1}\mathbf{K}_{n}$. Similar results have also been established in [39]. In
this critical regime, to better utilize the kernel matrix
$n^{-1}\mathbf{K}_{n},$ one needs to develop a high-dimensional nonlinear
signal and noise separation procedure.

Finally, in the current paper, we use an integral operator to learn the
geometric structure of the underlying manifold. To enhance our understanding
of the nonlinear structures, it is of interest to explore the relationship
between the eigenfunctions of our integral operator and those of the Laplace-
Beltrami operator. The main challenge is that, in the high-dimensional noisy
setting, the normalized graph Laplacian would no longer converge to the
Laplace-Beltrami operator [32], so that novel methodologies are needed to
facilitate the comparison. Moreover, in the current paper, we focus on an
RKHS-based method. It will be interesting to study the graph-based methods for
high-dimensional noisy datasets under the current setting.

###  C.2 Comparisons with Some Existing Works

In Section 2.2 of the our manuscript, we have briefly summarized the
differences between the current work and those of [48, 63, 1]. In what
follows, we provide some more details.

First, in [48], the author proposed a graph-cut based spectral clustering
algorithm (”SPUDS”), in which a major step was to obtain the spectral
embedding based on some graph Laplacian matrix. Specifically, the method
starts by constructing a Gaussian kernel matrix
$\mathbf{K}_{n}=(K(i,j))_{1\leq i,j\leq n}$ using the following bandwidth

$h_{n}=sn^{-1/(r+4)},$ (120)

where $s$ is some user-specified value which is bounded both above and away
from zero, and then obtains the eigendecomposition of the graph Laplacian
matrix ${\bf L}_{n}={\bf
I}-\mathbf{D}_{n}^{-1/2}\widehat{\mathbf{K}}_{n}\mathbf{D}_{n}^{-1/2}$, where
$\mathbf{D}=\text{diag}(\sum_{j\neq 1}K(1,j),...,\sum_{j\neq n}K(n,j))$, and
$\widehat{\mathbf{K}}_{n}=({K}(i,j)\cdot 1_{\\{i\neq j\\}})_{1\leq i,j\leq n}$
is the zero-diagonal kernel matrix. For some given initial value for the
number of clusters, K-means is then used along with the spectral embedding
(i.e., the eigenvectors of ${\bf L}_{n}$) to determine the cluster labels and
update the number of clusters iteratively using the normalized cut. In
general, both SPUDS and our proposed method rely on some kernel spectral
embeddings based on some kernel matrices, and require careful selection of the
bandwidths. However, there are many significant differences between the two
methods. First, [48] focuses on the noiseless data sets when clusters are
closely related to the modes of the density functions whereas we consider
high-dimensional noisy data sets. Second, the two methods focus on the
spectrum of two very different matrices: our method relies on the
eigendecomposition of the kernel matrix itself, whereas SPUDS relies on the
graph Laplacian matrix of the zero-diagonal kernel matrix. Third, the
bandwidth selection schemes are also very different. On the one hand, in [48],
$h_{n}$ in (120) is chosen to ensure the consistency of the normalized cut
(see Theorem 2 therein) rather than the kernel matrices or graph Laplacian
matrices. In fact, it is not clear from [48] that whether such a choice of
bandwidth ensures consistency and meaningful convergence of the kernel
matrices or graph Laplacian matrices. In contrast, our selected bandwidth is
capable of capturing the underlying signals and therefore guarantees the
convergence of the kernel matrices to some population integral operator. On
the other hand, our bandwidth selection procedure is data-driven and does not
require any prior knowledge. In contrast, (120) relies on the values of both
the (ambient) dimension $r$ which is in general unknown and some constant $s.$
To choose $s,$ for practical purpose, in Section 4 of [48], the author suggest

$h_{n}=1.2\sqrt{\bar{\lambda}_{r}}\bigg{(}\frac{4}{(2+r)n}\bigg{)}^{\frac{1}{r+4}},$
(121)

where $\bar{\lambda}_{r}$ is the average of the largest $r$ eigenvalues of the
covariance matrix of the data. However, in our setting even when $r$ is fixed,
$\bar{\lambda}_{r}$ usually diverges which violates the boundedness
requirement of $s$ in (120). Numerically, in Section 4.1.2 and Figure 1 of the
main manuscript, compared to the method in [48], we show that our proposed
method has better performance in terms of spectral clustering using two high-
dimensional noisy data sets.

Second, in [63], the authors studied the spectral embedding and clustering of
high-dimensional data under the Gaussian mixture model. The spectral embedding
$\hat{\bf Y}\in\mathbb{R}^{n\times k}$ was obtained based on the spectral
decomposition

${\bf
Y}=\sum_{i=1}^{\min\\{n,p\\}}\hat{\lambda}_{i}\hat{u}_{i}\hat{v}_{i}^{\top}\in\mathbb{R}^{n\times
p},\qquad\hat{\lambda}_{1}\geq...\geq\hat{\lambda}_{\min\\{p,n\\}},$

of the observed data matrix ${\bf Y}$, rather than some kernel matrices.
Consequently, they do not need to select a bandwidth. For clustering, they
used $\hat{\bf Y}=\hat{\bf U}\hat{\bf\Lambda}$, where
$\hat{\bf\Lambda}=\text{diag}(\hat{\lambda}_{1},...\hat{\lambda}_{k})$ and
$\hat{\bf U}=(\hat{u}_{1},...,\hat{u}_{k})$, with $k$ being the target
embedding dimension. K-means was then used on $\hat{\bf Y}$ to determine the
clusters. Compared with Step 3(ii) of our Algorithm 1, the embedding $\hat{\bf
Y}$ shared a similar form as our proposal. However, the linear embedding
$\hat{\bf Y}$ may not be as effective as our method in dealing with data sets
with nonlinear structures. For illustrations, in Section 4.1.2 of the main
manuscript, we compared the numerical performance of these methods under both
a Gaussian mixture model and a nonlinear nested sphere model. Figure 1 of the
main manuscript shows that, while the method of [63] performed equally well as
ours under the Gaussian mixture model, its performance was significantly worse
than our method under the nonlinear nested sphere model.

Finally, in [1], the authors developed a perturbation theory for a hollowed
version of PCA in general Hilbert spaces. Specifically, the model

$\mathbf{y}^{\prime}_{i}={\bf x}^{\prime}_{i}+{\bf
z}^{\prime}_{i}\in\mathcal{H},\qquad 1\leq i\leq n,$ (122)

was considered for some reproducing kernel Hilbert space (RKHS) $\mathcal{H}$
associated with the map $\phi:\mathcal{X}\to\mathcal{H}$ via the kernel
$K(\cdot,\cdot):\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ with
$K(\mathbf{x},\mathbf{y})=\langle\phi(\mathbf{x}),\phi(\mathbf{y})\rangle$. In
(122), $\\{\mathbf{x}^{\prime}_{i}\\}_{1\leq i\leq n}$ are the noiseless
signals, and $\\{{\bf z}^{\prime}_{i}\\}_{1\leq i\leq n}$ are sub-Gaussian
noise. The general idea, like our paper, was to learn the geometric
information of some underlying structure contained in the eigenfunctions of
some integral operator. The focus of [1] was on the perturbation analysis of
the eigenspaces associated to the noiseless infinite-dimensional self-adjoint
operator, or Gram matrix ${\bf H}_{n}^{*}={\bf X^{\prime}X^{\prime}}^{\top}$,
where

${\bf X}^{\prime}=\begin{bmatrix}{\bf x}^{\prime}_{1}\\\ {\bf
x}^{\prime}_{2}\\\ \vdots\\\ {\bf x}^{\prime}_{n}\end{bmatrix},$

such that ${\bf X}^{\prime}({\bf g})=(\langle{\bf x}^{\prime}_{1},{\bf
g}\rangle,...\langle{\bf x}^{\prime}_{n},{\bf g}\rangle).$ Furthermore, for
sub-Gaussian and Gaussian mixture models, the spectral embedding $\hat{\bf
Y}^{\prime}={\bf U}^{\prime}{\bf\Lambda^{\prime}}^{1/2}$ was proposed, where
${\bf\Lambda^{\prime}}=\text{diag}(\lambda^{\prime}_{1},...,\lambda^{\prime}_{k})$
and ${\bf U}^{\prime}=({\bf u}_{1},...,{\bf u}_{k})$ contain the $k$ leading
eigenvalues and eigenvectors of the zero-diagonal Gram matrix ${\bf
H}_{n}=\mathsf{H}({\bf Y}^{\prime}{\bf Y}^{\prime\top})$, with the operator
$\mathsf{H}(\cdot)$ zeroing out all diagonal entries of a square matrix.

In general, we agree that [1] is similar to our paper in the sense that both
works consider spectral embedding of noisy data and study the convergence of
some eigenvectors based on the observed samples to their counterparts
associated with the noiseless samples. In what follows, we compare and
contrast the two papers in the following aspects. First, the noiseless data of
[1] is assumed to be generated from some RKHS as in (122). In contrast, our
noiseless data is sampled from a nonlinear manifold model as in (1) and
Assumption 1. Second, [1] focused on the spectrum of the Gram matrices ${\bf
H}_{n}^{*}$ and ${\bf H}_{n}$ defined above, whereas our paper focused on the
distance kernel matrices $\mathbf{K}_{n}^{*}$ and $\mathbf{K}_{n}$ as in (3)
and (19). The main difference between the two types of matrices can be seen as
follows. On the one hand, construction of Gram matrices relies on specifying
the map $\phi$, whereas our kernel matrices rely on the choice of kernel
functions and a careful selection of the bandwidth. On the other hand, a key
idea of [1] was to zero out the diagonal entries of the Gram matrix ${\bf
Y}^{\prime}{\bf Y}^{\prime\top}$ as in [57, 40, 67], which led to better
spectral convergence rates. However, for the kernel matrices considered in our
manuscript, since the diagonal part is isotropic (i.e., $f(0)\mathbf{I}$ for
some given kernel function $f(\cdot)$), zeroing out the diagonals does not
change the eigenvectors, nor the convergence results. Third, in terms of
theoretical assumptions, on the one hand, the assumptions of [1] on the
signal-to-noise ratio are comparable to ours when
$\mathcal{H}=\mathbb{R}^{p}$. For example, under our model (Assumption 2) with
finite $r$, $\eta=1$ and $\beta=0$, the assumptions of [1] (cf. Assumptions
2.4 and 2.6 therein) imply that $\alpha>1$, which is the same as our paper. On
the other hand, for general nonlinear space $\mathcal{H}$, since the
assumptions of [1] are made on the Gram matrix ${\bf H}_{n}^{*}$, they can be
relatively less interpretable compared to ours, which are made directly on the
underlying manifold (cf. Assumptions 1 and 2) and the kernel functions (cf.
Assumption 3). Finally, in terms of theoretical results, both papers showed
spectral convergence of some matrices to their noiseless counterparts.
However, the eigenvector perturbation bounds were obtained under different
discrepancy measures (cf. Theorem 2.1 of [1] and our Theorems 2 and 4). In
particular, our results highlighted the impact of signal-to-noise ratio, the
underlying manifold structures and the choice of kernel functions on the final
rate of convergence, and, by connecting to the population integral operator,
yielded the proper interpretation of the eigenvectors with respect to the
underlying manifold structures. Such results have not been established in [1].

In addition, in terms of the proof technique, [1] relies on the leave-one-out
analysis, and uses the Davis-Kahan theorem for proving the eigenvector
perturbation bounds, whereas our proof of the eigenvector convergence relies
on the integral representation of the eigenvectors along with the resovlent
expansion. Since bandwidth plays an important role in our algorithm, we also
need to establish a novel concentration inequality (cf. Proposition 1) for our
bandwidth selection scheme. Numerically, since for given kernel matrices, the
proposal of [1] is equivalent to Step 3 of our algorithm, and the main
advantage of our method lies in the construction of kernel matrices and
bandwidth selection. For illustrations, in Section 4.1.2, we compared our
method with [1] using two mixture models. The result is summarized in Figure 1
of the manuscript which demonstrates the advantage of our proposal especially
when data sets are highly nonlinear.

###  C.3 Some additional remarks

In this subsection, we provide a few more remarks.

First, we make a comparison with Davis-Kahan theorem in terms of the results
in (29). In the literature, there exist different variants of Davis-Kahan
theorem. For definiteness and easy comparison, we focus on the one from the
popular monograph [90] (see Theorem 4.5.5 therein); that is

######  Theorem 16 (Davis Kahan Theorem).

Let $\mathbf{S}$ and $\mathbf{T}$ be two symmetric matrices with the same
dimension. For each fixed $i,$ denote

$\delta:=\min_{j:j\neq i}|\lambda_{i}(\mathbf{S})-\lambda_{j}(\mathbf{S})|.$
(123)

Assuming that $\delta\geq\tau>0$ for some constant $\tau>0,$ we have that

$|\langle\mathbf{u}_{i},\mathbf{v}_{i}\rangle^{2}-1|\leq\frac{4\|\mathbf{S}-\mathbf{T}\|^{2}}{\delta^{2}},$
(124)

where $\mathbf{u}_{i},\mathbf{v}_{i}$ are respectively the eigenvectors of
$\mathbf{S}$ and $\mathbf{T}.$

We see that in order to correctly apply Davis-Kahan and obtain (124), we need
that $\delta$ is bounded from below by some constant. In contrast, in our
equation (29), our bound reads $\psi_{n}/\mathsf{r}_{i}^{2}.$

We now make a comparison for both the denominators and numerators of the
bounds. For the $\mathsf{r}_{i}^{2}$ part, on the one hand, our
$\mathsf{r}_{i}$ only depends on three population eigenvalues, i.e.,
$\mathsf{r}_{i}=\min\\{\gamma_{i}-\gamma_{i+1},\gamma_{i-1}-\gamma_{i}\\}$. In
fact, $\delta$ used in Davis-Kahan theorem, i.e., (123), requires much more
information than our $\mathsf{r}_{i}$ (even though this condition may be
relaxed using the arguments in [100]). On the other hand and more importantly,
our $\mathsf{r}_{i}$ is allowed to decay to zero once equation (28) is
satisfied. For the $\psi_{n}$ part, it is weaker than the results
$\psi_{n}^{2}$ from Davis-Kahan which relies on stronger assumption on
$\mathsf{r}_{i}.$ Nevertheless, if $\delta$ is bounded from below by a
constant, we can improve our results using Davis-Kahan. To summarize, Davis-
Kahan may result in faster convergence rates but at the cost of much stronger
assumptions on the eigenvalues separation whereas our results are obtained
under general and weaker separation conditions. In the current paper, for the
purpose of generality and practical usability, we keep the current results in
(29).

Second, it is true that the left-hand side of (28) is large than $\psi_{n}.$
Therefore, (28) means that $\psi_{n}$ is much smaller than
$\mathsf{r}_{i}^{2}$ so that (123) implies consistent estimation. Sorry for
the confusion. Finally, the reason that the left-hand side of (28) is larger
than the ideal term $\psi_{n}$ is technical. Our proof relies on an integral
representation via the resolvent, such an assumption is needed to make sure
that we can find a contour which only contains $\lambda_{i}$ so that the
residual theorem can be applied; see the discussions between equations (B.1)
and (79).

##  Appendix D Some Results on Gaussian Integral Operators

In this section, we collect some results on integral operators. Consider a
special case when $\\{\bm{x}_{i}\\}$ are Gaussian. In such a setting, the
eigenvalues and eigenfunctions of the integral operator can be calculated
explicitly. We first state the results when $r=1.$ Let $\mathsf{p}(x)$ be the
density function of $\mathcal{N}(0,\sigma^{2})$ random variable. For the ease
of statement, we write the Gaussian kernel function that
$k(x,y)=e^{-\frac{\|x-y\|_{2}^{2}}{2h}},$ for some fixed bandwidth $h.$
Consequently, we can construct the following integral operator

$\mathcal{K}^{\sigma}_{1}f(x)=\int k(x,y)f(y)\mathsf{p}(y)\mathrm{d}y.$

Denote $\\{\gamma_{i}\\}$ and $\\{\phi_{i}(x)\\}$ as the sequence of the
eigenvalues and eigenfunctions of $\mathcal{K}^{\sigma}_{1}$ with

$\mathcal{K}_{1}^{\sigma}\phi_{i}(x)=\gamma_{i}\phi_{i}(x),\qquad\int\phi_{i}^{2}(x)\mathsf{p}(x)\mathrm{d}x=1.$

Note that the eigenfunctions are in fact defined in $\mathbb{R}$ via the
relation that

$\phi_{i}(x)=\frac{1}{\gamma_{i}}\int
k(x,y)\phi_{i}(y)\mathsf{p}(y)\mathrm{d}y,\ \text{if}\ \lambda_{i}>0.$

The following lemma provides the exact form for the eigenvalues and
eigenfunctions.

######  Lemma 17.

Denote $\beta=2\sigma^{2}/h.$ The eigenvalues and eigenfunctions of
$\mathcal{K}_{1}^{\sigma}$ are

$\gamma_{i}=\frac{\sqrt{2}}{\sqrt{1+\beta+\sqrt{1+2\beta}}}\left(\frac{\beta}{1+\beta+\sqrt{1+2\beta}}\right)^{i},$
(125)
$\phi_{i}(x)=\frac{(1+2\beta)^{1/8}}{\sqrt{2^{i}i!}}\exp\left(-\frac{x^{2}}{2\sigma^{2}}\frac{\sqrt{1+2\beta}-1}{2}\right)H_{i}\left(\left(\frac{1}{4}+\frac{\beta}{2}\right)^{1/4}\frac{x}{\sigma}\right),$
(126)

where $H_{i}(\cdot)$ is the $i$th order Hermite polynomial.

###### Proof.

See Proposition 1 of [79]. ∎

The above results can also be generalized to the $r$-dimensional setting when
$\mathsf{p}(\bm{x})$ is the density function of
$\mathcal{N}(\mathbf{0},\Sigma)$ where $\Sigma$ is a diagonal matrix so that

$\Sigma=\operatorname{diag}\\{\theta_{1},\cdots,\theta_{r}\\},\qquad\text{where
$\theta_{i}=\sigma_{i}^{2}$}$

Due to invariance of the kernel, the associated integral operator, denoted as
$\mathcal{K}_{r},$ has a direct sum decomposition
$\mathcal{K}_{r}=\bigoplus_{i=1}^{r}\mathcal{K}_{1}^{\sigma_{i}}.$
Analogously, the eigenvalues and eigenfunctions of $\mathcal{K}_{r}$ can also
be calculated explicitly as follows. Recall that for given two operators
$\mathcal{F}_{1}$ and $\mathcal{F}_{2},$ the spectrum of their direct sum
$\mathcal{F}_{1}\bigoplus\mathcal{F}_{2}$ consists of pairwise products
$\gamma_{i}(\mathcal{F}_{1})\gamma_{j}(\mathcal{F}_{2}).$ Moreover, the
eigenfunction of the product can be written into a multindex form so that
corresponding to $\gamma_{i}(\mathcal{F}_{1})\gamma_{j}(\mathcal{F}_{2})$, we
have

$\phi_{[i,j]}(x_{1},x_{2})=\phi_{i}(x_{1})\psi_{j}(x_{2}),$

where $\\{\phi_{i}(x)\\}$ and $\\{\psi_{j}(x)\\}$ are the eigenfunctions of
the operators $\mathcal{F}_{1}$ and $\mathcal{F}_{2},$ respectively.
Consequently, we can obtain the following results.

######  Lemma 18.

Denote $\beta_{i}=2\sigma^{2}_{i}/h.$ Then the eigenvalues and eigenfunctions
of $\mathcal{K}_{r}$ can be represented using the multindex over all
components, i.e.,

$\gamma_{[i_{1},i_{2},\cdots,i_{r}]}=\prod_{j=1}^{r}\gamma_{i_{j}}(\mathcal{K}_{1}^{\sigma_{j}}),$

and for $\bm{x}=(x_{1},x_{2},\cdots,x_{r})$, we have

$\phi_{[i_{1},i_{2},\cdots,i_{r}]}(\bm{x})=\prod_{j=1}^{r}\phi_{i_{j}}(\mathcal{K}_{1}^{\sigma_{j}})(x_{j}).$

###### Proof.

See Section 3 of [78]. ∎

##  Appendix E Tuning Parameter Selection and Additional Numerical Results

###  E.1 A Resampling Method for Selecting Percentile $\omega$ in (5)

As shown by our theory in Section 3, under our assumptions, $\omega$ can be
chosen as any constant between 0 and 1 to have the final embeddings that
achieve the same asymptotic behavior. Especially, our proof of Proposition 1
depends on the aspect ratio (2) fundamentally. In practice, to optimize the
empirical performance and improve automation of the method, we recommend using
a resampling approach as follows.

A resampling-based algorithm for determining the percentile parameter was
proposed in Algorithm 1 of [32], which has also been used in [28, 33]. The
method provides a choice of $\mathsf{s}$ using resampling method to
distinguish the larger outlier eigenvalues and bulk eigenvalues of the kernel
matrix $\mathbf{K}_{n}$, where the outlier eigenvalues stand for the signal
parts and bulk eigenvalues are for the noise part. A key observation
underlying this approach is that the bulk eigenvalues are close to each other
so that the ratios of two consecutive eigenvalues are close to one (see Remark
2.9 of [32] for more details). Given the choice of $\mathsf{s},$ the algorithm
is summarized below for self-completeness.

  1. 1.

For a pre-selected sequence of percentiles$\\{\omega\\}_{i=1}^{T},$ calculated
the associated bandwidths according to (5), denoted as
$\\{h_{i}\\}_{i=1}^{T}.$

  2. 2.

For each $1\leq i\leq T,$ calculate the eigenvalues of $\mathbf{K}_{n,i}$,
constructed using the bandwidth $h_{i}.$ Denote the eigenvalues of
$\mathbf{K}_{n,i}$ in the decreasing order as
$\\{\lambda_{k}^{(i)}\\}_{k=1}^{n}.$

  3. 3.

For some $\mathsf{s}>0,$ denote

$\mathsf{k}(\omega_{i}):=\max_{1\leq k\leq
n-1}\left\\{k\bigg{|}\frac{\lambda_{k}^{(i)}}{\lambda_{k+1}^{(i)}}\geq
1+\mathsf{s}\right\\}.$

Choose the percentile $\omega$ such that

$\omega=\max_{i}[\operatorname*{argmax}_{\omega_{i}}\mathsf{k}(\omega_{i})].$
(127)

In particular, the last step can be replaced by other criteria according to
users’ purpose. Here we choose the largest one for the purpose of robustness.
To demonstrate the empirical performance of the proposed percentile selection
algorithm, we consider the four examples in Section 4.1, and evaluate the
empirical percentile chosen by (127), across a variety of sample sizes. In
each setting, we used $\textsf{s}=1$. Figure E.1 shows the empirical
percentile $\omega$ selected by the resampling method in each example. Along
with Figures E.3 and E.5 below, these results indicate the usefulness and
robustness of the percentile selection algorithm.

Figure E.1: Plots of $\omega$ selected by the resampling method (127). Top row
from left to right: ”smiley face” and ”mammoth”; Bottom row from left to
right: ”Cassini oval” and ”torus.”

###  E.2 Additional Results from Simulation Studies

In Figure E.2, we illustrate the manifold structures considered in the
simulation studies. In Figure E.3, we show the convergence of kernel matrices.
In Figure E.4, we show the usefulness of our proposed population bandwidth. In
Figure E.5, we compare the rates of convergence under different kernel
functions. In particular, the top row of Figure E.5 shows that the convergence
rate under the Gaussian kernel is faster than that under the Laplacian kernel,
whereas the bottom row of Figure E.5 shows that convergence rate under the
Gaussian kernel is similar to that under the quadratic polynomial kernel.
These empirical results agree with our theory in Section 3.4.

Figure E.2: Top: ”smiley face,” ”mammoth,” and ”torus” (from left to right);
Bottom: ”Cassini oval” and its two-dimensional projections.

Figure E.3: Spectral convergence of kernel matrices under various kernel
functions and $\omega\in\\{0.05,0.25,0.5,0.75,0.95\\}$. Left: Gaussian kernel;
Middle: Laplacian kernel; Right: polynomial kernel. From top to bottom:
”smiley face,” ”mammoth,” ”Cassini oval,” and ”torus.”

Figure E.4: Comparison of eigenvalues of kernel matrices under various
bandwidths. Left: proposed bandwidth $\mathsf{h}=\mathsf{h}_{n}$ with
$\omega=0.5$; Middle: large bandwidth $\mathsf{h}=n^{5}$; Right: small
bandwidth $\mathsf{h}=p$. From top to bottom: ”smiley face,” ”mammoth,”
”Cassini oval,” and ”torus.”

Figure E.5: Top: ratios between the convergence rates under the Laplacian
kernel and under the Gaussian kernel. Bottom: ratios between the convergence
rates under the polynomial kernel and under the Gaussian kernel. From left to
right: ”smiley face,” ”mammoth,” ”Cassini oval” and ”torus.”

###  E.3 Two Additional Real Data Applications

In this subsection, we provide two additional real data analysis examples.

####  E.3.1 Cell Cycle Reconstruction

Our first example concerns the reconstruction of cell cycles from single-cell
RNA-Seq data. The cell cycle, or cell-division cycle, is the series of events
that take place in a cell that cause it to divide into two daughter
cells333https://en.wikipedia.org/wiki/Cell_cycle. Determining the cell cycle
stages of individual cells analyzed during development is important for
understanding its wide-ranging effects on cellular physiology and gene
expression profiles. Our dataset contains 288 mouse embryonic stem cells,
whose cell cycle stages were determined using flow cytometry sorting. As a
result, one-third (96) of the cells are in the G1 stage, one-third in the S
stage, and the rest in the G2M stage. The raw count data were preprocessed and
normalized using the same method as in Section 4.2, leading to a dataset
consisting of standardized expression levels of
$p\in\\{2500,3000,3500,4000\\}$ most variable genes for the 288 cells. We
apply our proposed method with a variety of $\omega\in\\{0.25,0.5,0.75\\}$.
Observing that the leading eigenvector is approximately a constant vector and
therefore non-informative, we consider a two-dimensional embedding with
$\Omega=\\{2,3\\}$, expected to capture the underlying circle manifold. As a
comparison, we also obtain two-dimensional embeddings based on PCA, MDS,
Laplacian eigenmap, LLE and DM, using functions implemented in the R package
dimRed under their default settings. To recover the cell cycle stages, or the
underlying position of the cells on the circle manifold, we project each
embedding to the two-dimensional unit circle, and then identify the cell
stages with their respective angles on the unit circle. We compare the
reconstruction performance of different methods based on their Kendall’s tau
distance to the true stages, up to a possible circular shift of the
reconstructed cycles.

The right panel of Figure E.6 shows that our proposed method has outstanding
performance compared to the other methods. In addition, Figure E.7 shows the
unit-circle projections of the two-dimensional embeddings obtained from each
method under $p=3000$. It can be seen that the projected cells reconstructed
by the proposed method were more separated and well-ordered according to their
true cycle stages.

Figure E.6: Left: comparison of five methods for cell ordering in Section 4.2.
Right: comparison of six methods for cell cycle reconstruction in Section
E.3.1. In both plots, the proposed method Algorithm 1 uses $\omega=0.5$.

Figure E.7: Comparison of reconstructed cell cycles based on six embedding
methods when $p=3000$. Red crosses: cells in G1 stage. Green dots: cells in S
stage. Black triangles: cells in G2M stage.

####  E.3.2 Clustering of Hand-Written Digits

Our last example concerns the MNIST444http://yann.lecun.com/exdb/mnist/
dataset containing images of hand-written digits. Specifically, we focus on
$n=3946$ images of hand-written digits ”2,” ”4,” ”6” and ”8,” among which
there are about 1000 images for each digit. As each image contains $28\times
28$ pixels, they can be treated as $784$-dimensional vectors. We apply six
different low-dimensional embedding methods, namely DM, kPCA as defined in
Section 2.2, Laplacian eigenmap, MDS, PCA and the proposed method, for a
variety of embedding dimensions n.dim $\in\\{5,10,15,20,25,30,35\\}$. In
particular, to ensure fairness in comparison, we have not included two-stage
algorithms such as tSNE and UMAP, which essentially take one of the above
kernel spectral embeddings as an initialization, and then refine the low-
dimensional embedding to make the cluster pattern more salient based on some
local-metric adjustment [66, 5, 62, 20].

Figure E.8: Comparison of six embedding methods for $n=3946$ samples from the
MNIST data. The Silhouette index indicates how much of the underlying cluster
pattern is preserved in the low-dimensional embeddings.

We evaluate the embedding quality of a method by calculating the average
Silhouette index [71] of the final low-dimensional embedding with respect to
the underlying true cluster membership. In general, a higher average
Silhouette index of a method indicates the underlying clusters are more
separate in the final embedding. Again, for the proposed method we use
$\omega=0.5$. Similar results are obtained for $\omega\in\\{0.25,0.75\\}$ as
in Section E.4. In Figure E.8, we found that our proposed method along with DM
has overall the best performance among the six methods. In particular,
although the standard kPCA differs from our proposed method only by an
additional mean shift (8), in all cases its embedding quality is clearly worse
than the latter, demonstrating the important distinction between kPCA and the
proposed method, and the potential advantage of the latter in applications.

###  E.4 Supplementary Figures to Section 4.2

Figure E.9 shows the scatter plots of the 2-dimensional embeddings based on
PCA and the proposed method (KEF stands for ”kernel eigenfunctions”) when
$p=3000$. The advantage of the proposed embedding over its linear counterpart,
such as informativeness and robustness to outliers, is visible and
significant.

Figure E.9: Comparison between principal component embeddings (PCs) and the
proposed kernel embeddings (KEFs) when $p=3000$, where the cells were labelled
and colored according to their actual time courses.

## References

  * [1] [author] Abbe, EmmanuelE., Fan, JianqingJ. and Wang, KaizhengK. (2022). An $\ell_{p}$ theory of PCA and spectral clustering. The Annals of Statistics 50 2359–2385.
  * [2] [author] Adeli, EhsanE., Wu, GuorongG., Saghafi, BehrouzB., An, LeL., Shi, FengF. and Shen, DinggangD. (2017). Kernel-based joint feature selection and max-margin classification for early diagnosis of Parkinson’s disease. Scientific Reports 7 1–14.
  * [3] [author] Amini, Arash AA. A. and Razaee, Zahra SZ. S. (2021). Concentration of kernel matrices with application to kernel spectral clustering. The Annals of Statistics 49 531–556.
  * [4] [author] Arias-Castro, EryE. and Pelletier, BrunoB. (2013). On the Convergence of Maximum Variance Unfolding. Journal of Machine Learning Research 14\.
  * [5] [author] Arora, SanjeevS., Hu, WeiW. and Kothari, Pravesh KP. K. (2018). An analysis of the t-sne algorithm for data visualization. In Conference On Learning Theory 1455–1462. PMLR.
  * [6] [author] Bao, ZhigangZ., Ding, XiucaiX., Wang, JingmingJ. and Wang, KeK. (2022). Statistical inference for principal components of spiked covariance matrices. The Annals of Statistics 50 1144–1169.
  * [7] [author] Bao, ZhigangZ., Ding, XiucaiX. and Wang, KeK. (2021). Singular vector and singular subspace distribution for the matrix denoising model. The Annals of Statistics 49 370–392.
  * [8] [author] Barp, AlessandroA., Oates, ChrisC., Porcu, EmilioE. and Girolami, MarkM. (2018). A Riemann-Stein kernel method. arXiv preprint arXiv:1810.04946.
  * [9] [author] Bayati, MohsenM. and Montanari, AndreaA. (2011). The Dynamics of Message Passing on Dense Graphs, with Applications to Compressed Sensing. IEEE Transactions on Information Theory 57 764-785.
  * [10] [author] Belkin, M.M. and Niyogi, P.P. (2003). Laplacian Eigenmaps for Dimensionality Reduction and Data Representation. Neural Computation 15 1373-1396.
  * [11] [author] Belkin, MikhailM. and Niyogi, ParthaP. (2007). Convergence of Laplacian eigenmaps. Advances in Neural Information Processing Systems 19 129\.
  * [12] [author] Berlinet, AlainA. and Thomas-Agnan, ChristineC. (2011). Reproducing Kernel Hilbert Spaces in Probability and Statistics. Springer Science & Business Media.
  * [13] [author] Billingsley, PatrickP. (2008). Probability and Measure. John Wiley & Sons.
  * [14] [author] Bishop, Christopher MC. M. (2006). Pattern Recognition and Machine Learning. Springer.
  * [15] [author] Blanchard, GillesG., Bousquet, OlivierO. and Zwald, LaurentL. (2007). Statistical properties of kernel principal component analysis. Machine Learning 66 259–294.
  * [16] [author] Boothby, William MW. M. (2003). An introduction to differentiable manifolds and Riemannian geometry, Revised second ed. Academic Press.
  * [17] [author] Bordenave, CharlesC. (2013). On Euclidean random matrices in high dimension. Electronic Communications in Probability 18 1–8.
  * [18] [author] Borg, IngwerI. and Groenen, Patrick JFP. J. (2005). Modern Multidimensional Scaling: Theory and Applications. Springer Science & Business Media.
  * [19] [author] Braun, Mikio L.M. L. (2006). Accurate Error Bounds for the Eigenvalues of the Kernel Matrix. Journal of Machine Learning Research 7 2303-2328.
  * [20] [author] Cai, T TonyT. T. and Ma, RongR. (2021). Theoretical Foundations of t-SNE for Visualizing High-Dimensional Clustered Data. arXiv preprint arXiv:2105.07536.
  * [21] [author] Calder, JeffJ., Trillos, Nicolas GarciaN. G. and Lewicka, MartaM. (2020). Lipschitz regularity of graph Laplacians on random data clouds. arXiv preprint arXiv:2007.06679.
  * [22] [author] Cannoodt, RobrechtR., Saelens, WouterW., Sichien, DorineD., Tavernier, SimonS., Janssens, SophieS., Guilliams, MartinM., Lambrecht, BartB., De Preter, KatleenK. and Saeys, YvanY. (2016). SCORPIUS improves trajectory inference and identifies novel modules in dendritic cell development. bioRxiv 079509\.
  * [23] [author] Chen, YuxinY., Chi, YuejieY., Fan, JianqingJ. and Ma, CongC. (2021). Spectral methods for data science: A statistical perspective. Foundations and Trends® in Machine Learning 14 566–806.
  * [24] [author] Cheng, Ming-YenM.-Y. and Wu, Hau-tiengH.-t. (2013). Local linear regression on manifolds and its geometric interpretation. Journal of the American Statistical Association 108 1421–1434.
  * [25] [author] Cheng, XiuyuanX. and Singer, AmitA. (2013). The spectrum of random inner-product kernel matrices. Random Matrices: Theory and Applications 2 1350010\.
  * [26] [author] Coifman, R. R.R. R. and Lafon, S.S. (2006). Diffusion maps. Applied and Computational Harmonic Analysis 21 5-30.
  * [27] [author] Cuturi, MarcoM. (2009). Positive definite kernels in machine learning. arXiv preprint arXiv:0911.5367.
  * [28] [author] Ding, XiucaiX. (2021). Spiked sample covariance matrices with possibly multiple bulk components. Random Matrices: Theory and Applications 10 2150014\.
  * [29] [author] Ding, XiucaiX. and Ma, RongR. (2022). Supplement to ”Learning Low-Dimensional Nonlinear Structures from High-Dimensional Noisy Data: An Integral Operator Approach”.
  * [30] [author] Ding, XiucaiX. and Wu, Hau-TiengH.-T. (2021). On the Spectral Property of Kernel-Based Sensor Fusion Algorithms of High Dimensional Data. IEEE Transactions on Information Theory 67 640-670.
  * [31] [author] Ding, XiucaiX. and Wu, Hau-TiengH.-T. (2021). How do kernel-based sensor fusion algorithms behave under high dimensional noise? arXiv preprint arXiv:2111.10940.
  * [32] [author] Ding, XiucaiX. and Wu, Hau-TiengH.-T. (2023). Impact of Signal-to-Noise Ratio and Bandwidth on Graph Laplacian Spectrum From High-Dimensional Noisy Point Cloud. IEEE Transactions on Information Theory 69 1899-1931.
  * [33] [author] Ding, XiucaiX. and Yang, FanF. (2021). Spiked separable covariance matrices and principal components. The Annals of Statistics 49 1113–1138.
  * [34] [author] Do, YenY. and Vu, VanV. (2013). The spectrum of random kernel matrices: universality results for rough and varying kernels. Random Matrices: Theory and Applications 2 1350005\.
  * [35] [author] Donoho, David LD. L. and Grimes, CarrieC. (2003). Hessian eigenmaps: Locally linear embedding techniques for high-dimensional data. Proceedings of the National Academy of Sciences 100 5591–5596.
  * [36] [author] Dunson, David BD. B., Wu, Hau-TiengH.-T. and Wu, NanN. (2020). Graph based Gaussian processes on restricted domains. arXiv preprint arXiv:2010.07242.
  * [37] [author] Dunson, David B.D. B., Wu, Hau-TiengH.-T. and Wu, NanN. (2021). Spectral convergence of graph Laplacian and heat kernel reconstruction in $L_{\infty}$ from random samples. Applied and Computational Harmonic Analysis 55 282-336.
  * [38] [author] El Karoui, NoureddineN. (2010). The spectrum of kernel random matrices. The Annals of Statistics 38 1–50.
  * [39] [author] El Karoui, NoureddineN. (2010). On information plus noise kernel random matrices. The Annals of Statistics 38 3191–3216.
  * [40] [author] El Karoui, N.N. and Wu, H. T.H. T. (2016). Connection graph Laplacian methods can be made robust to noise. The Annals of Statistics 44 346-372.
  * [41] [author] Fan, JianqingJ., Wang, WeichenW. and Zhong, YiqiaoY. (2018). An $\ell_{\infty}$ eigenvector perturbation bound and its application to robust covariance estimation. Journal of Machine Learning Research 18 1–42.
  * [42] [author] Fan, ZhouZ. and Montanari, AndreaA. (2019). The spectral norm of random inner-product kernel matrices. Probability Theory and Related Fields 173 27–85.
  * [43] [author] Fogel, FajwelF., d’Aspremont, AlexandreA. and Vojnovic, MilanM. (2014). Serialrank: Spectral ranking using seriation. Advances in Neural Information Processing Systems 27 900–908.
  * [44] [author] García Trillos, NicolásN., Gerlach, MoritzM., Hein, MatthiasM. and Slepcev, DejanD. (2020). Error estimates for spectral convergence of the graph Laplacian on random geometric graphs toward the Laplace-Beltrami operator. Foundations of Computational Mathematics 20 827–887.
  * [45] [author] Giné, EvaristE. and Koltchinskii, VladimirV. (2006). Empirical graph Laplacian approximation of Laplace–Beltrami operators: Large sample results. In High Dimensional Probability 238–259. Institute of Mathematical Statistics.
  * [46] [author] Guo, FanF., Yan, LiyingL., Guo, HongshanH., Li, LinL., Hu, BoqiangB., Zhao, YangyuY., Yong, JunJ., Hu, YuqiongY., Wang, XiaoyeX., Wei, YuanY. et al. (2015). The transcriptome and DNA methylome landscapes of human primordial germ cells. Cell 161 1437–1452.
  * [47] [author] Hofmann, ThomasT., Schölkopf, BernhardB. and Smola, Alexander JA. J. (2008). Kernel methods in machine learning. The Annals of Statistics 36 1171–1220.
  * [48] [author] Hofmeyr, David PD. P. (2019). Improving spectral clustering using the asymptotic value of the normalized cut. Journal of Computational and Graphical Statistics 28 980–992.
  * [49] [author] Horn, R. A.R. A. and Johnson, C. R.C. R. (2012). Matrix Analysis, 2nd ed. Cambridge University Press.
  * [50] [author] The Smithsonian Institute (2020). Mammuthus primigenius (blumbach). https://3d.si.edu/object/3d/mammuthus-primigenius-blumbach:341c96cd-f967-4540-8ed1-d3fc56d31f12.
  * [51] [author] Ji, ZhichengZ. and Ji, HongkaiH. (2016). TSCAN: Pseudo-time reconstruction and evaluation in single-cell RNA-seq analysis. Nucleic Acids Research 44 e117–e117.
  * [52] [author] Johnstone, Iain MI. M. (2001). On the distribution of the largest eigenvalue in principal components analysis. The Annals of Statistics 29 295–327.
  * [53] [author] Jolliffe, I. T.I. T. (2002). Principal Component Analysis. Springer Series in Statistics.
  * [54] [author] Jörgens, KonradK. (1982). Linear Integral Operators 13\. Pitman Advanced Publishing Program.
  * [55] [author] Kato, TosioT. (1982). A Short Introduction to Perturbation Theory for Linear Operators. Springer-Verlag, New York-Berlin.
  * [56] [author] Kobak, DmitryD. and Berens, PhilippP. (2019). The art of using t-SNE for single-cell transcriptomics. Nature Communications 10 1–14.
  * [57] [author] Koltchinskii, VladimirV. and Giné, EvaristE. (2000). Random matrix approximation of spectra of integral operators. Bernoulli 113–167.
  * [58] [author] Koltchinskii, VladimirV. and Giné, EvaristE. (2000). Random matrix approximation of spectra of integral operators. Bernoulli 6 113 – 167.
  * [59] [author] Kung, Sun YuanS. Y. (2014). Kernel Methods and Machine Learning. Cambridge University Press.
  * [60] [author] Lee, John AJ. A. and Verleysen, MichelM. (2007). Nonlinear Dimensionality Reduction. Springer.
  * [61] [author] Lee, John MJ. M. (2013). Introduction to smooth manifolds, Second ed. Springer.
  * [62] [author] Linderman, George CG. C. and Steinerberger, StefanS. (2019). Clustering with t-SNE, provably. SIAM Journal on Mathematics of Data Science 1 313–332.
  * [63] [author] Löffler, MatthiasM., Zhang, Anderson YA. Y. and Zhou, Harrison HH. H. (2021). Optimality of spectral clustering in the Gaussian mixture model. The Annals of Statistics 49 2506–2530.
  * [64] [author] Manton, Jonathan HJ. H. and Amblard, Pierre-OlivierP.-O. (2015). A primer on reproducing kernel hilbert spaces. Foundations and Trends® in Signal Processing 8 1–126.
  * [65] [author] Marčenko, Vladimir AV. A. and Pastur, Leonid AndreevichL. A. (1967). Distribution of eigenvalues for some sets of random matrices. Mathematics of the USSR-Sbornik 1 457\.
  * [66] [author] McInnes, LelandL., Healy, JohnJ. and Melville, JamesJ. (2018). Umap: Uniform manifold approximation and projection for dimension reduction. arXiv preprint arXiv:1802.03426.
  * [67] [author] Montanari, AndreaA. and Sun, NikeN. (2018). Spectral algorithms for tensor completion. Communications on Pure and Applied Mathematics 71 2381–2425.
  * [68] [author] Moon, Kevin RK. R., van Dijk, DavidD., Wang, ZhengZ., Gigante, ScottS., Burkhardt, Daniel BD. B., Chen, William SW. S., Yim, KristinaK., van den Elzen, AntoniaA., Hirn, Matthew JM. J., Coifman, Ronald RR. R. et al. (2019). Visualizing structure and transitions in high-dimensional biological data. Nature Biotechnology 37 1482–1492.
  * [69] [author] Nash, JohnJ. (1956). The imbedding problem for Riemannian manifolds. Annals of Mathematics 63 20–63.
  * [70] [author] Rosasco, LorenzoL., Belkin, MikhailM. and Vito, Ernesto DeE. D. (2010). On Learning with Integral Operators. Journal of Machine Learning Research 11 905-934.
  * [71] [author] Rousseeuw, Peter JP. J. (1987). Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. Journal of Computational and Applied Mathematics 20 53–65.
  * [72] [author] Roweis, S. T.S. T. and Saul, L. K.L. K. (2000). Nonlinear Dimensionality Reduction by Locally Linear Embedding. Science 290 2323-2326. 10.1126/science.290.5500.2323
  * [73] [author] Schölkopf, B.B. and Smola, A. J.A. J. (2002). Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond. Adaptive computation and machine learning. MIT Press.
  * [74] [author] Schölkopf, BernhardB., Smola, AlexanderA. and Müller, Klaus-RobertK.-R. (1998). Nonlinear component analysis as a kernel eigenvalue problem. Neural Computation 10 1299–1319.
  * [75] [author] Shawe-Taylor, JohnJ. and Cristianini, NelloN. (2004). Kernel Methods for Pattern Analysis. Cambridge University Press.
  * [76] [author] Shen, ChaoC., Lin, Yu-TingY.-T. and Wu, Hau-TiengH.-T. (2022). Robust and scalable manifold learning via landmark diffusion for long-term medical signal processing. Journal of Machine Learning Research 23 1–30.
  * [77] [author] Shen, ChaoC. and Wu, Hau-TiengH.-T. (2020). Scalability and robustness of spectral embedding: landmark diffusion is all you need. arXiv preprint arXiv:2001.00801.
  * [78] [author] Shi, TaoT., Belkin, MikhailM. and Yu, BinB. (2008). Data Spectroscopy: Learning Mixture Models Using Eigenspaces of Convolution Operators. In Proceedings of the 25th International Conference on Machine Learning. ICML ’08 936–943.
  * [79] [author] Shi, TaoT., Belkin, MikhailM. and Yu, BinB. (2009). Data spectroscopy: Eigenspaces of convolution operators and clustering. The Annals of Statistics 37 3960 – 3984.
  * [80] [author] Shiokawa, YukaY., Date, YasuhiroY. and Kikuchi, JunJ. (2018). Application of kernel principal component analysis and computational machine learning to exploration of metabolites strongly associated with diet. Scientific Reports 8 1–8.
  * [81] [author] Singer, A.A. (2006). From graph to manifold Laplacian: The convergence rate. Applied and Computational Harmonic Analysis 21 128-134.
  * [82] [author] Singer, AmitA. and Wu, H-TH.-T. (2012). Vector diffusion maps and the connection Laplacian. Communications on Pure and Applied Mathematics 65 1067–1144.
  * [83] [author] Singer, AmitA. and Wu, Hau-TiengH.-T. (2017). Spectral convergence of the connection Laplacian from random samples. Information and Inference: A Journal of the IMA 6 58–123.
  * [84] [author] Smale, StephenS. and Zhou, Ding-XuanD.-X. (2007). Learning Theory Estimates via Integral Operators and Their Approximations. Constructive Approximation 26 153-172.
  * [85] [author] Smale, SteveS. and Zhou, Ding-XuanD.-X. (2009). Geometry on probability spaces. Constructive Approximation 30 311–323.
  * [86] [author] Smith, AndrewA., Zha, HongyuanH. and Wu, Xiao-MingX.-M. (2008). Convergence and Rate of Convergence of a Manifold-Based Dimension Reduction Algorithm. In Advances in Neural Information Processing Systems 21\.
  * [87] [author] Tenenbaum, Joshua BJ. B., Silva, Vin deV. d. and Langford, John CJ. C. (2000). A global geometric framework for nonlinear dimensionality reduction. Science 290 2319–2323.
  * [88] [author] Van der Maaten, LaurensL. and Hinton, GeoffreyG. (2008). Visualizing data using t-SNE. Journal of Machine Learning Research 9\.
  * [89] [author] Van Der Maaten, LaurensL., Postma, EricE. and Van den Herik, JaapJ. (2009). Dimensionality Reduction: A Comparative Review. Journal of Machine Learning Research 10 13\.
  * [90] [author] Vershynin, RomanR. (2018). High-dimensional probability: An introduction with applications in data science 47\. Cambridge university press.
  * [91] [author] von Luxburg, UlrikeU. (2007). A tutorial on spectral clustering. Statistics and Computing 17 395–416.
  * [92] [author] von Luxburg, UlrikeU., Belkin, MikhailM. and Bousquet, OlivierO. (2008). Consistency of Spectral Clustering. The Annals of Statistics 36 555–586.
  * [93] [author] von Luxburg, UlrikeU., Bousquet, OlivierO. and Belkin, MikhailM. (2004). On the convergence of spectral clustering on random samples: the normalized case. In International Conference on Computational Learning Theory 457–471.
  * [94] [author] Wang, BoB., Zhu, JunjieJ., Pierson, EmmaE., Ramazzotti, DanieleD. and Batzoglou, SerafimS. (2017). Visualization and analysis of single-cell RNA-seq data by kernel-based similarity learning. Nature Methods 14 414–416.
  * [95] [author] Wang, YingfanY., Huang, HaiyangH., Rudin, CynthiaC. and Shaposhnik, YaronY. (2021). Understanding How Dimension Reduction Tools Work: An Empirical Approach to Deciphering t-SNE, UMAP, TriMap, and PaCMAP for Data Visualization. J. Mach. Learn. Res. 22 1–73.
  * [96] [author] Weinberger, Kilian QK. Q. and Saul, Lawrence KL. K. (2006). An introduction to nonlinear dimensionality reduction by maximum variance unfolding. In AAAI 6 1683–1686.
  * [97] [author] Williams, Christopher KIC. K. and Rasmussen, Carl EdwardC. E. (2006). Gaussian processes for machine learning. MIT press Cambridge, MA.
  * [98] [author] Wormell, Caroline LC. L. and Reich, SebastianS. (2021). Spectral convergence of diffusion maps: Improved error bounds and an alternative normalization. SIAM Journal on Numerical Analysis 59 1687–1734.
  * [99] [author] Wu, Hau-TiengH.-T. and Wu, NanN. (2018). Think globally, fit locally under the manifold setup: Asymptotic analysis of locally linear embedding. The Annals of Statistics 46 3805–3837.
  * [100] [author] Yu, YiY., Wang, TengyaoT. and Samworth, Richard JR. J. (2015). A useful variant of the Davis–Kahan theorem for statisticians. Biometrika 102 315–323.
  * [101] [author] Zhang, ZhenyueZ. and Zha, HongyuanH. (2004). Principal manifolds and nonlinear dimensionality reduction via tangent space alignment. SIAM Journal on Scientific Computing 26 313–338.
  * [102] [author] Zhou, ZhixinZ. and Amini, Arash A.A. A. (2019). Analysis of spectral clustering algorithms for community detection: the general bipartite setting. Journal of Machine Learning Research 20 1-47.


Concept | Active Neurons
---|---
|