# Nonparametric Estimation of Ordinary Differential Equations: Snake and
Stubble

Christof Schötz<EMAIL_ADDRESS> 0000-0003-3528-4544 Potsdam
Institute for Climate Impact Research  Technical University of Munich

###### Abstract

We study nonparametric estimation in dynamical systems described by ordinary
differential equations (ODEs). Specifically, we focus on estimating the
unknown function $f\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$ that governs the
system dynamics through the ODE $\dot{u}(t)=f(u(t))$, where observations
$Y_{j,i}=u_{j}(t_{j,i})+\varepsilon_{j,i}$ of solutions $u_{j}$ of the ODE are
made at times $t_{j,i}$ with independent noise $\varepsilon_{j,i}$.

We introduce two novel models – the Stubble model and the Snake model – to
mitigate the issue of observation location dependence on $f$, an inherent
difficulty in nonparametric estimation of ODE systems. In the Stubble model,
we observe many short solutions with initial conditions that adequately cover
the domain of interest. Here, we study an estimator based on multivariate
local polynomial regression and univariate polynomial interpolation. In the
Snake model we observe few long trajectories that traverse the domain on
interest. Here, we study an estimator that combines univariate local
polynomial estimation with multivariate polynomial interpolation.

For both models, we establish error bounds of order
$n^{-\frac{\beta}{2(\beta+1)+d}}$ for $\beta$-smooth functions $f$ in an
infinite dimensional function class of Hölder-type and establish minimax
optimality for the Stubble model in general and for the Snake model under some
conditions via comparison to lower bounds from parallel work.

Keywords: Ordinary Differential Equations, Nonparametric, Upper Bounds,
Optimal Rate of Convergence, Minimax

###### Contents

  1. 1 Introduction
    1. 1.1 Contributions
    2. 1.2 Relation to Other Models
    3. 1.3 Overview
  2. 2 Preliminaries
    1. 2.1 Ordinary Differential Equations
    2. 2.2 Formal Description of the General Model
    3. 2.3 Need for Restriction of the General Model
    4. 2.4 Standard Regression Model
  3. 3 The Stubble Model and Estimation of the Increment Map
    1. 3.1 Lipschitz
      1. 3.1.1 Model
      2. 3.1.2 Estimator
      3. 3.1.3 Result
      4. 3.1.4 Proof
    2. 3.2 General
      1. 3.2.1 Model
      2. 3.2.2 Estimation
      3. 3.2.3 Result
      4. 3.2.4 Proof
  4. 4 The Snake Model and Estimation of the Observed Solutions
    1. 4.1 Lipschitz
      1. 4.1.1 Model
      2. 4.1.2 Estimator
      3. 4.1.3 Result
      4. 4.1.4 Proof
    2. 4.2 General
      1. 4.2.1 Model
      2. 4.2.2 Multivariate Polynomial Interpolation
      3. 4.2.3 Estimator
      4. 4.2.4 Results
      5. 4.2.5 Proof
  5. A Derivatives
  6. B Local Polynomial Regression
    1. B.1 Model
    2. B.2 Estimator
    3. B.3 Result
    4. B.4 Proof
  7. C Polynomial Interpolation
    1. C.1 Univariate Polynomial Interpolation
    2. C.2 Univariate Polynomial Comparison
    3. C.3 Multivariate Polynomial Interpolation
  8. D Ordinary Differential Equations and Smoothness
    1. D.1 Difference Bound
    2. D.2 Smoothness of Solutions in Time
    3. D.3 Smoothness of Solutions in the Initial Conditions

##  1 Introduction

In statistics, we want to extract information from observations of the world.
The workings of the world are described in physics (and many other sciences)
with differential equations [Str24]. Hence, statistical models of dynamical
systems described by differential equations seem to be not only compelling but
almost inevitable for study.

One of the simplest non-trivial instance of such a model can be described in
the equations

$Y_{i}=u(t_{i})+\varepsilon_{i}\,,i=1,\dots,n\qquad\dot{u}(t)=f(u(t))\,,$ (1)

where $\dot{u}$ denotes the time-derivative, i.e.,
$(\mathrm{d}u)/(\mathrm{d}t)$. Here, $u\colon\mathbb{R}\to\mathbb{R}^{d}$ is
the solution to an ordinary differential equation (ODE) described by an
unknown function $f\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$. It is measured at
times $t_{1}\leq\dots\leq t_{n}$, which creates our observations $Y_{i}$ by
adding measurement noise $\varepsilon_{i}$ to the true state of the system
$u(t_{i})$. There are different objectives one could have for such a model:
reconstructing the solution $u$ in the observed time interval $[t_{1},t_{n}]$
(reanalysis, borrowing a term from climate modeling), predicting a future
state $u(t)$, $t>t_{n}$ (forecasting, as in weather forecast), or estimation
of $f$ (learning the dynamics of the dynamical system). Here, we are
interested in the estimation of $f$, which is fundamental also for the other
tasks as an estimate $\hat{f}$ allows create an estimate $\hat{u}$ by solving
$\dot{\hat{u}}(t)=\hat{f}(\hat{u}(t))$.

We want to estimate $f$ in a nonparametric model, i.e., we assume that $f$
belongs to an infinite dimensional class of smooth functions, but we do not
make any assumptions on its precise functional form. The parametric version of
this model, in which $f$ belongs to a finite dimensional class of functions
(e.g., polynomials of fixed degree) has been studied in [Bru08, QZ10, GK12,
DK15]. There, it is shown that the typical parametric $\sqrt{n}$-rate of
convergence can be achieved. An overview of the state of the art for
statistics of dynamical systems (mostly concerned with parametric settings)
can be found in [MMP15, RH17, Dat21].

In the nonparametric setting, different algorithms for estimation have been
proposed, e.g., [Hei+18, Che+18, GR21]. In [Mar+23], theoretical results for
learning dynamics nonparametrically are shown, but in a density estimation
context that is rather different from the regression-type model studied here.
To the best of the author’s knowledge, only [Lah+24] take a theoretical view
on the ODE regression problem (additionally to proposing a new algorithm based
on reproducing kernel Hilbert spaces). They show an error bound of order
$n^{-1/4}$ (root mean squared error) for the reanalysis problem [Lah+24,
Theorem 1]. This result cannot be optimal as standard nonparametric regression
for $u$ (ignoring the ODE constraints) yields an error of smaller order,
namely $n^{-\beta/(2\beta+1)}$, where $\beta$ describes the smoothness of $u$
(in the case of [Lah+24, Theorem 1], we have $\beta=3$).

Estimation in the nonparametric setting seem much harder than in the
parametric one: If $f$ is a smooth function with bounds on its derivative but
otherwise unrestricted, an observation $Y_{i}$ only yields information on
$f(x)$ for $x$ close to $u(t_{i})$, i.e., the information is local. In
contrast, each observation from an ODE where $f$ is known to be a polynomial
of degree at most $N$ carries information for all coefficients of said
polynomial, i.e., the information is global. What makes estimation in a
nonparametric ODE model more difficult than in a nonparametric regression
model (which also suffers from locality of information) is that fact that the
location of the observation $u(t_{i})$ is itself depends on the unknown $f$.

We introduce two models, called Snake model and Stubble model, which provide a
remedy of this problem. For each of the two models, we propose estimators and
study their rate of convergence. We find that in both instances, one obtains
an error bound (root mean squared error) of order

$n^{-\frac{\beta}{2(\beta+1)+d}}$ (2)

for $\beta$-smooth functions $f$ under optimal circumstances. This error bound
is minimax optimal as we have corresponding lower bounds in the parallel work
[SS24].

###  1.1 Contributions

We now give more details on the contributions of this work.

A general model. In general, we assume the true model function $f^{\star}$ to
be an element of Hölder-type smoothness class denoted as $\bar{\Sigma}^{d\to
d}(\beta,L_{{\left\llbracket 0,\beta\right\rrbracket}})$, which means that the
$k$-th derivative of $f$ for $k\in{\left\llbracket
0,\beta\right\rrbracket}:=\\{0,\dots,\beta\\}$ is bounded by
$L_{k}\in\mathbb{R}_{>0}$, i.e.,
$\sup_{x\in\mathbb{R}^{d}}\|D^{k}f(x)\|_{\mathsf{op}}\leq L_{k}$, where
$\|\cdot\|_{\mathsf{op}}$ is the operator norm. Then, for given initial
conditions $x_{1},\dots,x_{m}$, we observe the solution $t\mapsto
U(f^{\star},x_{j},t)$ of the ODE $\dot{u}(t)=f^{\star}(u(t))$ with
$U(f^{\star},x_{j},0)=x_{j}$ at time points $i\Delta\\!t$ with time step
$\Delta\\!t\in\mathbb{R}_{>0}$, i.e.,

$Y_{j,i}=U(f^{\star},x_{j},i\Delta\\!t)+\varepsilon_{j,i}\,,\
j\in{\left\llbracket 1,m\right\rrbracket},i\in{\left\llbracket
1,n_{j}\right\rrbracket}\,.$ (3)

In total, we have $n=\sum_{j=1}^{m}n_{j}$ observations. The noise variables
$\varepsilon_{j,i}$ are assumed to be independent, centered, and to have a
finite second moment. Because of the dependency of the location of the
observation $u(t_{j,i})$ on $f^{\star}$, in general, consistent estimation in
all of $[0,1]^{d}$ is impossible. Thus, we introduce the Stubble model and the
Snake model, which restrict the general model.

The Stubble model. We use the notation $\asymp$ and $\preccurlyeq$ to mean
asymptotically equal and asymptotically lower than (up to a positive
constant), respectively (see 2.3). In the Stubble model, we observe many
($m\asymp n$) short ($n_{j}\preccurlyeq 1$) solutions. We assume that their
initial conditions $x_{j}$ cover the domain of interest $[0,1]^{d}$ suitably.
For this model, we construct an estimator $\hat{f}$ based on a multivaraite
local polynomial estimator and an univariate polynomial interpolation, which
is similar to the Adams–Bashforth method for numerical solutions of ODEs
[But16, Chapter 24]. We obtain

$\mathbf{E}\\!\left[\left\|\hat{f}(x)-f^{\star}(x)\right\|_{2}^{2}\right]\preccurlyeq\left(\Delta\\!t^{2}n\right)^{-\frac{2\beta}{2\beta+d}}+\Delta\\!t^{2\beta}$
(4)

for all $x\in[0,1]^{d}$, see Corollary 3.13. The rate is shown to be minimax
optimal by comparing it to lower bounds in [SS24].

The Snake model. In the Snake model, we observe few ($m\preccurlyeq 1$) long
($n_{j}\asymp n$) solutions. We require the trajectories
$U(f^{\star},x_{j},[0,n_{j}\Delta\\!t])$ to cover the domain of interest
$[0,1]^{d}$ suitably. For this model, we construct an estimator $\hat{f}$
based on an univariate local polynomial estimator and a multivariate
polynomial interpolation. Then

$\sup_{x\in[0,1]^{d}}\left\|\hat{f}(x)-f^{\star}(x)\right\|_{2}^{2}\in\mathbf{O}_{\mathbf{P}}\\!\left(\delta^{2\beta}+\left(\Delta\\!t\log
n\right)^{\frac{2\beta}{2(\beta+1)+1}}\right)\,,$ (5)

see Corollary 4.17, where we want to view $\delta\in\mathbb{R}_{>0}$ for now
as the largest distance between a point $x\in[0,1]^{d}$ and its closest state
$U(f^{\star},x_{j},t)$, $t\in[0,T]$. This interpretation is true for the case
$\beta=1$, but if $\beta>1$, the definition is more complex and the result
more restrictive. By comparing it to lower bounds in [SS24], the rate is shown
to be minimax optimal if $\delta$ and $T$ are in a certain relation. For the
rate to be optimal, we essentially require observations to be rather dense in
time, and temporally distant parts of the trajectories
$U(f^{\star},x_{j},[0,T])$ to be distant enough in state space.

Connection and further results. Note the complementary nature of the two
models and their estimators. In spite of this, in the optimal setting
regarding $\Delta\\!t$ and $\delta$, we obtain an error bound of order (2) in
both models (up to a $\log$-factor when considering the sup-norm), see
Corollary 3.15 and Corollary 4.19. In both models, we start with the case
$\beta=1$, which allows for simple estimators and a gentle introduction of the
main ideas as well as slightly stronger and more specific results. See Theorem
3.2 and Corollary 3.4 for the Stubble model, and Theorem 4.4 and Corollary 4.6
for the Snake model. For the case of general smoothness $\beta\in\mathbb{N}$,
we first show results for estimation strategies that can be used with an
arbitrary regression estimator, Theorem 3.11 and Theorem 4.16, respectively.
If the chosen regression estimator achieves the minimax rate for a standard
nonparametric regression problem, the estimation strategies achive the rates
(4) and (5), respectively. In Corollary 3.13 and Corollary 4.17, we apply the
general results with the (minimax optimal) local polynomial estimator as
regression estimator. Furthermore, we study how the smoothness of $f$
influences the smoothness of $t\mapsto U(f,x,t)$ and the smoothness of
$x\mapsto U(f,x,t)$ in Corollary D.7 and Corollary D.12, respectively.

###  1.2 Relation to Other Models

Regression on the solutions. The ODE model can be viewed as a regression
model, with the ODE as a constraint on the regression function
$u^{\star}=U(f^{\star},x,\cdot)$ in $[0,T]$. Solutions have a smoothness
parameter $\tilde{\beta}=\beta+1$ if the model function is $\beta$-smooth.
Thus, we can solve the reanalysis problem, i.e., estimation of $u^{\star}$ in
$[0,T]$, with a the standard nonparametric rate of convergence for the squared
error

$\mathbf{E}\\!\left[\left\|\hat{u}(t)-u^{\star}(t)\right\|^{2}\right]\preccurlyeq\left(\frac{T}{n}\right)^{\frac{2(\beta+1)}{2(\beta+1)+1}}$
(6)

by ignoring the ODE-constraint and using a suitable estimator $\hat{u}$, see,
e.g., [Tsy09]. Furthermore, the derivative of the solution can similarly be
estimated with

$\mathbf{E}\\!\left[\left\|\hat{\dot{u}}(t)-\dot{u}^{\star}(t)\right\|^{2}\right]\preccurlyeq\left(\frac{T}{n}\right)^{\frac{2\beta}{2(\beta+1)+1}}\,.$
(7)

As $\dot{u}^{\star}(t)=f^{\star}(u^{\star}(t))$, this essentially yields an
estimate of $f^{\star}$ on the set $u^{\star}([0,T])$. This view is taken for
our estimator in the Snake model of section 4.

Minimax rate of nonparametric regression. For standard nonparametric
estimation of the $s$-th derivative of a $\tilde{\beta}$-smooth regression
function $g^{\star}\colon\mathbb{R}^{d}\to\mathbb{R}$ (or
$\mathbb{R}^{d}\to\mathbb{R}^{d}$), we have the minimax rate

$\mathbf{E}\\!\left[\left\|\widehat{D^{s}g}(x)-D^{s}g^{\star}(x)\right\|^{2}\right]\preccurlyeq
n^{-\frac{2(\tilde{\beta}-s)}{2\tilde{\beta}+d}}\,,$ (8)

see, e.g., appendix B. In the ODE model, we estimate a first derivative
($s=1$) in $d$-dimensions of the solutions that have smoothness
$\tilde{\beta}=\beta+1$. So, the rate of convergence we obtain in this work
(2) fits the formula of the minimax regression rate (8), even though the
models are quite different.

Errors in variables. Assume that time steps are constant, i.e.,
$t_{i+1}-t_{i}=\Delta\\!t$. Consider the so-called propagator $x\mapsto
U(f^{\star},x,\Delta\\!t)$. If $\Delta\\!t\to 0$, we can infer $f$ from it, as

$\frac{U(f^{\star},x,\Delta\\!t)-x}{\Delta\\!t}=\frac{U(f^{\star},x,\Delta\\!t)-U(f^{\star},x,0)}{\Delta\\!t}\xrightarrow{\Delta\\!t\to
0}f^{\star}(x)\,.$ (9)

Estimating the propagator from the data $(Y_{i-1},Y_{i})$ can be viewed as an
errors-in-variables regression problem, e.g., [MRS12]. In such a model one has
observation of the response variable $Z_{i}=g(x_{i})+\varepsilon_{i}$ and
assumes the predictors $x_{i}$ to be unknown but instead to be measured
noisily, e.g., via $X_{i}=x_{i}+\xi_{i}$ for a noise variable $\xi_{i}$. With
$(X_{i},Z_{i})=(Y_{i-1},Y_{i})$, estimation of the propagator in the ODE model
gives rise to an errors-in-variables regression problem. The Stubble model of
section 3 removes the problems coming with errors-in-variables, by assuming
known initial conditions and observing many solutions for only a short time.

Nonparametric time series analysis. The ODE model with the sequence $(Y_{i})$
can be view from the perspective of nonparametric time series analysis, where
discrete-time stochastic processes in one variable (time) are studied. See,
e.g., [FY03, Chapter 6]. On one hand, the focus on discrete time and often
only one state dimension, allows for a deep understanding of more complex time
evolutions than in an one-dimensional ODE model. On the other hand, it seems
that the structure of higher dimensional, time continuous ODEs is typically
not captured in such models.

Stochastic differential equations. The noise in stochastic differential
equation (SDE) models, e.g., [Bis08, CG20] can be described as system noise in
contrast to the measurement noise of the ODE model: It changes the state of
the system and solutions are stochastic processes. In the ODE model, the state
of the system is not influenced by the noise and solutions are deterministic.
In comparison to a model with observed solutions of a SDE in one variable
(time), estimation in the ODE model seems more difficult: The observation
$Z_{t}$ at time $t$ of the solution $(Z_{t})$ in the SDE model is the true
state of the system at that time. If the SDE is Markovian, then the
distribution of the future states $(Z_{t^{\prime}})_{t^{\prime}>t}$ only
depend on the now known state $Z_{t}$. In the ODE model, we do not know the
true state $u^{\star}(t)$, but only a noisy observation of it. Even though the
future evolution is deterministic, it can be hard to predict if the initial
state is not know perfectly (see the topic of chaotic systems, e.g., in
[Str24]).

Hidden markov models. As the true state of the system $u^{\star}(t)$ is hidden
from us and only noisily observed, the ODE model is connected to hidden Markov
models (HMMs), e.g., [BRR98]. In HMMs, on observes $Z_{i}$ which depends on
the unobserved random variable $X_{i}$, which in turn depends on $X_{i-1}$.
But $Z_{i}$ is independent of $X_{i-1}$ and $Z_{i-1}$ given $X_{i}$. If we set
$X_{i}=u^{\star}(t_{i})$, $Z_{i}=Y_{i}$. We can view the ODE model as a HMM.
In contrast to HMMs, the transition from $u^{\star}(t_{i})$ to
$u^{\star}(t_{i+1})$ is deterministic, and the time and state domains are
continuous.

###  1.3 Overview

The remaining article is structured as follows. In section 2, we introduce
basic concepts for the study of ordinary differential equations (section 2.1),
give a full formal description of the general model (section 2.2), argue why
this model is not useful without further restriction (section 2.3), and
introduce some notation for standard regression problems that allows us to
formulate estimation strategies referring to arbitrary regression estimators
(section 2.4). Sections 3 and 4 formally introduce the Stubble and Snake
model, respectively, and describe the estimators $\hat{f}$ of $f^{\star}$,
upper bounds on the error $\|\hat{f}(x)-f^{\star}(x)\|_{2}$ and their proofs.
Both sections are separated into two parts: The first one (section 3.1 and 4.1
respectively) concerns the Lipschitz-case ($\beta=1$). The second one (section
3.2 and 4.2 respectively) concerns the general case ($\beta\in\mathbb{N}$).
The appendix makes this article largely self-contained: We first state some
basics of multivariate derivatives and smoothness classes (appendix A). Then,
we recall classical results on the upper error bound for multivariate local
polynomial regression in appendix B. In appendix C, uni- and multivariate
polynomial interpolation are discussed. And finally, appendix D provides
smoothness results for the solutions of ODEs.

##  2 Preliminaries

In this section, we recall some basic concepts related to ODEs, we introduce a
general statistical model for observing solutions of ODEs, we explain why this
model requires further restriction in order to be a useful model, and we
introduce some terminology for a standard regression problem, which will later
allows us to construct ODE estimators in a generic fashion.

###  2.1 Ordinary Differential Equations

We introduce some basic terminology and fundamental properties concerning
ODEs. See also [Har02, Arn06].

######  Notation 2.1.

  1. (i)

Let $\mathbb{Z}$ be the set of integers and $\mathbb{R}$ the set of reals. For
$\mathbb{K}\in\\{\mathbb{Z},\mathbb{R}\\}$ and $a\in\mathbb{R}$, denote
$\mathbb{K}_{>a}=\\{x\in\mathbb{K}\mid x>a\\}$. Define $\mathbb{K}_{\geq a}$,
$\mathbb{K}_{<a}$, $\mathbb{K}_{\leq a}$ accordingly.

  2. (ii)

Set $\mathbb{N}:=\mathbb{Z}_{\geq 1},\mathbb{N}_{0}:=\mathbb{Z}_{\geq 0}$. Let
$a,b\in\mathbb{Z}$ with $a\leq b$. Set ${\left\llbracket
a,b\right\rrbracket}:=\mathbb{Z}_{\geq a}\cap\mathbb{Z}_{\leq b}$. Set
${\left\llbracket a\right\rrbracket}:={\left\llbracket 1,a\right\rrbracket}$.

Let $d\in\mathbb{N}$. Let $f\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$. For a
differentiable function $u\colon\mathbb{R}\to\mathbb{R}^{d}$, denote its
derivative as $\dot{u}\colon\mathbb{R}\to\mathbb{R}^{d}$. Then

$\dot{u}(t)=f(u(t))\qquad\text{for }t\in\mathbb{R}$ (10)

is an autonomous, first-order, ordinary differential equation. It is of first-
order, as (10) only depends on the first derivative of $u$. It is autonomous,
as the right-hand side term $f(u(t))$ only depends on $t$ via $u$. In
contrast, a non-autonomous, first-order ODE has the form
$\dot{u}(t)=g(t,u(t))$ for a function
$g\colon\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}^{d}$. Any differentiable
$u\colon\mathbb{R}\to\mathbb{R}^{d}$ that fulfills (10) is a solution to the
ODE. The domain of $u$ is called time. The codomain of $u$ is called state
space. A single element of the state space is a state. The image of $u$ is
called trajectory. We call $f$ the model function of the ODE. Let
$x\in\mathbb{R}^{d}$. Consider the requirement

$u(0)=x\,.$ (11)

We call $x$ the initial conditions. We call (10) together with (11) initial
value problem (IVP). If $u$ fulfills (10) and (11), it is a solution to the
IVP.

Assume that $f$ is (globally) Lipschitz continuous. Then the IVP (10), (11)
has a unique solution (Picard–Lindelöf theorem, also known as the
Cauchy–Lipschitz theorem). Denote this solution as
$U(f,x,\cdot)\colon\mathbb{R}\to\mathbb{R}^{d},t\mapsto U(f,x,t)$. The
function
$U(f,\cdot,\cdot)\colon\mathbb{R}^{d}\times\mathbb{R}\to\mathbb{R}^{d}$ is
called flow of the ODE (10). Note the semigroup property of the flow,

$\displaystyle U(f,x,0)$ $\displaystyle=x\,,$ (12) $\displaystyle U(f,x,s+t)$
$\displaystyle=U(f,U(f,x,s),t)$ (13)

for all $x\in\mathbb{R}^{d}$, $s,t\in\mathbb{R}$. For a given time step
$\Delta\\!t\in\mathbb{R}_{>0}$, the mapping $x\mapsto U(f,x,\Delta\\!t)$ is
called propagator. If $u$ is a solution to the ODE (10), then
$u(t+\Delta\\!t)=U(f,u(t),\Delta\\!t)$ for all $t\in\mathbb{R}$. Furthermore,
we call

$\Upsilon(f,\Delta\\!t,x):=U(f,x,\Delta\\!t)-x=U(f,x,\Delta\\!t)-U(f,x,0)$
(14)

the increment. If $u$ is a solution to the ODE (10), then
$u(t+\Delta\\!t)=u(t)+\Upsilon(f,\Delta\\!t,u(t))$ for all $t\in\mathbb{R}$.

The setting of autonomous, first-order ODEs is not a strong restriction:
Consider the $d$-dimensional, non-autonomous, $\ell$-th-order ODE

$v^{(\ell)}(t)=g(t,v(t),v^{(1)}(t),\dots,v^{(\ell-1)}(t))\,,$ (15)

where $v^{(k)}$ for $k\in{\left\llbracket 0,\ell\right\rrbracket}$ denotes the
$k$-th derivative of the $\ell$-times differentiable function
$v\colon\mathbb{R}\to\mathbb{R}^{d}$ and
$g\colon\mathbb{R}\times(\mathbb{R}^{d})^{\ell}\to\mathbb{R}^{d}$ is a
function. We represent the derivatives and the time variable with new
variables: $u_{k}=v^{(k-1)}$ for $k\in{\left\llbracket\ell\right\rrbracket}$
and $u_{\ell+1}(t)=t$. For $u(t)=(u_{1}(t),\dots,u_{\ell+1}(t))^{\\!\top}\\!$,
we obtain the ODE

$\begin{pmatrix}\dot{u}_{1}(t)\\\ \vdots\\\ \dot{u}_{\ell-1}(t)\\\
\dot{u}_{\ell}(t)\\\
\dot{u}_{\ell+1}(t)\end{pmatrix}=\begin{pmatrix}u_{2}(t)\\\ \vdots\\\
u_{\ell}(t)\\\ g(u_{\ell+1}(t),u_{1}(t),\dots,u_{\ell}(t))\\\
1\end{pmatrix}\,.$ (16)

It is of the form $\dot{u}(t)=f(u(t))$ for a suitably chosen
$f\colon\mathbb{R}^{\tilde{d}}\to\mathbb{R}^{\tilde{d}}$, where
$\tilde{d}:=d\ell+1$. Hence, (16) is a $\tilde{d}$-dimensional, autonomous,
first-order ODE. If $u$ is a solution to the ODE (16), then $u_{1}$ is a
solution of the ODE (15).

###  2.2 Formal Description of the General Model

We introduce a general model for observations from solutions of an ODE. It is
similar to the one in [Lah+24].

######  Notation 2.2.

  1. (i)

Let $a,b\in\mathbb{Z}$ with $a\leq b$. For $i\in{\left\llbracket
a,b\right\rrbracket}$, let $x_{i}$ be some object. Set $x_{{\left\llbracket
a,b\right\rrbracket}}:=(x_{i})_{i\in{\left\llbracket a,b\right\rrbracket}}$.

  2. (ii)

Let $d\in\mathbb{N}$. For $x\in\mathbb{R}^{d}$, let $\left\|x\right\|_{2}$
denote the Euclidean norm. Let $k\in\mathbb{N}$. For a matrix
$A\in\mathbb{R}^{k\times d}$, denote the operator norm as
$\left\|A\right\|_{\mathsf{op}}:=\sup_{v\in\mathbb{R}^{d},\|v\|_{2}=1}\|Av\|_{2}$.
Let $d_{\mathsf{x}},d_{\mathsf{y}}\in\mathbb{N}$. Let
$f\colon\mathbb{R}^{d_{\mathsf{x}}}\to\mathbb{R}^{d_{\mathsf{y}}}$. Denote the
sup norm of $f$ as
$|f|_{\infty}:=\sup_{x\in\mathbb{R}^{d_{\mathsf{x}}}}\left\|f(x)\right\|_{2}$.
For linear operators, like the derivative $Df(x)$, we abuse the notation
slightly and set
$|Df|_{\infty}:=\sup_{x\in\mathbb{R}^{d_{\mathsf{x}}}}\left\|Df(x)\right\|_{\mathsf{op}}$.
For more details on derivatives on their norms, see appendix A.

  3. (iii)

Let $d_{\mathsf{x}},d_{\mathsf{y}},\beta\in\mathbb{N}$. Denote the set of
$\beta$-times continuously differentiable functions
$f\colon\mathbb{R}^{d_{\mathsf{x}}}\to\mathbb{R}^{d_{\mathsf{y}}}$ as
$\mathcal{D}^{\beta}(\mathbb{R}^{d_{\mathsf{x}}},\mathbb{R}^{d_{\mathsf{y}}})$.
Let $L_{{\left\llbracket
0,\beta\right\rrbracket}}\subseteq\mathbb{R}_{>0}\cup\\{\infty\\}$. Denote by
$\bar{\Sigma}^{d_{\mathsf{x}}\to d_{\mathsf{y}}}(\beta,L_{{\left\llbracket
0,\beta\right\rrbracket}})\subseteq\mathcal{D}^{\beta}(\mathbb{R}^{d_{\mathsf{x}}},\mathbb{R}^{d_{\mathsf{y}}})$
the set of $\beta$-times continuously differentiable functions
$f\colon\mathbb{R}^{d_{\mathsf{x}}}\to\mathbb{R}^{d_{\mathsf{y}}}$ with
$\left|D^{k}f\right|_{\infty}\leq L_{k}$ for $k\in{\left\llbracket
0,\beta\right\rrbracket}$. More details on derivatives and the smoothness
class $\bar{\Sigma}^{d_{\mathsf{x}}\to
d_{\mathsf{y}}}(\beta,L_{{\left\llbracket 0,\beta\right\rrbracket}})$ are
described in appendix A.

Let $d\in\mathbb{N}$ be the dimension of the state space $\mathbb{R}^{d}$. Let
$\beta\in\mathbb{N}$ be the smoothness parameter. Let $L_{{\left\llbracket
0,\beta\right\rrbracket}}\subseteq\mathbb{R}_{>0}$ be the Lipschitz
parameters. Let $\mathcal{F}_{d,\beta}:=\bar{\Sigma}^{d\to
d}(\beta,L_{{\left\llbracket 0,\beta\right\rrbracket}})$ be the smoothness
class. Let the true model function be $f^{\star}\in\mathcal{F}_{d,\beta}$. Let
$m\in\mathbb{N}$ be the number of observed solutions and let
$x_{1},\dots,x_{m}\in\mathbb{R}^{d}$ be their initial conditions. For
$j\in{\left\llbracket m\right\rrbracket}$, let $n_{j}\in\mathbb{N}$ be the
number of observations for the $j$-th solution. Denote the total number of
observations as $n:=\sum_{j=1}^{m}n_{j}$. For $j\in{\left\llbracket
m\right\rrbracket}$, let $T_{j}\in\mathbb{R}_{>0}$ be the observation duration
of the $j$-th solution. For $j\in{\left\llbracket
m\right\rrbracket},i\in{\left\llbracket 0,n_{j}\right\rrbracket}$, let
$t_{j,i}\in\mathbb{R}_{\geq 0}$ be the observation times such that
$0=t_{j,0}\leq t_{j,1}\leq\dots\leq t_{j,n_{j}}=T_{j}$. Let
$\sigma\in\mathbb{R}_{\geq 0}$ be the variance bound of the noise. Let the
noise variables $(\epsilon_{j,i})_{j\in{\left\llbracket
m\right\rrbracket},i\in{\left\llbracket n_{j}\right\rrbracket}}$ be
independent $\mathbb{R}^{d}$-valued random variables such that
$\mathbf{E}[\epsilon_{j,i}]=0$ and
$\mathbf{E}[\|\epsilon_{j,i}\|_{2}^{2}]\leq\sigma^{2}$. For
$j\in{\left\llbracket m\right\rrbracket}$, $i\in{\left\llbracket
n_{j}\right\rrbracket}$, let the observations be

$Y_{j,i}:=U(f^{\star},x_{j},t_{j,i})+\epsilon_{j,i}\,.$ (17)

We observe $Y_{j,i}$ and know $x_{j}$ and $t_{j,i}$, but $f^{\star}$ is
unknown and to be estimated. We assume $d$, $\beta$, $L_{{\left\llbracket
0,\beta\right\rrbracket}}$, and $\sigma$ to be fixed and we are interested in
upper bounds when $n\to\infty$ for the mean squared error for estimators
$\hat{f}$ of $f^{\star}$ on the domain of interest $[0,1]^{d}$.

###  2.3 Need for Restriction of the General Model

We argue that the model introduced in section 2.2 is not suitable for
estimation of the model function in a fixed domain of interest and give first
informal descriptions of two possible remedies of the problem: the Snake model
and the Stubble model.

######  Notation 2.3.

  1. (i)

Let
$(a_{n})_{n\in\mathbb{N}},(b_{n})_{n\in\mathbb{N}}\subseteq\mathbb{R}_{>0}$.
Write

$\displaystyle a_{n}\preccurlyeq b_{n}$ $\displaystyle\ \text{ for
}\limsup_{n\to\infty}a_{n}/b_{n}<\infty\,,$ (18) $\displaystyle a_{n}\prec
b_{n}$ $\displaystyle\ \text{ for }\limsup_{n\to\infty}a_{n}/b_{n}=0\,,$ (19)
$\displaystyle a_{n}\asymp b_{n}$ $\displaystyle\ \text{ for
}a_{n}\preccurlyeq b_{n}\text{ and }b_{n}\preccurlyeq a_{n}\,.$ (20)

Define $\succ$ and $\succcurlyeq$ accordingly.

For the general model, consistent estimation is not possible when considering
the domain of interest $[0,1]^{d}$. We can easily construct settings, where we
will never make any observation in say $[-1/2,1/2]^{d}$, i.e.,
$u(t_{j,i})\not\in[-1/2,1/2]^{d}$ for all $j\in{\left\llbracket
m\right\rrbracket},i\in{\left\llbracket 0,n_{j}\right\rrbracket}$, no matter
how large $n$ is. As we are considering a nonparametric class of model
functions $f^{\star}\in\mathcal{F}_{d,\beta}$, we have for the maximal risk

$\sup_{f\in\mathcal{F}_{d,\beta}}\mathbf{E}_{f^{\star}}\\!\left[\left\|f^{\star}(0)-\hat{f}(0)\right\|_{2}^{2}\right]\succcurlyeq
1\,.$ (21)

This behavior is different from parametric ODE models (finite dimensional
class $\mathcal{F}_{d,\beta}$): In the nonparametric setting, observations
contain only local information (information about $f^{\star}(x)$ for $x$ close
to $u(t_{j,i})$), whereas in a parametric setting, observations typically
contain global information so that consistent estimation is possible under
mild restrictions [Bru08, QZ10, GK12, DK15].

In this work, we introduce two restrictions to the general nonparametric model
that make consistent estimation possible:

In the Stubble model (section 3), we require the initial conditions
$(x_{j})_{j\in{\left\llbracket m\right\rrbracket}}$ to suitably cover the
domain of interest $[0,1]^{d}$. By doing so, we ensure that we obtain
information from every part of the domain of interest. In this model, we
assume to observe many short trajectories, i.e., $m\succ 1$ and
$\max_{j\in{\left\llbracket m\right\rrbracket}}n_{j}\preccurlyeq 1$, which
looks like stubble.

In the Snake model (section 4), we directly require the trajectories
$U(f^{\star},x_{j},[0,T_{j}])$ to cover the domain of interest $[0,1]^{d}$. In
this model, we assume to observe few long trajectories, i.e.,
$\min_{j\in{\left\llbracket m\right\rrbracket}}n_{j}\succ 1$. In the extreme
case, we have $m=1$. Consistent estimation on all of $[0,1]^{d}$ with one
observed solution is only possible, if its trajectory covers $[0,1]^{d}$
suitably. Thus, this trajectory must not intersect itself, as otherwise it
would mean that the solution is periodic. This behavior of self-avoidance in a
bounded domain resembles the video game Snake [Don10, chapter 22].

###  2.4 Standard Regression Model

In this article, we present estimation strategies based on black-box
regression estimators. To formalize what we mean by black-box regression
estimators and to be able to talk about their properties, we introduce some
terminology in this section and give some examples at the end of the section.

######  Notation 2.4.

  1. (i)

Let $(X_{n})_{n\in\mathbb{N}}$ be a sequence of real-valued random variables.
Let $(a_{n})_{n\in\mathbb{N}}\subseteq\mathbb{R}_{>0}$. We use the big O in
probability notation: By $X_{n}\in\mathbf{O}_{\mathbf{P}}(a_{n})$, we mean
that for all $\epsilon\in\mathbb{R}_{>0}$, there are $B\in\mathbb{R}_{>0}$ and
$n_{0}\in\mathbb{N}$ such that

$\forall n\in\mathbb{N}_{\geq
n_{0}}\colon\mathbf{P}\\!\left(\left|\frac{X_{n}}{a_{n}}\right|>B\right)<\epsilon\,.$
(22)

  2. (ii)

Let $k,d\in\mathbb{N}$, $v=(v_{1},\dots,v_{d})\in\mathbb{R}^{d}$. Denote the
projection to the $k$-th dimension as $\Pi_{k}$, i.e., $\Pi_{k}v=v_{k}$.

Let $d_{\mathsf{x}},d_{\mathsf{y}}\in\mathbb{N}$ be the dimensions of the
predictor and the target, respectively. Let $\mathcal{F}$ be a set of
functions of the form
$\mathbb{R}^{d_{\mathsf{x}}}\to\mathbb{R}^{d_{\mathsf{y}}}$, the smoothness
class. Let $f^{\star}\in\mathcal{F}$ be the true regression function. Let
$n\in\mathbb{N}$ be the number of observations. Let $T\in\mathbb{R}_{>0}$ be
the extent of the domain of interest $[0,T]^{d_{\mathsf{x}}}$. For
$i\in{\left\llbracket n\right\rrbracket}$, let
$x_{i}\in[0,T]^{d_{\mathsf{x}}}$. Let $\sigma\in\mathbb{R}_{>0}$ the variance
bound of the noise. For $i\in{\left\llbracket n\right\rrbracket}$, let
$\varepsilon_{i}$ be independent $\mathbb{R}^{d_{\mathsf{y}}}$-valued random
variables with $\mathbf{E}[\epsilon_{i}]=0$ and
$\mathbf{E}[\|\epsilon_{i}\|_{2}^{2}]\leq\sigma^{2}$. We call
$\varepsilon_{i}$ the noise. For $i\in{\left\llbracket n\right\rrbracket}$,
set the observations as

$Y_{i}=f^{\star}(x_{i})+\varepsilon_{i}\,.$ (23)

The standard regression problem is to estimate $f^{\star}$ with an estimator
$\hat{f}$ given the data $\mathcal{D}:=(x_{i},Y_{i})_{i\in{\left\llbracket
n\right\rrbracket}}$. We denote this statistical model as
$\mathfrak{P}^{d_{\mathsf{x}}\to
d_{\mathsf{y}}}(\mathcal{F},T,(x_{i})_{i\in{\left\llbracket
n\right\rrbracket}},\sigma)$.

A regression estimator for the regression function $f^{\star}$ in
$\mathfrak{P}^{d_{\mathsf{x}}\to
d_{\mathsf{y}}}(\mathcal{F},T,(x_{i})_{i\in{\left\llbracket
n\right\rrbracket}},\sigma)$ is any measurable map

$\mathcal{E}^{d_{\mathsf{x}}\to
d_{\mathsf{y}}}(\cdot,\cdot)\colon(\mathbb{R}^{d_{\mathsf{x}}}\times\mathbb{R}^{d_{\mathsf{y}}})^{n}\times\mathbb{R}^{d_{\mathsf{x}}}\to\mathbb{R}^{d_{\mathsf{y}}}\,.$
(24)

Given the data $\mathcal{D}$, it estimates the regression function $f^{\star}$
as $x\mapsto\hat{f}(x):=\mathcal{E}^{d_{\mathsf{x}}\to
d_{\mathsf{y}}}(\mathcal{D},x)$. We define the term regression estimator for
the derivative $Df^{\star}$ of the regression function accordingly.

We say $\mathcal{E}^{d_{\mathsf{x}}\to d_{\mathsf{y}}}$ is the componentwise
estimator $\mathcal{E}^{d_{\mathsf{x}}\to 1}$, when we apply
$\mathcal{E}^{d_{\mathsf{x}}\to 1}$ in each target dimension:

$\mathcal{E}^{d_{\mathsf{x}}\to
d_{\mathsf{y}}}(\mathcal{D},\cdot)\colon\mathbb{R}^{d_{\mathsf{x}}}\to\mathbb{R}^{d_{\mathsf{y}}},\
x\mapsto\left(\mathcal{E}^{d_{\mathsf{x}}\to
1}((x_{i},\Pi_{1}Y_{i})_{i\in{\left\llbracket
n\right\rrbracket}},x),\dots,\mathcal{E}^{d_{\mathsf{x}}\to
1}((x_{i},\Pi_{d_{\mathsf{y}}}Y_{i})_{i\in{\left\llbracket
n\right\rrbracket}},x)\right)\,.$ (25)

The maximal risk for the model
$\mathfrak{P}_{n}=\mathfrak{P}^{d_{\mathsf{x}}\to
d_{\mathsf{y}}}(\mathcal{F},T,(x_{i})_{i\in{\left\llbracket
n\right\rrbracket}},\sigma)$ at a point $x_{0}\in[0,T]^{d_{\mathsf{x}}}$ of an
estimator $\hat{f}=\mathcal{E}^{d_{\mathsf{x}}\to
d_{\mathsf{y}}}(\mathcal{D},\cdot)$ is

$r(\hat{f},\mathfrak{P}_{n},x_{0}):=\sup_{\mathbf{P}\in\mathfrak{P}_{n}}\mathbf{E}_{\mathbf{P}}\\!\left[\left\|\hat{f}(x_{0})-f^{\star}_{\mathbf{P}}(x_{0})\right\|_{2}^{2}\right]\,,$
(26)

where $\mathbf{P}\in\mathfrak{P}_{n}$ is the probability measure that gives
rise to (23) with $f^{\star}=f^{\star}_{\mathbf{P}}$. We say that
$r^{\mathsf{pr},\mathsf{sup}}_{0}(\hat{f},\mathfrak{P}_{n})$ is a sup-norm
error bound in probability of an estimator
$\hat{f}=\mathcal{E}^{d_{\mathsf{x}}\to d_{\mathsf{y}}}(\mathcal{D},\cdot)$
if, for all $\mathbf{P}\in\mathfrak{P}_{n}$, for $n\to\infty$, we have

$\sup_{x\in[0,T]^{d_{\mathsf{x}}}}\left\|\hat{f}(x)-f^{\star}_{\mathbf{P}}(x)\right\|_{2}\in\mathbf{O}_{\mathbf{P}}\\!\left(r_{0}^{\mathsf{pr},\mathsf{sup}}(\hat{f},\mathfrak{P}_{n})\right)\,.$
(27)

We define the sup-norm error bound of the first derivative accordingly and
denote it by
$r_{1}^{\mathsf{pr},\mathsf{sup}}(\widehat{Df},\mathfrak{P}_{n})$.

An example of a regression estimator $\hat{f}$ suitable for the Hölder-type
smoothness class $\mathcal{F}=\Sigma^{d_{\mathsf{x}}\to
d_{\mathsf{y}}}(\beta,L)$ with $\beta,L\in\mathbb{R}_{>0}$ (see Definition
A.2), is the componentwise local polynomial estimator described in appendix B.
Suitable here means that it achieves the lowest values possible for
$r(\hat{f},\mathfrak{P}_{n},x_{0})$ for large enough $n\in\mathbb{N}$, up to a
constant.

Further examples of (nonparametric) regression estimators that may be used in
the general estimation strategies constructed in this article include
orthogonal series estimators [Tsy09, chapter 1.7], wavelets [Ant07], penalized
splines [Xia19], and neural networks [Sch20].

##  3 The Stubble Model and Estimation of the Increment Map

In the Stubble model, we observe many solutions to an ODE with different
initial conditions. For each solution, we have only a few observation. The
initial conditions are known and located so that they cover the domain of
interest.

In this section, we first present an explicit estimation procedure for a
Liptschitz-continuous class of model functions with optimal rate of
convergence of the mean squared error at a point. The estimation procedure is
based on a local constant estimator for the increment map and a linear
interpolation for the model function. In the second part of this section, we
generalize these results: For a general Hölder-smoothness class, we present an
estimation strategy based on an arbitrary multivariate regression estimator
for different increment maps and a univariate polynomial interpolation for the
model function. If the nonparametric estimator enjoys certain optimality
criteria with respect to the standard nonparametric regression problem, it
induces an optimal procedure for nonparametric ODE estimation.

###  3.1 Lipschitz

In this specific instance of the general ODE estimation model, we consider
Lipschitz-continuous functions $f$, initial conditions $x_{j}$ on a uniform
grid, and an estimation procedure for $f^{\star}$ that is based on the local
constant estimator for the increments $\Upsilon(f^{\star},\Delta\\!t,\cdot)$.

####  3.1.1 Model

The following is a restriction of the general model of section 2.2.

Let $d\in\mathbb{N}$. Set $\beta=1$. Let $L_{0},L_{1}\in\mathbb{R}_{>0}$. Set
$\mathcal{F}:=\bar{\Sigma}^{d\to d}(1,L_{0},L_{1})$. Let
$f^{\star}\in\mathcal{F}$ and $n_{0}\in\mathbb{N}$. Set $n:=n_{0}^{d}$. Let
$x_{i}\in[0,1]^{d}$ for $i\in{\left\llbracket n\right\rrbracket}$ form a
uniform grid in $[0,1]^{d}$, i.e.,

$\left\\{x_{1},\dots,x_{n}\right\\}=\left\\{\begin{pmatrix}\frac{k_{1}}{n_{0}},\dots,\frac{k_{d}}{n_{0}}\end{pmatrix}^{\\!\top}\\!\,\kern-1.2pt\mathopen{}\left|\vphantom{\begin{pmatrix}\frac{k_{1}}{n_{0}},\dots,\frac{k_{d}}{n_{0}}\end{pmatrix}^{\\!\top}\\!}\,k_{1},\dots,k_{d}\in{\left\llbracket
n_{0}\right\rrbracket}\right.\mathclose{}\kern-1.2pt\right\\}\,.$ (28)

Let $\sigma\in\mathbb{R}_{\geq 0}$. Let $\epsilon_{j}:=\epsilon_{j,1}$,
$j\in{\left\llbracket n\right\rrbracket}$ be independent
$\mathbb{R}^{d}$-valued random variables such that
$\mathbf{E}[\epsilon_{i}]=0$ and
$\mathbf{E}[\|\epsilon_{i}\|_{2}^{2}]\leq\sigma^{2}$. Let
$\Delta\\!t\in\mathbb{R}_{>0}$. Set

$Y_{j}:=Y_{j,1}:=U(f^{\star},x_{j},\Delta\\!t)+\epsilon_{j}\,.$ (29)

We observe $Y_{j}$ and know $x_{j}$ and $\Delta\\!t$, but $f^{\star}$ is
unknown and to be estimated. We assume $d$, $L_{0},L_{1}$, and $\sigma$ to be
fixed. We are interested in upper bounds for the mean squared error at a point
$x_{0}\in[0,1]^{d}$ depending on the asymptotics of $n$ and $\Delta\\!t$.

####  3.1.2 Estimator

For $f\in\mathcal{F}$, $\Delta\\!t\in\mathbb{R}_{\geq 0}$,
$x\in\mathbb{R}^{d}$, recall the definition of the increment map
$\Upsilon(f,\Delta\\!t,x)=U(f,x,\Delta\\!t)-x$ and set
$\iota^{\star}(x):=\Upsilon(f^{\star},\Delta\\!t,x)$.

Let $\hat{\iota}$ be the componentwise (see section 2.4) local constant
estimator (Nadaraya–Watson) of $\iota^{\star}$ with kernel $K$ using the data
$(x_{j},Y_{j})_{j\in{\left\llbracket n\right\rrbracket}}$ and a bandwidth of
optimal order. See appendix B with $\beta=1,\ell=0,s=0$ for details on the
local constant estimator.

For $x\in\mathbb{R}^{d}$, define the scaled increment estimator of $f^{\star}$
as

$\hat{f}(x):=\frac{\hat{\iota}(x)}{\Delta\\!t}\,.$ (30)

####  3.1.3 Result

######  Assumption 3.1.

  * •

StrictKernel: The kernel $K\colon\mathbb{R}_{\geq 0}\to\mathbb{R}$ is
Lipschitz-continuous, nonnegative, non-trivial ($K\neq 0$), and there is
$C_{\mathsf{ker}}\in\mathbb{R}_{>0}$ such that

$K(x)\leq C_{\mathsf{ker}}\mathds{1}_{[0,1]}(x)$ (31)

for all $x\in\mathbb{R}_{\geq 0}$.

######  Theorem 3.2.

Use the model of section 3.1.1 and the estimator of section 3.1.2. Assume
StrictKernel. Assume $\Delta\\!t\preccurlyeq 1$. Then

$\mathbf{E}\\!\left[\left\|\hat{f}(x_{0})-f^{\star}(x_{0})\right\|_{2}^{2}\right]\preccurlyeq\left(\Delta\\!t^{2}n\right)^{-\frac{2}{2+d}}+\Delta\\!t^{2}\,.$
(32)

######  Remark 3.3.

By [SS24, Corollary 3.11], the error rate in Theorem 3.2 is minimax optimal.

The following corollary minimizes the error bound over $\Delta\\!t$, i.e., it
shows the asymptotic behavior of $\Delta\\!t$ that allows the best estimates
of $f$ for the same amount of data.

######  Corollary 3.4.

Use the model of section 3.2.1 and the estimator of section 3.1.2. Assume
StrictKernel. Assume $\Delta\\!t\asymp n^{-\frac{1}{4+d}}$. Then

$\mathbf{E}\\!\left[\left\|\hat{f}(x_{0})-f^{\star}(x_{0})\right\|_{2}^{2}\right]\preccurlyeq
n^{-\frac{2}{4+d}}\,.$ (33)

####  3.1.4 Proof

######  Lemma 3.5.

Assume $f\in\mathcal{F}$. Then
$\Upsilon(f,\Delta\\!t,\cdot)\in\bar{\Sigma}^{d\to
d}(1,\tilde{L}_{0},\tilde{L}_{1})$ with

$\tilde{L}_{0}=\Delta\\!tL_{0}\,,\quad\tilde{L}_{1}=\exp(L_{1}\Delta\\!t)-1\,.$
(34)

###### Proof.

By the definition of $U$ and the fundamental theorem of calculus,

$U(f,x,t)=x+\int_{0}^{t}f(U(f,x,s))\mathrm{d}s\,.$ (35)

Plugging this into the definition of $\iota$ yields

$\iota(x)=U(f,x,\Delta\\!t)-U(f,x,0)=\int_{0}^{\Delta\\!t}f(U(f,x,s))\mathrm{d}s\,.$
(36)

Thus, as $f$ is uniformly bounded by $L_{0}$,

$\left|\iota(x)\right|\leq\int_{0}^{\Delta\\!t}\left|f(U(f,x,s))\right|\mathrm{d}s\leq\int_{0}^{\Delta\\!t}L_{0}\mathrm{d}s=\Delta\\!tL_{0}\,.$
(37)

This is the value of $\tilde{L}_{0}$. Lemma D.9 yields the value of
$\tilde{L}_{1}$ by

$\left|D\iota\right|_{\infty}\leq\exp(L_{1}t)-1\,.$ (38)

∎

######  Lemma 3.6.

Assume $\Delta\\!t\preccurlyeq 1$. Then, for all $x_{0}\in[0,1]^{d}$,

$\mathbf{E}\\!\left[\left\|\hat{\iota}(x_{0})-\iota^{\star}(x_{0})\right\|_{2}^{2}\right]\preccurlyeq\left(\Delta\\!t^{d}n^{-1}\right)^{\frac{2}{2+d}}\,.$
(39)

###### Proof.

According to Lemma 3.5, $\iota^{\star}$ is Lipschitz-continuous with constant
$\tilde{L}_{1}=\exp(L_{1}\Delta\\!t)-1$. As we consider $L_{1}$ to be constant
and $\Delta\\!t\preccurlyeq 1$, we obtain
$\tilde{L}_{1}\preccurlyeq\Delta\\!t$.

By Corollary B.4 and Proposition B.9, the local constant estimator
$\hat{\iota}$ with optimal bandwidth then fulfills

$\mathbf{E}\\!\left[\left\|\hat{\iota}(x_{0})-\iota^{\star}(x_{0})\right\|_{2}^{2}\right]\preccurlyeq\Delta\\!t^{\frac{2d}{2+d}}n^{-\frac{2}{2+d}}\,.$
(40)

∎

######  Lemma 3.7.

Let $\Delta\\!t\in\mathbb{R}_{>0}$. Assume $f\in\mathcal{F}$. Let
$\iota=\Upsilon(f,\Delta\\!t,\cdot)$. Then

$\sup_{x\in[0,1]^{d}}\left\|\frac{\iota(x)}{\Delta\\!t}-f(x)\right\|_{2}\leq\frac{1}{2}L_{0}L_{1}\Delta\\!t\,.$
(41)

###### Proof.

Let $x\in[0,1]^{d}$. With (36), as $f$ is $L_{1}$-Lipschitz,

$\displaystyle\left\|\frac{\iota(x)}{\Delta\\!t}-f(x)\right\|_{2}$
$\displaystyle=\left\|\frac{1}{\Delta\\!t}\int_{0}^{\Delta\\!t}f(U(f,x,s))-f(x)\mathrm{d}s\right\|_{2}$
(42)
$\displaystyle\leq\frac{L_{1}}{\Delta\\!t}\int_{0}^{\Delta\\!t}\left\|U(f,x,s)-x\right\|_{2}\mathrm{d}s\,.$
(43)

Using $U(f,x,s)=x+\int_{0}^{s}f(U(f,x,r))\mathrm{d}r$, we obtain

$\displaystyle\int_{0}^{\Delta\\!t}\left\|U(f,x,s)-x\right\|_{2}\mathrm{d}s$
$\displaystyle\leq\int_{0}^{\Delta\\!t}\int_{0}^{s}\left\|f(U(f,x,r))\right\|_{2}\mathrm{d}r\mathrm{d}s$
(44)
$\displaystyle\leq\int_{0}^{\Delta\\!t}\int_{0}^{s}L_{0}\mathrm{d}r\mathrm{d}s$
(45) $\displaystyle=\frac{1}{2}L_{0}\Delta\\!t^{2}\,.$ (46)

Combining the above inequalities yields

$\displaystyle\left\|\frac{\iota(x)}{\Delta\\!t}-f(x)\right\|_{2}$
$\displaystyle\leq\frac{1}{2}L_{0}L_{1}\Delta\\!t\,.$ (47)

As this inequality holds for all $x\in[0,1]^{d}$, we have finished the proof.
∎

###### Proof of Theorem 3.2.

Combine the bounds in Lemma 3.6 and Lemma 3.7 to obtain

$\displaystyle\mathbf{E}\\!\left[\left\|\hat{f}(x_{0})-f^{\star}(x_{0})\right\|_{2}^{2}\right]$
$\displaystyle=\mathbf{E}\\!\left[\left\|\frac{\hat{\iota}(x_{0})-\iota^{\star}(x_{0})}{\Delta\\!t}+\frac{\iota^{\star}(x_{0})}{\Delta\\!t}-f^{\star}(x_{0})\right\|_{2}^{2}\right]$
(48)
$\displaystyle\preccurlyeq\frac{1}{\Delta\\!t^{2}}\mathbf{E}\\!\left[\left\|\hat{\iota}(x_{0})-\iota^{\star}(x_{0})\right\|_{2}^{2}\right]+\left\|\frac{\iota^{\star}(x_{0})}{\Delta\\!t}-f^{\star}(x_{0})\right\|_{2}^{2}$
(49)
$\displaystyle\preccurlyeq\Delta\\!t^{-\frac{4}{2+d}}n^{-\frac{2}{2+d}}+\Delta\\!t^{2}\,.$
(50)

∎

###  3.2 General

We generalize the Lipschitz setting of section 3.1 in the following way: We
consider an arbitrary smoothness parameter $\beta\in\mathbb{N}$ of the model
function $f^{\star}$. The initial conditions are not restricted to a uniform
grid, but must still be uniform enough in some sense. We present an estimation
strategy which can be used with an arbitrary regression estimator. If this
estimator achieves the optimal rate of convergence for a standard
nonparametric regression problem, the resulting ODE estimator is also optimal
in a minimax sense.

######  Notation 3.8.

  1. (i)

Let $d\in\mathbb{N}$. Let $\alpha\in\mathbb{N}_{0}^{d}$.

     * •

Denote $\left|\alpha\right|:=\sum_{i=1}^{d}\alpha_{i}$,

     * •

Denote $\alpha!:=\prod_{i=1}^{d}\alpha_{i}!$,

     * •

Let $x\in\mathbb{R}^{d}$. Denote
$x^{\alpha}=\prod_{i=1}^{d}x_{i}^{\alpha_{i}}$.

  2. (ii)

For $\ell,d\in\mathbb{N}$, denote $\mathcal{P}_{d,\ell}$ be the set of all
polynomials of degree at most $\ell$ defined over $\mathbb{R}^{d}$,

$\mathcal{P}_{d,\ell}:=\left\\{p\colon\mathbb{R}^{d}\to\mathbb{R}\,\kern-1.2pt\mathopen{}\left|\vphantom{p\colon\mathbb{R}^{d}\to\mathbb{R}}\,p(x)=\sum_{\alpha\in\mathbb{N}_{0}^{d},\left|\alpha\right|\leq\ell}\beta_{\alpha}x^{\alpha}\,,\
\beta_{\alpha}\in\mathbb{R}\right.\mathclose{}\kern-1.2pt\right\\}\,.$ (51)

  3. (iii)

Let $A$ be a set. Denote the indicator function of $A$ as
$\mathds{1}_{\\!A}\\!\\!\left(\cdot\right)$.

####  3.2.1 Model

The following is a restriction of the general model of section 2.2.

Let $d\in\mathbb{N}$, $\beta\in\mathbb{N}$, and $L_{{\left\llbracket
0,\beta\right\rrbracket}}\subseteq\mathbb{R}_{>0}$. Set the smoothness class
as $\mathcal{F}_{d,\beta}:=\bar{\Sigma}^{d\to d}(\beta,L_{{\left\llbracket
0,\beta\right\rrbracket}})$. Let $f^{\star}\in\mathcal{F}_{d,\beta}$. Let
$m\in\mathbb{N}$ and $x_{1},\dots,x_{m}\in\mathbb{R}^{d}$. Let
$\Delta\\!t\in\mathbb{R}_{>0}$. Set $n_{j}=\beta$ for all
$j\in{\left\llbracket m\right\rrbracket}$. Set $t_{j,i}=t_{i}=i\Delta\\!t$ for
$i\in{\left\llbracket 0,\beta\right\rrbracket}$. Let
$\sigma\in\mathbb{R}_{\geq 0}$. Let $(\epsilon_{j,i})_{j\in{\left\llbracket
m\right\rrbracket},i\in{\left\llbracket\beta\right\rrbracket}}$ be independent
$\mathbb{R}^{d}$-valued random variables such that
$\mathbf{E}[\epsilon_{j,i}]=0$ and
$\mathbf{E}[\|\epsilon_{j,i}\|_{2}^{2}]\leq\sigma^{2}$. For
$j\in{\left\llbracket m\right\rrbracket}$,
$i\in{\left\llbracket\beta\right\rrbracket}$, let

$Y_{j,i}=U(f^{\star},x_{j},i\Delta\\!t)+\epsilon_{j,i}\,.$ (52)

We observe $Y_{j,i}$, and know $x_{j}$ and $\Delta\\!t$, but $f^{\star}$ is
unknown and to be estimated. We assume $d$, $L_{{\left\llbracket
0,\beta\right\rrbracket}}$, and $\sigma$ to be fixed and we are interested in
asymptotic upper bounds for the mean squared error at a point
$x_{0}\in[0,1]^{d}$ depending on $n$ and $\Delta\\!t$.

######  Remark 3.9.

  1. (i)

The results below, are the same (up to a constant) if we allow
$n_{j}\geq\beta$ as long as $\max_{j\in{\left\llbracket
m\right\rrbracket}}n_{j}\preccurlyeq 1$.

####  3.2.2 Estimation

For $i\in{\left\llbracket\beta\right\rrbracket}$, let
$\iota^{\star}_{i}=\Upsilon(f^{\star},t_{i},\cdot)$. Let
$L\in\mathbb{R}_{>0}$. Let $\mathcal{E}^{d\to d}$ be an arbitrary regression
estimator for the standard regression model $\mathfrak{P}:=\mathfrak{P}^{d\to
d}(\Sigma^{d\to d}(\beta,L),1,(x_{j})_{j\in{\left\llbracket
m\right\rrbracket}},\sigma)$, see section 2.4. For
$i\in{\left\llbracket\beta\right\rrbracket}$, estimate the increment maps
$\iota^{\star}_{i}$ by $\hat{\iota}_{i}:=\mathcal{E}^{d\to
d}((x_{j},Y_{j,i}-x_{j})_{j\in{\left\llbracket m\right\rrbracket}},\cdot)$.
For convenience, set $\hat{\iota}_{0}(x):=0$ for all $x\in\mathbb{R}^{d}$.
Denote the maximal risk of the regression estimator, see (26), as,

$r(d,\beta,L,\sigma,(x_{j})_{j\in{\left\llbracket
m\right\rrbracket}}):=\sup_{x\in[0,1]^{d}}r(\mathcal{E}^{d\to
d},\mathfrak{P},x)\,.$ (53)

For $x_{0}\in\mathbb{R}^{d}$, let
$\hat{p}(x_{0},\cdot):=(\hat{p}_{1}(x_{0},\cdot),\dots,\hat{p}_{d}(x_{0},\cdot))$,
where $\hat{p}_{k}(x_{0},\cdot)\in\mathcal{P}_{1,\beta}$ is the (univariate)
polynomial of at most degree $\beta$ that interpolates
$(t_{i},\Pi_{k}\hat{\iota}_{i}(x_{0}))_{i\in{\left\llbracket
0,\beta\right\rrbracket}}$. To estimate $f^{\star}$, we locally approximate
$U(f^{\star},x_{0},t)$ by $x_{0}+\hat{p}(x_{0},t)$, i.e., we consider
$\dot{\hat{p}}(x_{0},t)\approx f^{\star}(x_{0}+\hat{p}(x_{0},t))$. As
$f^{\star}(x_{0})=f^{\star}(x_{0}+\hat{p}(x_{0},0))$, this approach yields the
estimator

$\hat{f}(x_{0}):=\dot{\hat{p}}(x_{0},0)\,.$ (54)

######  Remark 3.10.

  1. (i)

The value of $L$ is to be specified later and depends on $\Delta\\!t$,
$\beta$, and $L_{{\left\llbracket 0,\beta\right\rrbracket}}$. Typically the
estimator $\mathcal{E}^{d\to d}$ does not depend on $L$, but the risk
$r(\mathcal{E}^{d\to d},\mathfrak{P},x)$ does.

  2. (ii)

If we set $\beta=1$ and use the componentwise local constant estimator for
$\mathcal{E}^{d\to d}$, we obtain the procedure of section 3.1.

  3. (iii)

The polynomial interpolation in the second step can be viewed as a
Adams–Bashforth method: a linear multistep method for solving ordinary
differential equations numerically, see [But16, Chapter 24].

####  3.2.3 Result

######  Theorem 3.11.

Use the model of section 3.2.1 and the estimator of section 3.2.2. Assume
$\Delta\\!t\preccurlyeq 1$. Then

$\mathbf{E}\\!\left[\left\|\hat{f}(x_{0})-f^{\star}(x_{0})\right\|_{2}^{2}\right]\preccurlyeq\Delta\\!t^{-2}r(d,\beta,C\Delta\\!t,\sigma,(x_{j})_{j\in{\left\llbracket
m\right\rrbracket}})+\Delta\\!t^{2\beta}\,.$ (55)

for some constant $C\in\mathbb{R}_{>0}$ depending only on $\beta$ and
$L_{{\left\llbracket 0,\beta\right\rrbracket}}$.

The minimax rate for the squared error in the standard regression problem with
smoothness class $\Sigma^{d\to d}(\beta,L)$ is
$L^{\frac{2d}{2\beta+d}}n^{-\frac{2\beta}{2\beta+d}}$ up to some constants. It
is achieved, for example, by the local polynomial estimator with optimal
bandwidth under some conditions on $x_{i}$, which are fulfilled by a uniform
grid. See Corollary B.4. Considering $m\asymp n$, we then have

$r(d,\beta,C\Delta\\!t,\sigma,(x_{j})_{j\in{\left\llbracket
m\right\rrbracket}})\preccurlyeq\Delta\\!t^{\frac{2d}{2\beta+d}}n^{-\frac{2\beta}{2\beta+d}}\,.$
(56)

The conditions for the optimal rate for the local polynomial estimator with
kernel $K$ of appendix B are as follows.

######  Assumption 3.12.

  * •

Eigenvalue: There is $C_{\mathsf{egv}}\in\mathbb{R}_{>0}$ such that

$\lambda_{\mathsf{min}}(B(x))^{-1}\leq C_{\mathsf{egv}}$ (57)

for all $x\in[0,T]^{d}$, where $B(x)$ is defined in (175).

  * •

Cover: There is $C_{\mathsf{cvr}}\in\mathbb{R}_{>0}$ such that

$\frac{1}{n}\sum_{i=1}^{n}\mathds{1}_{\\!\mathrm{B}^{d}(x,r)}\\!\\!\left(x_{i}\right)\leq\max\\!\left(\frac{1}{n},\,C_{\mathsf{cvr}}\left(\frac{r}{T}\right)^{d}\right)\,,$
(58)

for all $x\in[0,T]^{d}$, $r\in\mathbb{R}_{>0}$.

  * •

Kernel: The support of the kernel $K\colon\mathbb{R}_{\geq 0}\to\mathbb{R}$
fulfills $\operatorname{supp}(K)\subseteq[0,1]$. Furthermore, there is
$C_{\mathsf{ker}}\in\mathbb{R}_{>0}$ such that

$K(z)\leq C_{\mathsf{ker}}$ (59)

for all $z\in\mathbb{R}_{\geq 0}$.

The following corollary shows the result of Theorem 3.11 using a minimax
optimal estimator $\mathcal{E}^{d\to d}$.

######  Corollary 3.13.

Use the model of section 3.2.1. Use the estimator of section 3.2.2 with the
componentwise local polynomial estimator of degree $\ell:=\beta-1$ with kernel
$K$ and optimal bandwidth as the estimator $\mathcal{E}^{d\to d}$. Assume
$\Delta\\!t\preccurlyeq 1$. Assume Cover, Eigenvalue, and Kernel. Then, for
all $x_{0}\in[0,1]^{d}$,

$\mathbf{E}\\!\left[\left\|\hat{f}(x_{0})-f^{\star}(x_{0})\right\|_{2}^{2}\right]\preccurlyeq\left(\Delta\\!t^{2}n\right)^{-\frac{2\beta}{2\beta+d}}+\Delta\\!t^{2\beta}\,.$
(60)

######  Remark 3.14.

  1. (i)

Of course, the local polynomials can be replaced by any estimator that
achieves the same (optimal) rate of convergence.

  2. (ii)

By [SS24, Corollary 3.11], the error rate in Corollary 3.13 is minimax
optimal.

The following corollary minimizes the error bound over $\Delta\\!t$, i.e., it
shows the asymptotic behavior of $\Delta\\!t$ that allows the best estimates
of $f$ for the same amount of data.

######  Corollary 3.15.

Use the setting and assumptions of Corollary 3.13. Assume

$\Delta\\!t\asymp n^{-\frac{1}{2(\beta+1)+d}}\,.$ (61)

Then

$\mathbf{E}\\!\left[\left\|\hat{f}(x_{0})-f^{\star}(x_{0})\right\|_{2}^{2}\right]\preccurlyeq
n^{-\frac{2\beta}{2(\beta+1)+d}}\,.$ (62)

####  3.2.4 Proof

For $x_{0}\in\mathbb{R}^{d}$, let
$p^{\star}(x_{0},\cdot)\in\mathcal{P}_{1,\beta}$ be the (univariate)
polynomial of at most degree $\beta$ that interpolates
$(t_{i},\iota^{\star}_{i}(x_{0}))_{i\in{\left\llbracket
0,\beta\right\rrbracket}}$. By the definition of $\hat{f}$, we can write

$\hat{f}(x_{0})-f^{\star}(x_{0})=\dot{\hat{p}}(x_{0},0)-\dot{p}^{\star}(x_{0},0)+\dot{p}^{\star}(x_{0},0)-f^{\star}(x_{0})\,.$
(63)

Using this, we split the error into two parts,

$\mathbf{E}\\!\left[\left\|\hat{f}(x_{0})-f^{\star}(x_{0})\right\|_{2}^{2}\right]\leq
2\mathbf{E}\\!\left[\left\|\dot{\hat{p}}(x_{0},0)-\dot{p}^{\star}(x_{0},0)\right\|_{2}^{2}\right]+2\left\|\dot{p}^{\star}(x_{0},0)-f^{\star}(x_{0})\right\|_{2}^{2}\,.$
(64)

We start with the second part: The polynomial $t\mapsto p^{\star}(x_{0},t)$
interpolates the translated solution
$\bar{u}^{\star}(x_{0},t):=U(f^{\star},x_{0},t)-x_{0}$. By Corollary D.7, we
have $\left|D^{\beta+1}\bar{u}^{\star}(x_{0},\cdot)\right|_{\infty}\leq
C_{L}$, where $C_{L}$ depends only on $\beta$ and $L_{{\left\llbracket
0,\beta\right\rrbracket}}$. By the definition of $\bar{u}^{\star}(x_{0},t)$,
it fulfills $\dot{\bar{u}}^{\star}(x_{0},0)=f(x_{0})$. Thus, Lemma C.2
provides a constant $C_{\beta}$ depending only on $\beta$ such that

$\displaystyle\left\|\dot{p}^{\star}(x_{0},0)-f^{\star}(x_{0})\right\|_{2}$
$\displaystyle=\left\|\dot{p}^{\star}(x_{0},0)-\dot{\bar{u}}^{\star}(x_{0},0)\right\|_{2}$
(65) $\displaystyle\leq C_{\beta}C_{L}\Delta\\!t^{\beta}\,.$ (66)

In other words,

$\left\|\dot{p}^{\star}(x_{0},0)-f^{\star}(x_{0})\right\|_{2}^{2}\preccurlyeq\Delta\\!t^{2\beta}\,.$
(67)

Let us turn our attention to the first part of the right-hand side of (64): By
Corollary D.12, $\left|D^{\beta}\iota^{\star}_{i}\right|_{\infty}\leq
C\Delta\\!t$ for some constant $C$ depending only on $d,\beta$,
$L_{{\left\llbracket 0,\beta\right\rrbracket}}$ and the maximal $\Delta\\!t$
as $f^{\star}\in\bar{\Sigma}^{d\to d}(\beta,L_{{\left\llbracket
0,\beta\right\rrbracket}})$ and $\Delta\\!t\preccurlyeq 1$. As
$\mathcal{E}^{d\to d}$ achieves the rate of convergence
$r(d,\beta,L,\sigma,(x_{j})_{j\in{\left\llbracket m\right\rrbracket}})$ on
$\Sigma^{d\to d}(\beta,L)$, we obtain

$\mathbf{E}\\!\left[\left\|\hat{\iota}_{i}(x_{0})-\iota^{\star}_{i}(x_{0})\right\|_{2}^{2}\right]\preccurlyeq
r(d,\beta,C\Delta\\!t,\sigma,(x_{j})_{j\in{\left\llbracket
m\right\rrbracket}})\,.$ (68)

Applying Lemma C.3 with the bound (68) yields

$\mathbf{E}\\!\left[\left\|\dot{\hat{p}}(x_{0},0)-\dot{p}^{\star}(x_{0},0)\right\|_{2}^{2}\right]\preccurlyeq\Delta\\!t^{-2}r(d,\beta,C\Delta\\!t,\sigma,(x_{j})_{j\in{\left\llbracket
m\right\rrbracket}})\,.$ (69)

We combine (67) and (69) and get

$\mathbf{E}\\!\left[\left\|\hat{f}(x_{0})-f^{\star}(x_{0})\right\|_{2}^{2}\right]\preccurlyeq\Delta\\!t^{-2}r(d,\beta,C\Delta\\!t,\sigma,(x_{j})_{j\in{\left\llbracket
m\right\rrbracket}})+\Delta\\!t^{2\beta}\,.$ (70)

##  4 The Snake Model and Estimation of the Observed Solutions

In the Snake model, we observe one (or a few) solution(s) to an ODE. For each
solution, we have many observation. The associated estimation problem depends
strongly on how the trajectories of the solutions are located in the state
space.

In this section, we first present an explicit estimation procedure for a
Lipschitz-continuous class of model functions and an upper bound on the
maximal risk in sup norm. The estimation procedure is based on a local linear
estimator for the solutions and a nearest neighbor interpolation for the model
function. It is rate-optimal in some settings. In the second part of this
section, we generalize these results: For a general Hölder-smoothness class,
we present an estimation strategy based on a generic nonparametric estimator
for the solutions and a multivariate polynomial interpolation for the model
function. If the nonparametric estimator enjoys certain optimality criteria
with respect to the standard nonparametric regression problem, it induces a
procedure for nonparametric ODE estimation that is optimal in some settings.

###  4.1 Lipschitz

In this specific instance of the general ODE estimation model, we consider
Lipschitz-continuous model functions $f^{\star}$ and an estimation procedure
for $\hat{f}$ that is based on local linear estimators $\hat{u}$ and
$\hat{\dot{u}}$ for the observed solution $U(f^{\star},x_{1},\cdot)$ and its
derivative, and a nearest neighbor interpolation of
$\hat{u}(t)\mapsto\hat{\dot{u}}(t)$ for the model function. We obtain an upper
bound on the expected sup norm error. For this, we either require the
trajectories $U(f^{\star},x_{1},[0,T_{1}])$ to cover the fixed domain of
interest $[0,1]^{d}$ in a suitable way or we adapt our domain of interest to
these trajectories, i.e., the upper bound only holds close to the trajectory.

######  Notation 4.1.

  1. (i)

Let $d\in\mathbb{N}$, $r\geq 0$, and $z\in\mathbb{R}^{d}$. Denote the closed
ball with radius $r$ around $z$ as

$\mathrm{B}^{d}(z,r):=\left\\{x\in\mathbb{R}^{d}\,\kern-1.2pt\mathopen{}\left|\vphantom{x\in\mathbb{R}^{d}}\,\left\|x-z\right\|_{2}\leq
r\right.\mathclose{}\kern-1.2pt\right\\}\,.$ (71)

Let $\mathcal{Z}\subseteq\mathbb{R}^{d}$. Denote the closed ball with radius
$r$ around $\mathcal{Z}$ as

$\mathrm{B}^{d}(\mathcal{Z},r):=\bigcup_{z\in\mathcal{Z}}\mathrm{B}^{d}(z,r)\,.$
(72)

####  4.1.1 Model

The following is a restricting of the general model of section 2.2.

Let $d\in\mathbb{N}_{\geq 2}$ and $L_{0},L_{1}\in\mathbb{R}_{>0}$. Set
$\mathcal{F}:=\bar{\Sigma}^{d\to d}(1,L_{0},L_{1})$. Let
$f^{\star}\in\mathcal{F}$. Set $m=1$. Let $x_{1}\in\mathbb{R}^{d}$,
$n=n_{1}\in\mathbb{N}$, and $\Delta\\!t\in\mathbb{R}_{>0}$. Set
$t_{1,i}:=t_{i}:=i\Delta\\!t$ for $i={\left\llbracket 0,n\right\rrbracket}$
and $T:=T_{1}:=n\Delta\\!t$. Let $\sigma\in\mathbb{R}_{\geq 0}$. Let
$\epsilon_{i}:=\epsilon_{1,i}$, $i\in{\left\llbracket n\right\rrbracket}$ be
independent $\mathbb{R}^{d}$-valued random variables such that
$\mathbf{E}[\epsilon_{i}]=0$ and
$\mathbf{E}[\|\epsilon_{i}\|_{2}^{2}]\leq\sigma^{2}$. Set

$Y_{i}:=Y_{1,i}:=U(f^{\star},x_{1},t_{i})+\epsilon_{i}\qquad\text{for
}i\in{\left\llbracket n\right\rrbracket}\,.$ (73)

We observe $Y_{i}$ and know $x_{1}$ and $\Delta\\!t$, but $f^{\star}$ is
unknown and to be estimated. We assume $d$, $L_{0},L_{1}$, and $\sigma$ to be
fixed. For an estimator $\hat{f}$, we are interested in asymptotic upper
bounds for the mean squared sup-norm of $f^{\star}-\hat{f}$ in the domain of
interest $[0,1]^{d}$ depending on the asymptotics of $n$ and $\Delta\\!t$.

####  4.1.2 Estimator

Let $\hat{u}\colon[0,T]\to\mathbb{R}$ be the componentwise local linear
estimator (as described in Appendix B, (174) with $\ell=1$, $s=0$, $d=1$) of
$U(f^{\star},x_{1},\cdot)$ using the data $(t_{i},Y_{i})_{i\in{\left\llbracket
n\right\rrbracket}}$. Similarly, denote by $\hat{\dot{u}}$ the local linear
estimator of the first derivative of $U(f^{\star},x_{1},\cdot)$ ($\ell=1$,
$s=1$, $d=1$ in terms of Appendix B). We assume that the bandwidth $h$ is
chosen of optimal order for the error bound in sup-norm, see Corollary B.7.

Let $\hat{f}$ be the nearest neighbor interpolation of
$(\hat{u},\hat{\dot{u}})$: For $x\in\mathbb{R}^{d}$ and a continuous function
$u\colon[0,T]\to\mathbb{R}^{d}$, let

$t_{\mathsf{NN}}(u,T,x)\in\operatorname*{arg\,min}_{t\in[0,T]}\|u(t)-x\|_{2}$
(74)

with an arbitrary choice if the minimizer is not unique. Then, for
$x\in\mathbb{R}^{d}$, we define the estimator of $f^{\star}$ as

$\hat{f}(x):=\hat{\dot{u}}\\!\left(t_{\mathsf{NN}}(\hat{u},T,x)\right)\,.$
(75)

####  4.1.3 Result

For our results in sup-norm, we require the noise to be sub-Gaussian.

######  Definition 4.2 (Sub-Gaussian).

A real-valued random variable $Z$, is called sub-Gaussian if there is
$v\in\mathbb{R}_{>0}$ such that

$\mathbf{P}\\!\left(\left|Z\right|\geq t\right)\leq
2\exp\\!\left(-\frac{t^{2}}{2v}\right)$ (76)

for all $t\in\mathbb{R}_{>0}$. Then $v$ is called variance parameter.

######  Assumption 4.3.

  * •

SubGaussian: The noise $\varepsilon_{i}$ is sub-Gaussian with variance
parameter $\sigma^{2}$ in each component.

For $\mathcal{X}\subseteq\mathbb{R}^{d}$ and a function
$u\colon[0,T]\to\mathbb{R}^{d}$, define
$\delta_{\mathsf{NN}}(u,T,\mathcal{X})$ as smallest radius so that the ball
around $u([0,T])$ includes $\mathcal{X}$:

$\delta_{\mathsf{NN}}(u,T,\mathcal{X}):=\sup_{x\in\mathcal{X}}\inf_{t\in[0,T]}\left\|u(t)-x\right\|_{2}\,.$
(77)

The results come in two versions with respect to the domain of interest: We
either fix the domain of interest, e.g., $\mathcal{X}=[0,1]^{d}$ and make the
error bound depend on $\delta_{\mathsf{NN}}(u^{\star},T,\mathcal{X})$, or we
fix $\delta$ (or a sequence $\delta=\delta_{n}$) and adapt our domain of
interest to be $\mathcal{X}=\mathrm{B}^{d}(u^{\star}([0,T]),\delta)$.

######  Theorem 4.4.

Use the model of section 4.1.1 and the estimator of section 4.1.2. Assume
StrictKernel, SubGaussian. Assume

$\left(\frac{n}{\log(n)}\right)^{-\frac{1}{4}}\preccurlyeq T\preccurlyeq
n\log(n)^{\frac{1}{4}}\,.$ (78)

  1. (i)

Let $\delta\in\mathbb{R}_{\geq 0}$, potentially changing with $n$. Let
$\mathcal{X}:=\mathrm{B}^{d}(u^{\star}([0,T]),\delta)$. Then

$\mathbf{E}\\!\left[\sup_{x\in\mathcal{X}}\left\|\hat{f}(x)-f^{\star}(x)\right\|_{2}^{2}\right]\preccurlyeq\delta^{2}+\left(\frac{T\log
n}{n}\right)^{\frac{2}{5}}\,.$ (79)

  2. (ii)

Let $\delta=\delta_{\mathsf{NN}}(u^{\star},T,[0,1]^{d})$. Then

$\displaystyle\mathbf{E}\\!\left[\sup_{x\in[0,1]^{d}}\left\|\hat{f}(x)-f^{\star}(x)\right\|_{2}^{2}\right]\preccurlyeq\delta^{2}+\left(\frac{T\log
n}{n}\right)^{\frac{2}{5}}\,.$ (80)

######  Remark 4.5.

  1. (i)

Comparing the error bound in Theorem 4.4 (ii) with [SS24, Corollary 4.11], we
see that it is optimal if

$\delta^{2}\succcurlyeq\left(\frac{T\log n}{n}\right)^{\frac{1}{5}}\,.$ (81)

  2. (ii)

Note that [SS24, Corollary 4.11] uses $m\succ 1$ many solutions. But in [SS24,
Appendix E], it is argued that $m=1$ suffices.

  3. (iii)

One can think of (81) as the requirement of large distances between temporally
different parts of trajectory. If this is not the case, then the estimator
seems suboptimal: For the estimation of $\hat{u}$ and $\hat{\dot{u}}$, it
ignores information from observations that are distant in time but close in
state space. An estimator that is optimal for regimes where different
trajectory parts are not well-separated, must not ignore this information.

######  Corollary 4.6.

Use the model of section 4.1.1 and the estimator of section 4.1.2. Assume
StrictKernel, SubGaussian.

  1. (i)

Let $\delta\in\mathbb{R}_{\geq 0}$, potentially changing with $n$. Let
$\mathcal{X}:=\mathrm{B}_{\delta}(u^{\star}([0,T]))$. Assume
$T\preccurlyeq\left(\frac{n}{\log n}\right)^{\frac{d-1}{4+d}}$ and
$\delta\preccurlyeq\left(\frac{n}{\log n}\right)^{-\frac{1}{4+d}}$. Then

$\mathbf{E}\\!\left[\sup_{x\in\mathcal{X}}\left\|\hat{f}(x)-f^{\star}(x)\right\|_{2}^{2}\right]\preccurlyeq\left(\frac{n}{\log
n}\right)^{-\frac{2}{4+d}}\,.$ (82)

  2. (ii)

Let $\delta=\delta_{\mathsf{NN}}(u^{\star},T,[0,1]^{d})$. Assume
$T\preccurlyeq\left(\frac{n}{\log n}\right)^{\frac{d-1}{4+d}}$ and
$\delta\preccurlyeq\left(\frac{n}{\log n}\right)^{-\frac{1}{4+d}}$. Then

$\mathbf{E}\\!\left[\sup_{x\in[0,1]^{d}}\left\|\hat{f}(x)-f_{\star}(x)\right\|_{2}^{2}\right]\preccurlyeq\left(\frac{n}{\log
n}\right)^{-\frac{2}{4+d}}\,.$ (83)

######  Remark 4.7.

  1. (i)

The conditions on $T$ and $\delta$ ensure an optimal trade-off between the
ability to reconstruct $u^{\star}$ and coverage of $[0,1]^{d}$. They can be
fulfilled for certain $f$, but it is rather restrictive. Furthermore, to cover
the hypercube $[0,1]^{d}$ with $\mathrm{B}_{\delta}(u([0,T]))$, where $u$ has
a speed bounded by $L$, we require $T\succcurlyeq\delta^{-(d-1)}$. Hence, the
conditions on $\delta$ and $T$ in Theorem 4.4 are tight in the sense that one
variable determines the other up to a constant.

####  4.1.4 Proof

######  Lemma 4.8.

Let $\delta>0$. Let $x_{0}\in\mathrm{B}^{d}(u^{\star}([0,T]),\delta)$. Then

$\displaystyle\left\|\hat{f}(x_{0})-f^{\star}(x_{0})\right\|_{2}\leq
L_{1}\left(\delta+2\sup_{t\in[0,T]}\left\|\hat{u}(t)-u^{\star}(t)\right\|_{2}\right)+\sup_{t\in[0,T]}\left\|\hat{\dot{u}}(t)-\dot{u}^{\star}(t)\right\|_{2}\,.$
(84)

###### Proof.

Denote $u^{\star}:=U(f^{\star},x_{1},\cdot)$,
$\tau^{\star}:=t_{\mathsf{NN}}(u^{\star},T,x_{0})$, and
$\hat{\tau}:=t_{\mathsf{NN}}(\hat{u},T,x_{0})$. Using the triangle inequality
and the definition of $\hat{f}$ and $\dot{u}^{\star}$, we obtain

$\displaystyle\left\|\hat{f}(x_{0})-f^{\star}(x_{0})\right\|_{2}$
$\displaystyle\leq\left\|\hat{f}(x_{0})-f^{\star}(u^{\star}(\hat{\tau}))\right\|_{2}+\left\|f^{\star}(u^{\star}(\hat{\tau}))-f^{\star}(x_{0})\right\|_{2}$
(85)
$\displaystyle=\left\|\hat{\dot{u}}(\hat{\tau})-\dot{u}^{\star}(\hat{\tau})\right\|_{2}+\left\|f^{\star}(u^{\star}(\hat{\tau}))-f^{\star}(x_{0})\right\|_{2}\,.$
(86)

Let us consider the second term of the last line. As $f^{\star}$ is
$L_{1}$-Lipschitz, we have

$\left\|f^{\star}(u^{\star}(\hat{\tau}))-f^{\star}(x_{0})\right\|_{2}\leq
L_{1}\left\|u^{\star}(\hat{\tau})-x_{0}\right\|_{2}\,.$ (87)

We use the triangle inequality in
$\left\|x_{0}-u^{\star}(\hat{\tau})\right\|_{2}\leq\left\|x_{0}-\hat{u}(\hat{\tau})\right\|_{2}+\left\|\hat{u}(\hat{\tau})-u^{\star}(\hat{\tau})\right\|_{2}$
and bound

$\displaystyle\left\|x_{0}-\hat{u}(\hat{\tau})\right\|_{2}$
$\displaystyle\leq\left\|x_{0}-\hat{u}(\tau^{\star})\right\|_{2}$ (88)
$\displaystyle\leq\left\|x_{0}-u^{\star}(\tau^{\star})\right\|_{2}+\left\|u^{\star}(\tau^{\star})-\hat{u}(\tau^{\star})\right\|_{2}\,,$
(89)

where we used the minimizing property of $\hat{\tau}$ and again the triangle
inequality. Putting all previous bounds together, we obtain

$\displaystyle\left\|\hat{f}(x_{0})-f^{\star}(x_{0})\right\|_{2}$ (90)
$\displaystyle\leq
L_{1}\left(\left\|x_{0}-u^{\star}(\tau^{\star})\right\|_{2}+\left\|u^{\star}(\tau^{\star})-\hat{u}(\tau^{\star})\right\|_{2}+\left\|\hat{u}(\hat{\tau})-u^{\star}(\hat{\tau})\right\|_{2}\right)+\left\|\hat{\dot{u}}(\hat{\tau})-\dot{u}^{\star}(\hat{\tau})\right\|_{2}\,.$
(91)

The claim follows by taking supremum over $\tau^{\star},\hat{\tau}\in[0,T]$
and noting the definition of $x_{0}$ and $\tau^{\star}$. ∎

###### Proof of Theorem 4.4.

As $f^{\star}\in\bar{\Sigma}(1,L_{0},L_{1})$, we have
$u^{\star}\in\bar{\Sigma}(2,\infty,L_{0},\tilde{L}_{1})$, where
$\tilde{L}_{1}$ depends on $L_{0}$ and $L_{1}$, by Corollary D.7. We want to
apply Corollary B.7 to obtain bounds on the sup-norm of the error for
$\hat{u}$ and $\hat{\dot{u}}$. StrictKernel together with the uniform grid
$t_{i}=i\Delta\\!t$ implies Kernel, Eigenvalue, and Cover by Proposition B.9.
SubGaussian is assumed and (78) implies (184). Thus, Corollary B.7 yields

$\displaystyle\mathbf{E}\\!\left[\sup_{t\in[0,T]}\left\|\hat{u}(t)-u^{\star}(t)\right\|_{2}^{2}\right]$
$\displaystyle\preccurlyeq\left(\frac{T\log(n)}{n}\right)^{\frac{4}{5}}\,,$
(92)
$\displaystyle\mathbf{E}\\!\left[\sup_{t\in[0,T]}\left\|\hat{\dot{u}}(t)-\dot{u}^{\star}(t)\right\|_{2}^{2}\right]$
$\displaystyle\preccurlyeq\left(\frac{T\log(n)}{n}\right)^{\frac{2}{5}}\,.$
(93)

Together with Lemma 4.8, we obtain

$\displaystyle\mathbf{E}\\!\left[\sup_{x\in[0,1]^{d}}\left\|\hat{f}(x)-f^{\star}(x)\right\|_{2}^{2}\right]$
$\displaystyle\preccurlyeq\delta^{2}+\mathbf{E}\\!\left[\sup_{t\in[0,T]}\left\|\hat{u}(t)-u^{\star}(t)\right\|_{2}^{2}\right]+\mathbf{E}\\!\left[\sup_{t\in[0,T]}\left\|\hat{\dot{u}}(t)-\dot{u}^{\star}(t)\right\|_{2}^{2}\right]$
(94)
$\displaystyle\preccurlyeq\delta^{2}+\left(\frac{T\log(n)}{n}\right)^{\frac{4}{5}}+\left(\frac{T\log(n)}{n}\right)^{\frac{2}{5}}$
(95)
$\displaystyle\preccurlyeq\delta^{2}+\left(\frac{T\log(n)}{n}\right)^{\frac{2}{5}}\,.$
(96)

∎

###  4.2 General

We now want to consider estimation under higher order smoothness. We
generalize the estimation procedure of section 4.1.2 as follows: We replace
the local linear estimation of the solution $u^{\star}$ and its derivative
$\dot{u}^{\star}$ by arbitrary regression estimators $\mathcal{E}^{1\to d}$
and $\widetilde{\mathcal{E}}^{1\to d}$, respectively. E.g., this could be
local polynomial estimators of the appropriate degree. Furthermore, the
nearest neighbor interpolation is replaced by multivariate polynomial
interpolation of the appropriate degree.

######  Notation 4.9.

  1. (i)

Let $A\subseteq\mathbb{R}^{d}$. Denote the diameter of $A$ as

$\operatorname{\mathsf{diam}}(A):=\sup_{a,a^{\prime}\in
A}\left\|a-a^{\prime}\right\|_{2}\,.$ (97)

  2. (ii)

Let $A\subseteq\mathbb{R}^{d}$. Denote the interior of the convex hull of $A$
as

$\mathring{\mathsf{ch}}(A):=\left\\{\sum_{k=1}^{K}w_{i}a_{i}\,\kern-1.2pt\mathopen{}\left|\vphantom{\sum_{k=1}^{K}w_{i}a_{i}}\,K\in\mathbb{N},\
a_{1},\dots,a_{K}\in A,\
w_{1},\dots,w_{K}\in(0,1),\sum_{k=1}^{K}w_{i}=1\right.\mathclose{}\kern-1.2pt\right\\}\,.$
(98)

####  4.2.1 Model

The following is a restricting of the general model of section 2.2.

Let $d\in\mathbb{N}_{\geq 2}$, $\beta\in\mathbb{N}_{\geq 2}$, and
$L_{{\left\llbracket 0,\beta\right\rrbracket}}\subseteq\mathbb{R}_{>0}$. Set
$\mathcal{F}_{d,\beta}:=\bar{\Sigma}^{d\to d}(\beta,L_{{\left\llbracket
0,\beta\right\rrbracket}})$. Let $f^{\star}\in\mathcal{F}_{d,\beta}$. Set
$m=1$. Let $x_{1}\in\mathbb{R}^{d}$, $n=n_{1}\in\mathbb{N}$, and
$T:=T_{1}\in\mathbb{R}_{>0}$. Set $t_{1,i}:=t_{i}\in[0,T]$ for
$i\in{\left\llbracket 0,n\right\rrbracket}$ with $0=t_{0}\leq\dots\leq
t_{n}=T$. Let $\sigma\in\mathbb{R}_{\geq 0}$. Let
$\epsilon_{i}:=\epsilon_{1,i}$, $i\in{\left\llbracket n\right\rrbracket}$ be
independent $\mathbb{R}^{d}$-valued random variables such that
$\mathbf{E}[\epsilon_{i}]=0$ and
$\mathbf{E}[\|\epsilon_{i}\|_{2}^{2}]\leq\sigma^{2}$. Set

$Y_{i}:=Y_{1,i}:=U(f^{\star},x_{1},t_{i})+\epsilon_{i}\qquad\text{for
}i\in{\left\llbracket n\right\rrbracket}\,.$ (99)

We observe $Y_{i}$ and know $x_{1}$ and $t_{i}$, but $f^{\star}$ is unknown
and to be estimated. We assume $d$, $L_{{\left\llbracket
0,\beta\right\rrbracket}}$, and $\sigma$ to be fixed. For an estimator
$\hat{f}$, we are interested in asymptotic upper bounds for the mean squared
sup-norm of $f^{\star}-\hat{f}$ in the domain of interest $[0,1]^{d}$
depending on the asymptotics of $n$ and $T$.

######  Remark 4.10.

  1. (i)

It does not matter whether $x_{1}$ is known or not: Below, we derive upper
bounds on the estimation error for the potentially more difficult problem of
unknown $x_{1}$. The lower bounds from [SS24, Corollary 4.11] are formulated
with known $x_{1}$, which are then also lower bounds for the case of of
unknown $x_{1}$.

  2. (ii)

This model and following estimation procedure can easily be generalized to
larger $m$ as long as $\min_{j\in{\left\llbracket
m\right\rrbracket}}n_{j}\succ 1$ and the $t_{j,i}$ behave in an appropriate
way.

####  4.2.2 Multivariate Polynomial Interpolation

The estimation strategy for the model described in section 4.2.1 uses
multivariate polynomial interpolation. In this section, we introduce some
basics of polynomial interpolation and some further objects that are required
to describe the ODE estimator. Recall 3.8.

######  Definition 4.11 (Polynomial interpolation).

Let $\ell,d_{\mathsf{x}},d_{\mathsf{y}}\in\mathbb{N}$.

  1. (i)

Set

$N:=N_{{d_{\mathsf{x}}},\ell}:=\dim(\mathcal{P}_{d_{\mathsf{x}},\ell})=\begin{pmatrix}\ell+{d_{\mathsf{x}}}\\\
{d_{\mathsf{x}}}\end{pmatrix}\,.$ (100)

  2. (ii)

For $x\in\mathbb{R}^{d_{\mathsf{x}}}$, denote the vector of monomials of
degree at most $\ell$ in $d_{\mathsf{x}}$ dimensions as

$\psi(x):=\psi_{\ell}(x):=(x^{\alpha})_{|\alpha|\leq\ell}\in\mathbb{R}^{N}\,.$
(101)

For $M\in\mathbb{N}$,
$\mathbf{x}=(x_{1},\dots,x_{M})\in(\mathbb{R}^{d_{\mathsf{x}}})^{M}$, denote
$\Psi(\mathbf{x}):=\Psi_{\ell}(\mathbf{x}):=\left(\psi(x_{1}),\dots,\psi(x_{M})\right)^{\\!\top}\\!\in\mathbb{R}^{M\times
N}$.

  3. (iii)

For $x\in\mathbb{R}^{d_{\mathsf{x}}}$,
$\mathbf{x}\in(\mathbb{R}^{d_{\mathsf{x}}})^{N}$ with
$\Psi(\mathbf{x})\in\mathbb{R}^{N\times N}$,
$\mathbf{y}\in(\mathbb{R}^{d_{\mathsf{y}}})^{N}=\mathbb{R}^{N\times{d_{\mathsf{y}}}}$,
denote by $I:=I_{\ell}$ the componentwise polynomial interpolation,

$I_{\ell}(\mathbf{x},\mathbf{y},\cdot)\colon\mathbb{R}^{d_{\mathsf{x}}}\to\mathbb{R}^{d_{\mathsf{y}}},\,x\mapsto\begin{cases}\psi(x)^{\\!\top}\\!\Psi(\mathbf{x})^{-1}\mathbf{y}&\text{
if }\Psi(\mathbf{x})\text{ is invertible,}\\\ 0&\text{ otherwise.}\end{cases}$
(102)

######  Remark 4.12.

  1. (i)

Let $\mathbf{x}=(x_{1},\dots,x_{N})\in(\mathbb{R}^{d_{\mathsf{x}}})^{N}$ and
$\mathbf{y}=(y_{1},\dots,y_{N})\in(\mathbb{R}^{d_{\mathsf{y}}})^{N}$. The
function $I(\mathbf{x},\mathbf{y},\cdot)$ is an interpolation as
$I(\mathbf{x},\mathbf{y},x_{k})=y_{k}$. It is polynomial as
$x\mapsto\Pi_{k}I(\mathbf{x},\mathbf{y},x)\in\mathcal{P}_{d_{\mathsf{x}},\ell}$
for each $k\in{\left\llbracket d_{\mathsf{y}}\right\rrbracket}$. It is
componentwise as

$I(\mathbf{x},\mathbf{y},x)=\left(I\\!\left(\mathbf{x},\left(\Pi_{k}y_{1},\dots\Pi_{k}y_{N}\right),x\right)\right)_{k\in{\left\llbracket
d_{\mathsf{y}}\right\rrbracket}}\,.$ (103)

  2. (ii)

The definition of $I(\mathbf{x},\mathbf{y},\cdot)=0$ if $\Psi(\mathbf{x})$ is
not invertible is for convenience only. The value 0 is not of importance.

  3. (iii)

If $d_{\mathsf{x}}=1$, the matrix $\Psi(\mathbf{x})\in\mathbb{R}^{N\times N}$
is also called Vandermonde matrix. It is invertible if and only if all
$x_{k}\in\mathbb{R}$ are distinct. In the multivariate case, conditions for
invertibility of $\Psi(\mathbf{x})$ (and hence, unique existence of the
polynomial interpolation) are more complex, see [Olv06, Theorem 4.1].

######  Definition 4.13 (Normalization).

Let $d,N\in\mathbb{N}$. For $\mathbf{x}\in(\mathbb{R}^{d})^{N}$ with
$\operatorname{\mathsf{diam}}(\mathbf{x})>0$, let

$\eta_{\mathbf{x}}\colon\mathbb{R}^{d}\to\mathbb{R}^{d}\,,x\mapsto\frac{x-\frac{1}{N}\sum_{k=1}^{N}x_{k}}{\operatorname{\mathsf{diam}}(\mathbf{x})}\,.$
(104)

Furthermore, denote
$\eta_{\mathbf{x}}(\tilde{\mathbf{x}})=(\eta_{\mathbf{x}}(\tilde{x}_{1}),\dots,\eta_{\mathbf{x}}(\tilde{x}_{M}))$
for
$\tilde{\mathbf{x}}=(\tilde{x}_{1},\dots,\tilde{x}_{M})\in(\mathbb{R}^{d})^{M}$,
$M\in\mathbb{N}$.

######  Definition 4.14.

Let $d,\ell\in\mathbb{N}$, $N=\dim(\mathcal{P}_{d,\ell})$. Let
$\mathcal{X}\subseteq\mathbb{R}^{d}$ and $s,\delta\in\mathbb{R}_{\geq 0}$.

  1. (i)

Define $\mathcal{G}_{\ell,\delta,s}(\mathcal{X})$ as the set of all points
$x\in\mathbb{R}^{d}$ with following property: There is a set of points
$\mathbf{x}\in\mathcal{X}^{N}$ with $x\in\mathring{\mathsf{ch}}(\mathbf{x})$
such that
$\left\|\Psi(\eta_{\mathbf{x}}(\mathbf{x}))^{-1}\right\|_{\mathsf{op}}<s$ and
$\operatorname{\mathsf{diam}}(\mathbf{x})\leq\delta$.

  2. (ii)

Define $\mathcal{H}_{\ell,s}(\mathcal{X},x)$ as the set of all sets
$\mathbf{x}\in\mathcal{X}^{N}$ such that
$x\in\mathring{\mathsf{ch}}(\mathbf{x})$ and
$\left\|\Psi(\eta_{\mathbf{x}}(\mathbf{x}))^{-1}\right\|_{\mathsf{op}}\leq s$.

  3. (iii)

Let $\chi_{\ell,s}(\mathcal{X},x)$ be an element
$\mathbf{x}\in\mathcal{H}_{\ell,s}(\mathcal{X},x)$ such that

$\left\|\Psi(\eta_{\mathbf{x}_{0}}(\mathbf{x}_{0}))^{-1}\right\|_{\mathsf{op}}\operatorname{\mathsf{diam}}(\mathbf{x}_{0})^{\ell+1}\leq(1+\epsilon)\min_{\mathbf{x}\in\mathcal{H}_{\ell,s}(\mathcal{X},x)}\left\|\Psi(\eta_{\mathbf{x}}(\mathbf{x}))^{-1}\right\|_{\mathsf{op}}\operatorname{\mathsf{diam}}(\mathbf{x})^{\ell+1}$
(105)

for $\mathbf{x}_{0}=\chi_{\ell,s}(\mathcal{X},x)$ and a fixed arbitrarily
small $\epsilon\in\mathbb{R}_{>0}$.

####  4.2.3 Estimator

Let $L\in\mathbb{R}_{>0}$. Let $\mathcal{E}^{1\to d}$ and
$\widetilde{\mathcal{E}}^{1\to d}$ be arbitrary regression estimator for the
regression function and its derivative in the standard regression model
$\mathfrak{P}:=\mathfrak{P}^{1\to d}(\Sigma^{1\to
d}(\beta+1,L),T,(t_{i})_{i\in{\left\llbracket n\right\rrbracket}},\sigma)$,
see section 2.4. With these estimators, we estimate the observed ODE-solution
and its derivative,

$\hat{u}(t):=\mathcal{E}^{1\to d}\\!\left((t_{i},Y_{i})_{i\in{\left\llbracket
n\right\rrbracket}},t\right)\qquad\text{and}\qquad\hat{\dot{u}}(t):=\widetilde{\mathcal{E}}^{1\to
d}\\!\left((t_{i},Y_{i})_{i\in{\left\llbracket
n\right\rrbracket}},t\right)\,.$ (106)

Denote the maximal sup-norm risk in probability of the regression estimators,
see (27), as,

$\displaystyle\Gamma(n,\beta+1,L,T)$
$\displaystyle:=r^{\mathsf{pr},\mathsf{sup}}_{0}(\mathcal{E}^{1\to
d},\mathfrak{P})\,,$ (107) $\displaystyle\Lambda(n,\beta+1,L,T)$
$\displaystyle:=r^{\mathsf{pr},\mathsf{sup}}_{1}(\widetilde{\mathcal{E}}^{1\to
d},\mathfrak{P})\,.$ (108)

Set $\ell:=\beta-1$. Fix $s\in\mathbb{R}_{>0}$. Let $x\in\mathbb{R}^{d}$.
Assume $\chi_{\ell,s}(\hat{u}([0,T]),x)$ exists. Let
$\hat{\mathbf{x}}(x):=(\hat{u}(\hat{\tau}_{1}),\dots,\hat{u}(\hat{\tau}_{N}))$,
where $\hat{\tau}_{i}$ is chosen so that
$\mathbf{\hat{x}}(x)=\chi_{\ell,s}(\hat{u}([0,T]),x)$. Let
$\hat{\mathbf{y}}(x):=(\hat{\dot{u}}(\hat{\tau}_{1}),\dots,\hat{\dot{u}}(\hat{\tau}_{N}))$.
Define the estimator

$\hat{f}(x):=I_{\ell}(\mathbf{\hat{x}}(x),\hat{\mathbf{y}}(x),x)\,.$ (109)

######  Remark 4.15.

  1. (i)

The estimator is only defined for
$x\in\mathcal{G}_{\ell,\infty,s}(\hat{u}([0,T]))$, which is a subset of the
convex hull of $\hat{u}([0,T])$. If this is not the case, a fallback such as
nearest neighbor (see section 4.1) could be applied in practice.

  2. (ii)

The variable $s$ is a regularity parameter for the matrix inversion executed
in the multivariate polynomial interpolation. It will be assumed to fulfill
$s\asymp 1$ asymptotically.

####  4.2.4 Results

For $\mathcal{X}\subseteq\mathbb{R}^{d}$ and a function
$u\colon[0,T]\to\mathbb{R}^{d}$, define

$\delta_{\mathsf{NN}}^{\ell,s}(u,T,\mathcal{X}):=\inf\mathopen{}\left\\{\delta^{\prime}\geq
0\colon\mathcal{X}\subseteq\mathcal{G}_{\ell,\delta^{\prime},s}(u([0,T]))\right\\}\mathclose{}\,.$
(110)

######  Theorem 4.16.

Use the model of section 4.2.1 and the estimator of section 4.2.3. Let
$\delta\in\mathbb{R}_{\geq 0}$, potentially changing with $n$. Assume $T\prec
n$, $\Gamma\prec 1$, $\Gamma\preccurlyeq\delta$,
$\Gamma\preccurlyeq\Lambda\preccurlyeq 1$, $s\asymp 1$. Set
$L:=\beta!\sup_{k\in{\left\llbracket\beta\right\rrbracket}}{L_{k}}^{\beta+1}$.

  1. (i)

Set $\mathcal{X}:=\mathcal{G}_{\ell,\delta,s}(u^{\star}([0,T]))$. Then

$\sup_{x\in\mathcal{X}}\left\|\hat{f}(x)-f(x)\right\|_{2}^{2}\in\mathbf{O}_{\mathbf{P}}\\!\left(\delta^{2\beta}+\Lambda(n,\beta+1,L,T)^{2}\right)\,.$
(111)

  2. (ii)

Set $\delta:=\delta_{\mathsf{NN}}^{\ell,s}(u^{\star},T,[0,1]^{d})$. Then

$\sup_{x\in[0,1]^{d}}\left\|\hat{f}(x)-f(x)\right\|_{2}^{2}\in\mathbf{O}_{\mathbf{P}}\\!\left(\delta^{2\beta}+\Lambda(n,\beta+1,L,T)^{2}\right)\,.$
(112)

We want to use the componentwise local polynomial estimators of appendix B as
the estimators $\mathcal{E}^{1\to d}$ and $\widetilde{\mathcal{E}}^{1\to d}$
(appendix B: $d=1$, $\ell=\beta$, and $s=0$ and $s=1$, respectively). We
assume Kernel, Eigenvalue, and Cover for the observation times $t_{i}$ and
SubGaussian for the noise $\varepsilon_{i}$. Furthermore, we require (184).
Then Corollary B.7 yields

$\displaystyle\Gamma(n,\beta+1,L,T)$
$\displaystyle\preccurlyeq\left(\frac{T\log
n}{n}\right)^{\frac{\beta+1}{2(\beta+1)+1}}$ (113)
$\displaystyle\Lambda(n,\beta+1,L,T)$
$\displaystyle\preccurlyeq\left(\frac{T\log
n}{n}\right)^{\frac{\beta}{2(\beta+1)+1}}\,,$ (114)

assuming $L\asymp 1$.

######  Corollary 4.17.

Use the model of section 4.2.1 and the estimator of section 4.2.3 with the
componentwise local polynomial estimator of degree $\beta$ (see appendix B) as
$\mathcal{E}^{1\to d}$ and $\widetilde{\mathcal{E}}^{1\to d}$. Let
$\delta\in\mathbb{R}_{\geq 0}$, potentially changing with $n$ and assume

$\left(\frac{T\log
n}{n}\right)^{\frac{\beta+1}{2(\beta+1)+1}}\preccurlyeq\delta\,.$ (115)

Assume Kernel, Eigenvalue, Cover, where the symbols $d$ and $x_{i}$ in the
assumptions are $d=1$ and $x_{{\left\llbracket
1,n\right\rrbracket}}=t_{{\left\llbracket 1,n\right\rrbracket}}$. Assume
SubGaussian, $s\asymp 1$, and

$\left(\frac{\log(n)}{n}\right)^{\frac{1}{2(\beta+1)}}\preccurlyeq T\prec
n\,.$ (116)

  1. (i)

Set $\mathcal{X}:=\mathcal{G}_{\ell,\delta,s}(u^{\star}([0,T]))$. Then

$\sup_{x\in\mathcal{X}}\left\|\hat{f}(x)-f(x)\right\|_{2}^{2}\in\mathbf{O}_{\mathbf{P}}\\!\left(\delta^{2\beta}+\left(\frac{T\log
n}{n}\right)^{\frac{2\beta}{2(\beta+1)+1}}\right)\,.$ (117)

  2. (ii)

Set $\delta:=\delta_{\mathsf{NN}}^{\ell,s}(u^{\star},T,[0,1]^{d})$. Then

$\sup_{x\in[0,1]^{d}}\left\|\hat{f}(x)-f(x)\right\|_{2}^{2}\in\mathbf{O}_{\mathbf{P}}\\!\left(\delta^{2\beta}+\left(\frac{T\log
n}{n}\right)^{\frac{2\beta}{2(\beta+1)+1}}\right)\,.$ (118)

######  Remark 4.18.

  1. (i)

From [SS24, Corollary 4.11], we obtain

$\mathbf{E}\\!\left[\sup_{x\in[0,1]^{d}}\left\|\hat{f}(x)-f(x)\right\|_{2}^{2}\right]\succcurlyeq\delta^{2\beta}+\left(\delta^{d-1}\frac{T\log(n)}{n}\right)^{\frac{2\beta}{2(\beta+1)+d}}\,.$
(119)

if $\log(n)\asymp\log(\delta^{-(d-1)}nT^{-1})$ and

$\max\\!\left(\left(\frac{T}{n}\right)^{2\beta+d+1},T^{-1}\right)\preccurlyeq\delta^{d-1}\preccurlyeq\min\\!\left(1,\frac{n}{T}\right)\,.$
(120)

Thus, the rate given in (118) is minimax optimal if $\delta^{2\beta}$ is the
dominating term the error bound, i.e., if

$\delta\succcurlyeq\left(\frac{T\log n}{n}\right)^{\frac{1}{2(\beta+1)+1}}\,.$
(121)

This condition is fulfilled in the setting of Corollary 4.19 below. See Remark
4.5 for a discussion of the condition for $\beta=1$.

  2. (ii)

In section 4.1 can extrapolate beyond the convex hull of $\hat{u}([0,T])$ in
the case of $\beta=1$ using nearest neighbors. Whereas our general results are
only available for interpolation, i.e., for
$x\in\mathring{\mathsf{ch}}(\hat{u}([0,T]))$, to be precise, only for the even
more restrictive assumption
$x\in\mathcal{G}_{\ell,\delta,s}(u^{\star}([0,T]))$. Extrapolation seems also
possible using polynomials of degree $>0$, but the technicalities seem more
difficult.

  3. (iii)

In practice, finding $\chi_{\ell,s}$ may be computationally demanding.
Furthermore, using only results of interpolation may be inconveniently
restrictive. Thus, one may want to replace the polynomial interpolation step
by a polynomial regression step on the $k\geq N$ nearest neighbors of
$\hat{u}(\tau_{j})$ for some chosen time points
$\tau_{1},\dots,\tau_{J}\in[0,T]$. Another alternative to polynomials would be
Gaussian process interpolation (or regression) [RW06], potentially also
restricted to $k$ nearest neighbors for performance reasons.

The lowest upper bound on the error is obtained if the two terms in (118) are
balanced. Additionally, we require $T\delta^{d-1}\succcurlyeq 1$ to be able to
cover the domain of interest with
$\mathrm{B}^{d}(U(f^{\star},x_{1},[0,T]),\delta)$, which is necessary for
$[0,1]^{d}\subseteq\mathcal{G}_{\ell,\delta,s}(u^{\star}([0,T]))$.

######  Corollary 4.19.

Use the setting and assumptions of Corollary 4.17. Set Set
$\delta:=\delta_{\mathsf{NN}}^{\ell,s}(u^{\star},T,[0,1]^{d})$. Assume
$T\preccurlyeq\left(\frac{n}{\log n}\right)^{\frac{d-1}{2(\beta+1)+d}}$ and
$\delta\preccurlyeq\left(\frac{n}{\log n}\right)^{-\frac{1}{2(\beta+1)+d}}$.
Then

$\sup_{x\in[0,1]^{d}}\left\|\hat{f}(x_{0})-f(x_{0})\right\|_{2}^{2}\in\mathbf{O}_{\mathbf{P}}\\!\left(\left(\frac{n}{\log
n}\right)^{-\frac{2\beta}{2(\beta+1)+d}}\right)\,.$ (122)

####  4.2.5 Proof

As matrix inversion is continuous for the set of invertible matrices, we have
following lemma.

######  Lemma 4.20.

Let $d\in\mathbb{N}$. For $n\in\mathbb{N}$, let
$A_{n},B_{n}\in\mathbb{R}^{d\times d}$. Assume that $A_{n}$ is invertible for
$n$ large enough. Assume $\left\|A_{n}-B_{n}\right\|_{\mathsf{op}}\prec 1$ and
$\|A_{n}^{-1}\|_{\mathsf{op}}\preccurlyeq 1$. Then $B_{n}$ is invertible for
$n$ large enough and $\left\|A_{n}^{-1}-B_{n}^{-1}\right\|_{\mathsf{op}}\prec
1$.

######  Lemma 4.21.

Let $\delta,s>0$. Let
$x_{0}\in\mathcal{G}_{\beta,\delta,s}(u^{\star}([0,T]))$. Let

$\displaystyle\Gamma$
$\displaystyle:=\sup_{t\in[0,T]}\left\|\hat{u}(t)-u^{\star}(t)\right\|_{2}\,,$
$\displaystyle\Lambda$
$\displaystyle:=\sup_{t\in[0,T]}\left\|\hat{\dot{u}}(t)-\dot{u}^{\star}(t)\right\|_{2}\,.$
(123)

Assume $\Gamma\prec 1$, $\Gamma\preccurlyeq\delta$,
$\Gamma\preccurlyeq\Lambda\preccurlyeq 1$, $s\asymp 1$. Then $\hat{f}(x_{0})$
eventually (meaning for $n$ large enough) exists and

$\displaystyle\left\|\hat{f}(x_{0})-f^{\star}(x_{0})\right\|_{2}\preccurlyeq\Lambda+\delta^{\beta}\,.$
(124)

###### Proof.

To show that $\hat{f}(x_{0})$ eventually exists, we need to prove that
$\mathcal{H}_{\ell,s}(\hat{u}([0,T]),x)$ is eventually nonempty. As
$x_{0}\in\mathcal{G}_{\beta,\delta,s}(u^{\star}([0,T]))$, there is
$\mathbf{x}\in u^{\star}([0,T])^{N}$ with
$x\in\mathring{\mathsf{ch}}(\mathbf{x})$ such that
$\left\|\Psi(\eta_{\mathbf{x}}(\mathbf{x}))^{-1}\right\|_{\mathsf{op}}<s$. As
$\Gamma\prec 1$, the inequality in $s$ is strict, and
$\mathring{\mathsf{ch}}(\mathbf{x})$ is open, there must eventually be a
$\hat{\mathbf{x}}\in\hat{u}([0,T])^{N}$ with
$x\in\mathring{\mathsf{ch}}(\hat{\mathbf{x}})$ and
$\left\|\Psi(\eta_{\hat{\mathbf{x}}}(\hat{\mathbf{x}}))^{-1}\right\|_{\mathsf{op}}<s$.
Thus, $\hat{\mathbf{x}}(x_{0})$ exists and so does $\hat{f}(x_{0})$.

Let us establish some notation. Let
$\hat{\mathbf{x}}:=\hat{\mathbf{x}}(x_{0})$, Recall that
$\hat{\mathbf{x}}=(\hat{u}(\hat{\tau}_{1}),\dots,\hat{u}(\hat{\tau}_{N}))=\chi_{\ell,s}(\hat{u}([0,T]),x_{0})$.
Let $(u^{\star}(\tau_{1}),\dots,u^{\star}(\tau_{N}))=\mathbf{x}$ be any set
such that $x_{0}\in\mathring{\mathsf{ch}}(\mathbf{x})$,
$\|\Psi(\bar{\mathbf{x}})^{-1}\|_{\mathsf{op}}<s$, and
$\operatorname{\mathsf{diam}}(\mathbf{x})\leq\delta$, which exists as
$x_{0}\in\mathcal{G}_{\beta,\delta,s}(u^{\star}([0,T]))$. We summarize and add
some additional notation:

$\displaystyle\mathbf{x}$
$\displaystyle=(u^{\star}(\tau_{1}),\dots,u^{\star}(\tau_{N}))\,,$
$\displaystyle\mathbf{y}$
$\displaystyle:=(\dot{u}^{\star}(\tau_{1}),\dots,\dot{u}^{\star}(\tau_{N}))\,,$
(125) $\displaystyle\bar{\mathbf{x}}$
$\displaystyle:=(u^{\star}(\hat{\tau}_{1}),\dots,u^{\star}(\hat{\tau}_{N}))\,,$
$\displaystyle\bar{\mathbf{y}}$
$\displaystyle:=(\dot{u}^{\star}(\hat{\tau}_{1}),\dots,\dot{u}^{\star}(\hat{\tau}_{N}))\,,$
(126) $\displaystyle\tilde{\mathbf{x}}$
$\displaystyle:=(\hat{u}(\tau_{1}),\dots,\hat{u}(\tau_{N}))\,,$
$\displaystyle\tilde{\mathbf{y}}$
$\displaystyle:=(\hat{\dot{u}}(\tau_{1}),\dots,\hat{\dot{u}}(\tau_{N}))\,,$
(127) $\displaystyle\mathbf{x}$
$\displaystyle=(u^{\star}(\tau_{1}),\dots,u^{\star}(\tau_{N}))\,,$
$\displaystyle\mathbf{y}$
$\displaystyle:=(\dot{u}^{\star}(\tau_{1}),\dots,\dot{u}^{\star}(\tau_{N}))\,.$
(128)

Using the definition of $\hat{f}$ and the triangle inequality, we can split
the error into two parts,

$\left\|\hat{f}(x_{0})-f^{\star}(x_{0})\right\|_{2}\leq\left\|I(\mathbf{\hat{x}},\hat{\mathbf{y}},x_{0})-I(\mathbf{\hat{x}},\bar{\mathbf{y}},x_{0})\right\|_{2}+\left\|I(\mathbf{\hat{x}},\bar{\mathbf{y}},x_{0})-f^{\star}(x_{0})\right\|_{2}\,.$
(129)

For first term, with Lemma C.5, we obtain

$\displaystyle\left\|I(\mathbf{\hat{x}},\hat{\mathbf{y}},x_{0})-I(\mathbf{\hat{x}},\bar{\mathbf{y}},x_{0})\right\|_{2}$
$\displaystyle=\left\|\psi(x_{0})^{\\!\top}\\!\Psi(\mathbf{\hat{x}})^{-1}\left(\hat{\mathbf{y}}-\bar{\mathbf{y}}\right)\right\|_{2}$
(130)
$\displaystyle=\left\|\psi(\eta_{\mathbf{\hat{x}}}(x_{0}))^{\\!\top}\\!\Psi(\eta_{\mathbf{\hat{x}}}(\mathbf{\hat{x}}))^{-1}\left(\hat{\mathbf{y}}-\bar{\mathbf{y}}\right)\right\|_{2}$
(131)
$\displaystyle\leq\left\|\psi(\eta_{\mathbf{\hat{x}}}(x_{0}))\right\|_{2}\left\|\Psi(\eta_{\mathbf{\hat{x}}}(\mathbf{\hat{x}}))^{-1}\right\|_{\mathsf{op}}\left\|\hat{\mathbf{y}}-\bar{\mathbf{y}}\right\|_{\mathsf{op}}$
(132) $\displaystyle\leq c_{N}s\Lambda\,,$ (133)

where we used $\|\eta_{\mathbf{\hat{x}}}(x_{0})\|_{2}\leq 1$,
$\mathbf{\hat{x}}\in\mathcal{H}_{\ell,s}(\hat{u}([0,T]),x_{0})$, and the
definition of $\Lambda$ in the last step. Thus, we have

$\left\|I(\mathbf{\hat{x}},\hat{\mathbf{y}},x_{0})-I(\mathbf{\hat{x}},\bar{\mathbf{y}},x_{0})\right\|_{2}\preccurlyeq\Lambda$
(134)

For the second term, we apply Lemma C.6 and get

$\left\|I(\mathbf{\hat{x}},\bar{\mathbf{y}},x_{0})-f^{\star}(x_{0})\right\|_{2}\leq
C_{\beta}(\mathbf{\hat{x}})\operatorname{\mathsf{diam}}(\mathbf{\hat{x}})^{\beta}\,.$
(135)

With Lemma C.7 and the minimizing property of
$\mathbf{\hat{x}}=\chi_{\ell,s}(\hat{u}([0,T]),x_{0})$, we obtain

$\displaystyle
C_{\beta}(\mathbf{\hat{x}})\operatorname{\mathsf{diam}}(\mathbf{\hat{x}})^{\beta}$
$\displaystyle\leq\frac{N^{\frac{3}{2}}}{\beta!}\left\|\Psi(\eta_{\hat{\mathbf{x}}}(\mathbf{\hat{x}}))^{-1}\right\|_{\mathsf{op}}\operatorname{\mathsf{diam}}(\mathbf{\hat{x}})^{\beta}$
(136)
$\displaystyle\leq(1+\epsilon)\frac{N^{\frac{3}{2}}}{\beta!}\left\|\Psi(\eta_{\tilde{\mathbf{x}}}(\tilde{\mathbf{x}})))^{-1}\right\|_{\mathsf{op}}\operatorname{\mathsf{diam}}(\tilde{\mathbf{x}})^{\beta}\,.$
(137)

With $\|\hat{u}(\tau_{j})-u^{\star}(\tau_{j})\|_{2}\leq\Gamma$ and the
definition of $\mathbf{x}$, we have

$\operatorname{\mathsf{diam}}(\tilde{\mathbf{x}})\leq\operatorname{\mathsf{diam}}(\mathbf{x})+2\Gamma\preccurlyeq\delta+\Gamma\,.$
(138)

As $\left\|\tilde{\mathbf{x}}-{\mathbf{x}}\right\|_{\mathsf{op}}\prec 1$, we
have

$\displaystyle\left\|\Psi(\eta_{\tilde{\mathbf{x}}}(\tilde{\mathbf{x}})))^{-1}\right\|_{\mathsf{op}}$
$\displaystyle\preccurlyeq\left\|\Psi(\eta_{\mathbf{x}}(\tilde{\mathbf{x}})))^{-1}\right\|_{\mathsf{op}}$
(139)
$\displaystyle\preccurlyeq\left\|\Psi(\eta_{{\mathbf{x}}}({\mathbf{x}})))^{-1}\right\|_{\mathsf{op}}+\left\|\Psi(\eta_{{\mathbf{x}}}(\tilde{\mathbf{x}})))^{-1}-\Psi(\eta_{{\mathbf{x}}}(\mathbf{x})))^{-1}\right\|_{\mathsf{op}}\,.$
(140)

Using $\left\|\tilde{\mathbf{x}}-{\mathbf{x}}\right\|_{\mathsf{op}}\prec 1$
again, we obtain

$\|\Psi(\eta_{\mathbf{x}}(\tilde{\mathbf{x}}))-\Psi(\eta_{{\mathbf{x}}}({\mathbf{x}}))\|_{\mathsf{op}}\prec
1\,.$ (141)

Together with
$\|\Psi(\eta_{{\mathbf{x}}}({\mathbf{x}}))^{-1}\|_{\mathsf{op}}<s$, Lemma 4.20
then yields

$\left\|\Psi(\eta_{\mathbf{x}}(\hat{\mathbf{x}}))^{-1}-\Psi(\eta_{\mathbf{x}}(\bar{\mathbf{x}}))^{-1}\right\|_{\mathsf{op}}\prec
1\,.$ (142)

Putting the previous inequalities into (135), we obtain the following bound on
the second term:

$\left\|I(\mathbf{\hat{x}},{\mathbf{y}},x)-f^{\star}(x)\right\|_{2}\preccurlyeq\left(\delta+\Gamma\right)^{\beta}\,.$
(143)

The bounds (134) and (143) on the two terms of the right-hand side of (129)
together yield

$\left\|\hat{f}(x_{0})-f^{\star}(x_{0})\right\|_{2}\preccurlyeq\Lambda+\left(\delta+\Gamma\right)^{\beta}\preccurlyeq\Lambda+\delta^{\beta}$
(144)

as $\Gamma\preccurlyeq\delta$. ∎

###### Proof of Theorem 4.16.

Set $\tilde{\beta}=\beta+1$, $\tilde{L}_{0}=\infty$ and
$\tilde{L}_{k}=(k-1)!\sup_{\ell\in{\left\llbracket
0,\beta\right\rrbracket}}L_{\ell}^{k}$ for
$k\in{\left\llbracket\tilde{\beta}\right\rrbracket}$. As
$f\in\bar{\Sigma}^{d\to d}(\beta,L_{{\left\llbracket
0,\beta\right\rrbracket}})$, Corollary D.7 yields $u\in\bar{\Sigma}^{1\to
d}(\beta+1,\tilde{L}_{{\left\llbracket 0,\tilde{\beta}\right\rrbracket}})$.
Hence, we can use

$\Gamma=\Gamma(n,\tilde{\beta},\tilde{L}_{\tilde{\beta}},T)\qquad\text{and}\qquad\Lambda=\Lambda(n,\tilde{\beta},\tilde{L}_{\tilde{\beta}},T)$
(145)

in Lemma 4.21 to obtain

$\displaystyle\sup_{x\in\mathcal{X}}\left\|\hat{f}(x)-f(x)\right\|_{2}^{2}\in\mathbf{O}_{\mathbf{P}}\\!\left(\Lambda(n,\tilde{\beta},\tilde{L}_{\tilde{\beta}},T)^{2}+\delta^{2\beta}\right)\,.$
(146)

∎

##  Appendix A Derivatives

In this section, we introduce notation for derivatives in higher dimensions
and a Hölder-type smoothness class. Furthermore, some elementary analytical
results are shown for reference in the main proofs.

For any finite-dimensional $\mathbb{R}$-vector space $V$, we denote the
Euclidean norm as $\|x\|_{2}$ for $x\in V$. For $k\in\mathbb{N}$ and finite-
dimensional $\mathbb{R}$-vector spaces $V$ and $W$, let $\mathcal{L}_{k}(V,W)$
be the set of $k$-multilinear functions from $V^{k}$ to $W$. An element of
$\mathcal{L}_{k}(V,W)$ is called symmetric if it is invariant under
permutations of its $k$ arguments. A function $f\colon V\to W$ is
differentiable at $x\in V$ if there is $A\in\mathcal{L}_{1}(V,W)$ such that

$\lim_{\left\|v\right\|_{2}\to
0}\frac{\left\|f(x+v)-f(x)-A(v)\right\|_{2}}{\left\|v\right\|_{2}}=0\,.$ (147)

In this case we write $Df(x)=A$. The function $f$ is differentiable if it is
differentiable at all $x\in V$. Denote the set of differentiable functions
from $V$ to $W$ as $\mathcal{D}(V,W)$. From now on, assume $V=\mathbb{R}^{d}$
for a $d\in\mathbb{N}$. Define
$\mathbb{D}_{d}:=\left\\{v\in\mathbb{R}^{d}\,\kern-1.2pt\mathopen{}\left|\vphantom{v\in\mathbb{R}^{d}}\,\left\|v\right\|_{2}=1\right.\mathclose{}\kern-1.2pt\right\\}$.
The directional derivative of $f\in\mathcal{D}(V,W)$ at $x\in V$ in the
direction $v\in\mathbb{D}_{d}$ is

$D_{v}f(x):=Df(x)(v)\,.$ (148)

We set $D^{0}$ to the identity, i.e., $D^{0}f=f$. Let $k\in\mathbb{N}$. Define
the set of $k$-times differentiable functions $\mathcal{D}^{k}(V,W)$ and the
$k$-th derivative operator $D^{k}$ recursively: For $k=1$, we set
$\mathcal{D}^{1}(V,W):=\mathcal{D}(V,W)$ and $D^{1}:=D$. Let
$k\in\mathbb{N}_{\geq 2}$. A function $f\in\mathcal{D}^{k}(V,W)$ is $k$-times
differentiable at $x\in V$ if there is symmetric $A\in\mathcal{L}_{k}(V,W)$
such that

$\lim_{\left\|v\right\|_{2}\to
0}\sup_{\mathbf{v}\in\mathbb{D}_{d}^{k-1}}\frac{\left\|D^{k-1}_{\mathbf{v}}f(x+v)-D^{k-1}_{\mathbf{v}}f(x)-A((\mathbf{v},v))\right\|_{2}}{\left\|v\right\|_{2}}=0\,.$
(149)

In this case we write $D^{k}f(x)=A$. We may use square brackets for the
arguments of the linear operator to make formulas more clear, e.g.,
$Df(x)[v]:=Df(x)(v)$ The function $f$ is $k$-times differentiable if it is
$k$-times differentiable at all $x\in V$. Denote the set of $k$-times
differentiable functions from $V$ to $W$ as $\mathcal{D}^{k}(V,W)$. The $k$-th
directional derivative of $f\in\mathcal{D}^{k}(V,W)$ at $x\in V$ in the
directions $\mathbf{v}\in\mathbb{D}_{d}^{k}$ is

$D_{\mathbf{v}}f(x):=D^{k}f(x)(\mathbf{v})\,.$ (150)

We may indicate that differentiation takes place with respect to $x$ by
$D_{x,\mathbf{v}}g(x,y)=D_{\mathbf{v}}f_{y}(x)$, where $f_{y}(x)=g(x,y)$. Let
$p\in\mathbb{N}$. For $f\in\mathcal{D}^{k}(\mathbb{R}^{d},\mathbb{R}^{p})$ the
derivative operator acts component-wise, i.e., for $f=(f_{1},\dots,f_{p})$
with $f_{\ell}\in\mathcal{D}^{k}(\mathbb{R}^{d},\mathbb{R})$, we have

$D_{\mathbf{v}}f(x)=\begin{pmatrix}D_{\mathbf{v}}f_{1}(x)\\\ \vdots\\\
D_{\mathbf{v}}f_{p}(x)\end{pmatrix}\,,$ (151)

for $x\in\mathbb{R}^{d}$ and $\mathbf{v}\in\mathbb{D}_{d}^{k}$. The operator
norm for $A\in\mathcal{L}_{k}(V,W)$ is

$\left\|A\right\|_{\mathsf{op}}:=\sup_{\mathbf{v}\in\mathbb{D}_{d}^{k}}\left\|A(\mathbf{v})\right\|_{2}\,.$
(152)

Define the sup-norm of the $k$-th derivative as

$|D^{k}f|_{\infty}:=\sup_{x\in\mathbb{R}^{d}}\|D^{k}f(x)\|_{\mathsf{op}}=\sup_{x\in\mathbb{R}^{d}}\sup_{\mathbf{v}\in\mathbb{D}_{d}^{k}}\left\|D_{\mathbf{v}}f(x)\right\|_{2}\,.$
(153)

######  Notation A.1.

Let $x\in\mathbb{R}$. Denote the largest integer strictly smaller than $x$ as
$\llfloor x\rrfloor:=\max\mathbb{Z}_{<x}$.

######  Definition A.2 (Hölder-smoothness classes).

Let $d\in\mathbb{N}$, For $L\in\mathbb{R}_{>0}$, define

$\Sigma^{d\to
1}(0,L):=\left\\{f\colon\mathbb{R}^{d}\to\mathbb{R}\,\kern-1.2pt\mathopen{}\left|\vphantom{f\colon\mathbb{R}^{d}\to\mathbb{R}}\,\left|f\right|_{\infty}\leq
L\right.\mathclose{}\kern-1.2pt\right\\}\,.$ (154)

For $L\in\mathbb{R}_{>0}$, $\beta\in(0,1]$, define

$\Sigma^{d\to
1}(\beta,L):=\left\\{f\colon\mathbb{R}^{d}\to\mathbb{R}\,\kern-1.2pt\mathopen{}\left|\vphantom{f\colon\mathbb{R}^{d}\to\mathbb{R}}\,\left|f(x)-f(\tilde{x})\right|\leq
L\|x-\tilde{x}\|_{2}^{\beta}\right.\mathclose{}\kern-1.2pt\right\\}\,.$ (155)

Let $\beta\in\mathbb{R}_{>1}$ and $\ell:=\llfloor\beta\rrfloor$. Define

$\Sigma^{d\to
1}(\beta,L):=\left\\{f\in\mathcal{D}^{\ell}(\mathbb{R}^{d},\mathbb{R})\,\kern-1.2pt\mathopen{}\left|\vphantom{f\in\mathcal{D}^{\ell}(\mathbb{R}^{d},\mathbb{R})}\,\|D^{\ell}f(x)-D^{\ell}f(\tilde{x})\|_{\mathsf{op}}\leq
L\|x-\tilde{x}\|_{2}^{\beta-\ell}\right.\mathclose{}\kern-1.2pt\right\\}\,.$
(156)

For $\beta\in\mathbb{R}_{\geq 0}$, use the short notation

$\Sigma(\beta,L):=\Sigma^{1\to 1}(\beta,L)\,.$ (157)

Let $L_{0},\dots,L_{\ell}\in\mathbb{R}_{>0}\cup\\{\infty\\}$,
$L_{\beta}\in\mathbb{R}_{>0}$. Define

$\Sigma^{d\to 1}(\beta,L_{{\left\llbracket
0,\ell\right\rrbracket}},L_{\beta}):=\Sigma^{d\to
1}(\beta,L_{\beta})\cap\bigcap_{k=0}^{\ell}\Sigma^{d\to 1}(k,L_{k})\,.$ (158)

Let $d_{\mathsf{in}},d_{\mathsf{out}}\in\mathbb{N}$. Define

$\Sigma^{d_{\mathsf{in}}\to d_{\mathsf{out}}}(\beta,L_{{\left\llbracket
0,\ell\right\rrbracket}},L_{\beta}):=\left(\Sigma^{d_{\mathsf{in}}\to
1}(\beta,L_{{\left\llbracket
0,\ell\right\rrbracket}},L_{\beta})\right)^{d_{\mathsf{out}}}\,,$ (159)

where $f=(f_{1},\dots,f_{d_{\mathsf{out}}})\in\Sigma^{d_{\mathsf{in}}\to
d_{\mathsf{out}}}(\beta,L_{{\left\llbracket
0,\ell\right\rrbracket}},L_{\beta})$ is treated as a function

$f\colon\mathbb{R}^{d_{\mathsf{in}}}\to\mathbb{R}^{d_{\mathsf{out}}},x\mapsto(f_{1}(x),\dots,f_{d_{\mathsf{out}}}(x))^{\\!\top}\\!\,.$
(160)

If $\beta$ is an integer, we denote

$\bar{\Sigma}^{d_{\mathsf{in}}\to d_{\mathsf{out}}}(\beta,L_{{\left\llbracket
0,\beta\right\rrbracket}}):=\Sigma^{d_{\mathsf{in}}\to
d_{\mathsf{out}}}(\beta,L_{{\left\llbracket
0,\beta\right\rrbracket}})\cap\mathcal{D}^{\beta}(\mathbb{R}^{d_{\mathsf{in}}},\mathbb{R}^{d_{\mathsf{out}}})\,.$
(161)

A bound on the derivative yields the Lipschitz constant, as the next lemma
shows.

######  Lemma A.3.

Let $d\in\mathbb{N}$. Let $k\in\mathbb{N}$. Let
$f\in\mathcal{D}^{k}(\mathbb{R}^{d},\mathbb{R})$. Let $L\in\mathbb{R}_{>0}$.
Assume $|D^{k}f|_{\infty}\leq L$. Then $f\in\Sigma^{d\to 1}(k,L)$.

###### Proof of Lemma A.3.

Let $\mathbf{v}\in\mathbb{D}_{d}^{k-1}$. Then, by the mean value theorem,

$\left|D^{k-1}_{\mathbf{v}}f(x)-D^{k-1}_{\mathbf{v}}f(\tilde{x})\right|\leq\sup_{\tilde{v}\in\mathbb{D}_{d}}\left|D^{k}_{(\mathbf{v},\tilde{v})}f\right|_{\infty}\left\|x-\tilde{x}\right\|_{2}$
(162)

with $(\mathbf{v},\tilde{v})\in\mathbb{D}_{d}^{k}$. Thus,

$\left\|D^{k-1}f(x)-D^{k-1}f(\tilde{x})\right\|_{\mathsf{op}}\leq\left|D^{k}f\right|_{\infty}\left\|x-\tilde{x}\right\|_{2}\,.$
(163)

∎

The sub-multiplicativity of the operator norm in the following lemma is a
direct consequence of its definition.

######  Lemma A.4.

Let $k,d,p\in\mathbb{N}$. Let
$A\in\mathcal{L}_{k}(\mathbb{R}^{d},\mathbb{R}^{p})$ and
$v\in\mathbb{D}_{d}^{k}$. Then

$\left\|A(v_{1},\dots,v_{k})\right\|_{2}\leq\left\|A\right\|_{\mathsf{op}}\prod_{j=1}^{k}\left\|v_{j}\right\|_{2}$
(164)

For reference, we state the chain and product rule as well es Taylor’s
approximation theorem in the following lemmas.

######  Lemma A.5.

Let $d\in\mathbb{N}$. Let
$A\colon\mathbb{R}^{d}\to\mathcal{L}_{k}(\mathbb{R}^{d},\mathbb{R}^{d})$
symmetric. Let $v\in\mathbb{D}_{d}$.

  1. (i)

Product rule: Let $f_{\ell}\in\mathcal{D}(\mathbb{R}^{d},\mathbb{R}^{d})$ for
$\ell\in{\left\llbracket k\right\rrbracket}$. Then

$\displaystyle D_{v}\left(A(x)(f_{1}(x),\dots,f_{k}(x))\right)$ (165)
$\displaystyle=(D_{v}A(x))(f_{1}(x),\dots,f_{k}(x))\,+$ (166)
$\displaystyle\phantom{=}\
\sum_{\ell=1}^{k}A(x)\\!\left(f_{1}(x),\dots,f_{\ell-1}(x),D_{v}f_{\ell}(x),f_{\ell+1}(x),\dots,f_{k}(x)\right)\,.$
(167)

  2. (ii)

Chain rule: Let $f\in\mathcal{D}(\mathbb{R}^{d},\mathbb{R}^{d})$. Then

$D_{v}\left(A(f(x))(a_{1},\dots,a_{k})\right)=(DA)(f(x))(a_{1},\dots,a_{k},D_{v}f(x))\,.$
(168)

######  Notation A.6.

Let $d,m\in\mathbb{N}$. Let $x\in\mathbb{R}^{d}$. Denote
$\\{x\\}^{m}=(x,\dots,x)\in(\mathbb{R}^{d})^{m}$.

######  Lemma A.7 (Taylor theorem with Peano remainder).

Let $d\in\mathbb{N}$.

  1. (i)

Let $\ell\in\mathbb{N}_{0}$. Let
$f\in\mathcal{D}^{\ell}(\mathbb{R}^{d},\mathbb{R})$. There is
$\gamma\in[0,1]$, such that

$f(x)=\sum_{k=0}^{\ell-1}\frac{D^{k}f(x_{0})(\\{x-x_{0}\\}^{k})}{k!}+\frac{D^{\ell}f\\!\left(\tilde{x}\right)(\\{x-x_{0}\\}^{\ell})}{\ell!}$
(169)

with $\tilde{x}:=x_{0}+\gamma(x-x_{0})$.

  2. (ii)

Let $\beta,L\in\mathbb{R}_{>0}$. Set $\ell:=\llfloor\beta\rrfloor$. Let
$f\in\Sigma^{d\to 1}(\beta,L)$. Then

$\left|f(x)-\sum_{k=0}^{\ell}\frac{D^{k}f(x_{0})(\\{x-x_{0}\\}^{k})}{k!}\right|\leq
L\frac{\left\|x-x_{0}\right\|_{2}^{\beta}}{\ell!}\,.$ (170)

##  Appendix B Local Polynomial Regression

In order to make this article more self-contained, we provide here upper
bounds on the estimation error for a multivariate regression problem using
local polynomial estimators. We follow [Tsy09] and extend the proof there from
the domain $[0,1]$ to the domain $[0,T]^{d}$ and also discuss estimation of
derivatives. See also [RW94]. It is important for the results in the main part
of the article that the dependence on $T$ and on the Lipschitz constant $L$ of
the regression function are given explicitly.

###  B.1 Model

Let $d\in\mathbb{N}$. Let $\beta,L\in\mathbb{R}_{>0}$. Set
$\ell:=\llfloor\beta\rrfloor$. Let $f\in\Sigma^{d\to 1}(\beta,L)$. Let
$n\in\mathbb{N}$. Let $T\in\mathbb{R}_{>0}$. For $i\in{\left\llbracket
n\right\rrbracket}$, let $x_{i}\in[0,T]^{d}$. Let $\sigma\in\mathbb{R}_{>0}$.
For $i\in{\left\llbracket n\right\rrbracket}$, let $\varepsilon_{i}$ be
independent random variables with $\mathbf{E}[\epsilon_{i}]=0$ and
$\mathbf{E}[\epsilon_{i}^{2}]\leq\sigma^{2}$. For $i\in{\left\llbracket
n\right\rrbracket}$, set

$Y_{i}=f(x_{i})+\varepsilon_{i}\,.$ (171)

We observe $Y_{i}$ and know $x_{i}$. The function $f$ is unknown. Let
$s\in{\left\llbracket 0,\ell\right\rrbracket}$. Let
$\mathbf{v}\in\mathbb{D}_{d}^{s}$ be $s$-many directions. Our target is to
estimate the $s$-th directional derivative
$D_{\mathbf{v}}f\colon\mathbb{R}^{d}\to\mathbb{R}$ of $f$. For the local
polynomial estimator $\hat{g}$ defined below, we want to show upper bounds for
the mean squared error $\mathbf{E}[\|\hat{g}(x)-D_{\mathbf{v}}f(x)\|^{2}]$ at
a point $x\in[0,T]^{d}$ depending on $d,\beta,L,n,T,\sigma,s$.

###  B.2 Estimator

Recall 3.8. Set the number of scalar coefficients for polynomials
$\mathbb{R}^{d}\to\mathbb{R}$ of degree at most $\ell$ as

$N:=N_{d,\ell}:=\dim(\mathcal{P}_{d,\ell})=\begin{pmatrix}\ell+d\\\
d\end{pmatrix}\,.$ (172)

For $x\in\mathbb{R}^{d}$, denote
$\psi(x):=\psi_{\ell}(x):=((\alpha!)^{-1}x^{\alpha})_{\alpha\in\mathbb{N}_{0}^{d},|\alpha|\leq\ell}\in\mathbb{R}^{N}$.
Let $K\colon\mathbb{R}_{\geq 0}\to\mathbb{R}$. Let $h\in\mathbb{R}_{>0}$.
Denote the local polynomial mean squared error coefficients as

$\hat{\theta}(x)\in\operatorname*{arg\,min}_{\theta\in\mathbb{R}^{N}}\sum_{i=1}^{n}\left(Y_{i}-\theta^{\\!\top}\\!\psi\\!\left(\frac{x_{i}-x}{h}\right)\right)^{2}K\\!\left(\frac{\left\|x_{i}-x\right\|_{2}}{h}\right)\,.$
(173)

Then the local polynomial estimator of degree $\ell$ for the $s$-th derivative
in directions $\mathbf{v}\in\mathbb{D}_{d}^{s}$ is defined as

$\hat{g}(x):=h^{-s}D_{\mathbf{v}}\psi(0)^{\\!\top}\\!\hat{\theta}(x)\,.$ (174)

For $x\in\mathbb{R}^{d}$, define the symmetric matrix
$B(x)\in\mathbb{R}^{N\times N}$ and the vector $a(x)\in\mathbb{R}^{N}$ as

$\displaystyle B(x)$
$\displaystyle:=\frac{T^{d}}{nh^{d}}\sum_{i=1}^{n}\psi\\!\left(\frac{x_{i}-x}{h}\right)\psi\\!\left(\frac{x_{i}-x}{h}\right)^{\\!\top}\\!K\\!\left(\frac{\left\|x_{i}-x\right\|_{2}}{h}\right)\,,$
(175) $\displaystyle a(x)$
$\displaystyle:=\frac{T^{d}}{nh^{d}}\sum_{i=1}^{n}Y_{i}\psi\\!\left(\frac{x_{i}-x}{h}\right)K\\!\left(\frac{\left\|x_{i}-x\right\|_{2}}{h}\right)\,.$
(176)

######  Notation B.1.

Let $k\in\mathbb{N}$ and $A\in\mathbb{R}^{k\times k}$ be a symmetric matrix.
We write $\lambda_{\mathsf{min}}(A)$ for the smallest eigenvalue of $A$.

Assume $\lambda_{\mathsf{min}}(B(x))>0$. Then

$\hat{\theta}(x)=B(x)^{-1}a(x)\,.$ (177)

Furthermore, with the weights

$w_{i}(x):=w_{i,s}(x):=\frac{T^{d}}{nh^{d+s}}D_{\mathbf{v}}\psi(0)^{\\!\top}\\!B(x)^{-1}\psi\\!\left(\frac{x_{i}-x}{h}\right)K\\!\left(\frac{x_{i}-x}{h}\right)\,,$
(178)

we can write the estimator as

$\hat{g}(x)=h^{-s}D_{\mathbf{v}}\psi(0)^{\\!\top}\\!\hat{\theta}(x)=\sum_{i=1}^{n}w_{i}(x)Y_{i}\,.$
(179)

###  B.3 Result

######  Notation B.2.

All lower case $c$, with or without index, are elements of $\mathbb{R}_{>0}$
and universal insofar as they only depend on the variables written as index,
e.g., $c_{d,\beta}$ depends only on $\beta$ and $d$. In particular, a constant
$c$ with no index refers to a fixed positive number. Every occurrence of such
a variable may refer to a different value.

######  Theorem B.3.

Use the model of section B.1 and the estimator of section B.2. Assume Kernel,
Eigenvalue, and Cover. Assume $n\geq C_{\mathsf{cvr}}^{-1}h^{-d}T^{d}$. Then,
for all $x\in[0,T]^{d}$, we have

$\mathbf{E}\\!\left[\left(\hat{g}(x)-D_{\mathbf{v}}f(x)\right)^{2}\right]\leq
c_{\ell}\left(C_{\mathsf{egv}}^{2}C_{\mathsf{ker}}^{2}C_{\mathsf{cvr}}^{2}L^{2}h^{2(\beta-s)}+C_{\mathsf{egv}}^{2}C_{\mathsf{ker}}^{2}C_{\mathsf{cvr}}\frac{\sigma^{2}T^{d}}{nh^{d+2s}}\right)\,.$
(180)

######  Corollary B.4.

Use the model of section B.1 and the estimator of section B.2. Assume Kernel,
Eigenvalue, and Cover. Assume $n\geq n_{0}$, where $n_{0}\in\mathbb{N}$
depends only on
$d,\beta,C_{\mathsf{egv}},C_{\mathsf{ker}},C_{\mathsf{cvr}},\sigma,T,L$. Then
there are constants $C_{1},C_{2}\in\mathbb{R}_{>0}$, depending only on
$d,\beta,C_{\mathsf{egv}},C_{\mathsf{ker}},C_{\mathsf{cvr}}$, with the
following property: Set the bandwidth as

$h=C_{1}\left(\frac{\sigma^{2}T^{d}}{L^{2}n}\right)^{\frac{1}{2\beta+d}}\,.$
(181)

Then

$\mathbf{E}\\!\left[\left(\hat{g}(x)-D_{\mathbf{v}}f(x)\right)^{2}\right]\leq
C_{2}L^{\frac{2d+4s}{2\beta+d}}\left(\frac{\sigma^{2}T^{d}}{n}\right)^{\frac{2(\beta-s)}{2\beta+d}}\,.$
(182)

######  Remark B.5.

  1. (i)

The error rates obtained in Corollary B.4 are minimax optimal up to the
constant $C_{2}$, which can be seen by comparing them to lower bounds for
nonparametric regression, e.g., [SS24, Theorem 8].

  2. (ii)

Integrating the point-wise error over $[0,T]^{d}$ yields the minimax optimal
rate for the integrated error. For a bound on the error in sup-norm, we need
to do some extra work.

######  Theorem B.6.

Assume Kernel, Eigenvalue, and Cover. Assume SubGaussian. Assume $n\geq
C_{\mathsf{cvr}}^{-1}h^{-d}T^{d}$ Assume $T/h\geq 2$. Then, we have

$\mathbf{E}\\!\left[\left|\hat{g}(x)-D_{\mathbf{v}}f(x)\right|_{\infty}^{2}\right]\leq
c_{d,\ell,K}\left(\left(C_{\mathsf{egv}}C_{\mathsf{cvr}}Lh^{\beta-s}\right)^{2}+C_{\mathsf{cvr}}C_{\mathsf{egv}}^{2}\frac{\sigma^{2}T^{d}}{nh^{d+2s}}\log(T/h)\right)\,.$
(183)

######  Corollary B.7.

Assume Kernel, Eigenvalue, and Cover. Assume SubGaussian. There are
$C_{1},C_{2},C_{3},C_{4}\in\mathbb{R}_{>0}$ depending only on
$d,\beta,C_{\mathsf{egv}},C_{\mathsf{cvr}},K$ with the following property:
Assume

$C_{1}\left(\frac{\log(n)}{n}\right)^{\frac{d}{2\beta}}\leq
C_{2}\left(\frac{L^{2}}{\sigma^{2}}\right)^{\frac{d}{2\beta}}T^{d}\leq
n\log(n)^{\frac{d}{2\beta}}\,.$ (184)

Set the bandwidth as

$h=C_{3}\left(\frac{\sigma^{2}T^{d}\log(n)}{L^{2}n}\right)^{\frac{1}{2\beta+d}}\,.$
(185)

Then

$\mathbf{E}\\!\left[\left|\hat{g}-D_{\mathbf{v}}f\right|_{\infty}^{2}\right]\leq
C_{4}L^{\frac{2d+4s}{2\beta+d}}\left(\frac{\sigma^{2}T^{d}\log(n)}{n}\right)^{\frac{2(\beta-s)}{2\beta+d}}\,.$
(186)

######  Remark B.8.

  1. (i)

The error rates obtained in Corollary B.7 are minimax optimal up to the
constant $C_{4}$, which can be seen by comparing them to lower bounds for
nonparametric regression, e.g., [SS24, Theorem 8].

Following [Tsy09, Lemma 1.4, Lemma 1.5], we show that the uniform grid in
$[0,T]^{d}$ fulfills Cover, and also Eigenvalue if the kernel is not
degenerated.

######  Proposition B.9.

Let $n_{0}\in\mathbb{N}$. Set $n:=n_{0}^{d}$. Let $x_{i}\in[0,T]^{d}$ for
$i\in{\left\llbracket n\right\rrbracket}$ form a uniform grid in $[0,T]^{d}$,
i.e.,

$\left\\{x_{1},\dots,x_{n}\right\\}=\left\\{\begin{pmatrix}T\frac{k_{1}}{n_{0}},\dots,T\frac{k_{d}}{n_{0}}\end{pmatrix}^{\\!\top}\\!\,\kern-1.2pt\mathopen{}\left|\vphantom{\begin{pmatrix}T\frac{k_{1}}{n_{0}},\dots,T\frac{k_{d}}{n_{0}}\end{pmatrix}^{\\!\top}\\!}\,k_{1},\dots,k_{d}\in{\left\llbracket
n_{0}\right\rrbracket}\right.\mathclose{}\kern-1.2pt\right\\}\,.$ (187)

Assume StrictKernel with constants $C_{\mathsf{ker}}\in\mathbb{R}_{>0}$. Then
Kernel, Eigenvalue, and Cover are fulfilled with $C_{\mathsf{ker}}$,
$C_{\mathsf{egv}}=c_{d,\ell,K}(1+\frac{T}{n_{0}h})$, $C_{\mathsf{cvr}}=4^{d}$.

###  B.4 Proof

######  Lemma B.10.

Let $q\in\mathcal{P}_{d,\ell}$. Assume $\lambda_{\mathsf{min}}(B(x))>0$. Then

$\sum_{i=1}^{n}w_{i}(x)q(x_{i})=D_{\mathbf{v}}q(x)\,.$ (188)

###### Proof.

We set $Y_{i}$ to $q(x_{i})$ and consider the local polynomial fit. The
coefficients $\theta_{q}(x)$ in the fit are

$\theta_{q}(x)\in\operatorname*{arg\,min}_{\theta\in\mathbb{R}^{N}}\sum_{i=1}^{n}\left(q(x_{i})-\theta^{\\!\top}\\!\psi\\!\left(\frac{x_{i}-x}{h}\right)\right)^{2}K\\!\left(\frac{\left\|x_{i}-x\right\|_{2}}{h}\right)\,.$
(189)

Let $z\in\mathbb{R}^{d}$. Since $q(z)$ is a polynomial of degree at most
$\ell$ and $\lambda_{\mathsf{min}}(B(x))>0$, the fit of the degree $\ell$
polynomial with base $\psi((z-x)/h)$ retrieves the original function, i.e.,

$\displaystyle
q(z)=\theta_{q}(x)^{\\!\top}\\!\psi\\!\left(\frac{z-x}{h}\right)\,.$ (190)

Denote by $D_{z,\mathbf{v}}$ the $s$-th directional derivative in directions
$\mathbf{v}$ with respect to $z$. Then,

$D_{\mathbf{v}}q(z)=\theta_{q}(x)^{\\!\top}\\!D_{z,\mathbf{v}}\left(\psi\\!\left(\frac{z-x}{h}\right)\right)=h^{-s}\theta_{q}(x)^{\\!\top}\\!\left(D_{\mathbf{v}}\psi\right)\\!\left(\frac{z-x}{h}\right)\,.$
(191)

Thus,

$D_{\mathbf{v}}q(x)=h^{-s}\theta_{q}(x)^{\\!\top}\\!\left(D_{\mathbf{v}}\psi\right)\\!\left(0\right)=\sum_{i=1}^{n}w_{i}(x)q(x_{i})\,.$
(192)

∎

######  Lemma B.11.

Assume Kernel, Eigenvalue, and Cover. Assume $n\geq
C_{\mathsf{cvr}}^{-1}h^{-d}T^{d}$.

  1. (i)

If $\|x-x_{i}\|\geq h$, we have $w_{i}(x)=0$.

  2. (ii)

We have

$\sum_{i=1}^{n}\left|w_{i}(x)\right|\leq
c_{\ell}C_{\mathsf{egv}}C_{\mathsf{ker}}C_{\mathsf{cvr}}h^{-s}\,.$ (193)

  3. (iii)

We have

$\sum_{i=1}^{n}w_{i}(x)^{2}\leq
c_{\ell}C_{\mathsf{egv}}^{2}C_{\mathsf{ker}}^{2}C_{\mathsf{cvr}}\frac{T^{d}}{nh^{d+2s}}\,.$
(194)

###### Proof.

Using Kernel and Eigenvalue, we obtain from (178)

$\displaystyle nh^{d+s}T^{-d}\left|w_{i}(x)\right|$
$\displaystyle\leq\left\|D_{\mathbf{v}}\psi(0)\right\|_{2}\|B(x)^{-1}\|_{\mathsf{op}}\left\|\psi\\!\left(\frac{x_{i}-x}{h}\right)\right\|_{2}|K|_{\infty}\mathds{1}_{\mathrm{B}^{d}(x_{i},h)}$
(195) $\displaystyle\leq
c_{\ell}C_{\mathsf{egv}}C_{\mathsf{ker}}\mathds{1}_{\mathrm{B}^{d}(x_{i},h)}\,,$
(196)

which we use in each of the following parts.

  1. (i)

Trivial.

  2. (ii)

By Cover,

$\displaystyle\left(C_{\mathsf{egv}}C_{\mathsf{ker}}\right)^{-1}\sum_{i=1}^{n}\left|w_{i}(x)\right|$
$\displaystyle\leq
c_{\ell}T^{d}h^{-(d+s)}\frac{1}{n}\\#\left\\{i\in{\left\llbracket
n\right\rrbracket}\,\kern-1.2pt\mathopen{}\left|\vphantom{i\in{\left\llbracket
n\right\rrbracket}}\,x_{i}\in\mathrm{B}^{d}(x,h)\right.\mathclose{}\kern-1.2pt\right\\}$
(197) $\displaystyle\leq
c_{\ell}T^{d}h^{-(d+s)}\max\\!\left(\frac{1}{n},\,C_{\mathsf{cvr}}\left(\frac{h}{T}\right)^{d}\right)$
(198) $\displaystyle\leq c_{\ell}C_{\mathsf{cvr}}h^{-s}\,,$ (199)

as we assume $n\geq C_{\mathsf{cvr}}^{-1}h^{-d}T^{d}$.

  3. (iii)

With

$\left|w_{i}(x)\right|\leq
c_{\ell}C_{\mathsf{egv}}C_{\mathsf{ker}}\frac{T^{d}}{nh^{d+s}}$ (200)

and (ii), we obtain

$\displaystyle\sum_{i=1}^{n}\left|w_{i}(x)\right|^{2}$
$\displaystyle\leq\sup_{i\in{\left\llbracket
n\right\rrbracket}}\left|w_{i}(x)\right|\,\sum_{i=1}^{n}\left|w_{i}(x)\right|$
(201) $\displaystyle\leq
c_{\ell}C_{\mathsf{egv}}^{2}C_{\mathsf{ker}}^{2}C_{\mathsf{cvr}}\frac{T^{d}}{nh^{d+2s}}\,.$
(202)

∎

###### Proof of Theorem B.3.

Bias. By the model equation (171) and the linear representation of the
estimator (179), we have

$\mathbf{E}\\!\left[\hat{g}(x)\right]=\sum_{i=1}^{n}w_{i}(x)f(x_{i})\,.$ (203)

Using a Taylor expansion (Lemma A.7) of $f(x)$ at $x_{i}$ yields the existence
of $\gamma_{i}\in[0,1]$ such that

$f(x)=\sum_{k=0}^{\ell-1}\frac{D^{k}f(x)(\\{x_{i}-x\\}^{k})}{k!}+\frac{D^{\ell}f\\!\left(x+\gamma_{i}(x_{i}-x)\right)(\\{x_{i}-x\\}^{\ell})}{\ell!}\,.$
(204)

Denote the remainder term as

$R_{i}(x):=\frac{D^{\ell}f\\!\left(x+\gamma_{i}(x_{i}-x)\right)(\\{x_{i}-x\\}^{\ell})-D^{\ell}f\\!\left(x\right)(\\{x_{i}-x\\}^{\ell})}{\ell!}\,.$
(205)

By applying Lemma B.10 to

$q(x)=\sum_{k=0}^{\ell}\frac{D^{k}f(x)(\\{x_{i}-x\\}^{k})}{k!}\,,$ (206)

we obtain

$\displaystyle\sum_{i=1}^{n}w_{i}(x)f(x_{i})$
$\displaystyle=\sum_{i=1}^{n}w_{i}(x)\sum_{k=0}^{\ell}\frac{D^{k}f(x)(\\{x_{i}-x\\}^{k})}{k!}+\sum_{i=1}^{n}w_{i}(x)R_{i}(x)$
(207) $\displaystyle=D_{\mathbf{v}}f(x)+\sum_{i=1}^{n}w_{i}(x)R_{i}(x)\,.$
(208)

Hence, we have

$\displaystyle\mathbf{E}\\!\left[\hat{g}(x)\right]-D_{\mathbf{v}}f(x)=\sum_{i=1}^{n}w_{i}(x)R_{i}(x)\,.$
(209)

For bounding the remainder term, we note that $f\in\Sigma^{d\to 1}(\beta,L)$,
which yields

$\left|D^{\ell}f\\!\left(x+\gamma_{i}(x_{i}-x)\right)(\\{x_{i}-x\\}^{\ell})-D^{\ell}f\\!\left(x\right)(\\{x_{i}-x\\}^{\ell})\right|\leq
L\left\|x_{i}-x\right\|_{2}^{\beta}\,.$ (210)

Thus, with Lemma B.11 (i), we have

$\displaystyle\left|\mathbf{E}\\!\left[\hat{g}(x)\right]-D_{\mathbf{v}}f(x)\right|$
$\displaystyle\leq
L\sum_{i=1}^{n}\left|w_{i}(x)\right|\left\|x_{i}-x\right\|_{2}^{\beta}$ (211)
$\displaystyle\leq Lh^{\beta}\sum_{i=1}^{n}\left|w_{i}(x)\right|\,.$ (212)

Finally, by Lemma B.11 (ii), we obtain

$\left|\mathbf{E}\\!\left[\hat{g}(x)\right]-D_{\mathbf{v}}f(x)\right|\leq
c_{\ell}C_{\mathsf{egv}}C_{\mathsf{ker}}C_{\mathsf{cvr}}Lh^{\beta-s}\,.$ (213)

Variance. As the noise is independent and has a second moment bounded by
$\sigma^{2}$, we have

$\displaystyle\mathbf{E}\\!\left[\left(\hat{g}(x)-\mathbf{E}\\!\left[\hat{g}(x)\right]\right)^{2}\right]$
$\displaystyle=\mathbf{E}\\!\left[\left(\sum_{i=1}^{n}\varepsilon_{i}w_{i}(x)\right)^{2}\right]$
(214)
$\displaystyle=\sum_{i=1}^{n}w_{i}(x)^{2}\mathbf{E}\\!\left[\varepsilon_{i}^{2}\right]$
(215) $\displaystyle\leq\sigma^{2}\sum_{i=1}^{n}w_{i}(x)^{2}\,.$ (216)

Hence, Lemma B.11 (iii) yields

$\mathbf{E}\\!\left[\left(\hat{g}(x)-\mathbf{E}\\!\left[\hat{g}(x)\right]\right)^{2}\right]\leq
c_{\ell}C_{\mathsf{egv}}^{2}C_{\mathsf{ker}}^{2}C_{\mathsf{cvr}}\frac{\sigma^{2}T^{d}}{nh^{d+2s}}\,.$
(217)

Together. Putting (213) and (217) together, we obtain

$\displaystyle\mathbf{E}\\!\left[\left(\hat{g}(x)-D_{\mathbf{v}}f(x)\right)^{2}\right]$
$\displaystyle\leq
c_{\ell}\left(C_{\mathsf{egv}}^{2}C_{\mathsf{ker}}^{2}C_{\mathsf{cvr}}^{2}L^{2}h^{2(\beta-s)}+C_{\mathsf{egv}}^{2}C_{\mathsf{ker}}^{2}C_{\mathsf{cvr}}\frac{\sigma^{2}T^{d}}{nh^{d+2s}}\right)\,.$
(218)

∎

For the proof in sup-norm (Theorem B.6), we want an upper bound on the
supremum of a sub-gaussian process using the chaining technique. In order to
do this, we need some preparation:

######  Definition B.12 (Sub-Gaussian process).

Let $(\mathcal{X},\rho)$ be a metric space. A real-valued stochastic process
$(Z_{x})_{x\in\mathcal{X}}$ is sub-Gaussian if

$\mathbf{P}\\!\left(\left|Z_{x}-Z_{x^{\prime}}\right|\geq t\right)\leq
2\exp\\!\left(-\frac{t^{2}}{2\rho(x,x^{\prime})^{2}}\right)$ (219)

for all $t\in\mathbb{R}_{>0}$, $x,x^{\prime}\in\mathcal{X}$.

######  Lemma B.13.

Let $\epsilon_{i}$ be independent, centered, sub-Gaussian random variables
with variance parameters $v_{i}\in\mathbb{R}_{>0}$. Let $(\mathcal{X},\rho)$
be a metric space. Let $b_{i}\colon\mathcal{X}\to\mathbb{R}$. Define
$Z_{x}:=\sum_{i=1}^{n}b_{i}(x)\epsilon_{i}$. Then $Z_{x}$ is sub-Gaussian with
variance parameter $v_{x}:=\sum_{i=1}^{n}b_{i}(x)^{2}v_{i}$. Furthermore,
$(Z_{x})_{x\in\mathcal{X}}$ is sub-Gaussian with
$\rho(x,x^{\prime})^{2}:=\sum_{i=1}^{n}\left(b_{i}(x)-b_{i}(x^{\prime})\right)^{2}v_{i}$.

###### Proof.

The results stem from the well-known fact that a linear combination with non-
negative coefficients of independent sub-Gaussian variables is again sub-
Gaussian. It can be checked using the following equivalent condition for sub-
Gaussian: A random variable $X$ is sub-Gaussian with variance parameter $v$ if
and only if $\mathbf{E}[\exp(\lambda
X)]\leq\exp\\!\left(\frac{\lambda^{2}v}{2}\right)$ for all
$\lambda\in\mathbb{R}_{\geq 0}$. ∎

######  Definition B.14 (Covering Number).

Let $(\mathcal{X},\rho)$ be a metric space, $A\subseteq\mathcal{X}$, and
$r\in\mathbb{R}_{>0}$. The covering number $N(A,\rho,r)$ with balls
$\mathrm{B}^{\rho}(x,r):=\\{z\in\mathcal{X}\mid\rho(x,z)\leq r\\}$ of radius
$r$ is defined as

$N(A,\rho,r):=\min\left\\{n\in\mathbb{N}\,\kern-1.2pt\mathopen{}\left|\vphantom{n\in\mathbb{N}}\,\exists(x_{i})_{i\in{\left\llbracket
n\right\rrbracket}}\subseteq\mathcal{X}\colon
A\subseteq\bigcup_{i\in{\left\llbracket
n\right\rrbracket}}\mathrm{B}^{\rho}(x_{i},r)\right.\mathclose{}\kern-1.2pt\right\\}\,.$
(220)

######  Lemma B.15.

Let $d\in\mathbb{N}$, $a,h,T\in\mathbb{R}_{>0}$. Define the metric $\rho$ in
$\mathbb{R}^{d}$ as

$\rho(x,x^{\prime}):=a\min\\!\left(1\,,\
\frac{\left\|x-x^{\prime}\right\|_{2}}{h}\right)\,.$ (221)

Then

$\int_{0}^{\infty}\sqrt{\log(N([0,T]^{d},\rho,r))}\mathrm{d}r\leq
c_{d}a\max\\!\left(1,\sqrt{\log\\!\left(\frac{T}{h}\right)}\right)\,.$ (222)

###### Proof.

First consider $a=1$. For $r\leq 1$,

$N([0,T]^{d},\rho,r)\leq c_{d}T^{d}h^{-d}r^{-d}\,.$ (223)

For $r\geq 1$,

$N([0,T]^{d},\rho,r)=1\,.$ (224)

Thus,

$\displaystyle\int_{0}^{\infty}\sqrt{\log(N([0,T]^{d},\rho,r))}\mathrm{d}r$
$\displaystyle\leq\int_{0}^{1}\sqrt{\log(c_{d}T^{d}h^{-d}r^{-d})}\mathrm{d}r$
(225) $\displaystyle\leq
c_{d}+\sqrt{d\log(T/h)}+\sqrt{d}\int_{0}^{1}\sqrt{\log(1/r)}\mathrm{d}r\,.$
(226)

We have

$\int_{0}^{1}\sqrt{\log(1/r)}\mathrm{d}r=\int_{0}^{\infty}\sqrt{t}\exp(-t)\mathrm{d}t=\Gamma\\!\left(\frac{3}{2}\right)=\frac{1}{2}\sqrt{\pi}\,,$
(227)

where $\Gamma(\alpha)$ is the Gamma-function. All in all, we obtain

$\int_{0}^{\infty}\sqrt{\log(N([0,T]^{d},\rho,r))}\mathrm{d}r=c_{d}\max\\!\left(1\,,\,\sqrt{\log(T/h)}\right)\,.$
(228)

For arbitrary $a\in\mathbb{R}_{>0}$, and any metric space
$(\mathcal{X},\tilde{\rho})$, we note that

$\int_{0}^{\infty}\sqrt{\log(N(\mathcal{X},a\tilde{\rho},r))}\mathrm{d}r=\int_{0}^{\infty}\sqrt{\log(N(\mathcal{X},\tilde{\rho},r/a))}\mathrm{d}r=a\int_{0}^{\infty}\sqrt{\log(N(\mathcal{X},\tilde{\rho},s))}\mathrm{d}s\,.$
(229)

∎

######  Lemma B.16 (Chaining).

Let $(\mathcal{X},\rho)$ be a metric space. Let $(Z_{x})_{x\in\mathcal{X}}$ be
a centered, sub-Gaussian process. Let $x_{0}\in\mathcal{X}$. Let
$k\in\mathbb{N}$. Then

$\mathbf{E}\\!\left[\sup_{x\in\mathcal{X}}\left|Z_{x}-Z_{x_{0}}\right|^{k}\right]\leq
c_{k}\left(\int_{0}^{\infty}\sqrt{\log(N(\mathcal{X},\rho,r))}\mathrm{d}r\right)^{k}\,.$
(230)

###### Proof.

E.g., [Tal21, chapter 2]. To be more precise, combine the results and
definitions of Theorem 2.7.13, Definition 2.7.3, Exercise 2.7.6 (b), (2.40),
(2.37), and Exercise 2.3.8 (b) in that reference. ∎

###### Proof of Theorem B.6.

Bias. As

$\mathbf{E}\\!\left[\left|\hat{g}-D_{\mathbf{v}}f\right|_{\infty}^{2}\right]\leq
c\left(\mathbf{E}\\!\left[\sup_{x\in[0,T]^{d}}\left|\hat{g}(x)-\mathbf{E}\\!\left[\hat{g}(x)\right]\right|^{2}\right]+\sup_{x\in[0,T]^{d}}\left|\mathbf{E}\\!\left[\hat{g}(x)\right]-D_{\mathbf{v}}f(x)\right|^{2}\right)$
(231)

we can use the same bound of the bias term as in Theorem B.3, i.e., (213),

$\sup_{x\in[0,T]^{d}}\left|\mathbf{E}\\!\left[\hat{g}(x)\right]-D_{\mathbf{v}}f(x)\right|\leq
c_{\ell}C_{\mathsf{egv}}C_{\mathsf{ker}}C_{\mathsf{cvr}}Lh^{\beta-s}\,.$ (232)

Variance. For the variance term, we have to bound the expectation of the
supremum of an empirical process,

$\mathbf{E}\\!\left[\sup_{x\in[0,T]^{d}}\left|\hat{g}(x)-\mathbf{E}\\!\left[\hat{g}(x)\right]\right|^{2}\right]=\mathbf{E}\\!\left[\sup_{x\in[0,T]^{d}}\left|\sum_{i=1}^{n}w_{i}(x)\epsilon_{i}\right|^{2}\right]\,.$
(233)

Write

$w_{i}(x)=\eta^{\\!\top}\\!B(x)^{-1}b_{i}(x)$ (234)

with

$\displaystyle b_{i}(x)$
$\displaystyle:=\psi\\!\left(\frac{x_{i}-x}{h}\right)K\\!\left(\frac{\left\|x_{i}-x\right\|_{2}}{h}\right)\,,$
(235) $\displaystyle\eta$
$\displaystyle:=\frac{T^{d}}{nh^{d+s}}D_{\mathbf{v}}\psi(0)\,.$ (236)

Then

$\displaystyle\left|\sum_{i=1}^{n}w_{i}(x)\epsilon_{i}\right|$
$\displaystyle\leq\left|\eta^{\\!\top}\\!B(x)^{-1}\sum_{i=1}^{n}b_{i}(x)\epsilon_{i}\right|$
(237)
$\displaystyle\leq\|\eta\|\left\|B(x)^{-1}\right\|_{\mathsf{op}}\left\|\sum_{i=1}^{n}b_{i}(x)\epsilon_{i}\right\|_{2}$
(238) $\displaystyle\leq
c_{\ell}\frac{T^{d}}{nh^{d+s}}C_{\mathsf{egv}}\left\|\sum_{i=1}^{n}b_{i}(x)\epsilon_{i}\right\|_{2}\,.$
(239)

Hence, we have to bound the expectation of

$\sup_{x\in[0,T]^{d}}\left\|\sum_{i=1}^{n}b_{i}(x)\epsilon_{i}\right\|_{2}\leq\sum_{j=1}^{N}\sup_{x\in[0,T]^{d}}\left|\sum_{i=1}^{n}e_{j}^{\\!\top}\\!b_{i}(x)\epsilon_{i}\right|\,.$
(240)

Fix a $j\in{\left\llbracket N\right\rrbracket}$ and define
$Z_{i}(x)=e_{j}^{\\!\top}\\!b_{i}(x)\epsilon_{i}$. We want to find an upper
bound for

$\mathbf{E}\\!\left[\sup_{x\in[0,T]^{d}}\left|\sum_{i=1}^{n}Z_{i}(x)\right|^{2}\right]\,.$
(241)

Because of SubGaussian and Lemma B.13, $Z(x)=\sum_{i=1}^{n}Z_{i}(x)$ is a sub-
Gaussian process on $[0,T]^{d}$ with metric

$\rho(x,x^{\prime})^{2}:=\sigma^{2}\sum_{i=1}^{n}\left(e_{j}^{\\!\top}\\!\left(b_{i}(x)-b_{i}(x^{\prime})\right)\right)^{2}\,.$
(242)

As $z\mapsto\psi(z)K(z)$ is smooth, it is Lipschitz-continuous with a constant
$\tilde{L}\in\mathbb{R}_{>0}$ depending only on $\ell,d,K$. If
$\max(\left\|x_{i}-x\right\|_{2},\left\|x_{i}-x^{\prime}\right\|_{2})\leq h$,
then

$\left|e_{j}^{\\!\top}\\!\left(b_{i}(x)-b_{i}(x^{\prime})\right)\right|\leq
c\tilde{L}h^{-1}\left\|x-x^{\prime}\right\|_{2}\,.$ (243)

If $\min(\left\|x_{i}-x\right\|_{2},\left\|x_{i}-x^{\prime}\right\|_{2})\leq
h\leq\max(\left\|x_{i}-x\right\|_{2},\left\|x_{i}-x^{\prime}\right\|_{2})$,
then

$\left|e_{j}^{\\!\top}\\!\left(b_{i}(x)-b_{i}(x^{\prime})\right)\right|\leq
c\tilde{L}\,.$ (244)

If $\min(\left\|x_{i}-x\right\|_{2},\left\|x_{i}-x^{\prime}\right\|_{2})\geq
h$, then

$\left|e_{j}^{\\!\top}\\!\left(b_{i}(x)-b_{i}(x^{\prime})\right)\right|=0\,.$
(245)

Hence, counting the number of indices with $x_{i}$ in a ball of radius $h$
yields

$\rho(x,x^{\prime})\leq
c\sqrt{C_{\mathsf{cvr}}nh^{d}T^{-d}}\tilde{L}\sigma\min\\!\left(1,\,\frac{\left\|x-\tilde{x}\right\|_{2}}{h}\right)\,.$
(246)

Now Lemma B.16 and Lemma B.15 yield

$\mathbf{E}\\!\left[\sup_{x\in[0,T]^{d}}\left(Z(x)-Z(x_{0})\right)^{2}\right]\leq
c_{d,\ell,K}C_{\mathsf{cvr}}nh^{d}T^{-d}\sigma^{2}\log(T/h)$ (247)

for an arbitrary $x_{0}\in[0,T]^{d}$. Using independence of $\varepsilon_{i}$,
we have

$\mathbf{E}\\!\left[Z(x_{0})^{2}\right]\leq\sum_{i=1}^{n}(e_{j}^{\\!\top}\\!b_{i}(x))^{2}\mathbf{E}\\!\left[\epsilon_{i}^{2}\right]\leq
c_{d,\ell,K}C_{\mathsf{cvr}}nh^{d}T^{-d}\sigma^{2}\,.$ (248)

Hence, with $\log(T/h)\geq c$ as $T\geq 2h$, we get

$\mathbf{E}\\!\left[\sup_{x\in[0,T]^{d}}Z(x)^{2}\right]\leq
c_{d,\ell,K}C_{\mathsf{cvr}}nh^{d}T^{-d}\sigma^{2}\log(T/h)\,.$ (249)

We finally obtain a bound on the variance term (233),

$\displaystyle\mathbf{E}\\!\left[\sup_{x\in[0,T]^{d}}\left|\sum_{i=1}^{n}w_{i}(x)\epsilon_{i}\right|^{2}\right]$
$\displaystyle\leq
c_{d,\ell}\left(\frac{T^{d}}{nh^{d+s}}C_{\mathsf{egv}}\right)^{2}\mathbf{E}\\!\left[\sup_{x\in[0,T]^{d}}\left|\sum_{i=1}^{n}Z_{i}(x)\right|^{2}\right]$
(250) $\displaystyle\leq
c_{d,\ell,K}\left(\frac{T^{d}}{nh^{d+s}}C_{\mathsf{egv}}\right)^{2}C_{\mathsf{cvr}}nh^{d}T^{-d}\sigma^{2}\log(T/h)$
(251) $\displaystyle\leq
c_{d,\ell,K}C_{\mathsf{cvr}}C_{\mathsf{egv}}^{2}\frac{\sigma^{2}T^{d}}{nh^{d+2s}}\log(T/h)\,.$
(252)

Together. Combining the bias and variance term yields

$\mathbf{E}\\!\left[\left|\hat{g}-D_{\mathbf{v}}f\right|_{\infty}^{2}\right]\leq
c_{d,\ell,K}\left(\left(C_{\mathsf{egv}}C_{\mathsf{cvr}}Lh^{\beta-s}\right)^{2}+C_{\mathsf{cvr}}C_{\mathsf{egv}}^{2}\frac{\sigma^{2}T^{d}}{nh^{d+2s}}\log(T/h)^{2}\right)\,.$
(253)

∎

###### Proof of Proposition B.9.

StrictKernel directly implies Kernel with $C_{\mathsf{ker}}$.

For Cover, we note that $\mathrm{B}^{d}(z,r)$ is a subset of a hypercube with
side length $2r$. For each dimension, this hypercube can only contain
$2n_{0}r/T+1$ points separated by a distance $T/n_{0}$. Hence

$\frac{1}{n}\sum_{i=1}^{n}\mathds{1}_{\\!\mathrm{B}^{d}(x,r)}\\!\\!\left(x_{i}\right)\leq\begin{cases}\frac{1}{n}&\text{if
}2n_{0}r/T<1\,,\\\
\frac{1}{n}\left(2n_{0}r/T+1\right)^{d}&\text{otherwise}\,.\end{cases}$ (254)

If $2n_{0}r/T\geq 1$,

$\displaystyle\frac{1}{n}\left(2n_{0}r/T+1\right)^{d}\leq\frac{1}{n}\left(4n_{0}r/T\right)^{d}\leq
4^{d}\left(r/T\right)^{d}\,.$ (255)

Hence, we can set $C_{\mathsf{cvr}}=4^{d}$.

Now consider Eigenvalue. Let $v\in\mathbb{D}_{d}$ and $x_{0}\in[0,T]^{d}$. We
want to find a lower bound on $v^{\\!\top}\\!B(x_{0})v$. Let
$h\in\mathbb{R}_{>0}$ be the bandwidth. Denote $z(x)=(x-x_{0})/h$ and
$z_{i}:=z(x_{i})$. Denote the square with side length $2h$ around $x_{0}$ as
$\mathcal{S}(x_{0}):=\mathcal{S}_{h}(x_{0}):=\bigtimes_{k=1}^{d}[\Pi_{k}x_{0}-h,\Pi_{k}x_{0}+h]$.
Denote indices of the grid points in $\mathcal{S}(x_{0})$ as
$\mathcal{I}(x_{0}):=\\{i\in{\left\llbracket n\right\rrbracket}\mid
x_{i}\in\mathcal{S}(x_{0})\\}$. From the definition of $B(x_{0})$, (175), we
obtain

$\displaystyle v^{\\!\top}\\!B(x_{0})v$
$\displaystyle=\frac{T^{d}}{nh^{d}}\sum_{i=1}^{n}\left(v^{\\!\top}\\!\psi\\!\left(z_{i}\right)\right)^{2}K(\left\|z_{i}\right\|_{2})$
(256)
$\displaystyle=\frac{T^{d}}{nh^{d}}\sum_{i\in\mathcal{I}(x_{0})}\left(v^{\\!\top}\\!\psi\\!\left(z_{i}\right)\right)^{2}K(\left\|z_{i}\right\|_{2})\,.$
(257)

The $x_{i}$ form a grid of side length $T/n_{0}$ in $\mathcal{S}(x_{0})$,
which has a volume of $(2h)^{d}$. Assume $2h\geq T/n_{0}$. Then, the Riemann
sum approximation for an $L$-Lipschitz continuous function
$g\colon\mathcal{S}(x_{0})\to\mathbb{R}$ yields