Scalable Statistical Inference in Non-parametric Least Squares

To describe the structure of regression function $f$ in non-parametric regression model (1.1), we adopt the standard framework of a reproducing kernel Hilbert space (RKHS, [7, 27, 28]) by assuming $f^{\ast}=\mathop{\mathrm{argmin}}_{f}\mathbb{E}\big{[}(f(X)-Y)^{2}\big{]}$ to belong to an RKHS $\mathbb{H}$. Let $\mathbb{P}_{X}$ to denote the marginal distribution of the random design $X$, and $L^{2}(\mathbb{P}_{X})=\big{\{}f:\,\mathcal{X}\to\mathbb{R}\,\big{|}\,\int_{\mathcal{X}}f^{2}(X)\,\mathbb{P}_{X}(dx)<\infty\big{\}}$ to denote the space of all square-integrable functions over $\mathcal{X}$ with respect to $\mathbb{P}_{X}$. Briefly speaking, an RKHS is a Hilbert space $\mathbb{H}\subset L^{2}(\mathbb{P}_{X})$ of functions defined over a set $\mathcal{X}$, equipped with inner product $\langle\cdot,\,\cdot\rangle_{\mathbb{H}}$, so that for any $x\in\mathcal{X}$, the evaluation functional at $x$ defined by $L_{x}(f)=f(x)$ is a continuous linear functional on the RKHS. Uniquely associated with $\mathbb{H}$ is a positive-definite function $K:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$, called the reproducing kernel. The key property of the reproducing kernel is that it satisfies the reproducing property: the evaluation functional $L_{x}$ can be represented by the reproducing kernel function $K_{x}:\,=K(x,\,\cdot)$ so that $f(x)=L_{x}(f)=\langle K_{x},\,f\rangle_{\mathbb{H}}$. According to Mercer’s theorem [6], kernel function $K$ has the following spectral decomposition:

$$K(x,x^{\prime})=\sum_{j=1}^{\infty}\mu_{j}\,\phi_{j}(x)\,\phi_{j}(x^{\prime}),\,\,\,\,x,x^{\prime}\in\mathcal{X},$$

(2.1)

where the convergence is absolute and uniform on $\mathcal{X}\times\mathcal{X}$. Here, $\mu_{1}\geq\mu_{2}\geq\cdots\geq 0$ is the sequence of eigenvalues, and $\{\phi_{j}\}_{j=1}^{\infty}$ are the corresponding eigenfunctions forming an orthonormal basis in $L^{2}(\mathbb{P}_{X})$, with the following property: for any $j,k\in\mathbb{N}$,

$$\langle\phi_{j},\phi_{k}\rangle_{L^{2}(\mathcal{P}_{X})}=\delta_{jk}\quad\mbox{and}\quad\langle\phi_{j},\phi_{k}\rangle_{\mathbb{H}}=\delta_{jk}/\mu_{j},$$

where $\delta_{jk}=1$ if $j=k$ and $\delta_{jk}=0$ otherwise. Moreover, any $f\in\mathbb{H}$ can be decomposed into $f=\sum_{j=1}^{\infty}f_{j}\phi_{j}$ with $f_{j}=\langle f,\phi_{j}\rangle_{L_{2}(\mathbb{P}_{X})}$, and its RKHS norm can be computed via $\|f\|_{\mathbb{H}}^{2}=\sum_{j=1}^{\infty}\mu_{j}^{-1}f_{j}^{2}$.

We introduce some technical conditions on the reproducing kernel $K$ in terms of its spectral decomposition.

The uniform boundedness condition in Assumption A1 is common in the literature [29]. Assumption A2 assumes the kernel to have polynomially decaying eigenvalues. Assumption A1-A2 together also implies the kernel function is bounded as $\sup_{x}K(x,x)\leq c^{2}_{\phi}\sum_{k=1}^{\infty}k^{-\alpha}:=R^{2}$. One special class of kernels satisfying Assumptions A1-A2 is composed of translation-invariant kernels $K(t,s)=g(t-s)$ for some even function $g$ of period one. In fact, by utilizing the Fourier series expansion of the kernel function $g$, we observe that the eigenfunctions of the corresponding kernel matrix $K$ are trigonometric functions

$$\phi_{2k-1}(x)=\sin(\pi kx),\quad\phi_{2k}(x)=\cos(\pi kx),\quad k=1,2,\dots$$

on $\mathcal{X}=[0,1]$. It is easy to see that we can choose $c_{\phi}=1$ and $L=\pi$ to satisfy Assumption A1. Although we primarily consider kernels with eigenvalues that decay polynomially for the sake of clarity in this paper, it is worth mentioning that our theory extends to other kernel classes, such as squared exponential kernels and polynomial kernels [30].

To motivate functional SGD in RKHS, we first review SGD in Euclidean setting for minimizing the expected loss function $\mathcal{L}(\theta)=\mathbb{E}_{Z}[\ell(\theta;Z)]$, where $\theta\in\mathbb{R}^{d}$ is the parameter of interest, $\ell:\,\mathbb{R}^{d}\times\mathcal{Z}\to\mathbb{R}$ is the loss function and $Z$ denotes a generic random sample, e.g. $Z=(X,Y)$ in the non-parametric regression setting (1.1). By first-order Taylor’s expansion, one can locally approximate $\mathcal{L}(\theta+s)$ for any small deviation $s$ by $\mathcal{L}(\theta+s)\approx\mathcal{L}(\theta)+\langle\nabla\mathcal{L}(\theta),\,s\rangle$, where $\nabla\mathcal{L}(\theta)$ denotes the gradient (vector) of $\mathcal{L}(\cdot)$ evaluated at $\theta$. The gradient $\nabla\mathcal{L}(\theta)$ therefore encodes the (infinitesimal) steepest descent direction of $L$ at $\theta$, leading to the following gradient decent (GD) updating formula:

$\displaystyle\widehat{\theta}_{i}=\widehat{\theta}_{i-1}-\gamma_{i}\,\nabla L(\widehat{\theta}_{i-1}),\quad\mbox{for}\quad i=1,2,\ldots,$

starting from some initial value $\widehat{\theta}_{0}$, where $\gamma_{i}>0$ is the step size (also called learning rate) at iteration $i$. GD typically requires the computation of the full gradient $\nabla\mathcal{L}(\theta)$, which is unavailable due to the unknown data distribution of $Z$. In stochastic approximation, SGD takes a more efficient approach by using an unbiased estimate of the gradient as $G_{i}(\theta)=\nabla\ell(\theta,Z_{i})$ based on one sample $Z_{i}$ to substitute $\nabla\mathcal{L}(\theta)$ in the updating rule.

Accordingly, the SGD updating rule takes the form of

$\displaystyle\widehat{\theta}_{i}=\widehat{\theta}_{i-1}-\gamma_{i}\,G_{i}(\widehat{\theta}_{i-1}),\quad\mbox{for}\quad i=1,2,\ldots.$

Let us now extend the concept of SGD from minimizing an expected loss function in Euclidean space to minimizing an expected loss functional in function space. Here for concreteness, we develop SGD for minimizing the expected squared error loss $\mathcal{L}(f)=\mathbb{E}\big{[}(f(X)-Y)^{2}\big{]}$ over an RKHS $\mathbb{H}$ equipped with inner product $\langle\cdot,\cdot\rangle_{\mathbb{H}}$. Let us begin by extending the concept of the “gradient”. By identifying the gradient (operator) $\nabla L:\mathbb{H}\to\mathbb{H}$ of functional $\mathcal{L}(\cdot)$ as a steepest descent “direction” in $\mathbb{H}$ through the following first-order “Taylor expansion”

$\displaystyle\mathcal{L}(f)=\mathcal{L}(g)+\langle\nabla\mathcal{L}(g),\,f-g\rangle_{\mathbb{H}}+\mathcal{O}\big{(}\|f-g\|_{\mathbb{H}}^{2}\big{)},\ \ \mbox{as }f\to g,$

we obtain after some simple algebra that

	$\displaystyle\langle\nabla\mathcal{L}(g),\,f-g\rangle_{\mathbb{H}}+\mathcal{O}\big{(}\|f-g\|_{\mathbb{H}}^{2}\big{)}$	$\displaystyle\,=\mathcal{L}(f)-\mathcal{L}(g)$
		$\displaystyle\,=\mathbb{E}\big{[}\big{(}f(X)-g(X)\big{)}\cdot\big{(}g(X)-Y\big{)}\big{]}+\mathbb{E}\big{[}(f(X)-g(X))^{2}\big{]}.$

Now by using the reproducing property $h(x)=\langle h,\,K_{x}\rangle_{\mathbb{H}}$ for any $h\in\mathbb{H}$, we further obtain

$\displaystyle\langle\nabla\mathcal{L}(g),\,f-g\rangle_{\mathbb{H}}=\big{\langle}\mathbb{E}\big{[}\big{(}g(X)-Y\big{)}K_{X}\big{]},\,f-g\big{\rangle}_{\mathbb{H}}+\mathcal{O}\big{(}\|f-g\|_{\mathbb{H}}^{2}\big{)}.$

(2.2)

Here, we have used the fact that by Cauchy-Schwarz inequality,

$$(f(x)-g(x))^{2}=\langle f-g,\,K_{x}\rangle^{2}\leq\|f-g\|_{\mathbb{H}}^{2}\cdot\|K_{x}\|_{\mathbb{H}}^{2}=K(x,x)\,\|f-g\|_{\mathbb{H}}^{2}=\mathcal{O}\big{(}\|f-g\|_{\mathbb{H}}^{2}\big{)},$$

since Assumptions A1-A2 together with Mercer’s expansion (2.1) imply $K$ to be uniformly bounded, or $K(x,x)\leq c_{K}\sum_{j=1}^{\infty}\mu_{k}\leq C\sum_{j=1}^{\infty}j^{-\alpha}\leq C^{\prime}$, as long as $\alpha>1$. From equation (2.2), we can identify the gradient $\nabla\mathcal{L}(g)$ at $g\in\mathbb{H}$ as

$\displaystyle\nabla\mathcal{L}(g)=\mathbb{E}\big{[}\big{(}g(X)-Y\big{)}K_{X}\big{]}\in\mathbb{H}.$

Throughout the rest of the paper, we will refer to above $\nabla\mathcal{L}(g)$ as the RKHS gradient of functional $L$ at $g$.

Upon the arrival of the $i$th data point $(X_{i},Y_{i})$, we can form an unbiased estimator $G_{i}(g)$ of the RKHS gradient $\nabla\mathcal{L}(g)$ as $G_{i}(g)=\big{(}g(X_{i})-Y_{i}\big{)}K_{X_{i}}$. This leads to the following SGD in RKHS for solving non-parametric least squares: for a given initial estimate $\widehat{f}_{0}$, the SGD recursively updates the estimate of $f$ upon the arrival of each data point as

$$\widehat{f}_{i}=\widehat{f}_{i-1}-\gamma_{i}\,G_{i}(\widehat{f}_{i-1})=\widehat{f}_{i-1}+\gamma_{i}\big{(}Y_{i}-\widehat{f}_{i-1}(X_{i})\big{)}K_{X_{i}},\quad\mbox{for }i=1,2,\ldots.$$

(2.3)

By utilizing the reproducing property, the above iterative updating formula can be rewritten as

$\displaystyle\widehat{f}_{i}=\widehat{f}_{i-1}+\gamma_{i}\,\big{(}Y_{i}-\langle\widehat{f}_{i-1},K_{X_{i}}\rangle_{\mathbb{H}}\big{)}\,K_{X_{i}}=(I-\gamma_{i}\,K_{X_{i}}\otimes K_{X_{i}})\,\widehat{f}_{i-1}+\gamma_{i}\,Y_{i}\,K_{X_{i}},$

(2.4)

where $I$ denotes the identity map on $\mathbb{H}$, and $\otimes$ is the tensor product operator defined through $g\otimes h(f)=\langle f,h\rangle_{\mathbb{H}}\,g$ for all $g,h,f\in\mathbb{H}$. Formula (2.3) is more straightforward to use for practical implementation, while formula (2.4) provides a more tractable expression that will facilitate our theoretical analysis. Following [2] and [31], we consider the so-called Polyak averaging scheme to further improve the estimation accuracy by averaging over the entire updating trajectory, i.e. we use $\bar{f}_{n}=n^{-1}\sum_{i=1}^{n}\,\widehat{f}_{i}$ as the final functional SGD estimator of $f$ based on a dataset of sample size $n$. Note that this averaged estimator can be efficiently computed without storing all past estimators by using the recursively updating formula $\bar{f}_{i}=(1-i^{-1})\,\bar{f}_{i-1}\,+\,i^{-1}\,\widehat{f}_{i}$ for $i=1,\dots,n$ on the fly. We will refer to the above SGD as functional SGD in order to differentiate it from the SGD in Euclidean space, and $\bar{f}_{n}$ as the functional SGD estimator (using $n$ samples). Throughout the remainder of the paper, we consider a zero initialization, $\widehat{f}_{0}=0$, without loss of generality.

In functional SGD with total sample size (time horizon) $n$, the only adjustable component is the step size scheme $\{\gamma_{i}:\,i=1,2,\ldots,n\}$, which is crucial for achieving fast convergence and accurate estimations (c.f. Remark 3.1). We examine two common schemes [15, 32]: (1) constant step size scheme where $\gamma_{i}\equiv\gamma=\gamma(n)$ only depends on the total sample size $n$; (2) non-constant step size scheme where $\gamma_{i}=i^{-\xi}$ decays polynomially in $i$ for $i=1,2,\ldots,n$ and some $\xi>0$. While the constant step scheme is more amenable to theoretical analysis, it suffers from two notable drawbacks: (1) it assumes prior knowledge of the sample size $n$, which is typically unavailable in streaming data scenarios, and (2) the optimal estimation error is only achieved at the $n$-th iteration, leading to suboptimal performance before that time point. In contrast, the non-constant step size scheme, despite significantly complicating our theoretical analysis, overcomes the aforementioned limitations and leads to a truly online algorithm that achieves rate-optimal estimation at any intermediate time point (c.f. Theorem 3.1). Due to this characteristic, we will also refer to the non-constant step size scheme as the online scheme.

Although functional SGD operates in the infinite-dimensional RKHS, it can be implemented using a finite-dimensional representation enabled by the kernel trick. Concretely, upon the arrival of the $i$-th observation $(X_{i},Y_{i})$, we can express the time-$i$ intermediate estimator $\widehat{f}_{i}$ as $\widehat{f}_{i}=\sum_{j=1}^{i}\widehat{\beta}_{j}\,K_{X_{j}}$ due to equation (2.3) and the zero initialization $\widehat{f}^{\ast}=0$ condition, where only the last entry $\widehat{\beta}_{i}$ in the coefficient vector $(\widehat{\beta}_{1},\,\widehat{\beta}_{2},\,\dots,\,\widehat{\beta}_{i})^{\top}$ needs to be updated,

$\displaystyle\widehat{\beta}_{i}=\gamma_{i}\,\big{(}Y_{i}-\widehat{f}_{i-1}(X_{i})\big{)}=\gamma_{i}\,Y_{i}-\gamma_{i}\sum_{j=1}^{i-1}\widehat{\beta}_{j}\,K(X_{j},\,X_{i}).$

Note that the computational complexity at time $i$ is $\mathcal{O}(i)$ for $i=1,2,\ldots,n$. Correspondingly, the functional SGD estimator at time $i$ can be computed through $\bar{f}_{i}=(1-i^{-1})\,\bar{f}_{i-1}+i^{-1}\widehat{f}_{i}=\sum_{j=1}^{i}\bar{\beta}_{j}\,K_{X_{j}}$, where (can be proved by induction)

$\displaystyle\bar{\beta}_{j}=\Big{(}1-\frac{j-1}{i}\Big{)}\,\widehat{\beta}_{j},\quad\mbox{for}\quad j=1,2,\ldots,i.$

Consequently, the overall time complexity of the resulting algorithm is $\mathcal{O}(n^{2})$, and the space complexity is $\mathcal{O}(n)$.

Our objective is to develop online inference for the non-parametric regression function $f^{\ast}$ based on the functional SGD estimator $\bar{f}_{n}$. Specifically, we aim to construct level-$\beta$ pointwise confidence intervals (local inference) $CI_{n}(x;\,\beta)=[U_{n}(x;\,\beta),\,V_{n}(x;\,\beta)]$ for $f^{\ast}(x)$, where $x\in\mathcal{X}$, and a level-$\beta$ simultaneous confidence band (global inference) $CB_{n}(\beta)=\big{\{}g:\,\mathcal{X}\to\mathbb{R}\,\big{|}\,g(x)\in[\bar{f}_{n}(x)-b_{n}(\beta),\,\bar{f}_{n}(x)+b_{n}(\beta)],\ \forall x\in\mathcal{X}\big{\}}$ for $f^{\ast}$. We require these intervals and band to be asymptotically valid, meaning that the coverage probabilities, i.e., the probabilities of $f^{\ast}(x)$ or $f^{\ast}$ falling within $CI_{n}(x;\,\beta)$ or $CB_{n}(\beta)$ respectively, are close to their nominal level $\beta$. Mathematically, this means $\mathbb{P}[f^{\ast}(x)\in CI_{n}(x;\,\beta)]=\beta+o(1)$ and $\mathbb{P}[f^{\ast}\in CB_{n}(\beta)]=\beta+o(1)$ as $n\to\infty$.

The coverage probability analysis of these intervals and band requires us to examine and prove the distributional convergence of two random quantities (with appropriate rescaling) based on the functional SGD estimator $\bar{f}_{n}$: the pointwise difference $\bar{f}_{n}(x)-f^{\ast}(x)$ for $x\in\mathcal{X}$ and the supremum norm $\|\bar{f}_{n}-f^{\ast}\|_{\infty}$ of $\bar{f}_{n}-f^{\ast}$. In particular, the appropriate rescaling choice determines a precise convergence rate of $\bar{f}$ towards $f^{\ast}$ under the supremum norm metric. The characterization of the convergence rate of a non-parametric regression estimator under the supremum norm metric is a challenging and important problem in its own right. We note that the distribution of the supremum norm $\|\bar{f}_{n}-f^{\ast}\|_{\infty}$ after a proper re-scaling behaves like the supreme norm of a Gaussian process in the asymptotic sense, which is not practically feasible to estimate. Therefore, for inference purposes, it is not necessary to explicitly characterize this distributional limit; instead, we will prove a bootstrap consistency by showing that the Kolmogorov distance between the sampling distributions of this supremum norm and its bootstrapping counterpart converges to zero as $n\to\infty$.

In our theoretical development to address these problems, we will utilize a recursive expansion of the functional SGD updating formula to construct a higher-order expansion of $\bar{f}_{n}$ under the $\|\cdot\|_{\infty}$ norm metric. Building upon this expansion, we will establish in Section 3 the distributional convergence of the two aforementioned random quantities and characterize their limiting distributions with an explicit representation of the limiting variance for $\bar{f}_{n}(x)-f^{\ast}(x)$ in the large-sample setting. However, these limiting distributions and variances depend on the spectral decomposition of the kernel $K$, the marginal distribution of the design variable $X$, and the unknown noise variance $\sigma^{2}$, which are either inaccessible or computationally expensive to evaluate in an online learning scenario. To overcome this challenge, we will propose a scalable bootstrap-based inference method in Section 4, enabling efficient online inference for $f^{\ast}$.

We begin by decomposing the functional SGD update of $\widehat{f}_{n}-f^{\ast}$ into two leading recursive formulas and a higher-order remainder term. This decomposition allows us to distinguish between the deterministic term responsible for the estimation bias and the stochastic fluctuation term contributing to the estimation variance. Concretely, we obtain the following by plugging $Y_{i}=f^{\ast}(X_{i})+\epsilon_{i}$ into the recursive updating formula (2.4),

$$\widehat{f}_{i}-f^{\ast}=(I-\gamma_{i}\,K_{X_{i}}\otimes K_{X_{i}})\,(\widehat{f}_{i-1}-f^{\ast})+\gamma_{i}\,\epsilon_{i}\,K_{X_{i}}.$$

(3.1)

Let $\Sigma:\,=\mathbb{E}[K_{X_{1}}\otimes K_{X_{1}}]:\,\mathbb{H}\to\mathbb{H}$ denote the population-level covariance operator, so that for any $f$, $g\in\mathbb{H}$ we have $\langle f,\,\Sigma\,g\rangle_{\mathbb{H}}=\mathbb{E}[f(X_{1})\,g(X_{1})]$. Now we recursively define the leading bias term through

$$\eta_{0}^{bias,0}=\widehat{f}_{0}-f^{\ast}=-f^{\ast}\quad\mbox{and}\quad\eta_{i}^{bias,0}=(I-\gamma_{i}\,\Sigma)\,\eta_{i-1}^{bias,0}\quad\mbox{for}\quad i=1,2,\ldots$$

(3.2)

that collects the leading deterministic component in (3.1); and the leading noise term through

$\displaystyle\eta_{0}^{noise,0}=0\quad\mbox{and}\quad\eta_{i}^{noise,0}=(I-\gamma_{i}\,\Sigma)\,\eta^{noise,0}_{i-1}+\gamma_{i}\,\epsilon_{i}\,K_{X_{i}}\quad\mbox{for}\quad i=1,2,\ldots$

(3.3)

that collects the leading stochastic fluctuation component in (3.1); so that we have the following decomposition for the recursion:

$$\widehat{f}_{i}-f^{\ast}=\underbrace{\eta_{i}^{bias,0}}_{\text{leading bias}}+\underbrace{\eta_{i}^{noise,0}}_{\text{leading noise}}+\ \ \underbrace{\big{(}\widehat{f}_{i}-f^{\ast}-\eta_{i}^{bias,0}-\eta_{i}^{noise,0}\big{)}}_{\text{remainder term}}\quad\mbox{for}\quad i=1,2,\ldots.$$

(3.4)

Correspondingly, we define $\bar{\eta}_{i}^{bias,0}=i^{-1}\sum_{j=1}^{i}\eta_{j}^{bias,0}$ and $\bar{\eta}_{i}^{noise,0}=i^{-1}\sum_{j=1}^{i}\eta_{j}^{noise,0}$ as the leading bias and noise terms, respectively, in the functional SGD estimator (after averaging). The following Theorem 3.1 presents finite-sample bounds for the two leading terms and the remainder term associated with $\bar{f}_{n}$ under the supreme norm metric. The results indicate that the remainder term is of strictly higher order (in terms of dependence on $n$) compared to the two leading terms, validating the term “leading” for them.

As we discussed in Section 2.3, for inference purposes, it is not necessary to explicitly characterize distributional convergence limit of the supremum norm $\|\bar{f}_{n}-f^{\ast}\|_{\infty}$; instead, we will prove a bootstrap consistency by showing that the Kolmogorov distance between the sampling distributions of this supremum norm and its bootstrapping counterpart converges to zero as $n\to\infty$. However, the pointwise convergence limit of $\bar{f}_{n}(z_{0})-f^{\ast}(z_{0})$ for fixed $z_{0}\in[0,1]$ has an easy characterization. Therefore, we present the pointwise convergence limit and use it to discuss the impact of online estimation in the non-parametric regression model in the following subsection.

According to Theorem 3.1, the large-sample behavior of the functional SGD estimator $\bar{f}_{n}$ is completely determined by the two leading processes: bias term and noise term. According to (3.2), under the constant step size $\gamma$, the leading bias term has an explicit expression as

	$\displaystyle\bar{\eta}_{n}^{bias,0}(x)=$	$\displaystyle\frac{1}{n}\gamma^{-1}\Sigma^{-1}\,(I-\gamma\Sigma)\,[I-(I-\gamma\Sigma)^{n}\,]f^{\ast}(x)$		(3.5)
	$\displaystyle=$	$\displaystyle\frac{1}{\sqrt{\gamma}n}\sum_{k=1}^{n}\sum_{\nu=1}^{\infty}\langle f^{\ast},\phi_{\nu}\rangle_{L_{2}}\mu_{\ell}^{-1/2}(1-\gamma\mu_{\nu})^{k}(\gamma\mu_{\nu})^{1/2}\phi_{\nu}(x),\quad\forall x\in\mathcal{X},$

and the leading noise term is

	$\displaystyle\bar{\eta}_{n}^{noise,0}(x)=$	$\displaystyle\,\frac{1}{n}\sum_{k=1}^{n}\Sigma^{-1}\big{[}I-(I-\gamma\Sigma)^{n+1-k}\big{]}\,K(X_{k},\,x)\,\epsilon_{k}$		(3.6)
	$\displaystyle=$	$\displaystyle\,\frac{1}{n}\sum_{k=1}^{n}\ \epsilon_{k}\,\cdot\,\underbrace{\bigg{\{}\sum_{\nu=1}^{\infty}\big{[}1-(1-\gamma\mu_{\nu})^{n+1-k}\big{]}\,\phi_{\nu}(X_{k})\,\phi_{\nu}(x)\bigg{\}}}_{\Omega_{n,k}(x)},\quad\forall x\in\mathcal{X}.$

For each fixed $z_{0}\in\mathcal{X}$, conditioning on the design $\{X_{i}\}_{i=1}^{n}$, the leading noise term $\bar{\eta}_{n}^{noise,0}(z_{0})$ is a weighted average of $n$ independent and centered normally distributed random variables. This representation enables us to identify the limiting distribution of $\bar{\eta}_{n}^{noise,0}(z_{0})$ (this subsection) and conduct local inference (i.e. pointwise confidence intervals) by a bootstrap method (next section). Under Assumption A2, the weight $\Omega_{n,k}(z_{0})$ associated with the $k$-th observation pair $(X_{k},\,Y_{k})$ is of order $\sum_{\nu=1}^{\infty}\big{[}1-(1-\gamma\mu_{\nu})^{n+1-k}\big{]}\asymp\big{[}(n+1-k)\gamma\big{]}^{1/\alpha}$, which decreases in $k$. This diminishing impact trend is inherent to online learning, as later observations tend to have a smaller influence compared to earlier observations. This characteristic is radically different from offline estimation settings, where all observations contribute equally to the final estimator, and will change the asymptotic variance (i.e., the $\sigma^{2}_{z_{0}}$ in Theorem 3.2).

Furthermore, the entire leading noise process $\bar{\eta}_{n}^{noise,0}(\cdot)$ can be viewed as a weighted and non-identically distributed empirical process indexed by the spatial location. This characterization enables us to conduct global inference (i.e. simultaneous confidence band) for non-parametric online learning by borrowing and extending the recent developments [26, 34, 35] on Gaussian approximation and multiplier bootstraps for suprema of (equally-weighted and identically distributed) empirical processes, which will be the main focus of next section.

In the following Theorem 3.2, we prove, by analyzing the leading noise term $\bar{\eta}_{n}^{noise,0}$, a finite-sample upper bound on the Kolmogorov distance between the sampling distribution of $\bar{f}_{n}(z_{0})-f^{\ast}(z_{0})$ and the distribution of a standard normal random variable (i.e. supreme norm between the two cumulative distributions) for any $z_{0}\in\mathcal{X}$.