Proof of Lemma 1. Following the idea of [19], take \(\Omega\) as an arbitrary bounded open set such that \(B_{x}(c_0)\subseteq \Omega\) for all \(x\in\mathcal{X}\) for some \(c_0>0\) and we can assume that \(f\in H^{1}(\Omega,g)\) and \(\|f\|_{H^{1}(\Omega,g)}\lesssim \|f\|_{H^{1}(\mathcal{X},g)}\) without loss of generality due to the existence of an extension operator \(E:H^{1}(\mathcal{X},g)\rightarrow H^{1}(\Omega,g)\) such that \(Ef\) satisfies these properties, see Theorem 1 in Chapter 5.4 in [42]. Also, since \(C^{\infty}(\Omega)\) is dense in \(H^{1}(\Omega,g)\) and the integral in Lemma 1 is continuous in \(H^{1}(\Omega,g)\), we can assume \(f_{g}:=f/g^{r}\in C^{\infty}(\Omega)\) so that \[\begin{align} f_{g}(x')-f_{g}(x)=\int_{0}^{1}\nabla f_{g}(x+t(x'-x))^{T}(x'-x)dt. \end{align}\] Then, we have by symmetry in the first step: \[\begin{align} \label{s1decomp} &2\langle L_{w,\epsilon}f,f\rangle_{g^{p-r}} \nonumber\\&=\frac{1}{\epsilon^{d+2}}\int_{\mathcal{X}}\int_{\mathcal{X}}\frac{\eta\left(\frac{\|x-y\|}{\epsilon}\right)}{g(x)^{1-q/2}g(y)^{1-q/2}}\left|\frac{f(x)}{g(x)^{r}}-\frac{f(y)}{g(y)^{r}}\right|^{2}g(x)g(y)dxdy\nonumber\\&=\frac{1}{\epsilon^{d+2}}\int_{\mathcal{X}}\int_{\mathcal{X}}\frac{\eta\left(\frac{\|x-y\|}{\epsilon}\right)}{g(x)^{1-q/2}g(y)^{1-q/2}}\left(\int_{0}^{1}\nabla f_{g}(y+t(x-y))^{T}(x-y)dt\right)^{2}g(x)g(y)dxdy\nonumber\\&\le \frac{1}{\epsilon^{d+2}}\int_{\mathcal{X}}\int_{\mathcal{X}}\int_{0}^{1}\frac{\eta\left(\frac{\|x-y\|}{\epsilon}\right)}{g(x)^{1-q/2}g(y)^{1-q/2}}\left(\nabla f_{g}(y+t(x-y))^{T}(x-y)\right)^{2}g(x)g(y)dtdxdy\nonumber\\&\le \int_{\mathcal{X}}\int_{B_{\mathbf{0}}(1)}\int_{0}^{1}\left(\nabla f_{g}(y+\epsilon tz)^{T}z\right)^{2}\frac{\eta(\|z\|)}{g(y+\epsilon z)^{1-q/2}g(y)^{1-q/2}}g(y+\epsilon z)g(y)dtdzdy,\nonumber\\&\qquad \text{with}\;(x-y)/\epsilon=z\nonumber\\&\lesssim \int_{\mathcal{X}}\int_{B_{\mathbf{0}}(1)}\int_{0}^{1}\left(\nabla f_{g}(y+\epsilon tz)^{T}z\right)^{2}\eta(\|z\|)g(y+\epsilon tz)^{q}dtdzdy\nonumber\\&\le \int_{\Omega}\int_{B_{\mathbf{0}}(1)}\int_{0}^{1}\left(\nabla f_{g}(\tilde{y})^{T}z\right)^{2}\eta(\|z\|)g(\tilde{y})^{q}dtdzd\tilde{y},\quad \tilde{y}=y+\epsilon tz\in \Omega. \end{align}\tag{8}\] Since we have \(\left(\nabla f_{g}(\tilde{y})^{T}z\right)^{2}=\left(\sum_{i=1}^{d}(\nabla f_{g}(\tilde{y}))^{(i)}z^{(i)}\right)^{2}\) and \(\eta(\|z\|)\) is invariant with respect to the rotation, it yields that \[\begin{align} \label{s1decomppart} \int_{B_{\mathbf{0}}(1)}\left(\nabla f_{g}(\tilde{y})^{T}z\right)^{2}\eta(\|z\|)dz&=\sum_{i,j=1}^{d}(\nabla f_{g}(\tilde{y}))^{(i)}(\nabla f_{g}(\tilde{y}))^{(j)}\int_{B_{\mathbf{0}}(1)}z^{(i)}z^{(j)}\eta(\|z\|)dz\nonumber\\&=\sum_{i=1}^{d}\left((\nabla f_{g}(\tilde{y}))^{(i)}\right)^{2}\int_{B_{\mathbf{0}}(1)}\left(z^{(i)}\right)^{2}\eta(\|z\|)dz\nonumber\\&=\sigma_{1}\left\|\nabla\left(\frac{f(\tilde{y})}{g(\tilde{y})^{r}}\right)\right\|^{2}. \end{align}\tag{9}\] Plugging 9 in 8 , we conclude \[\begin{align} 2\langle L_{w,\epsilon}f,f\rangle_{g^{p-r}}\lesssim \sigma_1M^{2}. \end{align}\] This finishes the proof. ◻

Proof. By Cauchy-Schwarz inequality, we have \[\begin{align} |L_{w,\epsilon}f(x)|^2&=\frac{1}{\epsilon^{2(d+2)}}\left(\int_{\mathcal{U}}g(x)^{1-p}\frac{\eta\left(\frac{\|x-z\|}{\epsilon}\right)}{g(x)^{1-q/2}g(z)^{1-q/2}}(g(x)^{-r}f(x)-g(z)^{-r}f(z))g(z)dz\right)^{2}\\&\lesssim \frac{1}{\epsilon^{2(d+2)}}g(x)^{2(1-p)}\int_{\mathcal{U}}\frac{\eta\left(\frac{\|x-z\|}{\epsilon}\right)}{g(x)^{1-q/2}g(z)^{1-q/2}}(f_{g}(x)-f_{g}(z))^2dz\cdot\int_{\mathcal{X}}\frac{\eta\left(\frac{\|x-z\|}{\epsilon}\right)}{g(x)^{1-q/2}g(z)^{1-q/2}}dz\\&\lesssim \frac{2\sigma_0}{\epsilon^{4+d}}g(x)^{2(q-p)}\int_{\mathcal{U}}\eta\left(\frac{\|x-z\|}{\epsilon}\right)(|f_{g}(x)|^2+|f_{g}(z)|^2)dz. \end{align}\] Then, we have \[\begin{align} \|L_{w,\epsilon}f\|_{L^{2}(\mathcal{U},g^{p+r})}^{2}&=\int_{\mathcal{U}}g^{p-r}(x)|L_{w,\epsilon}f(x)|^{2}dx\\&\lesssim \frac{2}{\epsilon^{4+d}}\int_{\mathcal{U}}\int_{\mathcal{U}}g(x)^{2(q-p)+p-r}(x)\eta\left(\frac{\|x-z\|}{\epsilon}\right)(|f_{g}(x)|^2+|f_{g}(z)|^2)dzdx\\&\lesssim \frac{2}{\epsilon^{4+d}}\int_{\mathcal{U}}\int_{\mathcal{U}}\eta\left(\frac{\|x-z\|}{\epsilon}\right)(|f_{g}(x)|^2+|f_{g}(z)|^2)dzdx\\&\lesssim \frac{4}{\epsilon^{4+d}}\int_{\mathcal{U}}\int_{\mathcal{U}}\eta\left(\frac{\|x-z\|}{\epsilon}\right)|f_{g}(x)|^2dzdx\\&\le \frac{4}{\epsilon^{4}}\int_{\mathcal{U}}|g(x)^{p+r}(x)f_{g}(x)|^2dx\\&\lesssim \frac{4}{\epsilon^{4}}\|f_{g}\|_{L^{2}(U,g^{p+r})}^{2} \end{align}\] ◻

Proof. Note that \[\begin{align} E_{w,\epsilon}(f;\mathcal{U})&=\frac{1}{\epsilon^{d+2}}\int_{\mathcal{U}}\int_{\mathcal{U}}(g(x)^{-r}f(x)-g(z)^{-r}f(z))^2\frac{\eta\left(\frac{\|x-z\|}{\epsilon}\right)}{g(x)^{1-q/2}g(z)^{1-q/2}}g(x)g(z)dxdz\\&\le \frac{2}{\epsilon^{d+2}}\int_{\mathcal{U}}\int_{\mathcal{U}}(|g(x)^{-r}f(x)|^2+|g(z)^{-r}f(z)|^2)\frac{\eta\left(\frac{\|x-z\|}{\epsilon}\right)}{g(x)^{1-q/2}g(z)^{1-q/2}}g(x)g(z)dxdz\\&= \frac{4}{\epsilon^{d+2}}\int_{\mathcal{U}}\int_{\mathcal{U}}|g(x)^{-r}f(x)|^2\eta\left(\frac{\|x-z\|}{\epsilon}\right)g(x)^{q/2}g(z)^{q/2}dxdz\\&\lesssim \frac{4}{\epsilon^{d+2}}\int_{\mathcal{U}}\int_{\mathcal{U}}|g(x)^{-r}f(x)|^2\eta\left(\frac{\|x-z\|}{\epsilon}\right)g(x)^{p+r}dxdz\\&\lesssim \frac{4}{\epsilon^{2}}\int_{\mathcal{U}}|g(x)^{-r}f(x)|^2g(x)^{p+r}dx. \end{align}\] ◻

Proof. Without loss of generality, we assume both \(g\) and \(f\) are \(C^{\infty}(\mathcal{X})\) due to the fact that \(C^{\infty}(\mathcal{X})\) is dense in both \(H^{s}(\mathcal{X},g)\) and \(C^{s-1}(\mathcal{X})\) and the norm in the statements is continuous with respect to \(\|\cdot\|_{H^{s}(\mathcal{X},g)}\) and \(\|\cdot\|_{C^{s-1}(\mathcal{X})}\).

Actually, we claim the following stronger result: for \(t<s/2\) and every \(x\in \mathcal{X}_{t\epsilon}\), \[\begin{align} \label{inter} L_{w,\epsilon}^{t}f(x)=\sigma_1^{t}\mathcal{L}_{w}f(x)+ \sum_{j=1}^{\lfloor (s-1)/2\rfloor-t}r_{2(j+t)}(x)\epsilon^{2j}+r_{s}(x)\epsilon^{s-2t}, \end{align}\tag{10}\] for some functions \(r_{j}\) such that \[\begin{align} \label{rjbound} \|r_j\|_{H^{s-j}(\mathcal{X}_{t\epsilon},g)}\le C\|g\|_{C^{s-1}(\mathcal{X})}^{t}M. \end{align}\tag{11}\] Note that the dependence of the functions \(r_{j}\) on \(t\) is suppressed in the notation.

The key idea underlying the proof of (10 ) is to consider the following Taylor expansion. For an \(s\)-times differentiable function \(F:\mathcal{X}\rightarrow \mathbb{R}\) and \(x\in\mathcal{X}\), define the following operator \(d_{x}^{s}\): \[\begin{align} (d_{x}^{s}F)(z):=\sum_{|\alpha|=s}D^{\alpha}F(x)z^{\alpha}. \end{align}\] Also, define \(d^{s}F:=\sum_{|\alpha|=s}D^{\alpha}F\). Then, for \(\phi\in C^{s}(\mathcal{X})\) and some \(h>0\), \(z\in\mathcal{X}_{h}\), \(x\in B_{z}(h)\), the Taylor expansion at \(z\) is given as: \[\begin{align} \phi(x)=\phi(z)+\sum_{j=1}^{s-1}\frac{1}{j!}(d_{x}^{j}\phi)(x-z)+R_{s}(x,z;\phi). \end{align}\] Here, we note that \((d_{x}^{j}\phi)(z)\) is a polynomial of degree \(j\) and we have for any \(y\in\mathbb{R}\): \[\begin{align} (d_{x}^{j}\phi)(yz)=y^j (d_{x}^{j}\phi)(z). \end{align}\] The remainder term \(R_{j}(x,z;\phi)\) is \[\begin{align} R_{j}(x,z;\phi):=\frac{1}{(j-1)!}\int_{0}^{1}(1-\theta)^{j-1}(d_{z+\theta(x-z)}^{j}\phi)(x-z)d\theta, \end{align}\] such that for any \(x^{*}\in B_{\mathbf{0}}(1)\), \[\begin{align} \underset{x\in\mathcal{X}_{h}}{\sup}|R_{j}(x,x+hx^{*};\phi)|\le Ch^{j}\|\phi\|_{C^{j}(\mathcal{X})}, \end{align}\] and \[\begin{align} \int_{\mathcal{X}_{h}}|R_{j}(z+\theta hx,z;\phi)|^{2}dz\le h^{2j}\int_{\mathcal{X}_{h}}\int_{0}^{1}|(d_{z+\theta hx}^{j}\phi)(z)|^{2}d\theta dz\le h^{2j}\|d^{j}\phi\|_{L^{2}(\mathcal{X})}^{2}. \end{align}\]

Now, we apply the above Taylor expansion on the function \(f_{g}(x):=f(x)/g(x)^{r}\) up to order \(s\) and the function \(g^{q/2}(x)\) up to order \(S\) in \(L_{w,\epsilon}f(x)\), where \(S=1\) if \(s=1\) and otherwise \(S=s-1\) and obtain: \[\begin{align} L_{w,\epsilon}f(x)&=\frac{1}{\epsilon^{d+2}}\sum_{j_1=1}^{s-1}\sum_{j_2=0}^{S-1}\frac{1}{j_1!j_2!}\int_{\mathcal{X}}g(x)^{q/2-p}\eta\left(\frac{\|x-z\|}{\epsilon}\right)(d_{x}^{j_1}f_{g})(x-z)(d_{x}^{j_2}g^{q/2})(z-x)dz\\&+\frac{1}{\epsilon^{d+2}}\sum_{j=1}^{s-1}\frac{1}{j!}\int_{\mathcal{X}}g(x)^{q/2-p}\eta\left(\frac{\|x-z\|}{\epsilon}\right)(d_{x}^{j_1}f_{g})(x-z)R_{S}(x,z;g^{q/2})dz\\&+\frac{1}{\epsilon^{d+2}}\int_{\mathcal{X}}g(x)^{p/2-p}\eta\left(\frac{\|x-z\|}{\epsilon}\right)R_{s}(x,z;f_g)g(z)^{q/2}dz. \end{align}\] Now, with the transformation \(y=(z-x)/\epsilon\), we have \[\begin{align} L_{w,\epsilon}f(x)&=-\frac{1}{\epsilon^{2}}\sum_{j_1=1}^{s-1}\sum_{j_2=0}^{S-1}\frac{\epsilon^{j_1+j_2}}{j_1!j_2!}\int_{B_{\mathbf{0}}(1)}g(x)^{q/2-p}\eta\left(\|y\|\right)(d_{x}^{j_1}f_{g})(y)(d_{x}^{j_2}g^{q/2})(y)dy\\&-\frac{1}{\epsilon^{2}}\sum_{j=1}^{s-1}\frac{\epsilon^{j}}{j!}\int_{B_{\mathbf{0}}(1)}g(x)^{q/2-p}\eta\left(\|y\|\right)(d_{x}^{j_1}f_{g})(y)R_{S}(x,\epsilon y+x;g^{q/2})dy\\&+\frac{1}{\epsilon^{2}}\int_{B_{\mathbf{0}}(1)}g(x)^{p/2-p}\eta\left(\|y\|\right)R_{s}(x,\epsilon y+x;f_g)g(\epsilon y+x)^{q/2}dy\\&=:L_{1}(x)+L_{2}(x)+L_{3}(x). \end{align}\] We will now prove 10 by induction on \(t,\) and throughout this proof, with a slight abuse of notation, the functions \(r_{j}\) in 10 may vary from line to line depending on \(t\) at the induction step but they will always satisfy the condition 11 as we are only interested in the bounds.

Firstly, we start with \(L_{1}(x)\). If \(s=1\), we can see \(L_{1}(x)=0\). Therefore, in the following, we only focus on \(s\ge 2\). Now, we define \[\begin{align} l_{j_1,j_2}(x):=\int_{B_{\mathbf{0}}(1)}g(x)^{q/2-p}\eta(\|y\|)(d_{x}^{j_1}f_{g})(y)(d_{x}^{j_2}g^{q/2})(y)dy, \end{align}\] such that \[\begin{align} L_{1}(x)=-\frac{1}{\epsilon^{2}}\sum_{j_1=1}^{s-1}\sum_{j_2=0}^{S-1}\frac{\epsilon^{j_1+j_2}}{j_1!j_2!}l_{j_1,j_2}(x). \end{align}\] Since \((d_{x}^{j}f_{g})(y)\) is a polynomial of degree \(j\), \(l_{j_1,j_2}\) actually depends on the sum \(j_1+j_2\) and \(d_{x}^{j_1}d_{x}^{j_2}\) is an order \(j_1+j_2\) multivariate monomial. Therefore, when \(j_1+j_2\) is odd, we have \[\begin{align} l_{j_1,j_2}(x)=0. \end{align}\] Then, when \(s=2\), we have \(j_1+j_2=1\) and \(L_{1}(x)=0\). As for \(s\ge 3\), we notice that the lowest order term of \(L_{11}(x)\) is from \(j_1+j_2=2\), which means either \(j_1=1,j_2=1\) or \(j_1=2,j_2=0\). We have \[\begin{align} l_{1,1}(x)&=\int_{B_{\mathbf{0}}(1)}g(x)^{q/2-p}\eta(\|y\|)(d_{x}^{1}f_{g})(y)(d_{x}^{1}g^{q/2})(y)dy\\&=\sum_{i_1=1,i_2=1}^{d}g(x)^{q/2-p}(Df_{g}(x))^{(i_1)}(Dg^{q/2}(x))^{(i_{2})}\int_{B_{\mathbf{0}}(1)}\|y\|^{2}\eta(\|y\|)dy, \end{align}\] and \[\begin{align} \frac{1}{2}l_{2,0}(x)&=\frac{1}{2}\int_{B_{\mathbf{0}}(1)}g(x)^{q/2-p}\eta(\|y\|)(d_{x}^{2}f_{g})(y)(d_{x}^{2}g^{q/2})(y)dy\\&\\&=\frac{1}{2}\sum_{i=1}^{d}g(x)^{q/2-p}((Df_{g}(x))^{(i)})^{2}g(x)^{q/2}\int_{B_{\mathbf{0}}(1)}\|y\|^{2}\eta(\|y\|)dy. \end{align}\] Therefore, we have by definition: \[\begin{align} \mathcal{L}_{w}f(x)=-\frac{1}{2g(x)^{p}}\left(\nabla g(x)^{q}\cdot \nabla\left(\frac{f(x)}{g(x)^{r}}\right)+g(x)^{q}\Delta\left(\frac{f(x)}{g(x)^{r}}\right)\right), \end{align}\] and \[\begin{align} -(l_{1,1}(x)+\frac{1}{2}l_{2,0}(x))=\sigma_1\mathcal{L}_{w}f(x). \end{align}\] This is exactly the leading term. We remark here that in [1], the negative sign is missing, which does not actually give the Laplacian operator by the leading term. Now, it remains to bound the higher order terms with \(j_1+j_2>2\). We will show that \[\begin{align} L_{1}(x)=\sigma_1\mathcal{L}_{w}+\sum_{j=1}^{\lfloor (s-1)/2\rfloor-1}r_{2(j+1)}(x)\epsilon^{2j}+r_{s}(x)\epsilon^{s-2}. \end{align}\] It suffices to show for \(j_1+j_2>2\), \(l_{j_1,j_2}\) satisfies 11 for \(j=\min\{j_1+j_2-2,s-2\}\). Through the multi-index notation, we write that \[\begin{align} l_{j_1,j_2}(x)=g(x)^{q/2-p}\sum_{|\alpha_1|=j_1,|\alpha_2|=j_2}D^{\alpha_1}f_g(x)D^{\alpha_2}g^{q/2}(x)\int_{B_{\mathbf{0}}(1)}y^{\alpha_1}y^{\alpha_2}\eta(\|y\|)dy, \end{align}\] where \(|\int_{B_{\mathbf{0}}(1)}y^{\alpha_1}y^{\alpha_2}\eta(\|y\|)dy|<\infty\) for all \(\alpha_1,\alpha_2\). Then, by Hölder’s inequality, we have for \(|\alpha_1|=j_1\), \(|\alpha_{2}|=j_2\), \[\begin{align} \|g(x)^{q/2-p}D^{\alpha_1}f_g D^{\alpha_2}g^{q/2}\|_{H^{s-(j+2)}(\mathcal{X},g)}&\lesssim \|D^{\alpha_1}f_g \|_{H^{s-(j+2)}(\mathcal{X},g)}\|g(x)^{q/2-p}D^{\alpha_2}g^{q/2}\|_{C^{s-(j+2)}(\mathcal{X})}\\&\lesssim \|D^{\alpha_1}f_g \|_{H^{s-j_1}(\mathcal{X},g)}\|D^{\alpha_2}g^{q/2}\|_{C^{s-(j_2+1)}(\mathcal{X})}\\&\le M\|g\|_{C^{s-1}}. \end{align}\] Summing over all \(|\alpha_1|=j_1\) and \(|\alpha_2|=j_2\), we obtain that \(l_{j_1,j_2}\) satisfies 11 .

Next, as for \(L_2(x)\), note that if \(s=1\), \(L_{2}(x)=0\). We want to show that for \(s\ge 2\), \[\begin{align} \|L_{2}\|_{L^{2}(\mathcal{X}_{\epsilon},g^{p+r})}\le C\epsilon^{s-2}M\|g\|_{C^{s-1}(\mathcal{X})}. \end{align}\] Clearly, if \(s=1\), \(L_2(x)=0\). Now, for \(s\ge 2\), we have \(S=s-1\) and since \(|R_{s-1}(x,x+\epsilon x^{*})|\le C\epsilon^{s-1}\|g\|_{C^{s-1}(\mathcal{X})}\) for any \(x^{*}\in B_{\mathbf{0}}(1)\) and \(d_{x}^{j}(\cdot)\) is a \(j\)-homogeneous function, we have \[\begin{align} |L_{2}(x)|&\le \sum_{j=1}^{s-1}\frac{\epsilon^{j-2}}{j!}\int_{B_{\mathbf{0}}(1)}g(x)^{q/2-p}\eta\left(\|y\|\right)|(d_{x}^{j_1}f_{g})(y)|\cdot|R_{S}(x,\epsilon y+x;g^{q/2})|dy\\&\le C\epsilon^{s-2}\|g\|_{C^{s-1}(\mathcal{X})}\sum_{j=1}^{s-1}\frac{1}{j!}\int_{B_{\mathbf{0}}(1)}g(x)^{q/2-p}\eta\left(\|y\|\right)|(d_{x}^{j_1}f_{g})(y)|dy. \end{align}\] Moreover, we have by Cauchy–Schwarz inequality, \[\begin{align} &\int_{\mathcal{X}_{\epsilon}}g(x)^{p+r}\left(\int_{B_{\mathbf{0}}(1)}g(x)^{q/2-p}\eta\left(\|y\|\right)|(d_{x}^{j_1}f_{g})(y)|dy\right)^{2}dx\\&\le \int_{\mathcal{X}_{\epsilon}}g(x)^{q-p+r}\left(\int_{B_{\mathbf{0}}(1)}\eta\left(\|y\|\right)|(d_{x}^{j_1}f_{g})(y)|^2dy\right)\left(\int_{B_{\mathbf{0}}(1)}\eta(\|y\|)dy\right)dx\\&\le \sigma_0 \int_{B_{\mathbf{0}}(1)}\int_{\mathcal{X}_{\epsilon}}g(x)^{q-p+r}\eta(\|y\|)((d^{j}f_g)(x))^{2}dxdy\\&\lesssim \sigma_0^{2}\int_{\mathcal{X}_{\epsilon}}g(x)^{p+r}((d^{j}f_g)(x))^{2}dx\\&=\sigma_{0}^{2}\|d^{j}f_g\|_{L^{2}(\mathcal{X}_{\epsilon},g^{p+r})}^{2}, \end{align}\] where in the last step, we use the fact that \(|d_{x}^{j}f(y)|\le |d^{j}f(x)|\) for all \(y\in B_{\mathbf{0}}(1)\). Therefore, it yields that \[\begin{align} &\int_{\mathcal{X}_{\epsilon}}g(x)^{p+r}|L_{2}(x)|^{2}dx\\&\le C\left(\epsilon^{s-2}\|g\|_{C^{s-1}(\mathcal{X})}\right)^{2}\sum_{j=1}^{s-1}\int_{\mathcal{X}_{\epsilon}}g(x)^{p+r}\left(\frac{1}{j!}\int_{B_{\mathbf{0}}(1)}g(x)^{q/2-p}\eta\left(\|y\|\right)|(d_{x}^{j_1}f_{g})(y)|dy\right)^{2}dx\\&\le C\left(\epsilon^{s-2}\|g\|_{C^{s-1}(\mathcal{X})}\right)^{2}\sum_{j=1}^{s-1}\|d^{j}f_{g}\|_{L^{2}(\mathcal{X}_{\epsilon},g^{p+r})}^{2}. \end{align}\] We obtain the desired bound.

Finally, similar to \(L_{2}(x)\), we obtain the same bound for \(L_3(x).\) Combining the obtained bounds for \(L_1(x) - L_3(x),\) we obtain 10 for \(t=1\).

Now, we perform the induction step. Assuming the bound 10 holds up to some \(t<s/2\), we want to show it also holds for \(t+1\), with \(t+1 < s/2\). For convenience, we introduce the following notation: for any \(1\le j\le l\le s\), denote by \(r_{j,l}(x)=r_{(s-l)+j}(x)\). Note again that the functions \(r_{j,l}\) implicitly depend on \(t\) at the induction step thus they may vary in the below arguments from line to line. For a function \(r\in H^{l}(\mathcal{X}_{t\epsilon},g;C\|g\|_{C^{s-1}(\mathcal{X})}^{t}M)\) for some \(l\le s\), if \(l\le 2\), we have by the inductive hypothesis that for any \(x\in\mathcal{X}_{(t+1)\epsilon}\), \[\begin{align} L_{w,\epsilon}r(x)=r_{l}^{l}(x)\epsilon^{l-2}. \end{align}\] On the other hand, if \(2<l\le s\), then by the inductive hypothesis, it holds that for any \(x\in\mathcal{X}_{(t+1)\epsilon}\), \[\begin{align} \label{inductive1} L_{w,\epsilon}r(x)=\sigma_{1}\mathcal{L}_{w}r(x)+\sum_{j=1}^{\lfloor (l-1)/2\rfloor-1}r_{2j+2,l}(x)\epsilon^{2j}+r_{l,l}(x)\epsilon^{l-2}. \end{align}\tag{12}\] Then, we have \[\begin{align} \label{inductive2} L_{w,\epsilon}^{t+1}f(x)&=(L_{w,\epsilon}\circ L_{w,\epsilon}^{t}f)(x)\nonumber\\&=\sigma_1^{t}L_{w,\epsilon}\mathcal{L}_{w}^{t}f(x)+\sum_{j=1}^{\lfloor (s-1)/2\rfloor-k}L_{w,\epsilon}r_{2(j+t)}(x)\epsilon^{2j}+L_{w,\epsilon}r_{s}(x)\epsilon^{s-2t}. \end{align}\tag{13}\] In the following, we will bound these terms on the right-hand side individually. First of all, since \(\mathcal{L}_{w}^{t}f\in H^{s-2t}(\mathcal{X},g;C\|g\|_{C^{s-1}(\mathcal{X})}^{t}M)\), applying 12 yields \[\begin{align} \label{t4311} L_{w,\epsilon}\mathcal{L}_{w}^{t}f(x)&=\sigma_1 \mathcal{L}_{w}^{t+1}f(x)+\sum_{j=1}^{\lfloor (s-2t-1)/2\rfloor-1}r_{2j+2,s-2t}(x)\epsilon^{2j}+r_{s-2t,s-2t}(x)\epsilon^{s-2t-2}\nonumber\\&=\sigma_1 \mathcal{L}_{w}^{t+1}f(x)+\sum_{j=1}^{\lfloor (s-1)/2\rfloor-(t+1)}r_{2(t+1+j)}(x)\epsilon^{2j}+r_{s}(x)\epsilon^{s-2(t+1)}, \end{align}\tag{14}\] where we apply the fact mentioned before that \(r_{j,l}(x)=r_{(s-l)+j}(x)\).

Next, suppose \(j<\lfloor (s-1)/2\rfloor-t\). We apply 12 and obtain \[\begin{align} L_{w,\epsilon}r_{2(j+t)}(x)&=\sigma_1 \mathcal{L}_{w}r_{2(j+t)}(x)+\sum_{i=1}^{\lfloor (s-2j-2t-1)/2\rfloor-1}r_{2i+2,s-2(j+t)}(x)\epsilon^{2i}\\&\quad+r_{s-2(j+t),s-2(j+t)}(x)\epsilon^{s-2(j+t)-2}\\&=r_{2(j+t+1)}(x)+\sum_{i=1}^{\lfloor (s-1)/2\rfloor-(j+t+1)}r_{2(i+j+t+1)}(x)\epsilon^{2i}+r_{s}(x)\epsilon^{s-2(j+t+1)}, \end{align}\] where we use the fact that \(r_{j,l}(x)=r_{(s-l)+j}(x)\) and \(\sigma_1\mathcal{L}_w r_{2(j+t)}(x)=r_{2,s-2(j+t)}(x)=r_{2(j+t+1)}(x)\). Therefore, we have \[\begin{align} \label{t4312} L_{w,\epsilon}r_{2(j+t)}(x)\epsilon^{2j}=r_{2(j+t+1)}(x)\epsilon^{2j}+\sum_{m=1}^{\lfloor (s-1)/2\rfloor-(k+1)}r_{2(m+t+1)}(x)\epsilon^{2m}+r_{s}(x)\epsilon^{s-2(k+1)}, \end{align}\tag{15}\] where the last equality is by changing the variable \(m=i+j\). Moreover, when \(j=\lfloor (s-1)/2\rfloor-t\), we have \(2(j+t)=2\lfloor (s-1)/2\rfloor\) and we simply calculate that \[\begin{align} \label{t4313} L_{w,\epsilon}r_{2(j+t)}(x)\epsilon^{2j}=r_{s-2(j+t)}^{s-2(j+t)}(x)\epsilon^{s-2(j+k)}\epsilon^{2j}=r_{s}(x)\epsilon^{s-2(k+1)}. \end{align}\tag{16}\] Finally, according to 12 , we have \[\begin{align} \label{t4314} L_{w,\epsilon}r_{s}(x)\epsilon^{s-2t}=r_{s}(x)\epsilon^{s-2(t+1)}. \end{align}\tag{17}\] Combining 14 17 with 13 , we obtain the proof for \(t+1\). ◻

Proof. Note that according to Lemma 2, we have \[\begin{align} \|L_{w,\epsilon}^{t}f\|_{L^{2}(\partial_{t\epsilon} \mathcal{X},g^{p+r})}^{2}\lesssim \frac{1}{\epsilon^{4}} \|L_{w,\epsilon}^{t-1}f\|_{L^{2}(\partial_{t\epsilon} \mathcal{X},g^{p+r})}^{2}\lesssim\ldots\lesssim \frac{1}{\epsilon^{4t}}\|f_{g}\|_{L^{2}(\partial_{t\epsilon} \mathcal{X},g^{p+r})}^{2}. \end{align}\] Therefore, it suffices to show for all \(\epsilon<c\), \[\begin{align} \label{boundarybound} \|f_{g}\|_{L^{2}(\partial_{t\epsilon} \mathcal{X},g^{p+r})}^{2}\lesssim\epsilon^{2s}\|f\|_{H^{s}(\mathcal{X},g)}^{2}. \end{align}\tag{18}\] In order to deal with \(f_{g}\) near the boundary, we will take a similar procedure used in [1] and [53] as follows. With loss of generality, we take \(t=1\) as one can view \(\epsilon<c/t\) for proving for the general case.

Step I: Local patch. We assume that for some \(c_0>0\) and a Lipschitz mapping \(\phi:\mathbb{R}^{d-1}\rightarrow [-c_0,c_0]\) and since \(f\in H^{s}_{0}(\mathcal{X},g;M)\), without loss of generality, we can assume that \(f_{g}\in C_{c}^{\infty}(U_{\psi}(c_0))\) with \[\begin{align} U_{\psi}(c_0):=\{y\in Q(0,c_0):\psi(y^{(-d)})\le y^{(d)}\}, \end{align}\] where \(Q(0,c_0)\) is the \(d\)-dimensional hypercube of side length \(c_0\) centered at \(\mathbf{0}\). Now, following step 1 in [1] by replacing \(f\) as \(f_{g}\), we have \[\begin{align} |f_{g}(y)|^2&\lesssim \epsilon^{2(s-1)}\left(\int_{\psi(y^{(-d)})}^{y^{(d)}}|(D^sf_{g}(y^{(-d)},z))^{(d)}|dz\right)^{2}\\&\lesssim \epsilon^{2s-1}\int_{\psi(y^{(-d)})}^{y^{(d)}}|(D^sf_{g}(y^{(-d)},z))^{(d)}|^{2}dz. \end{align}\] Then, we obtain: \[\begin{align} \label{step1:1} \int_{V_{\psi}(\epsilon)}g(y)^{p+r}|f_{g}(y)|^2dy&\lesssim \int_{Q_{d-1}(c_0)}\int_{\psi(y^{(-d)})}^{\psi(y^{(-d)})+\epsilon}|f_{g}(y^{(-d)},y^{(d)})|^{2}dy^{(d)}dy^{(-d)}\nonumber\\&\lesssim \epsilon^{2s-1}\int_{Q_{d-1}(c_0)}\int_{\psi(y^{(-d)})}^{\psi(y^{(-d)})+\epsilon}\int_{\psi(y^{(-d)})}^{y^{(d)}}|(D^sf_{g}(y^{(-d)},z))^{(d)}|^{2}dzdy^{(d)}dy^{(-d)}, \end{align}\tag{19}\] where \(Q_{d-1}(0,c_0)\) is the \(d\)-1 dimensional hypercube of side length \(c_0\) centered at \(\mathbf{0}\). Also, by changing the integration order, it yields that \[\begin{align} \label{step1:2} \int_{\psi(y^{(-d)})}^{\psi(y^{(-d)})+\epsilon}\int_{\psi(y^{(-d)})}^{y^{(d)}}|(D^sf_{g}(y^{(-d)},z))^{(d)}|^{2}dzdy^{(d)}&\lesssim \epsilon\int_{\psi(y^{(-d)})}^{\psi(y^{(-d)})+\epsilon}|(D^sf_{g}(y^{(-d)},z))^{(d)}|^{2}dz\nonumber\\&\lesssim \epsilon\int_{\psi(y^{(-d)})}^{c_0}|(D^sf_{g}(y^{(-d)},z))^{(d)}|^{2}dz. \end{align}\tag{20}\] Combining 19 and 20 , we obtain: \[\begin{align} \int_{V_{\psi}(\epsilon)}g(y)^{p+r}|f_{g}(y)|^2dy&\lesssim \epsilon^{2s}\int_{Q_{d-1}(c_0)}\int_{\psi(y^{(-d)})}^{c_0}g(y^{(-d)},z)^{q}|(D^sf_{g}(y^{(-d)},z))^{(d)}|^{2}dzdy^{(-d)}\\&\lesssim \epsilon^{2s}\|f\|_{H^{s}(U_{\psi}(c_0),g)}^{2}. \end{align}\]

Step 2: Rigid motion of local patch Now suppose at a point \(x_0\in \partial \mathcal{X}\), there exits a rigid motion \(T:\mathbb{R}^{d}\rightarrow \mathbb{R}^{d}\) such that \(T(x_0)=0\), and a number \(C_0\) such that we have all \(C_0\epsilon\le c_0\), \[\begin{align} \label{localrigid} T(Q_{T}(x_0,c_0)\cap \partial_{\epsilon} \mathcal{X})\subseteq V_{\psi}(C_0\epsilon)\quad \text{and}\quad T(Q_{T}(x_0,c_0)\cap \mathcal{X})=U_{\psi}(c_0), \end{align}\tag{21}\] where \(Q_{T}(x_0,c_0)\) is a hypercube in \(\mathbb{R}^{d}\) of side length \(c_0\) centered at \(x_0\) (not necessarily coordinate-axis-aligned). Let \(v_{g}(y):=f_{g}(T^{-1}(y))\) and \(v(y):=f(T^{-1}(y))\) for all \(y\in U_{\psi}(c_0)\). Then, if \(f_{g}\in C_{c}^{\infty}(\mathcal{X})\), we have \(v_{g}\in C_{c}^{\infty}(U_{\psi}(c_0))\) such that \(\|v_{g}\|_{H^{s}(U_{\psi}(c_0))}^{2}=\|f_g\|_{H^{s}(Q_{T}(x_0,c_0))\cap \mathcal{X}}^{2}\). Therefore, according to Step 1, we have \[\begin{align} \int_{V_{\psi}(C_0\epsilon)}g(x)^{p+r}|v_g(y)|^2dy\lesssim \epsilon^{2s}\|v\|_{H^{s}(U_{\psi}(c_0)),g}^{2}. \end{align}\] Then, it yields that \[\begin{align} &\int_{Q_{T}(x_0,c_0)\cap \partial_{\epsilon} \mathcal{X}}g^{p+r}(x)|f_g(x)|^{2}dx\\&=\int_{T(Q_{T}(x_0,c_0)\cap \partial_{\epsilon} \mathcal{X})}g^{p+r}(y)|v_g(y)|^{2}dx\\&\lesssim \int_{V_{\psi}(C_0\epsilon)}g(x)^{p+r}|v_g(y)|^2dy\\&\lesssim \epsilon^{2s}\|v\|_{H^{s}(U_{\psi}(c_0),g)}^{2}\\&\lesssim \epsilon^{2s}\|f\|_{H^{s}(Q_{T}(x_0,c_0)\cap \partial_{\epsilon} \mathcal{X},g)}^{2}\lesssim \epsilon^{2s}\|f\|_{H^{s}(\mathcal{X},g)}^{2}. \end{align}\]

Step 3: Lipschitz domain. Now we arrive at the last step where we shall deal with the case: \(\mathcal{X}\) is assumed to be an open, bounded subset of \(\mathbb{R}^{d}\) with Lipschitz boundary. Again, following the procedure in [1]. In this case, for every \(x_0\in \partial\mathcal{X}\), there exists a rigid motion \(T_{x_0}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}\) such that \(T_{x_0}(x_0)=0\), a number \(c_0(x_0)\), a Lipshitz mapping \(\psi_{x_0}:\mathbb{R}^{d-1}\rightarrow [-c_0(x_0),c_0(x_0)]\) and a number \(C_{0}(x_0)\) satisfying for all \(C_0(x_0)\epsilon\le c_0(x_0)\), 21 holds for replacing \(c_0,C_0,T,\psi\) by \(c_0(x_0),C_0(x_0),T_{x_0},\psi_{x_0}\) respectively. Therefore, by Step 2, we have \[\begin{align} \int_{Q_{T_{x_0}}(x_0,c_0(x_0))\cap \partial_{\epsilon} \mathcal{X}}g^{p+r}(x)|f_g(x)|^{2}dx\lesssim_{x_0} \epsilon^{2s}\|f\|_{H^{s}(\mathcal{X},g)}^{2}. \end{align}\] Although the constant in the last bound depends on \(x_0\), by compactness assumption, there exists a finite subset (denoted by \(x_{0,1},\ldots,x_{0,N}\)) of the collection of hypercubes \(\{Q_{T_{x_0}}(x_0,c_0(x_0)/2):x_{0}\in\partial\mathcal{X}\}\) which covers \(\partial\mathcal{X}\). Then, by taking the minimum of all constants with respect to \(x_{0,1},\ldots,x_{0,N}\), we can conclude that \[\begin{align} \partial_{\epsilon}(\mathcal{X})\subseteq \bigcup_{i=1}^{N}Q_{T_{x_{0,i}}}(x_{0,i},c_{0}(x_{0,i})). \end{align}\] Consequently, we have \[\begin{align} \int_{\partial \mathcal{X}}g^{p+r}(x)|f_{g}(x)|^{2}dx\lesssim \sum_{i=1}^{N}\int_{Q_{T_{x_{0,i}}}(x_{0,i},c_0(x_{0,i}))\cap \partial_{\epsilon} \mathcal{X}}g^{p+r}(x)|f_g(x)|^{2}dx\lesssim \epsilon^{2s}\|f\|_{H^{s}(\mathcal{X},g)}^{2}. \end{align}\] Therefore, we proved the desired result 18 . ◻

Proof of Lemma 6. Note here that we fix the assumption that \(f\in H^{1}(\mathcal{X},g;M)\) besides \(f\in H_{0}^{s}(\mathcal{X},g;M)\), which is missing in the statement of [1]. In general, it is not true for \(\mathcal{X}\neq\mathbb{R}^{d}\) that \(H_{0}^{1}(\mathcal{X},g;M)=H^{1}(\mathcal{X},g;M)\). Based on Lemma 1, we will prove Lemma 6 in a recursive way for \(s>1\). Recall that \(L_{w,\epsilon}\) is self-adjoint with respect to the weighted inner product, meaning \(\langle L_{w,\epsilon}f_1,f_2\rangle_{g^{p-r}}=\langle f_1,L_{w,\epsilon}f_2\rangle_{g^{p-r}}\) for any \(f_1/g^r,f_2/g^r\in L^{2}(\mathcal{X},g^{p+r})\). Also recall the definition of the Dirichlet energy given in Lemma 3, which can be stated as \(E_{w,\epsilon}(f,\mathcal{X}) = 2\langle L_{w,\epsilon}f,f\rangle_{g^{p-r}}\).

Following the procedure in [1] ³, when \(s=2t+1\) for \(t\ge 1\), by using self-adjointness, we have \[\begin{align} \langle L_{w,\epsilon}^{s}f,f\rangle_{g^{p-r}}=\langle L_{w,\epsilon}^{t+1}f,L_{w,\epsilon}^{t}f\rangle_{g^{p-r}}=\frac{1}{2}E_{w,\epsilon}(L_{w,\epsilon}^{t}f,\mathcal{X}). \end{align}\] We divide the Dirichlet energy into two parts: \[\begin{align} &E_{w,\epsilon}(L_{w,\epsilon}^{t}f,\mathcal{X})\\&=\frac{1}{\epsilon^{d+2}}\int_{\mathcal{X}_{t\epsilon}}\int_{\mathcal{X}_{t\epsilon}}(g(x)^{-r}f(x)-g(z)^{-r}f(z))^2\frac{\eta\left(\frac{\|x-z\|}{\epsilon}\right)}{g(x)^{1-q/2}g(z)^{1-q/2}}g(x)g(z)dxdz\\&\quad+\frac{1}{\epsilon^{d+2}}\int_{\partial_{t\epsilon} \mathcal{X}}\int_{\partial_{t\epsilon} \mathcal{X}}(g(x)^{-r}f(x)-g(z)^{-r}f(z))^2\frac{\eta\left(\frac{\|x-z\|}{\epsilon}\right)}{g(x)^{1-q/2}g(z)^{1-q/2}}g(x)g(z)dxdz\\&=:E_{w,\epsilon}(L_{w,\epsilon}^{t}f,\mathcal{X}_{t\epsilon})+E_{w,\epsilon}(L_{w,\epsilon}^{t}f,\partial_{t\epsilon} \mathcal{X}), \end{align}\] where \(\mathcal{X}_{t\epsilon}\) and \(\partial \mathcal{X}_{t\epsilon}\) have been introduced right before Lemma 4 (\(\partial \mathcal{X}_{t\epsilon} \subset \mathcal{X}\) consists of points \(t\epsilon\)-close to the boundary of \(\mathcal{X}\), and \(\mathcal{X}_{t\epsilon} = \mathcal{X} \setminus \partial \mathcal{X}_{t\epsilon}\)). By Jensen’s inequality, we have \[\begin{align} E_{w,\epsilon}(L_{w,\epsilon}^{t}f,\mathcal{X}_{t\epsilon})&\le 3\sigma_{1}^{2t}E_{w,\epsilon}\left(\sigma_{1}^{t}L_{w}^{t}f,\mathcal{X}_{t\epsilon}\right)\\&\quad+\frac{6}{\epsilon^{d+2}}\int_{\mathcal{X}_{t\epsilon}}\int_{\mathcal{X}_{t\epsilon}}\left(g(x)^{-r}L_{w,\epsilon}^{t}f(x)-g(z)^{-r}\sigma_{1}^{t}L_{w}^{t}f(z)\right)^2\\&\qquad\qquad\qquad\frac{\eta\left(\frac{\|x-z\|}{\epsilon}\right)}{g(x)^{1-q/2}g(z)^{1-q/2}}g(x)g(z)dxdz. \end{align}\] By definition 5 , we have \(\mathcal{L}_{w}^{t}f\in H^{1}(\mathcal{X},g;C\|g\|_{C^{s-1}(\mathcal{X})}^{t}M)\) for some constant \(C>0\), an application of Lemma 1 shows \(E_{w,\epsilon}\left(\sigma_{1}^{t}L_{w}^{t}f,\mathcal{X}_{t\epsilon}\right)\lesssim M^{2}\). We then focus on the second term on the right-hand side of the above inequality. According to Lemma 4, we obtain: \[\begin{align} &\frac{1}{\epsilon^{d+2}}\int_{\mathcal{X}_{t\epsilon}}\int_{\mathcal{X}_{t\epsilon}}\left(g(x)^{-r}L_{w,\epsilon}^{t}f(x)-g(z)^{-r}\sigma_{1}^{t}\mathcal{L}_{w}^{t}f(z)\right)^2\frac{\eta\left(\frac{\|x-z\|}{\epsilon}\right)}{g(x)^{1-q/2}g(z)^{1-q/2}}g(x)g(z)dxdz\\&\lesssim \frac{1}{\epsilon^{d+2}}\int_{\mathcal{X}_{t\epsilon}}\int_{\mathcal{X}_{t\epsilon}}g(x)^{p+r}\left(g(x)^{-r}L_{w,\epsilon}^{t}f(x)-g(x)^{-r}\sigma_{1}^{t}\mathcal{L}_{w}^{t}f(x)\right)^2\eta\left(\frac{\|x-z\|}{\epsilon}\right)dxdz\\&\lesssim \frac{1}{\epsilon^{2}}\int_{\mathcal{X}_{t\epsilon}}g(x)^{p-r}\left(L_{w,\epsilon}^{t}f(x)-\sigma_{1}^{t}\mathcal{L}_{w}^{t}f(x)\right)^2dx\\&\lesssim M^{2}. \end{align}\] Furthermore, near the boundary, according to Lemma 3 and Lemma 5, it yields that \[\begin{align} E_{w,\epsilon}(L_{w,\epsilon}^{t}f,\partial_{t\epsilon} \mathcal{X})\lesssim \frac{1}{\epsilon^{2}}\|L_{w,\epsilon}^{t}f\|_{L^{2}(\partial_{t\epsilon}\mathcal{X},g^{p+r})}\lesssim M^2. \end{align}\]

Putting all pieces above together, we obtain the proof for the case when \(s\) is odd and \(t:=(s-1)/2\). Similar arguments can be applied to the case when \(s\) is even and \(t:=(s-2)/2\). Therefore, combining all above together, we obtain for all integer \(s>1\): \[\begin{align} \langle L_{w,\epsilon}^{s}f,f \rangle_{g^{p-r}}\lesssim M^{2}. \end{align}\] ◻

Proof of Theorem 1. By Cauchy-Schwarz inequality, we have: for all \(s\in\mathbb{N}_{+}\): \[\begin{align} \|\hat{f}-f\|_{w,n}^2\le 2(\|\mathbb{E}\hat{f}-f|_{w,n}^2+\|\hat{f}-\mathbb{E}\hat{f}\|_{w,n}^2). \end{align}\] Then, according to PCR-WLE algorithm in Section 2.2, we obtain \[\begin{align} \label{decomp:bias} \|\mathbb{E}\hat{f}-f\|_{w,n}^{2}=\sum_{k=K+1}^{n}\langle v_k,f \rangle_{w,n}^2\le \frac{\langle L_{w,n,\epsilon}^{s}f,f\rangle_{w,n}}{\lambda_{K+1}^{s}}, \end{align}\tag{22}\] and \[\begin{align} \|\hat{f}-\mathbb{E}\hat{f}\|_{w,n}^2=\sum_{k=1}^{K}\langle v_k,\varepsilon\rangle_{w,n}^{2}. \end{align}\] Since \(\langle v_k,\varepsilon\rangle_{w,n}\) is normally distributed with \(0\) mean and variance: \[\begin{align} \label{var1} \mathop{\mathrm{Var}}\langle v_k,\varepsilon\rangle_{w,n}=\frac{1}{n^2}v_k^TD^{\frac{2(p-1-r)}{q-1}}v_k, \end{align}\tag{23}\] where \(\langle v_k,v_k\rangle_{w,n}=\frac{1}{n}v_{k}^{T}D^{\frac{p-1-r}{q-1}}v_k=1\). Note that \(\langle v_k/\sqrt{n},v_k/\sqrt{n}\rangle_{g^{p-r}}=1\), then we have \[\begin{align} \label{eigenD} \underset{v_{k}/\sqrt{n}\in\mathbb{R}^{n}}{\min}~\frac{1}{n}v_k^{T}D^{\frac{p-1-r}{q-1}}D^{\frac{p-1-r}{q-1}}v_{k} \end{align}\tag{24}\] is the smallest eigenvalue of the matrix \(D^{\frac{p-1-r}{q-1}}\) with respect to the inner product \(\langle\cdot,\cdot\rangle_{g^{p-r}}\). As \(D\) is a diagonal matrix with the \((i,i)\)-element as \(d_i\), according to Section 6.1, it is bounded from below, say by a constant \(C>0\), almost surly for \(n\) large enough. Then, combining 23 and 24 , we have: \[\begin{align} \|\hat{f}-\mathbb{E}\hat{f}\|_{w,n}^2=\frac{1}{n}\sum_{k=1}^{K}(\sqrt{n}\langle v_k,\epsilon\rangle_{w,n})^{2}, \end{align}\] with \(\sqrt{n}\langle v_k,\epsilon\rangle_{w,n}\) being normal with mean \(0\) and variance \[\begin{align} \mathop{\mathrm{Var}}(\sqrt{n} \langle v_k,\epsilon\rangle_{w,n})\ge C>0. \end{align}\] According to an exponential inequality for chi-square distributions from [54], we obtain: \[\begin{align} \label{decomp:variance} \mathbb{P}\left(\|\hat{f}-\mathbb{E}\hat{f}\|_{w,n}^2\ge \frac{CK}{n}+2\frac{\sqrt{K}}{n}\sqrt{t}+2\frac{t}{n}\right)\le e^{-t}. \end{align}\tag{25}\] With 22 and 25 , it yields \[\begin{align} \label{biasvariancedecomp} \|\hat{f}-f\|_{w,n}^{2}\le \frac{\langle L_{w,n,\epsilon}^{s}f,f\rangle_{w,n}}{\lambda_{K+1}^{s}}+\frac{CK}{n}, \end{align}\tag{26}\] with probability at least \(1-e^{-K}\) if \(1\le K\le n\). Then, it remains to bound the empirical weighted Sobolev seminorm \(\langle L_{w,n,\epsilon}^{s}f,f\rangle_{w,n}\) and the graph weighted Laplacian eigenvalue \(\lambda_{K+1}^{s}\).

We will first focus on \(\langle L_{w,n,\epsilon}^{s}f,f\rangle_{w,n}\) for \(s=1\). By definition 4 , we have by symmetry: \[\begin{align} \label{disexpress} \mathbb{E}\langle L_{w,n,\epsilon}f,f\rangle_{w,n}=\frac{1}{2}\mathbb{E}\left(\frac{1}{\epsilon^{d+2}}|d_{i}^{-\frac{r}{q-1}}f(X_i)-d_{j}^{-\frac{r}{q-1}}f(X_j)|^{2}d_{i}^{\frac{1-p}{q-1}}\frac{\eta\left(\frac{\|X_i-X_{j}\|}{\epsilon}\right)}{\tilde{d}_{i}^{1-q/2}\tilde{d}_{j}^{1-q/2}}\right). \end{align}\tag{27}\] We would like to point out here that the normalization factor \(\epsilon^{-(d+2)}\) is motivated by the fact that a factor of \(\epsilon^{-d}\) is needed to scale \(\eta\left(\frac{\|X_i-X_j\|}{\epsilon}\right)\) and the remaining factor, \(\epsilon^{-2},\) stabilized the squared differences of \(d_{i}^{-\frac{r}{q-1}}\) under the expectation.

According to Section 6.1 and by conditioning on \(X_i\) and the law of iterated expectation, we have for \(n\) large enough, \[\begin{align} \label{distononlocal} {!}{\left|\mathbb{E}\left(\frac{1}{\epsilon^{d+2}}|d_{i}^{-\frac{r}{q-1}}f(X_i)-d_{j}^{-\frac{r}{q-1}}f(X_j)|^{2}d_{i}^{\frac{1-p}{q-1}}\frac{\eta\left(\frac{\|X_i-X_{j}\|}{\epsilon}\right)}{\tilde{d}_{i}^{1-q/2}\tilde{d}_{j}^{1-q/2}}\right)-2\langle L_{w,\epsilon}f,f\rangle_{g^{p-r}}\right| \lesssim \Delta(n,\epsilon,\eta,g)+\epsilon,} \end{align}\tag{28}\] where \[\begin{align} \Delta(n,\epsilon,\eta,g): =\frac{1}{n}g_{\max}+\frac{\eta(0)}{n\epsilon^{d}}+\frac{n-1}{n}\Big(\sqrt{\frac{|\log \epsilon|}{n\epsilon^{d}}}+\epsilon\Big) \to 0\quad\text{as }\;n\to \infty. \end{align}\]

Combining 27 , 28 and Lemma 1, we obtain: \[\begin{align} \mathbb{E}\langle L_{w,n,\epsilon}f,f\rangle_{w,n}\lesssim M^{2}+\Delta(n,\epsilon,\eta,g)+\epsilon. \end{align}\] Consequently, by Markov’s inequality, we have: for any \(\delta\in(0,1)\), \[\begin{align} \label{markovs611} \langle L_{w,n,\epsilon}f,f\rangle_{w,n}\lesssim \frac{1}{\delta}\left(M^{2}+\Delta(n,\epsilon,\eta,g)+\epsilon\right), \end{align}\tag{29}\] with probability at least \(1-\delta\). Note that the above bound on the expected weighted Sobolev seminorm generalizes the results in [1] to the weighted Laplacians by some properties of KDE.

Next, we proceed to the higher order case when \(s>1\) for \(\langle L_{w,n,\epsilon}^{s}f,f\rangle_{w,n}\). We define the following difference operator: \[\begin{align} D_{j}f(x)=(d_{\cdot}^{-\frac{r}{q-1}}f(x)-d_{j}^{-\frac{r}{q-1}}f(X_j))d_{\cdot}^{\frac{1-p}{q-1}}w_{\cdot,j}^{\epsilon}, \end{align}\] where \(d_{\cdot}\) and \(w_{\cdot,j}^{\epsilon}\) are defined by replacing \(X_i\) by \(x\) in both \(d_{i}\) and \(w_{i,j}^{\epsilon}.\) Furthermore, let \(D_{\mathbf{j}}f(x):=(D_{j_1}f\circ\ldots\circ D_{j_s}f)(x)\), where \(\boldsymbol{j} = (j_1,\ldots,j_s) \in [n]^s := \{1,\ldots,n\}^s\). Denote by \((n)^{s}\) the sub-collection of vectors in \([n]^s\) with no repeated indices and let by \(i\mathbf{j}\mathrel{\vcenter{:}}= (i,j_1,\ldots,j_{s}).\)

Following the idea of [1], we decompose the weighted Sobolev seminorm into a U-statistic, which is an unbiased estimator of the non-local Sobolev seminorm \(\langle L_{w,\epsilon}^{s}f,f\rangle_{g^{p-r}}\), and a pure bias term: \[\begin{align} \label{I1I2} \langle L_{w,n,\epsilon}^{s}f,f\rangle_{w,n}&=\frac{1}{n}\sum_{i=1}^{n}d_{i}^{\frac{p-1-r}{q-1}}L_{w,n,\epsilon}^{s}f(X_i)\cdot f(X_i)\nonumber\\&=\frac{1}{n\epsilon^{2s}}\sum_{i\mathbf{j}\in (n)^{s+1}}d_{i}^{\frac{p-1-r}{q-1}}D_{\mathbf{j}}f(X_i)\cdot f(X_i)\nonumber\\&\qquad+\frac{1}{n\epsilon^{2s}}\sum_{i\mathbf{j}\in [n]^{s+1}\backslash(n)^{s+1}}d_{i}^{\frac{p-1-r}{q-1}}D_{\mathbf{j}}f(X_i)\cdot f(X_i)\nonumber\\&=:I_1+I_2. \end{align}\tag{30}\]

Note that there are errors in [1] when bounding both \(\mathbb{E}I_{1}\) and \(\mathbb{E}I_{2}\). Specifically, in [1], there should not be a \(\delta\) appearing in Equation D.4 by Markov’s inequality and the power of \(\epsilon\) should be \(2s+d\). Although their final result is correct, we will fix these errors in the following proof. Now, determined by whether all \(i\mathbf{j}\) are distinct, the empirical weighted Sobolev seminorm can be divided into two parts, \(I_1\) and \(I_2\). The first one involves all distinct indices where we make approximation by the so-called non-local weighted sobolev norm \(\langle L_{w,\epsilon}^{s}f,f\rangle_{g^{p-r}}\); the second part focuses on the case where not all \(i\mathbf{j}\) are distinct and use the fact that it is related to a connected subgraph.

As for \(I_1\) from 30 , we have \[\begin{align} \mathbb{E}I_1&=\frac{1}{n\epsilon^{2s}}\frac{n!}{(n-s-1)!}\mathbb{E}\left(d_{i}^{\frac{p-1-r}{q-1}}D_{\mathbf{j}}f(X_i)\cdot f(X_i)\right)\\&=\frac{1}{n\epsilon^{2s}}\frac{n!}{(n-s-1)!}\mathbb{E}\langle D_{\mathbf{j}}f(X_i),f(X_i)\rangle_{g^{p-r}}, \end{align}\] where the operator \(D_{\mathbf{j}}\) is iterated for \(s\) different times due to the fact that \(i\mathbf{j}\) are all distinct. For each iteration, say \(s=1\), we have \[\begin{align} \label{alldiffers1} &\quad\mathbb{E}\langle D_{j}f(X_i),f(X_i)\rangle_{g^{p-r}}\nonumber\\&=\frac{\epsilon^{2}}{2n}\mathbb{E}\left(\frac{1}{\epsilon^{d+2}}|d_{i}^{-\frac{r}{q-1}}f(X_i)-d_{j}^{-\frac{r}{q-1}}f(X_j)|^{2}d_{i}^{\frac{1-p}{q-1}}\frac{\eta\left(\frac{\|X_i-X_{j}\|}{\epsilon}\right)}{\tilde{d}_{i}^{1-q/2}\tilde{d}_{j}^{1-q/2}}\right). \end{align}\tag{31}\] Then, plugging 28 in 31 , we obtain \[\begin{align} \left|\mathbb{E}\langle D_{j}f(X_i),f(X_i)\rangle_{g^{p-r}}-\frac{\epsilon^{2}}{n}\langle L_{w,\epsilon}f,f\rangle_{g^{p-r}}\right|\lesssim \frac{\epsilon^{2}}{2n}(\Delta(n,\epsilon,\eta,g)+\epsilon). \end{align}\] After \(s\) times iteration, it yields that \[\begin{align} \left|\mathbb{E}\langle D_{\mathbf{j}}f(X_i),f(X_i)\rangle_{g^{p-r}}-\frac{\epsilon^{2s}}{n^{s}}\langle L_{w,\epsilon}^{s}f,f\rangle_{g^{p-r}}\right|\lesssim \frac{\epsilon^{2s}}{2^{s}n^{s}}(\Delta(n,\epsilon,\eta,g)+\epsilon). \end{align}\] Putting all above results back in \(\mathbb{E}I_1\), we conclude that for \(n\) large enough, \[\begin{align} \label{beforestirling} \left|\mathbb{E}I_1-\frac{n!}{n^{s+1}(n-s-1)!}\langle L_{w,\epsilon}^{s}f,f \rangle_{g^{p-r}}\right|\lesssim \frac{n!}{n^{s+1}(n-s-1)!}(\Delta(n,\epsilon,\eta,g)+\epsilon). \end{align}\tag{32}\] The Stirling’s formula shows \[\begin{align} \underset{n\rightarrow\infty}{\lim}~\frac{n!}{n^{s+1}(n-s-1)!}=1. \end{align}\] Therefore, by 32 , we have for \(n\) large enough, \[\begin{align} \mathbb{E}I_{1}\lesssim \langle L_{w,\epsilon}^{s}f,f \rangle_{g^{p-r}}+(\Delta(n,\epsilon,\eta,g)+\epsilon). \end{align}\] According to Lemma 6, it yields that \[\begin{align} \label{ExpectI1} \mathbb{E}I_{1}\lesssim M^{2}+(\Delta(n,\epsilon,\eta,g)+\epsilon). \end{align}\tag{33}\]

We next shift our attention to \(I_{2}\) in 30 : \[\begin{align} \frac{1}{n\epsilon^{2s}}\sum_{i\mathbf{j}\in [n]^{s+1}\backslash(n)^{s+1}}d_{i}^{\frac{p-1-r}{q-1}}D_{\mathbf{j}}f(X_i)\cdot (f(X_i)-f(X_{j_1})). \end{align}\]

For \(i\mathbf{j}\) not all distinctive, if they contains a total of \((k+1)\) distinct indices for example for \(1\le k\le s-1\), we have by symmetry: \[\begin{align} \sum_{i\mathbf{j}\in [n]^{s+1}\backslash(n)^{s+1}}d_{i}^{\frac{p-1-r}{q-1}}D_{\mathbf{j}}f(X_i)\cdot f(X_i)=\frac{1}{2}\cdot \sum_{i\mathbf{j}\in [n]^{s+1}\backslash(n)^{s+1}}d_{i}^{\frac{p-1-r}{q-1}}D_{\mathbf{j}}f(X_i)\cdot (f(X_i)-f(X_{j_1})). \end{align}\] Observe that in order for \[\begin{align} d_{i}^{\frac{p-1-r}{q-1}}|D_{\mathbf{j}}f(X_i)|\cdot |f(X_i)-f(X_{j_1})| \end{align}\] to be non-zero, it must be the case that the graph \(G_{n,\epsilon}(X_{i\mathbf{j}})\) which is the subgraph induced by the vertices \(X_i,X_{j_1},\ldots,X_{j_{s}}\) is complete. Since we have: \[\begin{align} D_{ij}f(x)&=D_{i}(D_{j}f(x))\\&=D_{i}\left((d_{\cdot}^{-\frac{r}{q-1}}f(x)-d_{j}^{-\frac{r}{q-1}}f(X_j))d_{\cdot}^{\frac{1-p}{q-1}}w_{\cdot,j}^{\epsilon}\right)\\&=(d_{\cdot}^{-\frac{r}{q-1}}D_j f(x)-d_{i}^{-\frac{r}{q-1}}D_{j}f(X_i))d_{\cdot}^{\frac{1-p}{q-1}}w_{\cdot,i}^{\epsilon}, \end{align}\] then \[\begin{align} |D_{j_{1}j_{2}}f(X_{i})|\le \left(d_{i}^{-\frac{r}{q-1}}|D_{j_2}f(X_{i})|+d_{j_1}^{-\frac{r}{q-1}}|D_{j_2}f(X_{j_1})|\right)d_{i}^{\frac{1-p}{q-1}}w_{i,j_1}^{\epsilon}. \end{align}\] Repeating the above computation and by induction, it yields that for \(s\ge 2\), \[\begin{align} |D_{\mathbf{j}}f(X_i)|\le (s-1)d_{\max/\min}^{-\frac{(s-1)r}{q-1}}d_{\max/\min}^{\frac{(s-1)(1-p)}{q-1}}\left(w_{\max}^{\epsilon}\right)^{s-1}\sum_{j\in i\mathbf{j}\backslash\{j_{s}\}}|D_{j_{s}}f(X_{j})|, \end{align}\] where \(d_{\max}:=\underset{i=1,\ldots,n}{\max}~d_{i}\), \(d_{\min}:=\underset{i=1,\ldots,n}{\min}~d_{i}\), \(w_{\max}:=\underset{i,j=1,\ldots,n}{\max}~w_{i,j}\) and \(d_{\max/\min}\) means it is \(d_{\max}\) if \(-(s-1)r/(d-1)\) (respectively \((s-1)(1-p)/(q-1)\)) are positive and it is \(d_{\min}\) otherwise.

According to Section 6.1, we have for \(n\) large enough, \(d_{\max}\) is bounded from above and \(d_{\min}\) is bounded from below a.s. and \[\begin{align} w_{\max}\lesssim \frac{1}{n\epsilon^{d}}, \end{align}\] almost surely.

Consequently, it yields that \[\begin{align} \label{iteration} &d_{i}^{\frac{p-1-r}{q-1}}|D_{\mathbf{j}}f(X_i)|\cdot |f(X_i)-f(X_{j_1})|\nonumber\\&=d_{i}^{\frac{p-1-r}{q-1}}|D_{\mathbf{j}}f(X_i)|\cdot |f(X_i)-f(X_{j_1})|\cdot\mathbf{1}_{\{G_{n,\epsilon}(X_{i\mathbf{j}})\;\text{is connected}\}}\nonumber\\&\lesssim \frac{1}{(n\epsilon^{d})^{s-1}}\sum_{j\in i\mathbf{j}\backslash\{j_s\}}\left(d_{i}^{\frac{p-1-r}{q-1}}|D_{j_{s}}f(X_{j})|\cdot|f(X_i)-f(X_{j_1})|\cdot\mathbf{1}_{\{G_{n,\epsilon}(X_{i\mathbf{j}})\;\text{is connected}\}}\right)\nonumber\\&=\frac{\epsilon^{2}}{n^{s}\epsilon^{d(s-1)}}\sum_{j\in i\mathbf{j}\backslash\{j_s\}}\Bigg(\frac{1}{\epsilon^{d+2}}d_{i}^{\frac{p-1-r}{q-1}}|d_{j}^{-\frac{r}{q-1}}f(X_j)-d_{j_s}^{-\frac{r}{q-1}}f(X_{j_{s}})|d_{j}^{\frac{1-p}{q-1}}\frac{\eta\left(\frac{\|X_j-X_{j_s}\|}{\epsilon}\right)}{\tilde{d}_{j}^{1-q/2}\tilde{d}_{j_{s}}^{1-q/2}}\nonumber\\&\qquad|f(X_i)-f(X_{j_1})| \mathbf{1}_{\{G_{n,\epsilon}(X_{i\mathbf{j}})\;\text{is connected}\}}\Bigg), \end{align}\tag{34}\] where we again assign \(\epsilon^{d+2}\) as a normalization factor into the expectation as 27 .

Now, note that for \(j=i\) in the summand on the right-hand side of 34 , we have according to Section 6.1: \[\begin{align} \label{j61i} &\mathbb{E}\Bigg(\frac{1}{\epsilon^{d+2}}d_{i}^{-\frac{r}{q-1}}|d_{i}^{-\frac{r}{q-1}}f(X_i)-d_{j_s}^{-\frac{r}{q-1}}f(X_{j_{s}})|\frac{\eta\left(\frac{\|X_i-X_{j_s}\|}{\epsilon}\right)}{\tilde{d}_{i}^{1-q/2}\tilde{d}_{j_{s}}^{1-q/2}}|f(X_i)-f(X_{j_1})|\nonumber\\&\qquad \mathbf{1}_{\{G_{n,\epsilon}(X_{i\mathbf{j}})\;\text{is connected}\}}\Bigg)\nonumber\\&\lesssim\mathbb{E}\left(\left(\frac{1}{\epsilon^{d+2}}|d_{i}^{-\frac{r}{q-1}}f(X_i)-d_{j_s}^{-\frac{r}{q-1}}f(X_{j_{s}})|\frac{\eta\left(\frac{\|X_i-X_{j_s}\|}{\epsilon}\right)}{\tilde{d}_{i}^{1-q/2}\tilde{d}_{j_{s}}^{1-q/2}}|d_{i}^{-\frac{r}{q-1}}f(X_i)-d_{j_{1}}^{-\frac{r}{q-1}}f(X_{j_{1}})|\right.\right.\nonumber\\&\qquad+(\Delta(n,\epsilon,\eta,g)+\epsilon)\Bigg)\mathbf{1}_{\{G_{n,\epsilon}(X_{i\mathbf{j}})\;\text{is connected}\}}\Bigg)\nonumber\\&\lesssim \mathbb{E}\left(\left(\frac{1}{\epsilon^{d+2}}|d_{i}^{-\frac{r}{q-1}}f(X_i)-d_{j_s}^{-\frac{r}{q-1}}f(X_{j_{s}})|^{2}\frac{\eta\left(\frac{\|X_i-X_{j_s}\|}{\epsilon}\right)}{\tilde{d}_{i}^{1-q/2}\tilde{d}_{j_{s}}^{1-q/2}}+(\Delta(n,\epsilon,\eta,g)+\epsilon)\right)\right.\nonumber\\&\qquad\mathbf{1}_{\{G_{n,\epsilon}(X_{i\mathbf{j}})\;\text{is connected}\}}\Bigg), \end{align}\tag{35}\] where the last inequality is by Cauchy–Schwarz inequality and \(X_1,\ldots,X_{n}\) being i.i.d. data. Then, by integrating out all indices in \(\mathbf{j}\) not equal to \(i\) or \(j_{s}\), it yields that \[\begin{align} \label{k431} &\mathbb{E}\left(\left(\frac{1}{\epsilon^{d+2}}|d_{i}^{-\frac{r}{q-1}}f(X_i)-d_{j_s}^{-\frac{r}{q-1}}f(X_{j_{s}})|^{2}\frac{\eta\left(\frac{\|X_i-X_{j_s}\|}{\epsilon}\right)}{\tilde{d}_{i}^{1-q/2}\tilde{d}_{j_{s}}^{1-q/2}}+(\Delta(n,\epsilon,\eta,g)+\epsilon)\right)\right.\nonumber\\&\qquad\mathbf{1}_{\{G_{n,\epsilon}(X_{i\mathbf{j}})\;\text{is connected}\}}\Bigg)\nonumber\\&\lesssim \left(C\epsilon^{d} g_{\max}V_{d}\right)^{k-1}\mathbb{E}\left(\left(\frac{1}{\epsilon^{d+2}}|d_{i}^{-\frac{r}{q-1}}f(X_i)-d_{j_s}^{-\frac{r}{q-1}}f(X_{j_{s}})|^{2}\frac{\eta\left(\frac{\|X_i-X_{j_s}\|}{\epsilon}\right)}{\tilde{d}_{i}^{1-q/2}\tilde{d}_{j_{s}}^{1-q/2}}\right.\right.\nonumber\\&\qquad+(\Delta(n,\epsilon,\eta,g)+\epsilon)\Bigg)\Bigg). \end{align}\tag{36}\]

Therefore, according to 35 , 36 , 28 and Lemma 1, we obtain \[\begin{align} \label{jifinal} &\mathbb{E}\Bigg(\frac{1}{\epsilon^{d+2}}d_{i}^{-\frac{r}{q-1}}|d_{i}^{-\frac{r}{q-1}}f(X_i)-d_{j_s}^{-\frac{r}{q-1}}f(X_{j_{s}})|\frac{\eta\left(\frac{\|X_i-X_{j_s}\|}{\epsilon}\right)}{\tilde{d}_{i}^{1-q/2}\tilde{d}_{j_{s}}^{1-q/2}}|f(X_i)-f(X_{j_1})|\nonumber\\&\qquad \mathbf{1}_{\{G_{n,\epsilon}(X_{i\mathbf{j}})\;\text{is connected}\}}\Bigg)\nonumber\\&\lesssim \epsilon^{d(k-1)}\left(M^{2}+\Delta(n,\epsilon,\eta,g)+\epsilon\right). \end{align}\tag{37}\] Applying a similar approach to all \(j\neq j_{s}\) and plugging 37 in 34 and 30 , we have \[\begin{align} \mathbb{E}I_{2}&\lesssim \frac{1}{n\epsilon^{2s}}\frac{1}{n^{s}\epsilon^{d(s-1)}}\sum_{k=1}^{s-1}\epsilon^{d(k-1)}\left(M^{2}+\Delta(n,\epsilon,\eta,g)\right)n^{k431}\\&\lesssim \frac{\epsilon^{2}}{n\epsilon^{2s}}\left(M^{2}+\Delta(n,\epsilon,\eta,g)+\epsilon\right)\sum_{k=1}^{s-1}\frac{(n\epsilon^{d})^{k}}{(n\epsilon^{d})^{s}}n. \end{align}\] Note that the above sum is bounded from above when \(k=s-1\) by the assumption \(n\epsilon^{d}\ge 1\). Finally, we conclude that \[\begin{align} \label{ExpectI2} \mathbb{E}I_{2}\lesssim \frac{\epsilon^{2}}{n\epsilon^{2s+d}}\left(M^{2}+\Delta(n,\epsilon,\eta,g)+\epsilon\right). \end{align}\tag{38}\]

Finally, combining 30 , 33 and 38 , we obtain: \[\begin{align} \mathbb{E}\langle L_{w,n,\epsilon}^{s}f,f\rangle_{w,n}&\lesssim M^{2}+(\Delta(n,\epsilon,\eta,g)+\epsilon)+\frac{\epsilon^{2}}{n\epsilon^{2s+d}}\left(M^{2}+\Delta(n,\epsilon,\eta,g)+\epsilon\right)\\&\lesssim M^{2}+(\Delta(n,\epsilon,\eta,g)+\epsilon), \end{align}\] where the last step is by the assumption that \(\epsilon\gtrsim n^{-1/(2(s-1)+d)}\). By Markov’s inequality, we have for any \(\delta\in (0,1)\), \[\begin{align} \label{markovs621} \langle L_{w,n,\epsilon}^{s}f,f\rangle_{w,n}\lesssim \frac{1}{\delta}\left(M^{2}+(\Delta(n,\epsilon,\eta,g)+\epsilon\right), \end{align}\tag{39}\] with probability at least \(1-2\delta\). This bound can be considered as a higher order variant of 29 for \(s>1\).

Now, recall the bound 26 . We have bounded the empirical weighted Sobolev seminorm by 29 and 39 . It remains to bound the eigenvalues \(\lambda_{K+1}\).

According to Lemma 7, we have: \[\begin{align} \label{eigenboundreadytouse} \lambda_{k}=\lambda_{k}(L_{w,n,\epsilon})\gtrsim \lambda_{k}(\mathcal{L}_{w}) \wedge \epsilon^{2},\;\text{for all}\;2\le k\le n, \end{align}\tag{40}\] with probability at least \(1-Cne^{-cn\epsilon^{d}}\) for some constants \(C,c>0\).

For \(s=1\), combining 26 , 29 and 40 , we have with probability at least \(1-\delta-Cne^{-cn\epsilon^{d}}-e^{-K}\) and \(n\) large enough: \[\begin{align} \|\hat{f}-f\|_{w,n}^{2}\lesssim \frac{M^{2}}{\delta\left(\lambda_{K+1}(\mathcal{L}_{w}) \wedge \epsilon^{2}\right)}+\frac{K}{n}. \end{align}\] Furthermore, based on the assumption \(\epsilon\lesssim K^{-1/d}\) and Proposition 6, the above inequality becomes: \[\begin{align} \label{tradeoffres} \|\hat{f}-f\|_{w,n}^{2}\lesssim \frac{M^{2}}{\delta}(K+1)^{-2/d}+\frac{K}{n}. \end{align}\tag{41}\] By balancing the two terms on the right-hand side, we pick \(K=\lfloor M^2n \rfloor^{d/(2+d)}\). Then, it yields that \[\begin{align} \label{finallyminimaxopt} \|\hat{f}-f\|_{w,n}^{2}\lesssim \frac{1}{\delta}M^2(M^2n)^{-2/(2+d)}. \end{align}\tag{42}\] If \(M^2<n^{-1}\), we can take \(K=1\) and obtain from 41 that: \[\begin{align} \|\hat{f}-f\|_{w,n}^{2}\lesssim \frac{1}{n\delta}. \end{align}\] If \(M>n^{1/d}\), we take \(K=n\) and in this case, we actually have \(\hat{f}(X_i)=Y_i\) for \(i=1,\ldots,n\) and \[\begin{align} \|\hat{f}-f\|_{w,n}^{2}=\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}^{2}\lesssim C, \end{align}\] with probability at least \(1-e^{-n}\) for some constant \(C\). Combining all above cases depending on choices of \(K\), it yields that bound in Theorem 1 .

For \(s>1\), the proof follows in a similar way by considering 39 instead of 29 . ◻

Proof of Theorem 2. Recall the construction of the estimator based on Lepski’s procedure: \(\hat{f}_{\mathsf{adapt}}=\hat{f}_{\hat{s},\hat{M}}\) with \(\hat{s},\hat{M}\) given in Section 3.3. Let the event \(\mathcal{E}_{j}\) be that \(\hat{s}=s_{j}\) and suppose \(s=s_{i}\) for the true smooth parameter.

First of all, it suffices to consider \(M\in \mathcal{D}\) by realizing that if \(M\in (M_{j-1},M_{j})\), then \(f\in H^{s}(\mathcal{X},g;M)\) with \(H^{s}(\mathcal{X},g;M_{j-1})\subset H^{s}(\mathcal{X},g;M)\subset H^{s}(\mathcal{X},g;M_{j})\). Now, we also suppose \(M=M_i\) correspondingly and consider bounding the sum: \[\begin{align} \sum_{j=1}^{N_l}\left(\|\hat{f}_{s_j}-f\|_{w,n}^{2}M_i^{-2}(M_i^2 n/\log n )^{2s_i/(2s_i+d)}\mathbf{1}_{\mathcal{E}_{j}}\right), \end{align}\] conditional on the event that the sample points \(X_1,\ldots,X_n\) satisfy 26 and 40 with \(K=\lfloor M_i^2n \rfloor^{d/(2s_i+d)}\). These two statements hold with probability at least \(1-Cne^{-Cn\epsilon^{d}}-e^{-\lfloor M_i^2n \rfloor^{d/(2s_i+d)}}\). As we will see, the fact that this sum does not explode, relies on the fact that the probabilities of the sets \({\mathcal{E}}_j\) get small as \(n \to \infty\).

First, note that by Cauchy-Schwarz inequality, we have \[\begin{align} &\sum_{j61i}^{N_l}\left(\|\hat{f}_{s_j}-f\|_{w,n}^{2}M_i^{-2}(M_i^2 n/\log n )^{2s_i/(2s_i+d)}\mathbf{1}_{\mathcal{E}_{j}}\right)\\&\le \sum_{j61i}^{N_l}\left(\|\hat{f}_{s_j}-\hat{f}_{s_i}+\hat{f}_{s_j}-f\|_{w,n}^{2}M_i^{-2}(M_i^2 n/\log n )^{2s_i/(2s_i+d)}\mathbf{1}_{\mathcal{E}_{j}}\right)\\&\le \sum_{j61i}^{N_l}\left(2c_0^2 \mathbf{1}_{\mathcal{E}_{j}}+2\left(\|\hat{f}_{s_i}-f\|_{w,n}^{2}M_i^{-2}(M_i^2 n/\log n)^{2s_i/(2s_i+d)}\mathbf{1}_{\mathcal{E}_{j}}\right)\right)\\&\le 2c_0^2+2\left(\|\hat{f}_{s_i}-f\|_{w,n}^{2}M_i^{-2}(M_i^2 n/\log n)^{2s_i/(2s_i+d)}\right). \end{align}\] Therefore, according to Theorem 1, we have: for any \(\delta\in(0,1)\), \[\begin{align} \sum_{j61i}^{N_l}\left(\|\hat{f}_{s_j}-f\|_{w,n}^{2}M_i^{-2}(M_i^2 n/\log n )^{2s_i/(2s_i+d)}\mathbf{1}_{\mathcal{E}_{j}}\right)\lesssim \frac{1}{\delta}, \end{align}\] with probability at least \(1-\delta \log^{-2s_i/(2s_i+d)}n-Cne^{-Cn\epsilon^{d}}-e^{-\lfloor M_i^2n \rfloor^{d/(2s_i+d)}}\).

Next, we consider the other part when \(j<i\): \[\begin{align} \label{sumj60i1} \sum_{j=1}^{i-1}\left(\|\hat{f}_{s_j}-f\|_{w,n}^{2}M_i^{-2}(M_i^2 n/\log n )^{2s_i/(2s_i+d)}\mathbf{1}_{\mathcal{E}_{j}}\right). \end{align}\tag{43}\]

By the definition, on the event \(\mathcal{E}_{j}\), there exists \(s'\in \mathcal{B}\) with \(s'<s_i\) such that \(\|\hat{f}_{s_i}-\hat{f}_{s'}\|_{w,n}>c_0 M'^{-2}(M'^2 n/\log n)^{-s'/(2s'+d)}\). This means \(\|\hat{f}_{s_i}-\hat{f}_{s'}\|_{w,n}^2 M'^{-2}(M'^2 n/\log n)^{2s'/(2s'+d)}>c_0^{2}\). By triangle inequality, this implies we have either \(\|\hat{f}_{s_i}-f\|_{w,n}^2 M'^{-2}(M'^2 n/\log n)^{2s'/(2s'+d)}>c_0^2/4\) or \(\|\hat{f}_{s'}-f\|_{w,n}^2 M'^{-2}(M'^2 n/\log n)^{2s'/(2s'+d)}>c_0^2/4\). Then, we have \[\begin{align} \label{PmalE} \mathbb{P}(\mathcal{E}_{j})&\le \sum_{l=1}^{i-1}\left(\mathbb{P}\left(\|\hat{f}_{s_i}-f\|_{w,n}^2 M_l^{-2}(M_l^2 n/\log n)^{2s_l/(2s_l+d)}>c_0^2/4\right)\right.\nonumber\\&\quad\left.+\mathbb{P}\left(\|\hat{f}_{s_l}-f\|_{w,n}^2 M_l^{-2}(M_l^2 n/\log n)^{2s_l/(2s_l+d)}>c_0^2/4\right)\right). \end{align}\tag{44}\] Since \(l<i\), we have \(f\in H^{s_i}(\mathcal{X},g;M_l)\subset H^{s_{l}}(\mathcal{X},g;M_l)\) for all \(l<i\). Therefore, it suffices to focus on the concentration inequality of \(\hat{f}_{s_l}\) to \(f\), i.e., bounding \[\begin{align} \label{j601part1} \mathbb{P}\left(\|\hat{f}_{s_l}-f\|_{w,n}^2 M_l^{-2}(M_l^2 n/\log n)^{2s_l/(2s_l+d)}>c_0^2/4\right). \end{align}\tag{45}\]

Note that the key problem here is the rate of convergence of \(\|\hat{f}_{s_j}-f\|_{w,n}^{2}\) in 43 does not match the rate \((n/\log n)^{2s_i/(2s_i+d)}\) given there. However, this can be dealt with by controlling the probability of the event \(\mathcal{E}_{j}\). The strategy here is we need a better concentration inequality than what has been proven previously as 39 otherwise the probability of the event \(\mathcal{E}_{j}\) will not decay to \(0\). Observe that the concentration 39 : for \(n\) large enough and with probability smaller than \(1-2\delta\), \[\begin{align} \langle L_{w,n,\epsilon}^{s}f,f\rangle_{w,n}\lesssim \delta^{-1} M^2, \end{align}\] is from the application of Markov’s inequality with \[\begin{align} \mathbb{E}\langle L_{w,n,\epsilon}^{s}f,f\rangle_{w,n}\lesssim M^2, \end{align}\] for \(n\) large enough. While bounding the first moment gives a concentration inequality with probability \(1-2\delta\), establishing a higher moment bound, e.g. the second moment, would result in a better concentration inequality with higher probability similar to [19], which fits in our proof technique.

Starting with \(s=1\) and similar to 28 , we have: for \(n\) large enough, \[\begin{align} \label{varianceconcs611} &\mathop{\mathrm{Var}}\langle L_{w,n,\epsilon}f,f \rangle_{w,n}\nonumber\\\lesssim& \mathop{\mathrm{Var}}\left(\frac{1}{2}\frac{1}{n^2\epsilon^{d+2}}\sum_{i,j=1}^{n}(g(X_i)^{-r}f(X_i)-g(X_j)^{-r}f(X_j))^2g(X_{i})^{1-p}\frac{\eta\left(\frac{\|X_i-X_j\|}{\epsilon}\right)}{g(X_i)^{1-q/2}g(X_j)^{1-q/2}}\right). \end{align}\tag{46}\] For \(i,j\in 1,\ldots,n\), let \[V_{ij}:=(g(X_i)^{-r}f(X_i)-g(X_j)^{-r}f(X_j))^2g(X_{i})^{1-p}\frac{\eta\left(\frac{\|X_i-X_j\|}{\epsilon}\right)}{g(X_i)^{1-q/2}g(X_j)^{1-q/2}}.\]

We have: \[\begin{align} &\mathop{\mathrm{Var}}\left(\sum_{i,j=1}^{n}(g(X_i)^{-r}f(X_i)-g(X_j)^{-r}f(X_j))^2g(X_{i})^{1-p}\frac{\eta\left(\frac{\|X_i-X_j\|}{\epsilon}\right)}{g(X_i)^{1-q/2}g(X_j)^{1-q/2}}\right)\\&=\sum_{i,j=1}^{n}\sum_{l,m=1}^{n}\text{Cov}(V_{ij},V_{lm}). \end{align}\] Now, consider the following four scenarios depending on the cardinality of \(\{i,j,l,m\}\).

• If \(|\{i,j,l,m\}|=4\), since \(V_{ij}\) and \(V_{lm}\) are independent, we have \(\text{Cov}(V_{ij},V_{lm})=0\).

• If \(|\{i,j,l,m\}|=3\), without loss of generality, say \(i=l\), we have by Lipschitz condition, \[\begin{align} \text{Cov}(V_{ij},V_{im})&\le \mathbb{E}[V_{ij}V_{im}]\\&\lesssim \epsilon^{2d+4}M^{4}. \end{align}\]

• If \(|\{i,j,l,m\}|=2\), without loss of generality, say \(i=l\) and \(j=m\), similarly, we obtain \[\begin{align} \text{Cov}(V_{ij},V_{ij})&\le \mathbb{E}V_{ij}^2\\&\lesssim \epsilon^{d+4}M^{4}. \end{align}\]

• If \(|\{i,j,l,m\}|=1\), we have \(V_{ij}=V_{lm}=0\).

Plugging the above results in 46 , it yields that for \(n\) large enough, \[\begin{align} \mathop{\mathrm{Var}}\langle L_{w,n,\epsilon}f,f \rangle_{w,n}\lesssim \frac{1}{4n^{4}\epsilon^{2d+4}}\left(n^{3}\epsilon^{2d+4}M^{4}+n^2\epsilon^{d+4}M^{4}\right)\lesssim n^{-1}M^4, \end{align}\] where the last step follows by the assumption that \(n\epsilon^{d}\ge 1\). Then, by Markov’s inequality, we obtain: for any \(\delta\in (0,1)\), \[\begin{align} \label{markovvariances611} \mathbb{P}\left(\left|\langle L_{w,n,\epsilon}f,f \rangle_{w,n}-\mathbb{E} \langle L_{w,n,\epsilon}f,f \rangle_{w,n}\right|\ge \frac{1}{\delta} M^{2}\right)\lesssim \frac{\delta^{2}}{n}. \end{align}\tag{47}\] Combining 47 and 29 , we conclude that for \(n\) large enough, \[\begin{align} \langle L_{w,n,\epsilon}f,f \rangle_{w,n}\lesssim \frac{1}{\delta}M^2 \end{align}\] holds with probability not less than \(1-\frac{\delta}{n^2}\). Furthermore, following a similar argument in Lemma 6, one can show the above high-probability bound also holds for the case \(s>1\). Thus, under the additional Lipschitz assumption that \(|f_{g}(x)-f_{g}(x')|\le M\|x-x'\|\), we establish a better bound for the empirical weighted Sobolev seminorm: for all \(s\in\mathbb{N}_{+}\) and \(n\) large enough, \[\begin{align} \langle L_{w,n,\epsilon}^s f,f \rangle_{w,n}\lesssim \frac{1}{\delta}M^2, \end{align}\] with probability at least \(1-C\frac{\delta^2}{n}\).

Conditional on the event that the sample points \(X_1,\ldots,X_n\) satisfy 26 and 40 with \(K=\lfloor M^2n \rfloor^{d/(2s+d)}\), following the proof of Theorem 1 to obtain 42 by using the better concentration inequality we derived above instead, we have for \(n\) large enough, \[\begin{align} \|\hat{f}-f\|_{w,n}^{2}\lesssim \frac{1}{\delta}M^2(M^2n)^{-2s/(2s+d)}, \end{align}\] with probability at least \(1-C\delta^2 n^{-1}-Cne^{-Cn\epsilon^{d}}-e^{-\lfloor M^2n \rfloor^{d/(2s+d)}}\) under the minimax optimal setting for \(M\).

Now, returning to our mission 45 , by setting \(\delta^{-1}=c_{0}^2/4\cdot \log^{2s_l/(2s_l+d)} n\), we have: \[\begin{align} &\mathbb{P}\left(\|\hat{f}_{s_l}-f\|_{w,n}^2 M_l^{-2}(M_l^2 n/\log n)^{2s_l/(2s_l+d)}>c_0^2/4\right)\\&\le 16C c_0^{-4} n^{-1}\log^{-2s_l/(2s_l+d)} n+Cne^{-Cn\epsilon^{d}}+e^{-\lfloor M_{\min}^2 n \rfloor^{d/(2s+d)}}. \end{align}\] With 44 , we obtain: \[\begin{align} \mathbb{P}(\mathcal{E}_{j})\le 16C c_0^{-4} n^{-1}\log^{1-2s_{\min}/(2s_{\min}+d)} n+Cn e^{-Cn\epsilon^{d}}\log n +e^{-\lfloor M_{\min}^2 n \rfloor^{d/(2s+d)}}\log n . \end{align}\] Combining the above result with 43 and noting that on \(\mathcal{E}_{j}^{c}\), \(\mathbf{1}_{\mathcal{E}_{j}}=0\), it yields that \[\begin{align} \sum_{j=1}^{i-1}\left(\|\hat{f}_{s_j}-f\|_{w,n}^{2}M_i^{-2}(M_i^2 n/\log n )^{2s_i/(2s_i+d)}\mathbf{1}_{\mathcal{E}_{j}}\right)\lesssim \frac{1}{\delta}, \end{align}\] with probability at least \[1-\delta \log^{-2s_i/(2s_i+d)}n-16C c_0^{-4} n^{-1}\log^{2-2s_{\min}/(2s_{\min}+d)} n-Cn e^{-Cn\epsilon^{d}}\log^2 n -e^{-\lfloor M_{\min}^2 n \rfloor^{d/(2s+d)}}\log^2 n.\] ◻