Random matrices based schemes for stable and robust nonparametric and
functional regression estimators.

Asma Ben Sabera and Abderrazek Karouia 111Corresponding author: Abderrazek
Karoui
Email addresses<EMAIL_ADDRESS>(A. Karoui),
<EMAIL_ADDRESS>(A. Ben Saber)
This work was supported in part by the DGRST research grant LR21ES10 and the
PHC-Utique research project 20G1503.

a University of Carthage, Department of Mathematics, Faculty of Sciences of
Bizerte, Tunisia.

Abstract— In the first part of this work, we develop a novel scheme for
solving nonparametric regression problems. That is the approximation of
possibly low regular and noised functions from the knowledge of their
approximate values given at some random points. Our proposed scheme is based
on the use of the pseudo-inverse of a random projection matrix, combined with
some specific properties of the Jacobi polynomials system, as well as some
properties of positive definite random matrices. This scheme has the
advantages to be stable, robust, accurate and fairly fast in terms of
execution time. In particular, we provide an $L_{2}$ as well as an
$L_{2}-$risk errors of our proposed nonparametric regression estimator.
Moreover and unlike most of the existing nonparametric regression estimators,
no extra regularization step is required by our proposed estimator. Although,
this estimator is initially designed to work with random sampling set of uni-
variate i.i.d. random variables following a Beta distribution, we show that it
is still works for a wide range of sampling distribution laws. Moreover, we
briefly describe how our estimator can be adapted in order to handle the
multivariate case of random sampling sets. In the second part of this work, we
extend the random pseudo-inverse scheme technique to build a stable and
accurate estimator for solving linear functional regression (LFR) problems. A
dyadic decomposition approach is used to construct this last stable estimator
for the LFR problem. Alaso, we give an $L_{2}-$risk error of our proposed LFR
estimator. Finally, the performance of the two proposed estimators are
illustrated by various numerical simulations. In particular, a real dataset is
used to illustrate the performance of our nonparametric regression estimator.

Keywords: Nonparametric regression, linear functional regression, Jacobi
polynomials, eigenvalues, random matrix, random pseudo-inverse.

##  1 Introduction

In the first part of this work, we consider the following nonparametric
regression model:

$Y_{i}=f(X_{i})+\varepsilon_{i}=\widetilde{f}(X_{i}),\quad 1\leq i\leq n,$ (1)

where $(X_{i})_{1\leq i\leq n}$ are random variables (or inputs) with
distribution $\rho_{X}$ and the noise terms $(\varepsilon_{i})_{1\leq i\leq
n}$ are i.i.d. real-valued centered random variables with variance
$\sigma^{2}.$ For a given complete metric space $\mathcal{X},$ a measurable
subset $\mathcal{Y}$ of $\mathbb{R}$ and a training set
$\\{(X_{i},Y_{i}),1\leq i\leq n\\}$ in $\mathcal{X}\times\mathcal{Y},$ we are
interested in the construction of an estimator $\widehat{f}$ such that for
$x\in\mathcal{X},$ $\widehat{f}(x)$ is a good estimate for the corresponding
regression function $f.$ Recall that if the $(X_{i},Y_{i})$ are drawn from a
joint probability measure $\rho$ on $\mathcal{X}\times\mathcal{Y},$ then the
regression function $f$ is given by

$f(x)=\int_{\mathcal{Y}}yd\rho(y|x)=\mathbb{E}(Y|X=x),\quad x\in\mathcal{X},$

For simplicity, we will assume that the $X_{i}$ take values in the interval
$I=[-1,1].$ Nonetheless, our results are easily adapted to a more general
compact interval of $\mathbb{R}.$ Also, we show how our nonparametric
regression estimator can be adapted in order to handle the more general case
of random sampling points $X_{i}\in\mathbb{R}^{d},\,d\geq 2.$

Note that the kernel ridge regression (KRR) scheme is a widely used scheme for
solving the regression problem modeled by (1). For more details on the KRR
based estimator, the reader is referred for example to [17, 18]. Other
estimators for nonparametric regression problem have been proposed in the
literature. To cite but a few, F. Comte and her co-authors have proposed
estimators based on projections associated with some classical orthogonal
polynomials [7]. Also, in the literature, there exists another class of
convolution kernels with variable bandwidths based robust estimators for the
nonparametric regression problem, see for example [2, 21].

Our proposed estimator for solving problem (1) is based on combining the
Random Sampling Consensus (RANSAC) iterative algorithm with a random
pseudoinverse based scheme for nonparametric regression. That is at each
iteration of the RANSAC algorithm, a stable nonparametric regression estimator
is constructed by using the pseudoinverse of a random projection matrix,
associated with Jacobi orthogonal polynomials. Note that the RANSAC procedure,
see for example [15] consists in randomly selecting $m<n$ measurements from
model (1). From these $m$ training data points $(X_{i},Y_{i}),$ construct an
estimator and compute the mean squared prediction error estimation. This
procedure is iterated till an estimator with high consensus is obtained. We
should mention that in the literature, pseudoinverse scheme has been
successfully used in a wide range of applications from different scientific
area, including machine learning, see [13] and the references therein. The
usual way of computing the pseudoinverse of an $n\times m$ real coefficients
matrix is to combine the Moore-Penrose inverse and the SVD techniques. It is
well known that the time complexity of this last procedure is
$O\big{(}\min(nm^{2},nm^{2})\big{)}.$ Nonetheless in [13], the authors have
proposed an accurate and Fast PesudoInverse algorithm for sparse feature
matrices.

For the sake of simplicity, the emphasis will be on the construction of a
stable estimator for problem (1) by using a random training set
$\\{(X_{i},Y_{i}),\,\,1\leq i\leq n\\}.$ More precisely, for real parameters
${\displaystyle\alpha,\beta\geq-\frac{1}{2},}$ we consider the weight function
$\omega_{\alpha,\beta}(x)=(1-x)^{\alpha}(1+x)^{\beta}$ and the associated
weighted $L^{2}_{\omega}(I)-$space. This last space is given by the set of
real valued functions $f$ that are measurable on $I$ and satisfying
${\displaystyle\|f\|_{\omega}^{2}=\int_{I}|f(x)|^{2}\omega_{\alpha,\beta}(x)\,dx<+\infty.}$
For a positive integer $n\in\mathbb{N},$ let $\\{X_{i},\,1\leq i\leq n\\}$ be
a sampling set of i.i.d random variables following the beta distribution
$B(\alpha+1,\beta+1)$ over $I.$ Note that this last condition is not mandatory
and it is used in order to simplify the different theoretical results related
to our nonparametric regression estimator. In fact, we show how this estimator
can be adapted in order to handle random sampling set drawn from a more
general probability distribution law or even from an unknown probability
distribution.

Our nonparametric regression estimator is briefly described as follows. We use
the notation $\boldsymbol{u}^{\prime}$ for the transpose of $\boldsymbol{u},$
then for a positive integer $N\in\mathbb{N}$ with $N+1\leq n,$
$\widehat{f}_{n,N}$ is given in terms of the first $N+1$ normalized Jacobi
polynomials $\widetilde{P}^{(\alpha,\beta)}_{k}(x),\ k=0,\ldots,N.$ More
precisely, we have

$\widehat{f}_{n,N}(x)=\sum_{k=1}^{N+1}\widehat{c}_{k}\widetilde{P}^{(\alpha,\beta)}_{k-1}(x),\quad
x\in I.$ (2)

Here, the expansion coefficients vector
$\widehat{\boldsymbol{C}}=(\widehat{c}_{1},\ldots,\widehat{c}_{N+1})^{\prime}$
is given by the over-determined system

$B_{N}\cdot\widehat{\boldsymbol{C}}=\Big{[}\frac{1}{\sqrt{n}}\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{j})\Big{]}_{j,k}\cdot\widehat{\boldsymbol{C}}=\mathbf{Z}_{n},$
(3)

where $\mathbf{Z}_{n}=\frac{1}{\sqrt{n}}(Y_{1},\ldots,Y_{n})^{\prime}.$ We
refer to the $n\times(N+1)$ matrix $B_{N}$ as the random projection matrix.
One of the main result of this work is to prove that with high probability,
the $N+1$ dimensional positive definite matrix $A_{N}=B_{N}^{\prime}B_{N}$ is
invertible and well conditioned. More precisely, we show that with high
probability, the $2-$norm condition number of $A_{N},$ denoted by
$\kappa_{2}(A_{N})$ is bounded by a fairly small constant. This ensures the
stability of the estimator $\widehat{f}_{n,N}.$ Consequently, the expansion
coefficients vector $\widehat{\boldsymbol{C}}$ is simply given by

$\widehat{\boldsymbol{C}}=\big{(}B_{N}^{\prime}B_{N}\big{)}^{-1}B_{N}^{\prime}\mathbf{Z}_{n}.$
(4)

Also, we provide an upper bound for the regression error
$\|f-\widehat{f}_{n,N}\|_{\omega},$ that holds with high probability. In
particular, we provide more precise regression errors in the special cases
where $f$ has a Lipschitzian $p-$th order derivative or $f$ is a
$c-$bandlimited function, for some $c>0.$ Moreover, we provide a bound for the
$L^{2}-$risk error
${\displaystyle\mathbb{E}\Big{[}\|f-\widehat{f}_{n,N}\|^{2}_{\omega}\Big{]}}.$

In the second part of this work, we push forward the random pseudo-inverse
based scheme for solving a linear functional regression (LFR) problem, that
was introduced in [4]. This is done by combining the random pseudo-inverse
scheme with a dyadic decomposition technique. This combination yields a stable
and a highly accurate estimator for solving the LFR problem. Recall that for a
compact interval $J,$ the LFR model is given as follows

$Y_{i}=\int_{J}X_{i}(s)\,\beta_{0}(s)\,ds+\varepsilon_{i},\quad i=1,\ldots,n.$
(5)

Here, the $X_{i}(\cdot)\in L^{2}(J)$ are random functional predictors, the
$\varepsilon_{i}$ are i.i.d. centered white noise independent of the
$X_{i}(\cdot)$ and $\beta_{0}(\cdot)\in L^{2}(J)$ is the unknown slope
function to be recovered. As usual, we assume that the $X_{i}(\cdot)$ are
given by their series expansions with respect to an orthonormal set of
$L^{2}(J).$ That is for an orthonormal family
${\displaystyle\\{\varphi_{k}(\cdot),\,k\geq 1\\}}$ of $L^{2}(J),$ we have
${\displaystyle X_{i}(s)=\sum_{k\geq 1}\xi_{k}Z_{i,k}\varphi_{k}(s),}$ where
the $Z_{i,k}$ are i.i.d. centered random variables with variance
$\sigma_{Z}^{2}$ and $(\xi_{k})_{k\geq 1}$ is a deterministic sequence of
$\mathbb{R}\setminus\\{0\\}.$ Our proposed stable random pseudo-inverse based
estimator $\widehat{\beta}_{n,N}(\cdot)$ of
${\displaystyle\beta_{0}(\cdot)=\sum_{k\geq 1}c_{k}\varphi_{k}(\cdot)}$ is
briefly described as follows. For an integer $N\geq 2,$ we let
$K_{N}=\Big{[}\frac{\log N}{\log 2}\Big{]}+1$ and we consider the dyadic
decomposition of the set $[[1,N]]=\\{1,\ldots,N\\}$ into $K_{N}$ subsets
$I_{k}=[[2^{k-1}+1,\min(2^{k},N)]],$ $1\leq k\leq K_{N}.$ On each subset
$I_{k},$ we consider the random matrices $G_{k},F_{k}$ given by

$G_{k}=\displaystyle\left[\frac{1}{n}\sum_{i=1}^{n}\xi_{j}\xi_{l}Z_{i,j}Z_{i,l}\right]_{j,l\in
I_{k}}=F_{k}^{\prime}F_{k},\quad
F_{k}=\frac{1}{\sqrt{n}}\big{[}\xi_{j}Z_{i,j}\big{]}_{\begin{subarray}{c}\scriptstyle
1\leq i\leq N\\\ \scriptstyle j\in I_{k}\end{subarray}}.$

Under the previous notation, our proposed estimator for the LFR problem (5) is
given by

$\widehat{\beta}_{n,N}(s)=\sum_{k=1}^{K_{N}}\sum_{j\in
I_{k}}\widehat{c}_{n,j}\varphi_{j}(s),\quad\widehat{\boldsymbol{c}}_{n,k}=[\widehat{c}_{n,j}]_{j\in
I_{k}}^{\prime}=G_{k}^{-1}\Big{(}F_{k}^{\prime}\frac{1}{\sqrt{n}}[Y_{i}^{k}]^{\prime}_{1\leq
i\leq n}\Big{)},\quad 1\leq k\leq K_{N}.$

Note that the $\boldsymbol{Y}^{k}=[Y_{i}^{k}]_{1\leq i\leq n}$ are obtained
from the output vector $\boldsymbol{Y}=[Y_{i}]_{1\leq i\leq n}$ by
substituting in (5), $X_{i}(s)$ with its projection ${X_{i}^{k}(s)=\sum_{j\in
I_{k}}\xi_{j}Z_{i,j}\varphi_{j}(s)}.$ In particular, we provide an upper bound
for the mean squared estimation error
$\|\widehat{\beta}_{n,N}(\cdot)-\beta_{0}(\cdot)\|_{2}^{2},$ that holds with
high probability. This upper error bound involves the quantity
${\displaystyle\sum_{k=1}^{K_{N}}\kappa_{2}(G_{k}).}$ In practice, the
deterministic sequence $\boldsymbol{\xi}=(\xi_{j})_{j}$ satisfies a decay
condition of the type $\xi_{j}\asymp j^{-s}$ for some $s>0.$ Under this last
assumption, it is easy to check that our combined dyadic decomposition and
random pseudo-inverse scheme is stable in the sense that
${\displaystyle\sum_{k=1}^{K_{N}}\kappa_{2}(G_{k})\lesssim 2^{s}\log N,}$
which is a relatively small cumulative condition number. This ensures the
stability of our proposed LFR estimator $\widehat{\beta}_{n,N}(\cdot).$ Also,
we give an estimate for the $L^{2}-$risk error
$\mathbb{E}\Big{[}\|\beta_{0}(\cdot)-\widetilde{\beta}_{N,M}(\cdot)\|^{2}_{2}\Big{]},$
where $\widetilde{\beta}_{N,M}(\cdot)$ is a truncated version of the estimator
$\widehat{\beta}_{n,N}(\cdot).$

This work is organized as follows. In section 2, we give some mathematical
preliminaries that are frequently used in this work. In section 3, we study
our random pseudo-inverse based estimator for the stable and robust
approximation of the regression function $f$ given by (1). In section 4, we
give an error analysis of our combined dyadic decomposition, random pseudo-
inverse based estimator for the LFR problem (5). In section 5, we give some
numerical simulations that illustrate the different results of this work. In
particular, we illustrate the performance of our nonparametric regression by
applying it on a real dataset given by the daily Covid-19 cases of some
world’s countries.

##  2 Mathematical preliminaries

We first recall that for two real numbers $\alpha,\beta>-1,$ the classical
Jacobi polynomials $P^{(\alpha,\beta)}_{k}$ are defined for $x\in[-1,1],$ by
the following Rodrigues formula,

$P^{(\alpha,\beta)}_{k}(x)=\frac{(-1)^{k}}{2^{k}k!}\frac{1}{\omega_{\alpha,\beta}(x)}\frac{\mbox{d}^{k}}{\mbox{d}x^{k}}\Big{(}\omega_{\alpha,\beta}(x)(1-x^{2})^{k}\Big{)},\quad\omega_{\alpha,\beta}(x)=(1-x)^{\alpha}(1+x)^{\beta},$
(6)

with

$P^{(\alpha,\beta)}_{k}(1)={k+\max(\alpha,\beta)\choose{k}}=\frac{\Gamma(k+\max(\alpha,\beta)+1)}{k!\,\Gamma(\max(\alpha,\beta)+1)}.$

Here, $\Gamma(\cdot)$ denotes the usual Gamma function. In the sequel, we let
$\widetilde{P}^{(\alpha,\beta)}_{k}$ denote the normalized Jacobi polynomial
of degree $k,$ given by

$\widetilde{P}^{(\alpha,\beta)}_{k}(x)=\frac{1}{\sqrt{h_{k}^{\alpha,\beta}}}P^{(\alpha,\beta)}_{k}(x),\quad
h_{k}^{\alpha,\beta}=\frac{2^{\alpha+\beta+1}\Gamma(k+\alpha+1)\Gamma(k+\beta+1)}{k!(2k+\alpha+\beta+1)\Gamma(k+\alpha+\beta+1)}.$
(7)

In this case, we have

$\|\widetilde{P}^{(\alpha,\beta)}_{k}\|^{2}_{\omega}=\int_{-1}^{1}(\widetilde{P}^{(\alpha,\beta)}_{k}(y))^{2}\omega_{\alpha,\beta}(y)\,dy=1.$
(8)

The following lemma provides us with some useful upper bounds for the
normalized Jacobi polynomials $\widetilde{P}^{(\alpha,\beta)}_{k},$ for $k\geq
2.$

######  Lemma 1.

For $\alpha,\beta\geq-\frac{1}{2},$ let $\mu=\max(\alpha,\beta),$
$c_{\alpha,\beta}=\frac{\alpha+\beta+1}{2}$ and let

$\eta_{\alpha,\beta}=\frac{\exp\big{(}\frac{2\max(\mu,0)}{12}+\frac{\max(\mu^{2}+\alpha\beta,0)}{8}\big{)}}{2^{(\alpha+\beta)/2}\Gamma(\mu+1)}$
(9)

Then, we have

$\max_{x\in[-1,1]}|\widetilde{P}^{(\alpha,\beta)}_{k}(x)|\leq\eta_{\alpha,\beta}\gamma_{k}^{\alpha,\beta},\qquad\gamma_{k}^{\alpha,\beta}=k^{\mu}\sqrt{k+c_{\alpha,\beta}},\quad\forall\,k\geq
2.$ (10)

Proof: It is well known that under the previous condition on $\alpha,\beta,$
we have

$\max_{x\in[-1,1]}|P^{(\alpha,\beta)}_{k}(x)|=\max\big{(}|P^{(\alpha,\beta)}_{k}(1)|,|P^{(\alpha,\beta)}_{k}(-1)|\big{)}=\frac{\Gamma(k+\mu+1)}{k!\,\Gamma(\mu+1)},\quad\mu=\max(\alpha,\beta).$
(11)

On the other hand, from [20], for an integer $n\geq 2$ and for two real
numbers $a,b>-1,$ we have

$\frac{\Gamma(n+a)}{\Gamma(n+b)}\leq C_{n}^{a,b}n^{a-b},\quad
C_{n}^{a,b}=\exp\Big{(}\frac{a-b}{2(n+b-1)}+\frac{1}{12(n+a-1)}+\frac{(a-1)(b-1)}{n}\Big{)}.$
(12)

By combining (7), (11) and (12), one gets

$\max_{x\in[-1,1]}|\widetilde{P}^{(\alpha,\beta)}_{k}(x)|\leq\frac{\sqrt{C_{k}^{\mu+1,1}C_{k}^{\alpha+\beta+1,\alpha+\beta-\mu+1}}}{\Gamma(\mu+1)2^{(\alpha+\beta)/2}}\,\,k^{\mu}\,\,\sqrt{k+\frac{(\alpha+\beta+1)}{2}}.$

But from (12), it can be easily checked that for $k\geq 2,$ we have

$C_{k}^{\mu+1,1}C_{k}^{\alpha+\beta+1,\alpha+\beta-\mu+1}\leq\eta_{\alpha,\beta}.$

By combining the previous two inequalities, we get

$\max_{x\in[-1,1]}|\widetilde{P}^{(\alpha,\beta)}_{k}(x)|\leq\eta_{\alpha,\beta}\gamma_{k}^{\alpha,\beta}=\eta_{\alpha,\beta}k^{\mu}\,\,\sqrt{k+\frac{(\alpha+\beta+1)}{2}},\quad
k\geq 2.$ (13)

Next, we recall from [11], the following Gershgorin circle theorem that is
used for the estimation of $2-$norm condition numbers of various random
matrices used in this work.

######  Theorem 1 (Gershgorin circle theorem).

Let $A=[a_{ij}]$ be a complex $n\times n$ matrix. For $1\leq i\leq n,$ let
${\displaystyle R_{i}=\sum_{j\neq i}|a_{ij}|.}$ Then, every eigenvalue of $A$
lies within at least one of the discs $D(a_{ii},R_{i}).$

The following McDiarmid’s concentration inequality is also used in the
spectral analysis of the random matrices of the next section.

######  Theorem 2 (McDiarmid’s inequality).

Let $X_{1},\ldots,X_{n}$ be independent random variables with respective
ranges in metric spaces $\boldsymbol{\chi}_{1},\ldots,\boldsymbol{\chi}_{n}.$
Let
$f:\boldsymbol{\chi}_{1}\times\cdots\times\boldsymbol{\chi}_{n}\rightarrow\mathbb{C}$
be a function satisfying the bounded differences property

$\Big{|}f(x_{1},\ldots,x_{i-1},x_{i},x_{i+1},\ldots,x_{n})-f(x_{1},\ldots,x_{i-1},x^{\prime}_{i},x_{i+1},\ldots,x_{n})\Big{|}\leq
c_{i},\quad 1\leq i\leq n.$

Then, for any $\epsilon>0,$ we have

$\mathbb{P}\Big{(}f(X_{1},\ldots,X_{n})-\mathbb{E}\big{[}f(X_{1},\ldots,X_{n})\big{]}\leq\epsilon\Big{)}\geq
1-\exp\left(\frac{2\epsilon^{2}}{\sum_{i=1}^{n}c_{i}^{2}}\right).$

##  3  Random projection on Jacobi polynomials system and nonparametric
regression estimator

The main part of this section is devoted to the study of the stability and
error analysis of our nonparametric regression estimator in the framework of a
sampling set drawn from a Beta distribution. Then, at the end of this section,
we show how our estimator can be adapted in order to handle random sampling
sets drawn from more general distribution laws or even from an unknown
distribution law.

We first consider two real numbers $\alpha,\beta\geq-\frac{1}{2}$ and two
positive integers $n\geq N+1$, then we let $\\{X_{i},\,1\leq i\leq n\\}$ be a
sampling set of i.i.d random variables following the beta
$B(\alpha+1,\alpha+1)$ distribution over the interval $I.$ Then, under the
assumption that the regression function $f$ (solution of the nonparametric
regression problem (1)) is well approximated by its projection
$f_{N}=\pi_{N}f$ over the subspace
$\mathcal{H}_{N}=\mbox{Span}\\{\widetilde{P}^{(\alpha,\beta)}_{k},\,0\leq
k\leq N\\},$ our random pseudo-inverse based estimator $\widehat{f}_{n,N}$ of
$f\in L^{2}(I)$ is given by

$\widehat{f}_{n,N}(x)=\sum_{k=1}^{N+1}\widehat{c}_{k}\widetilde{P}^{(\alpha,\beta)}_{k-1}(x),\quad
x\in I.$ (14)

We require that this estimator is stable and robust. This estimator is given
as the Least square solution of the over-determined system

$\sum_{k=1}^{N+1}c_{k}\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{j})=Y_{j},\quad
j=1,\ldots,n.$ (15)

Note that if the regression function $f\in L^{2}_{\omega}(I),$ then from the
uniqueness of its expansion coefficients with respect to the orthonormal basis
$\\{\widetilde{P}^{(\alpha,\beta)}_{k}(\cdot),\,k\geq 0\\},$ for any integer
$N\in\mathbb{N},$ its first $N$ coefficients $c_{i}(f)$ do not depend on the
specific values of $y_{j}=f(X_{j}),$ at a specific random sampling set
$\\{X_{j},j=1,\ldots,n\\}.$ These coefficients are uniquely given by

$c_{k}(f)=<\widetilde{f},\widetilde{P}^{(\alpha,\beta)}_{k-1}>=\int_{I}\widetilde{f}(x)\widetilde{P}^{(\alpha,\beta)}_{k-1}(x)\omega_{\alpha,\beta}(x)\,dx,\quad
1\leq k\leq N+1.$

Note that by re-scaling by a factor ${\frac{1}{\sqrt{n}},}$ the over-
determined system (15) is written as

$B_{N}\cdot\boldsymbol{C}=\boldsymbol{Y},\quad
B_{N}=\Big{[}\frac{1}{\sqrt{n}}\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{j})\Big{]}_{1\leq
j\leq n,1\leq k\leq N+1},\quad\boldsymbol{Y}=[Y_{j}]^{\prime}_{1\leq j\leq
n}.$ (16)

We show that under some conditions on $f$ as well as on the integers $n,N,$
the least square solution of the previous system is a good approximation of
the set of the first $N+1$ expansion coefficients of $f$. For this purpose, we
define an $(N+1)\times(N+1)$ random Jacobi projection matrix $A_{N}$ by

$A_{N}=B_{N}^{\prime}B_{N}=\left[\frac{1}{n}\sum_{j=1}^{n}\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{j})\widetilde{P}^{(\alpha,\beta)}_{l-1}(X_{j})\right]_{1\leq
k,l\leq N+1}=\sum_{j=1}^{n}\mathbf{D}_{j},$ (17)

Note that the random matrix $\mathbf{D}_{j}$ is also positive semidefinite. In
fact, $\mathbf{D}_{j}=b_{j}^{\prime}b_{j},$ where $b_{j}$ is the
$1\times(N+1)$ matrix given by ${\displaystyle
b_{j}=\frac{1}{\sqrt{n}}\Big{[}\widetilde{P}^{(\alpha,\beta)}_{0}(X_{j})\cdots\widetilde{P}^{(\alpha,\beta)}_{N}(X_{j})\Big{]}}.$
Consequently, for any $X\in\mathbb{R}^{N+1},$ we have

$X^{\prime}\cdot\mathbf{D}_{j}X=X^{\prime}b_{j}^{\prime}\cdot b_{j}X\geq 0.$

That is $\mathbf{D}_{j}$ is a positive semi-definite matrix. An important
property of the matrix $A_{N}$ is that under some conditions on the parameters
$\alpha,\beta,n,N,$ this matrix is well conditioned. To prove this first main
result, we need the following proposition that provides us with an upper and a
lower bound for $\lambda_{\max}(A_{N})$ and $\lambda_{\min}(A_{N}),$
respectively.

######  Proposition 1.

Under the notation of Lemma 1, for any two positive integers $n>N\geq 2,$
satisfying

$0.63\,n>m_{\alpha,\beta}^{2}\cdot(N+1)^{2\mu+2}\,\log(N+1),$ (18)

where ${\displaystyle
m_{\alpha,\beta}^{2}=\frac{1+\frac{1}{2}\sqrt{\frac{c_{\alpha,\beta}}{2}}}{\mu+3/2}\eta^{2}_{\alpha,\beta}}$
and $\eta_{\alpha,\beta}$ is as given by (9). Then, we have

$\mathbb{E}\big{(}\lambda_{\max}(A_{N})\big{)}\leq
1.72+m_{\alpha,\beta}^{2}\frac{(N+1)^{2\mu+2}\,\log(N+1)}{n}$ (19)

and

$\mathbb{E}\big{(}\lambda_{\min}(A_{N})\big{)}\geq
0.63-m_{\alpha,\beta}^{2}\frac{(N+1)^{2\mu+2}\,\log(N+1)}{n}.$ (20)

Proof: We first note that from (7), we have

$\displaystyle\max_{x\in[-1,1]}|\widetilde{P}^{(\alpha,\beta)}_{0}(x)|$
$\displaystyle=$
$\displaystyle\big{(}2^{\alpha+\beta+1}B(\alpha+1,\beta+1)\big{)}^{-1/2}$
$\displaystyle\max_{x\in[-1,1]}|\widetilde{P}^{(\alpha,\beta)}_{1}(x)|$
$\displaystyle=$
$\displaystyle(1+\mu)\sqrt{\frac{3+\alpha+\beta}{(\alpha+1)(\beta+1)}}\big{(}2^{\alpha+\beta+1}B(\alpha+1,\beta+1)\big{)}^{-1/2}.$

On the other hand, if $\gamma_{k}^{\alpha,\beta}$ is as given by (10), then it
is easy to check that for $\mu\geq-\frac{1}{2}$ and for
$\alpha,\beta\geq-\frac{1}{2},$ the function $x\rightarrow
x^{\mu}\sqrt{x+c_{\alpha,\beta}}$ is increasing on $[2,+\infty).$
Consequently, we have

$\displaystyle\sum_{k=2}^{N}\gamma_{k}^{\alpha,\beta}$ $\displaystyle\leq$
$\displaystyle\int_{2}^{N+1}x^{\mu}\sqrt{x+c_{\alpha,\beta}}\,dx\leq\int_{2}^{N+1}x^{\mu+\frac{1}{2}}\Big{(}1+\frac{1}{2}\sqrt{\frac{c_{\alpha,\beta}}{2}}\Big{)}\,dx$
$\displaystyle\leq$
$\displaystyle\frac{\Big{(}1+\frac{1}{2}\sqrt{\frac{c_{\alpha,\beta}}{2}}\Big{)}}{\mu+3/2}(N+1)^{\mu+3/2}.$

Hence, by using (10) and straightforward computation, one gets

$\max_{x\in[-1,1]}|\widetilde{P}^{(\alpha,\beta)}_{0}(x)|+\max_{x\in[-1,1]}|\widetilde{P}^{(\alpha,\beta)}_{1}(x)|+\cdots+\max_{x\in[-1,1]}|\widetilde{P}^{(\alpha,\beta)}_{N}(x)|\leq\eta_{\alpha,\beta}\int_{N}^{N+1}x^{\mu}\sqrt{x+c_{\alpha,\beta}}\,dx.$

The previous two inequalities give us

$\displaystyle\max_{x\in[-1,1]}|\widetilde{P}^{(\alpha,\beta)}_{N}(x)|\sum_{k=1}^{N+1}\Big{(}\max_{x\in[-1,1]}|\widetilde{P}^{(\alpha,\beta)}_{k-1}(x)|\Big{)}\leq\eta^{2}_{\alpha,\beta}N^{\mu}\sqrt{N+c_{\alpha,\beta}}\frac{\Big{(}1+\frac{1}{2}\sqrt{\frac{c_{\alpha,\beta}}{2}}\Big{)}}{\mu+3/2}(N+1)^{\mu+3/2}$
(21)
$\displaystyle\qquad\qquad\quad\qquad\qquad\qquad\leq\frac{\Big{(}1+\frac{1}{2}\sqrt{\frac{c_{\alpha,\beta}}{2}}\Big{)}}{\mu+3/2}\eta_{\alpha,\beta}^{2}(N+1)^{2\mu+2}=m_{\alpha,\beta}^{2}(N+1)^{2\mu+2}.$

Next, let $\mathbf{D}_{j}$ be the positive semi-definite matrix given by (17).
Since the random samples $X_{i}$ follow the Beta-distribution on $[-1,1]$ with
parameters $(\alpha+1,\beta+1),$ then we have
$\mathbb{E}\left(\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{j})\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{j})\right)=\delta_{kl}.$
That is
${\displaystyle\mathbb{E}\Big{(}\sum_{j=1}^{n}\mathbf{D}_{j}\Big{)}=I_{N+1},}$
where $I_{N+1}$ is the identity matrix of dimension $N+1.$ Hence, we have

$\Lambda_{\min}=\lambda_{\min}\Big{(}\mathbb{E}\big{(}\sum_{j=1}^{n}\mathbf{D}_{j}\big{)}\Big{)}=1,\quad\Lambda_{\max}=\lambda_{\max}\Big{(}\mathbb{E}\big{(}\sum_{j=1}^{n}\mathbf{D}_{j}\big{)}\Big{)}=1.$

Also, from Gershgorin circle theorem, we have

$\lambda_{\max}(\mathbf{D}_{j})\leq\frac{1}{n}\max_{x\in[-1,1]}|\widetilde{P}^{(\alpha,\beta)}_{N}(x)|\sum_{k=1}^{N+1}\Big{(}\max_{x\in[-1,1]}|\widetilde{P}^{(\alpha,\beta)}_{k-1}(x)|\Big{)}\leq\frac{1}{n}m_{\alpha,\beta}^{2}(N+1)^{2\mu+2}=\frac{L_{N}(\alpha,\beta)}{n}.$
(22)

Moreover, from [19], we have

$\mathbb{E}\Big{(}\lambda_{\min}(\sum_{j=1}^{n}\mathbf{D}_{j}\big{)}\Big{)}\geq
0.63\Lambda_{\min}-\frac{L_{N}(\alpha,\beta)}{n}\log(N+1),\quad\mathbb{E}\Big{(}\lambda_{\max}\big{(}\sum_{j=1}^{n}\mathbf{D}_{j}\big{)}\Big{)}\leq
1.72\Lambda_{\max}+\frac{L_{N}(\alpha,\beta)}{n}\log(N+1).$ (23)

This concludes the proof of the proposition.

The following theorem provides us with an upper bound for the actual condition
number $\kappa_{2}(A_{N}).$

######  Theorem 3.

Under the previous notation and the same conditions as proposition 1, for any
$0<\delta<1,$ we have

$\kappa_{2}(A_{N})\leq\frac{1.72+m_{\alpha,\beta}^{2}(N+1)^{2\mu+2}\,\Big{(}\frac{\log(N+1)}{n}+\sqrt{\frac{2}{n}\log\big{(}\frac{2}{\delta}\big{)}}\Big{)}}{0.63-m_{\alpha,\beta}^{2}(N+1)^{2\mu+2}\,\Big{(}\frac{\log(N+1)}{n}+\sqrt{\frac{2}{n}\log\big{(}\frac{2}{\delta}\big{)}}\Big{)}}$
(24)

with probability at least $(1-\delta)^{2}.$

Proof: Given $x=(x_{1},\cdots,x_{n})\in[-1,1]^{n}$, we consider the
$(N+1)\times(N+1)$ matrix with entries
$\frac{1}{n}\sum_{j=1}^{n}\widetilde{P}^{(\alpha,\beta)}_{k-1}(x_{j})\widetilde{P}^{(\alpha,\beta)}_{l-1}(x_{j})$.
We use the notation $A_{N}(x)$ for such a matrix. We denote by
$\lambda_{k}(A_{N}(x))$, $1\leq k\leq N+1,$ the eigenvalues of $A_{N}(x)$
arranged in decreasing order. For $1\leq k\leq N+1,$ we use McDiarmid’s
concentration inequality for the $n-$variate mapping

$x\mapsto\phi_{k}(x)=\lambda_{k}(A_{N}(x))$

We prove that the previous mapping $\phi_{k}$ satisfies the bounded
differences assumption. That is, when only one of the $n$ coordinates differs
between $x$ and $x^{\prime}$, then

$|\phi_{k}(x)-\phi_{k}(x^{\prime})|\leq\frac{2}{n}m_{\alpha,\beta}^{2}(N+1)^{2\mu+2}=2\frac{L_{N}(\alpha,\beta)}{n}.$
(25)

Let us take this for granted for now and we will prove it later on. It follows
from McDiarmid’s inequality that

$\mathbb{P}\Big{(}\phi_{k}(X_{1},\cdots,X_{n})-\mathbb{E}\phi_{k}(X_{1},\cdots,X_{n})\geq\varepsilon\Big{)}\leq\exp\frac{-n\varepsilon^{2}}{2m_{\alpha,\beta}^{4}(N+1)^{4\mu+4}}=\exp\frac{-n\varepsilon^{2}}{2L_{N}^{2}(\alpha,\beta)}$
(26)

and

$\mathbb{P}\big{(}\phi_{k}(X_{1},\cdots,X_{n})-\mathbb{E}\phi_{k}(X_{1},\cdots,X_{n})\leq-\varepsilon\big{)}\leq\exp\frac{-n\varepsilon^{2}}{2L_{N}^{2}(\alpha,\beta)}.$
(27)

For the value of $k=1$, (26) gives us

$\mathbb{P}\big{(}\lambda_{\max}(A_{N}))-\mathbb{E}(\lambda_{\max}(A_{N}))\leq\varepsilon\big{)}\geq
1-\exp\frac{-n\varepsilon^{2}}{2L_{N}^{2}(\alpha,\beta)}$ (28)

Moreover, for the value of $k=N+1$, (27) gives us

$\mathbb{P}\big{(}\lambda_{\min}(A_{N})-\mathbb{E}(\lambda_{\min}(A_{N}))\geq-\varepsilon\big{)}\geq
1-\exp\frac{-n\varepsilon^{2}}{2L_{N}^{2}(\alpha,\beta)}$ (29)

By combining (23) and (28), one gets

$\lambda_{\max}(A_{N})\leq
1.72+\frac{L_{N}(\alpha,\beta)}{n}\log(N+1)+\varepsilon$ (30)

with probability at least
$1-\exp\frac{-n\varepsilon^{2}}{2L_{N}^{2}(\alpha,\beta)}$. Here,
$L_{N}(\alpha,\beta)$ is given by (22). Also, by combining (23) and (29), one
gets

$\lambda_{\min}(A_{N})\geq\big{(}0.63-\frac{L_{N}(\alpha,\beta)}{n}\log(N+1)\big{)}-\varepsilon$
(31)

with probability at least
$1-\exp\frac{-n\varepsilon^{2}}{2L_{N}^{2}(\alpha,\beta)}$. Hence, by
combining (30) and (31) , one gets

$\kappa_{2}(A_{N})\leq\frac{1.72+\frac{L_{N}(\alpha,\beta)}{n}\log(N+1))+\varepsilon}{0.63-\frac{L_{N}(\alpha,\beta)}{n}\log(N+1)-\varepsilon}$

with probability at least
${\displaystyle\Big{(}1-\exp\frac{-n\varepsilon^{2}}{2L_{N}^{2}(\alpha,\beta)}\Big{)}^{2}}$.
Finally to get (24), it suffices to let
${\displaystyle\delta=\exp\left(-\frac{n\varepsilon^{2}}{2L_{N}^{2}(\alpha,\beta)}\right)}.$

Let us now prove the bounded differences condition (25). We have

$\displaystyle|\phi_{k}(x)-\phi_{k}(x^{\prime})|$ $\displaystyle=$
$\displaystyle|\phi_{k}(x_{1},\ldots,x_{i-1},x_{i},x_{i+1}\cdots,x_{n})-\phi_{k}(x_{1},\ldots,x_{i-1},x_{i}^{\prime},x_{i+1}\cdots,x_{n})|$
$\displaystyle=$
$\displaystyle\lambda_{k}(A_{N}(x))-\lambda_{k}(A_{N}(x^{\prime}))$

Let $E$ be the $(N+1)\times(N+1)$ matrix with entries

$E_{k,l}=\frac{1}{n}\left(\widetilde{P}^{(\alpha,\beta)}_{k-1}(x_{i})\widetilde{P}^{(\alpha,\beta)}_{l-1}(x_{i})-\widetilde{P}^{(\alpha,\beta)}_{k-1}(x^{\prime}_{i})\widetilde{P}^{(\alpha,\beta)}_{l-1}(x^{\prime}_{i})\right),$

so that $A_{N}(x)=A_{N}(x^{\prime})+E$. From Weyl’s perturbation theorem of
the spectrum of a perturbed Hermitian matrix, see for example [5], we have

$\lambda_{k}(A_{N}(x^{\prime}))+\lambda_{N+1}(E)\leq\lambda_{k}(A_{N}(x^{\prime})+E)\leq\lambda_{k}(A_{N}(x^{\prime}))+\lambda_{1}(E),\quad
k=1,\cdots,N+1.$

That is

$|\lambda_{k}(A_{N}(x^{\prime}))-\lambda_{k}(A_{N}(x))|\leq||E||.$

Moreover, from Gershgorin circle theorem, we have, for $1\leq j\leq N+1$

$|\lambda_{j}(E)-E_{j,j}|\leq\sum_{\begin{subarray}{c}1\leq p\leq N+1\\\ p\neq
j\end{subarray}}|E_{j,p}|.$

Hence, one gets

$\displaystyle|\lambda_{j}(E)|$ $\displaystyle\leq$
$\displaystyle|E_{j,j}|+\sum_{\begin{subarray}{c}1\leq p\leq N+1\\\ p\neq
j\end{subarray}}|E_{j,p}|=\sum_{1\leq p\leq N+1}|E_{j,p}|$ $\displaystyle\leq$
$\displaystyle\frac{1}{n}\sum_{1\leq p\leq
N+1}\big{|}\widetilde{P}^{(\alpha,\beta)}_{j-1}(X_{i})\widetilde{P}^{(\alpha,\beta)}_{p-1}(X_{i})|+|\widetilde{P}^{(\alpha,\beta)}_{j-1}(X^{\prime}_{i})\widetilde{P}^{(\alpha,\beta)}_{p-1}(X^{\prime}_{i})\big{|}$
$\displaystyle\leq$
$\displaystyle\frac{2}{n}\max_{x\in[-1,1]}|\widetilde{P}^{(\alpha,\beta)}_{N}(x)|\sum_{k=1}^{N+1}\Big{(}\max_{x\in[-1,1]}|\widetilde{P}^{(\alpha,\beta)}_{k}(x)|\Big{)}\leq\frac{2}{n}m_{\alpha,\beta}^{2}(N+1)^{2\mu+2}=2\frac{L_{N}(\alpha,\beta)}{n}.$

This last inequality follows from (21). Since $E$ is Hermitian, then
$||E||=\displaystyle\max_{\begin{subarray}{c}1\leq j\leq
N+1\end{subarray}}|\lambda_{j}|$ which concludes the proof.

######  Remark 1.

For the special convenient values of $\alpha=\beta=-\frac{1}{2},$ the estimate
of $\kappa_{2}(A_{N})$ given by the previous theorem can be further improved.
In fact, in this case, the normalized Jacobi polynomials are reduced to the
normalized Chebyshev polynomials $\widetilde{T}_{n}$ defined on $[-1,1]$ by

$\widetilde{T}_{n}(x)=\frac{2}{\pi}T_{n}(x),\quad
T_{n}(cos\theta)=\cos(n\theta),\quad\theta\in[0,\pi].$

It is easy to see that in this case, the optimal value of ${\displaystyle
m_{-1/2,-1/2}=\frac{2}{\pi}.}$ Consequently, with high probability, we have

$\kappa_{2}(A_{N})\lesssim\frac{1.72+\frac{2}{\pi}\frac{(N+1)\,\log(N+1)}{n}}{0.63-\frac{2}{\pi}\frac{(N+1)\,\log(N+1)}{n}}.$
(32)

######  Remark 2.

Note that the upper bound of $\kappa_{2}(A_{N})$ given by (24) of the previous
theorem is the essential condition for the validity of our nonparametric
regression estimator. Numerical results given in the last section of this
work, indicate that the quantity $\kappa_{2}(A_{N})$ remains bounded by a
fairly convenient positive constant when the i.i.d. random sampling points
follow a more general probability law. That is our proposed nonparametric
regression estimator can be used in the framework of a more general random
sampling set $\\{X_{i},\,1\leq i\leq n\\}.$

Next, to study the regression estimation error in the
$\|\cdot\|_{\omega}-$norm, we need the following technical lemma.

######  Lemma 2.

Let $g\in L^{2}_{\omega}(I)$ be a bounded function. Under the previous
notations, for any $\alpha\geq-\frac{1}{2}$ and any $0<\delta<1,$ we have with
probability at least $1-\delta,$

$\left|\frac{\gamma_{\alpha,\beta}}{n}\sum_{j=1}^{n}\big{(}g(X_{j})\big{)}^{2}-\|g\|_{\omega}^{2}\right|\leq\frac{\gamma_{\alpha,\beta}}{\sqrt{n}}\sqrt{\log\big{(}\frac{2}{\delta}\big{)}}\|g\|_{\infty}^{2}.$
(33)

where $\gamma_{\alpha,\beta}=B(\alpha+1,\beta+1)2^{\alpha+\beta+1}.$

Proof: It suffices to use McDiarmid’s inequality with the real valued
$n-$variate function,

$\phi:I^{n}\rightarrow\mathbb{R},\,\,(x_{1},\ldots,x_{n})\rightarrow\frac{\gamma_{\alpha,\beta}}{n}\sum_{i=1}^{n}\big{(}g(x_{i})\big{)}^{2}.$

Then, for any $x_{1},\ldots,x_{i-1},x^{\prime}_{i},x_{i+1},\ldots,x_{n}\in I,$
we have

$\displaystyle|\phi(x_{1},\ldots,x_{i-1},x_{i},x_{i+1},\ldots,x_{n})-\phi(x_{1},\ldots,x_{i-1},x^{\prime}_{i},x_{i+1},\ldots,x_{n})|$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad=\frac{\gamma_{\alpha,\beta}}{n}\Big{|}\big{(}g(x_{i})\big{)}^{2}-\big{(}g(x^{\prime}_{i})\big{)}^{2}\Big{|}\leq\frac{\gamma_{\alpha,\beta}}{n}\|g\|_{\infty}^{2}=d_{i}.$

Consequently, we have
${\displaystyle\sum_{i=1}^{n}d_{i}^{2}=\frac{(\gamma_{\alpha,\beta})^{2}}{n}\|g\|_{\infty}^{4}.}$
On the other hand, we have

$\mathbb{E}\left(\frac{\gamma_{\alpha,\beta}}{n}\sum_{i=1}^{n}\big{(}g(X_{i})\big{)}^{2}\right)=\mathbb{E}\Big{(}\gamma_{\alpha,\beta}\big{(}g(X)\big{)}^{2}\Big{)}=\int_{I}|g(x)|^{2}\,\omega_{\alpha,\beta}(x)\,dx=\|g\|_{\omega}^{2}.$

The following theorem provides us with an upper bound in the
$\|\cdot\|_{\omega}-$norm for the regression error in terms of
$\|\pi_{N}f\|_{\omega}$ and $\|\pi_{N}f\|_{\infty}.$ Here, $\pi_{N}f$ is the
orthogonal projection of the regression function over the subspace
$\mathcal{H}_{N}=\mbox{Span}\\{\widetilde{P}^{(\alpha,\beta)}_{k}(\cdot),\,0\leq
k\leq N\\}.$

######  Theorem 4.

Under the previous notations and the hypotheses of Theorem 3, let
${\displaystyle\boldsymbol{\eta}_{n}=\max_{1\leq i\leq n}|\varepsilon_{i}|}$,
then for any $0<\delta<1,$ we have with probability at least $(1-\delta)^{2},$

$\|f-\widehat{f}_{n,N}\|_{\omega}\leq\|f-\pi_{N}f\|_{\omega}+\sqrt{2\kappa_{2}(A_{N})}\frac{{\Big{(}\frac{1}{n}\log\big{(}\frac{2}{\delta}\big{)}\Big{)}}^{\frac{1}{4}}\|f-\pi_{N}f\|_{\infty}+\frac{\|f-\pi_{N}f\|_{\omega}}{\sqrt{\gamma_{\alpha,\beta}}}+\boldsymbol{\eta}_{n}}{\frac{1}{\sqrt{\gamma_{\alpha,\beta}}}-{\Big{(}\frac{1}{n}\log\big{(}\frac{2}{\delta}\big{)}\Big{)}}^{\frac{1}{4}}\frac{\|\pi_{N}f\|_{\infty}}{\|\pi_{N}f\|_{\omega}}}$
(34)

Proof: Since for any $x\in I,$
$\pi_{N}f(x)=\sum_{k=1}^{N+1}c_{k}(f)\widetilde{P}^{(\alpha,\beta)}_{k-1}(x),$
where
$c_{k}(f)=\int_{I}f(t)\widetilde{P}^{(\alpha,\beta)}_{k-1}(t)\omega_{\alpha,\beta}(t)\,dt,$
then we have

$B_{N}\cdot\boldsymbol{C}=\Big{[}\frac{1}{\sqrt{n}}\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{j})\Big{]}_{j,k}\cdot\boldsymbol{C}=\mathbf{P}.$
(35)

where $\boldsymbol{C}=\big{(}c_{0}(f),\ldots,c_{N}(f)\big{)}^{\prime},$ and
$\mathbf{P}=\frac{1}{\sqrt{n}}\big{(}\pi_{N}f(X_{1}),\ldots,\pi_{N}f(X_{n})\big{)}^{\prime}.$
The estimator $\widehat{f}_{n,N}$ is given by
${\displaystyle\widehat{f}_{n,N}(x)=\sum_{k=1}^{N+1}\widehat{c}_{k}\widetilde{P}^{(\alpha,\beta)}_{k-1}(x),\,x\in
I}$ where
$\widehat{\boldsymbol{C}}=\big{(}\widehat{c}_{0},\ldots,\widehat{c}_{N}\big{)}^{\prime}$
is the least square solution of the over-determined system

$B_{N}\cdot\widehat{\boldsymbol{C}}=\Big{[}\frac{1}{\sqrt{n}}\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{j})\Big{]}_{j,k}\cdot\widehat{\boldsymbol{C}}=\mathbf{Z}_{n}.$
(36)

This can be written as

$B_{N}\cdot\widehat{\boldsymbol{C}}=\mathbf{P}+\mathbf{\Delta P}.$ (37)

where $\mathbf{\Delta
P}=\frac{1}{\sqrt{n}}\big{(}(f-\pi_{N}f)(X_{1})+\varepsilon_{1},\ldots,(f-\pi_{N}f)(X_{n})+\varepsilon_{n}\big{)}^{\prime}.$
Consequently, the least square solution of (36) is a perturbation of the least
square solution of (35). In this case, we have (see for example [6])

$\frac{\|\boldsymbol{C}-\widehat{\boldsymbol{C}}\|_{\ell_{2}}^{2}}{\|\boldsymbol{C}\|_{\ell_{2}}^{2}}\leq\kappa_{2}(A_{N})\frac{\|\mathbf{\Delta
P}\|_{\ell_{2}}^{2}}{\|\mathbf{P}\|_{\ell_{2}}^{2}}.$

On the other hand, since

${(f-\pi_{N}f)(X_{j})+\varepsilon_{j}}^{2}\leq
2[{(f-\pi_{N}f)(X_{j})}^{2}+{\varepsilon_{j}}^{2}],\quad\forall\;1\leq j\leq
n,$

then we have

$\|\mathbf{\Delta P}\|_{\ell_{2}}^{2}\leq
2(\|\mathbf{\widetilde{\Delta}P}\|_{\ell_{2}}^{2}+\boldsymbol{\eta}_{n}^{2})$

where $\boldsymbol{\eta}_{n}=\displaystyle\max_{\begin{subarray}{c}1\leq j\leq
n\end{subarray}}|\varepsilon_{j}|$ and
$\mathbf{{\widetilde{\Delta}P}}=\frac{1}{\sqrt{n}}\big{(}(f-\pi_{N}f)(X_{1}),\ldots,(f-\pi_{N}f)(X_{n})\big{)}^{\prime}.$
Hence, one gets

$\frac{\|\boldsymbol{C}-\widehat{\boldsymbol{C}}\|_{\ell_{2}}^{2}}{\|\boldsymbol{C}\|_{\ell_{2}}^{2}}\leq
2\kappa_{2}(A_{N})\frac{\|\mathbf{\widetilde{\Delta}P}\|_{\ell_{2}}^{2}+\boldsymbol{\eta}_{n}^{2}}{\|\mathbf{P}\|_{\ell_{2}}^{2}}.$
(38)

In order to estimate $\|f-\widehat{f}_{n,N}\|_{\omega},$ one needs to estimate
$\|\boldsymbol{C}-\widehat{\boldsymbol{C}}\|_{\ell_{2}}^{2}$ which requires
the estimation of
$\|\mathbf{\widetilde{\Delta}P}\|_{\ell_{2}}^{2}+\boldsymbol{\eta}_{n}^{2}$
and ${\|\mathbf{P}\|_{\ell_{2}}^{2}}$. Taking $g=f-\pi_{N}f$ in Lemma 2 gives,
with probability at least $(1-\delta)^{2},$

$\left|\frac{\gamma_{\alpha,\beta}}{n}\sum_{j=1}^{n}\big{(}(f-\pi_{N}f)(X_{j})\big{)}^{2}-\|f-\pi_{N}f\|_{\omega}^{2}\right|\leq\frac{\gamma_{\alpha,\beta}}{\sqrt{n}}\sqrt{\log\big{(}\frac{2}{\delta}\big{)}}\|f-\pi_{N}f\|_{\infty}^{2}.$

Hence, with probability at least $(1-\delta)^{2},$ we have

$\frac{1}{n}\sum_{j=1}^{n}\big{(}(f-\pi_{N}f)(X_{j})\big{)}^{2}\leq\sqrt{\frac{\log\big{(}\frac{2}{\delta}\big{)}}{n}}\|f-\pi_{N}f\|_{\infty}^{2}+\frac{\|f-\pi_{N}f\|_{\omega}^{2}}{\gamma_{\alpha,\beta}}.$

Equivalently, we have

$\|\mathbf{\widetilde{\Delta}P}\|_{\ell_{2}}^{2}+\boldsymbol{\eta}_{n}^{2}\leq\sqrt{\frac{\log\big{(}\frac{2}{\delta}\big{)}}{n}}\|f-\pi_{N}f\|_{\infty}^{2}+\frac{\|f-\pi_{N}f\|_{\omega}^{2}}{\gamma_{\alpha,\beta}}+\boldsymbol{\eta}_{n}^{2}$

with probability at least $(1-\delta)^{2}.$ In the same way, taking
$g=\pi_{N}f$ in Lemma 2, one gets with probability at least $(1-\delta)^{2},$

$\frac{1}{n}\sum_{j=1}^{n}\big{(}(\pi_{N}f)(X_{j})\big{)}^{2}\geq\frac{\|\pi_{N}f\|_{\omega}^{2}}{\gamma_{\alpha,\beta}}-\sqrt{\frac{\log\big{(}\frac{2}{\delta}\big{)}}{n}}\|\pi_{N}f\|_{\infty}^{2}$

or equivalently

${\|\mathbf{P}\|_{\ell_{2}}^{2}\geq\frac{\|\pi_{N}f\|_{\omega}^{2}}{\gamma_{\alpha,\beta}}-\sqrt{\frac{\log\big{(}\frac{2}{\delta}\big{)}}{n}}\|\pi_{N}f\|_{\infty}^{2}}.$

Since by Parseval’s equality, we have
$\|\boldsymbol{C}\|_{\ell_{2}}^{2}=\|\pi_{N}f\|_{\omega}^{2}$ and by using
(38), one gets

$\|\boldsymbol{C}-\widehat{\boldsymbol{C}}\|_{\ell_{2}}^{2}\leq
2\kappa_{2}(A_{N})\frac{\sqrt{\frac{\log\big{(}\frac{2}{\delta}\big{)}}{n}}\|f-\pi_{N}f\|_{\infty}^{2}+\frac{\|f-\pi_{N}f\|_{\omega}^{2}}{\gamma_{\alpha,\beta}}+\boldsymbol{\eta}_{n}^{2}}{\frac{1}{\gamma_{\alpha,\beta}}-\sqrt{\frac{\log\big{(}\frac{2}{\delta}\big{)}}{n}}\frac{\|\pi_{N}f\|_{\infty}^{2}}{\|\pi_{N}f\|_{\omega}^{2}}}$

with probability at least $(1-\delta)^{2}.$ To conclude for the proof, it
suffices to note that

$\|f-\widehat{f}_{n,N}\|_{\omega}\leq\|f-\pi_{N}f\|_{\omega}+\|\pi_{N}f-\widehat{f}_{n,N}\|_{\omega}$

and use the fact that
$\|\pi_{N}f-\widehat{f}_{n,N}\|_{\omega}=\|\boldsymbol{C}-\widehat{\boldsymbol{C}}\|_{\ell_{2}}.$

The following corollary provides us with more explicit estimation error for
the estimator $\widehat{f}_{n,N}$ in the special case where the regression
function has a Lipschitzian $p-$th derivative.

######  Corollary 1.

Assume that $f$ has $p$ continuous derivatives on $I$ and its $p-$th
derivative $f^{(p)}\in Lip(\gamma),0<\gamma<1.$ Assume that $\|f\|_{\omega}$
and $\|f\|_{\infty}$ are bounded away from zero. Moreover, we assume that
$p+\gamma\geq\max(\mu+\frac{1}{2},\frac{1}{2}-\tau),\,\mu=\max(\alpha,\beta),\,\tau=\min(\alpha,\beta).$
Then under the previous notation and assumption, for any $0<\delta<1$ and for
sufficiently large values of $n,N,$ we have with high probability,

$\|f-\widehat{f}_{n,N}\|_{\omega}\lesssim\frac{\log
N}{N^{p+\gamma}}+\sqrt{2\kappa_{2}(A_{N})}\frac{\frac{\log
N}{N^{p+\gamma-\mu-\frac{1}{2}}}+\boldsymbol{\eta}_{n}}{\frac{1}{\sqrt{\gamma_{\alpha,\beta}}}-2{\Big{(}\frac{1}{n}\log\big{(}\frac{2}{\delta}\big{)}\Big{)}}^{\frac{1}{4}}\frac{\|f\|_{\infty}}{\|f\|_{\omega}}}.$
(39)

Proof: We recall that under the hypotheses of the corollary, we have see [14],

$\|f-\pi_{N}f\|_{\infty}\leq c_{1}\frac{\log
N}{N^{p+\gamma-\mu-\frac{1}{2}}},\quad\|f-\pi_{N}f\|_{\omega}\leq
c_{2}\frac{\log N}{N^{p+\gamma}}.$

Here, $c_{1},c_{2}$ are two positive constants. Hence, for sufficiently large
values of $n,N$ and by using the previous inequalities together with
inequality (34), one gets the desired error estimate (39).

Next, assume that for some $c>0,$ the regression function is the restriction
to $I$ of a $c-$bandlimited function which we denote also by $\overline{f}.$
That is $\overline{f}\in L^{2}(\mathbb{R})$ with Fourier transform supported
on the compact interval $[-c,c].$ Then a more explicit estimation error of the
estimator $\widehat{f}_{n,N}$ is given in this case by the following
corollary.

######  Corollary 2.

Assume that for some real $c>0,$ $f$ is the restriction to $I$ of a
$c-$bandlimited function $\overline{f}.$ Then, for
$\beta=\alpha\geq-\frac{1}{2},$ any $0<\delta<1$ and for sufficiently large
values of $n,N,$ we have with high probability,

$\|f-\widehat{f}_{n,N}\|_{\omega}\lesssim\eta_{N}^{(\alpha)}\|\overline{f}{f}\|_{\omega}+\sqrt{2\kappa_{2}(A_{N})}\frac{{\Big{(}\frac{1}{n}\log\big{(}\frac{2}{\delta}\big{)}\Big{)}}^{\frac{1}{4}}+\frac{\eta_{N}^{(\alpha)}\|\tilde{f}\|_{\omega}}{\sqrt{\gamma_{\alpha,\beta}}}+\boldsymbol{\eta}_{n}}{\frac{1}{\sqrt{\gamma_{\alpha,\beta}}}-2{\Big{(}\frac{1}{n}\log\big{(}\frac{2}{\delta}\big{)}\Big{)}}^{\frac{1}{4}}\frac{\|f\|_{\infty}}{\|f\|_{\omega}}},$
(40)

where
$\eta_{N}^{(\alpha)}=c^{-1/2}C_{\alpha}\Big{(}\frac{ec}{2N+2}\Big{)}^{N+2},$
for some constant $C_{\alpha}.$

Proof: In [12], it has been shown that for such a function $f,$ we have

$\|f-\pi_{N}f\|_{\omega}\leq\eta_{N}^{(\alpha)}\|\overline{f}\|_{L^{2}(\mathbb{R})},\quad\|f-\pi_{N}f\|_{\infty}\leq\gamma_{N}^{(\alpha)}\|\overline{f}\|_{L^{2}(\mathbb{R})},$

where
$\eta_{N}^{(\alpha)}=c^{-1/2}C_{\alpha}\Big{(}\frac{ec}{2N+2}\Big{)}^{N+2},\,\gamma_{N}^{(\alpha)}=c^{\alpha}C_{\alpha}\Big{(}\frac{ec}{2N+2}\Big{)}^{N+3/2-\alpha},$
for some constant $C_{\alpha}.$ Hence, by combining the previous two
inequalities and (34), one gets the desired error estimate (40).

######  Remark 3.

Although the quantity ${\displaystyle\boldsymbol{\eta}_{n}=\max_{1\leq i\leq
n}|\eta_{i}|},$ given in Theorem 4 and the previous two corollaries depends on
$n,$ in practice and with high probability, the different $\eta_{i}$ are
uniformly bounded independently of $n.$ This is the case for example when the
$\eta_{i}$ are i.i.d. copies of the largely used centered Gaussian random
variable with variance $\sigma^{2}.$ In this case, for any fixed $k_{0}>0$ and
for any $i\in\mathbb{N},$ we have $|\eta_{i}|\leq k_{0}\sigma$ with
probability at least ${\displaystyle
1-\mbox{erf}\Big{(}\frac{k_{0}}{\sqrt{2}}\Big{)}\approx
1-\frac{e^{-k_{0}^{2}/2}}{k_{0}\sqrt{\pi/2}}}.$ This last quantity is very
close to $1$ even for reasonable small values of $k_{0}.$ Consequently, the
results of Theorem 4 and its two corollaries still hold with high probability,
with the quantity $\boldsymbol{\eta}_{n}$ replaced by a uniform constant
$\boldsymbol{\eta}.$

Next, in order to have an $L_{2}$-risk error, we need to define a truncated
version of our estimator under some conditions. Let $r>0$ and fix a constant
$0<c<0.63$ and consider two positive integers $n,N$ satisfying the inequality

$0.63-m_{\alpha,\beta}^{2}\frac{(N+1)^{2\mu+2}\,\log(N+1)}{n}-n^{-r}\geq c>0$
(41)

Also, assume that the regression function is almost everywhere bounded by a
constant $M$, that is

$|f(x)|\leq M,\quad\mbox{a.e.}x\in I.$ (42)

Let $\tilde{f}_{N,M}$ be the truncated version of the estimate
$\widehat{f}_{n,N}$ given by

$\tilde{f}_{N,M}(x)=\mbox{sign}(\widehat{f}_{n,N}(x))\min(M,|\widehat{f}_{n,N}(x)|)$
(43)

Under the usual assumption that the added random noise
$(\varepsilon_{i})_{1\leq i\leq n}$ are i.i.d. real-valued centered random
variables with variance $\sigma^{2}$ and independent from the $X_{i},$ we have
the following theorem that provides us with an estimate of the $L_{2}$-risk
error of the estimator $\tilde{f}_{N,M}$. The proof of this theorem is partly
inspired from the techniques developed in [22].

######  Theorem 5.

Under the previous notation and hypotheses, we have

$\mathbb{E}\Big{[}\|f-\tilde{f}_{N,M}\|_{\omega}^{2}\Big{]}\leq\frac{1}{\gamma_{\alpha,\beta}c^{2}}\left(\frac{N}{n}\sigma^{2}+\frac{m_{\alpha,\beta}^{2}(N+1)^{2\mu+2}}{n}\|f-\pi_{N}f\|_{\omega}^{2}\right)+\|f-\pi_{N}f\|_{\omega}^{2}+4M^{2}\gamma_{\alpha,\beta}n^{-r}.$
(44)

Proof: Recall that from (31), we have for any $\varepsilon>0$

$\lambda_{\min}(A_{N})\geq\big{(}0.63-m_{\alpha,\beta}^{2}\frac{(N+1)^{2\mu+2}}{n}\log(N+1)\big{)}-\varepsilon$

with probability at least ${\displaystyle
1-\exp\Big{(}\frac{-n\varepsilon^{2}}{2m_{\alpha,\beta}^{4}(N+1)^{4\mu+4}}\Big{)}}$.
In particular for $\varepsilon=\varepsilon_{r,n,N}=\sqrt{\frac{2r\log
n}{n}}(N+1)^{\mu+1}$, and under condition (41), one gets

$\mathbb{P}(\lambda_{\min}(A_{N})\geq c)\geq 1-n^{-r}$ (45)

As it is done in [22], we consider the two measurable subsets of $I^{n}$
denoted by $\Omega_{+}$ and $\Omega_{-}$, where $\Omega_{+}$ is the set of all
possible draw $(X_{1},\cdots,X_{n})$ giving $\lambda_{\min}(A_{N})\geq c$,
while $\Omega_{-}$ is the set of possible draw giving
$\lambda_{\min}(A_{N})<c$. Let $d\boldsymbol{\rho}$ be the probability measure
on $I^{n}$, given by the tensor product of unidimensional probability measure
associated to the $B(\alpha+1,\beta+1)$ distribution, that it is

$d\boldsymbol{\rho}=\otimes^{n}d\rho(x_{i}),\qquad
d\rho(x_{i})=\frac{1}{\gamma_{\alpha,\beta}}\omega_{\alpha,\beta}(x_{i})\mathbf{1}_{I}(x_{i})\,dx_{i},\qquad\gamma_{\alpha,\beta}=B(\alpha+1,\beta+1)2^{\alpha+\beta+1}.$
(46)

Then, we have

$\int_{\Omega_{-}}d\boldsymbol{\rho}=\mathbb{P}\big{\\{}(X_{1},\cdots,X_{n})\in
I^{n};\lambda_{\min}(A_{N})<c\big{\\}}\leq n^{-r}.$ (47)

Next, by using (43), it is easy to see that the truncated estimator
$\tilde{f}_{N,M}$ satisfies

$|f(x)-\tilde{f}_{N,M}|\leq|f(x)-\widehat{f}_{n,N}(x)|,\quad\forall\,x\in I$
(48)

and

$|f(x)-\tilde{f}_{N,M}(x)|\leq|f(x)|+|\tilde{f}_{N,M}|\leq
2M,\quad\forall\,x\in I.$ (49)

Hence, we have

$\mathbb{E}\big{(}\|f-\tilde{f}_{N,M}\|_{\omega}^{2}\big{)}=\int_{\Omega_{+}}\|f-\tilde{f}_{N,M}\|_{\omega}^{2}d\boldsymbol{\rho}+\int_{\Omega_{-}}\|f-\tilde{f}_{N,M}\|_{\omega}^{2}d\boldsymbol{\rho}$
(50)

By using (49), one gets

$\int_{\Omega_{-}}\|f-\tilde{f}_{N,M}\|_{\omega}^{2}d\boldsymbol{\rho}\leq
4M^{2}\gamma_{\alpha,\beta}n^{-r}$ (51)

To bound the first quantity of the right hand side of (50), we use (48) to get

$\displaystyle\int_{\Omega_{+}}\|f-\tilde{f}_{N,M}\|_{\omega}^{2}d\boldsymbol{\rho}$
$\displaystyle\leq$
$\displaystyle\int_{\Omega_{+}}\|f-\widehat{f}_{n,N}\|_{\omega}^{2}d\boldsymbol{\rho}$
(52) $\displaystyle\leq$
$\displaystyle\int_{\Omega_{+}}\|f-\pi_{N}f\|_{\omega}^{2}d\boldsymbol{\rho}+\int_{\Omega_{+}}\|\pi_{N}f-\widehat{f}_{n,N}\|_{\omega}^{2}d\boldsymbol{\rho}$

By using Parseval’s equality and our previous notation, we have, on
$\Omega_{+}$

$\displaystyle\|\pi_{N}f-\widehat{f}_{n,N}\|_{\omega}^{2}$ $\displaystyle=$
$\displaystyle\|\boldsymbol{C}-\widehat{\boldsymbol{C}}\|_{\ell_{2}}^{2}\leq\|A_{N}^{-1}\|_{2}^{2}\|B_{N}^{\prime}\mathbf{\Delta
P}\|_{\ell_{2}}^{2}$ $\displaystyle\leq$
$\displaystyle\frac{1}{c^{2}}\|B_{N}^{\prime}\mathbf{\Delta
P}\|_{\ell_{2}}^{2}$

That is

$\mathbb{E}\Big{[}\|\pi_{N}f-\widehat{f}_{n,N}\|_{\omega}^{2}\Big{]}\leq\frac{1}{c^{2}}\mathbb{E}\Big{[}\|B_{N}^{\prime}\mathbf{\Delta
P}\|_{\ell_{2}}^{2}\Big{]}.$

But

$\|B_{N}^{\prime}\mathbf{\Delta
P}\|_{\ell_{2}}^{2}=\frac{1}{n^{2}}\sum_{k=1}^{N+1}\sum_{j,l=1}^{n}\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{j})\Big{(}g_{N}(X_{j})+\varepsilon_{j}\Big{)}\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{l})\Big{(}g_{N}(X_{l})+\varepsilon_{l}\Big{)}$

where
$g_{N}=(f-\pi_{N}f)\quad\bot\quad\widetilde{P}^{(\alpha,\beta)}_{k-1},\quad
1\leq k\leq N+1.$ Using the independence of the $X_{j}$’s as well as the fact
that the $\varepsilon_{j}$’s are independent of the $X_{j}$’s, together with
the facts that $\mathbb{E}(\varepsilon_{j})=0$ and
$\mathbb{E}\big{[}\varepsilon_{j}^{2}\big{]}=\sigma^{2}$, one gets

$\mathbb{E}\Big{[}\|B_{N}^{\prime}\mathbf{\Delta
P}\|_{\ell_{2}}^{2}\Big{]}=\frac{1}{n^{2}}\sum_{k=1}^{N+1}\sum_{j=1}^{n}\mathbb{E}\Big{[}\varepsilon_{j}^{2}\big{(}\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{j})\big{)}^{2}\Big{]}+\frac{1}{n^{2}}\sum_{k=1}^{N+1}\mathbb{E}\Big{[}\sum_{j=1}^{n}\big{(}\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{j})\big{)}^{2}\big{(}g_{N}(X_{j})\big{)}^{2}\Big{]}.$
(53)

Since
${\displaystyle\mathbb{E}\Big{[}\varepsilon_{j}^{2}\big{(}\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{j})\big{)}^{2}\Big{]}=\mathbb{E}\big{[}\varepsilon_{j}^{2}\big{]}\mathbb{E}\Big{[}\big{(}\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{j})\big{)}^{2}\Big{]}=\sigma^{2}\frac{1}{\gamma_{\alpha,\beta}}}$
and since
${\displaystyle\mathbb{E}\Big{[}(g_{N}(X_{j}))^{2}\Big{]}=\frac{1}{\gamma_{\alpha,\beta}}\|g_{N}\|_{\omega}^{2}}$,
together with the fact that
${\displaystyle\sum_{k=1}^{N+1}(\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{j}))^{2}\leq
m_{\alpha,\beta}^{2}(N+1)^{2\mu+2}}$, then we have

$\frac{1}{n^{2}}\sum_{k=1}^{N+1}\mathbb{E}\Big{[}\sum_{j=1}^{n}\big{(}\widetilde{P}^{(\alpha,\beta)}_{k-1}(X_{j})\big{)}^{2}\big{(}g_{N}(X_{j})\big{)}^{2}\Big{]}\leq\frac{m_{\alpha,\beta}^{2}(N+1)^{2\mu+2}}{n\gamma_{\alpha,\beta}}\|g_{N}\|_{\omega}^{2}.$
(54)

Finally, by combining the previous equalities and inequalities, together with
(50)–(53), one gets the desired $L_{2}-$risk error of the estimator
$\tilde{f}_{N,M}.$

Next, we briefly describe how our estimator $\widehat{f}_{n,N}$ can be used in
the framework of a general random sampling set $\\{X_{i},\,1\leq i\leq n\\},$
drawn from a distribution law, other than the Beta law with parameters
$\alpha+1,\beta+1.$ There is two cases to consider. In the first case, we
assume that the $X_{i}$ are i.i.d. copies of a random variable with known
cumulative distribution function (CDF) $F_{X}(\cdot).$ In this case, it is
well known that the transformed random variables $x_{i}=F_{X}(X_{i})$ are
i.i.d. copies of a random variables following the uniform law over $(0,1).$
Moreover, it is well known that the CDF associated with the Beta distribution
with parameters $\alpha+1,\beta+1$ is given by $I_{x}(\alpha+1,\beta+1):$ the
regularized incomplete Beta function. Since this last function is invertible,
then the transformed random sampling points, given by

$\tau_{i}=I_{x}^{-1}(\alpha+1,\beta+1)\big{(}F_{X}(X_{i})\big{)},\quad 1\leq
i\leq n$ (55)

follow the $Beta(\alpha+1,\beta+1)$ distribution. For the interesting special
values of $\alpha=\beta=-\frac{1}{2},$ see Remark 1, the inverse of the
$I_{x}(1/2,1/2)$ is simply given by

$I_{x}^{-1}(1/2,1/2)(t)=\frac{1}{2}\Big{(}1+\sin\big{(}\pi
t-\frac{\pi}{2}\big{)}\Big{)},\quad t\in[0,1].$ (56)

In the second case where the CDF $F_{X}$ is not known, this is the case where
the sampling distribution law is unknown, then one can replace in formula
(55), the CDF $F_{X}$ by an estimator of this later. There is a fairly rich
literature devoted to build estimators of the CDF of unknown probability laws,
see for example [8] and the references therein. Note that to get accurate
approximate values of the outputs $\widetilde{Y}_{i}$ at the transformed
random sampling points $\tau_{i},$ one may use an interpolation scheme based
on the observed $Y_{i}$ at the original random sampling points $X_{i}.$ In
example 1 of the numerical simulations section, we provide some numerical
results that illustrate the stability of the estimator $\widehat{f}_{n,N}$
when the random sampling set is drawn from the standard normal distribution.

Finally, we briefly describe our random pseudoinverse scheme based estimator
can be generalized to the multivariate case, with random sampling sets
$\\{X_{i},1\leq i\leq n\\}\subset\mathbb{R}^{d}.$ In fact, we only need to
replace each Jacobi polynomial $\widetilde{P}^{(\alpha,\beta)},$ $0\leq k\leq
N$ by its tensor product $d-$dimensional version

$\Psi_{\boldsymbol{m}}^{\alpha,\beta}(\boldsymbol{x})=\prod_{j=1}^{d}\widetilde{P}^{(\alpha,\beta)}_{m_{j}}(x_{j}),\quad\boldsymbol{x}=(x_{1},\ldots,x_{d})\in
I^{d},\quad\boldsymbol{m}=(m_{1},\ldots,m_{d})\in\\{0,1,\ldots,N\\}^{d}.$

Note that
$\\{\Psi_{\boldsymbol{m}}^{\alpha,\beta},\,\boldsymbol{m}\in\mathbb{N}_{0}^{d}\\}$
is an orthonormal basis of $L^{2}(I^{d},\boldsymbol{\omega}_{\alpha,\beta}),$
where
${\displaystyle\boldsymbol{\omega}_{\alpha,\beta}(\boldsymbol{x})=\prod_{j=1}^{d}\omega_{\alpha,\beta}(x_{j})}.$
For the sake of simplicity, we restricted ourselves to the case $d=1.$ The
extension of this scheme to higher dimensions is the subject of a future work.

##  4 Random pseudo-inverse based scheme for functional regression estimator

In this paragraph, we first recall that for a compact interval $J,$ the LFR
problem is given by

$Y_{i}=\int_{J}X_{i}(s)\,\beta_{0}(s)\,ds+\varepsilon_{i},\quad i=1,\ldots,n,$
(57)

where the $X_{i}(\cdot)\in L^{2}(J)$ are random functional predictors, given
by ${\displaystyle X_{i}(s)=\sum_{k=1}^{\infty}\xi_{k}Z_{i,k}\varphi_{k}(s).}$
Here, the ${\displaystyle\varphi_{k}(\cdot)}$ is an orthonormal family of
$L^{2}(J),$ the $Z_{i,k}$ are i.i.d centered random variables with variance
$\sigma_{Z}^{2}$ and $(\xi_{k})_{1\leq k\leq N}$ is a real valued
deterministic sequence. The $\varepsilon_{i}$ are i.i.d centered white noise
independent of the $X_{i}(\cdot)$ and $\beta_{0}(\cdot)\in L^{2}(J)$ is the
unknown slope function to be recovered. We use a partition of the set
$[[1,N]]$ given by $K_{N}=\Big{[}\frac{\log N}{\log 2}\Big{]}+1$ subsets
$I_{k}=[[2^{k-1}+1,\min(2^{k},N)]].$ On each subset $I_{k},$ we consider the
random matrix $G_{k}$ given by

$G_{k}=\displaystyle\left[\frac{1}{n}\sum_{i=1}^{n}\xi_{j}\xi_{l}Z_{i,j}Z_{i,l}\right]_{j,l\in
I_{k}}$ (58)

Then, we have

$\widehat{\beta}_{n,N}(s)=\sum_{k=1}^{K_{N}}\sum_{j\in
I_{k}}\widehat{c}_{n,j}\varphi_{j}(s),\quad\widehat{\boldsymbol{c}}_{n,k}=[\widehat{c}_{n,j}]_{j\in
I_{k}}^{\prime},$ (59)

where

$\widehat{\mathbf{c}}_{n,k}=G_{k}^{-1}(F_{k}^{\prime}\frac{1}{\sqrt{n}}[Y_{i}^{k}]^{\prime}_{1\leq
i\leq n}),\quad 1\leq k\leq K_{N}.$ (60)

Here, the $\boldsymbol{Y}^{k}=[Y_{i}^{k}]_{1\leq i\leq n}$ are obtained from
$\boldsymbol{Y}=[Y_{i}]_{1\leq i\leq n}$ by substituting in (5), $X_{i}(s)$
with its projection ${X_{i}^{k}(s)=\sum_{j\in
I_{k}}\xi_{j}Z_{i,j}\varphi_{j}(s)}.$ Consequently, the model (57) is
substituted with the following model

$Y^{k}_{i}=\int_{J}X^{k}_{i}(s)\,\beta_{0}(s)\,ds+\varepsilon^{k}_{i},\quad
i=1,\ldots,n,\quad 1\leq k\leq K_{N},$ (61)

where the i.i.d. noises $\varepsilon^{k}_{i}$ are centered and
$\mathbb{E}\Big{[}\big{(}\varepsilon^{k}_{i}\big{)}^{2}\Big{]}=\sigma_{k}^{2}.$
Moreover, we assume that the $\varepsilon^{k}_{i}$ are independent of the
$Z_{ij}.$ Then we show that the estimator $\widehat{\beta}_{n,N}(\cdot),$
given by (59) and (60) is an accurate and stable estimator for the slope
function $\beta_{0}(\cdot),$ solution of the previous LFR problem. The
following theorem provides us with a bound for the $2-$condition number of
$G_{k},$ as well as a bound for the estimation error
$\|\beta_{0}(\cdot)-\widehat{\beta}_{n,N}(\cdot)\|^{2}_{2}.$ For this purpose,
we use the function
${\displaystyle\boldsymbol{\xi}(\cdot)=\sum_{j=1}^{N}\xi_{j}\varphi_{j}(\cdot)}.$
Moreover, for $1\leq k\leq K_{N},$ we use the notation
$g^{k}(\cdot)=\pi_{I_{k}}g(\cdot),$ the orthogonal projection of $g(\cdot)$
over $\mbox{Span}\\{\varphi_{j}(\cdot),\,j\in I_{k}\\}.$

######  Theorem 6.

Assume that ${\displaystyle\max_{1\leq j\leq N}{|Z_{ij}|}\leq M}$ with
probability at least $(1-\delta_{N}).$ Then, for any $\eta>0,$ we have with
probability at least
$1-2\Big{(}\exp\big{(}\frac{-n\eta^{2}}{2M_{\boldsymbol{\xi}}}\big{)}+\delta_{N}\Big{)}$

$\kappa_{2}(G_{k})\leq\frac{1.72\max_{\begin{subarray}{c}j\in
I_{k}\end{subarray}}\sigma_{Z}^{2}\xi_{j}^{2}+\frac{M_{\boldsymbol{\xi}}}{n}\log(2^{k-1})+\eta}{0.63\min_{\begin{subarray}{c}j\in
I_{k}\end{subarray}}\sigma_{Z}^{2}\xi_{j}^{2}-\frac{M_{\boldsymbol{\xi}}}{n}\log(2^{k-1})-\eta};\quad
1\leq k\leq K_{N}.$ (62)

Here $M_{\boldsymbol{\xi}}=M^{2}\max_{\begin{subarray}{c}1\leq j\leq
N\end{subarray}}{|\xi_{j}|}\|\boldsymbol{\xi}\|_{\ell_{1}}.$ Moreover, we have

$\|\widehat{\beta}_{n,N}(\cdot)-\beta_{0}(\cdot)\|_{2}^{2}\leq\frac{1}{n}\Big{(}\sum_{k=1}^{K_{N}}\kappa_{2}(G_{k})\Big{)}\max_{1\leq
k\leq
K_{N}}\frac{\|\beta^{k}_{0}(\cdot)\|_{2}^{2}\|(\varepsilon^{k}_{i})_{i}\|^{2}_{\ell_{2}}}{\sigma_{Z}^{2}\sum_{j\in
I_{k}}\xi_{j}^{2}c_{j}^{2}-\eta}.$ (63)

with probability at least
$\Big{(}1-\exp\big{(}-2n\frac{\eta^{2}}{M^{4}\|\boldsymbol{\xi}(\cdot)\|_{\ell_{2}}^{4}\|\beta_{0}(\cdot)\|_{L_{2}}^{4}}\big{)}\Big{)}^{K_{N}}.$

Proof: Since the i.i.d. random variables $Z_{i,j}$ are centred with variances
$\sigma_{Z}^{2}$, then

$\mathbb{E}(G_{k})=\begin{bmatrix}\sigma_{Z}^{2}\xi_{2^{k-1}+1}^{2}&&\\\
&\ddots\\\ &&\sigma_{Z}^{2}\xi_{2^{k}}\end{bmatrix}.$

Consequently, the minimum and the maximum eigenvalue of $\mathbb{E}(G_{k})$
are given by

$\Lambda_{\min}=\lambda_{\min}(\mathbb{E}(G_{k}))=\sigma_{Z}^{2}\displaystyle\min_{\begin{subarray}{c}j\in
I_{k}\end{subarray}}{\xi_{j}^{2}},\quad\Lambda_{\max}=\lambda_{\max}(\mathbb{E}(G_{k}))=\sigma_{Z}^{2}\displaystyle\max_{\begin{subarray}{c}j\in
I_{k}\end{subarray}}{\xi_{j}^{2}}.$

On the other hand, the random matrix $G_{k}$ is written in the following form

$G_{k}=\frac{1}{n}\sum_{i=1}^{n}D_{i,k},\quad
D_{i,k}=\left[\sum_{i=1}^{n}\xi_{j}\xi_{l}Z_{i,j}Z_{i,l}\right]_{j,l\in
I_{k}}.$

Note that each matrix $D_{i,k}$ is positive semi definite. This follows from
the fact that for any $\boldsymbol{x}\in\mathbb{R}^{2^{k-1}}$, we have

$\boldsymbol{x^{\prime}}D_{i,k}\boldsymbol{x}=\boldsymbol{x^{\prime}}B_{i,k}B_{i,k}\boldsymbol{x}\geq
0,\quad B_{i,k}=\frac{1}{\sqrt{n}}[\xi_{j}Z_{i,j}]_{j\in I_{k}}.$

By using Gershgorin circle theorem, one gets

$\begin{array}[]{lll}\lambda_{\max}(D_{i,k})&\leq&\frac{1}{n}\max_{\begin{subarray}{c}j\in
I_{k}\end{subarray}}{|\xi_{j}Z_{i,j}|}\displaystyle\sum_{l=2^{k-1}+1}^{2^{k}}|\xi_{l}Z_{i,l}|\\\
&\leq&\frac{M^{2}}{n}\max_{\begin{subarray}{c}j\in
I_{k}\end{subarray}}{|\xi_{j}|}\displaystyle\sum_{l=2^{k-1}+1}^{2^{k}}|\xi_{l}|\leq\frac{M_{\boldsymbol{\xi}}}{n}\end{array}$

Hence, we have

$\mathbb{E}(\lambda_{\min}(G_{k}))\geq 0.63\min_{\begin{subarray}{c}j\in
I_{k}\end{subarray}}\sigma_{Z}^{2}\xi_{j}^{2}-\frac{M_{\boldsymbol{\xi}}}{n}\log(2^{k-1})$

and

$\mathbb{E}(\lambda_{\max}(G_{k}))\leq 1.72\max_{\begin{subarray}{c}j\in
I_{k}\end{subarray}}\sigma_{Z}^{2}\xi_{j}^{2}+\frac{M_{\boldsymbol{\xi}}}{n}\log(2^{k-1})$

Then, we use McDiarmid’s concentration inequality to conclude that for any
$\eta>0$, we have

$\lambda_{\min}(G_{k})\geq\mathbb{E}(\lambda_{\min}(G_{k}))-\eta,\quad\lambda_{\max}(G_{k})\leq\mathbb{E}(\lambda_{\max}(G_{k}))+\eta$

with probability at least $1-\exp\frac{-n\eta^{2}}{2M_{\boldsymbol{\xi}}}$. By
combining the previous two inequalities, one gets (62).

Next to get the estimation error given by (63), we note that on each subset
$I_{k}$, the associated exact expansion coefficients vector of
$\beta_{0}(\cdot)$ is given by $\mathbf{c}^{k}=[c_{j}]_{j\in I_{k}}.$
Similarly, we define the expansion coefficients restricted to $I_{k}$ of the
estimator $\widehat{\beta}_{n,N},$ that we denote by
$\widehat{\mathbf{c}}_{n,k}.$ These two coefficients vectors are given by

${\mathbf{c}^{k}}^{\prime}=G_{k}^{-1}\Big{(}F_{k}^{\prime}\frac{1}{\sqrt{n}}[Y_{i}^{k}-\varepsilon_{i}^{k}]^{\prime}_{1\leq
i\leq
n}\Big{)}=G_{k}^{-1}\Big{(}F_{k}^{\prime}\frac{1}{\sqrt{n}}[\tilde{Y_{i}^{k}}]^{\prime}_{1\leq
i\leq
n}\Big{)},\quad\widehat{\mathbf{c}}_{n,k}=G_{k}^{-1}\Big{(}F_{k}^{\prime}\frac{1}{\sqrt{n}}[Y_{i}^{k}]^{\prime}_{1\leq
i\leq n}\Big{)},$ (64)

where,

$F_{k}=\frac{1}{\sqrt{n}}\big{[}\xi_{j}Z_{i,j}\big{]}_{1\leq i\leq n,j\in
I_{k}},\quad 1\leq i\leq n,\quad 1\leq k\leq K_{N}.$ (65)

Then we have the following perturbation result for the pseudo-inverse least
square solution of an over-determined system,

$\frac{\|\widehat{\mathbf{c}}_{n,k}-\mathbf{c}^{k}\|_{\ell_{2}}^{2}}{\|\mathbf{c}^{k}\|_{\ell_{2}}^{2}}\leq\kappa_{2}(G_{k})\frac{\|\mathbf{\widetilde{Y}}^{k}-\mathbf{Y}^{k}\|_{\ell_{2}}^{2}}{\|\mathbf{\widetilde{Y}}^{k}\|_{\ell_{2}}^{2}}=\kappa_{2}(G_{k})\frac{\frac{1}{n}\|(\varepsilon_{i}^{k})_{i}\|_{\ell_{2}}^{2}}{\frac{1}{n}\|(\widetilde{Y_{i}}^{k})_{i}\|_{\ell_{2}}^{2}}.$
(66)

On the other hand, since
$\Big{(}\widetilde{Y_{i}}^{k}\Big{)}^{2}=\Big{(}\sum_{j\in
I_{k}}\xi_{j}Z_{i,j}c_{j}\Big{)}^{2},$ then we have

$0\leq\big{(}\widetilde{Y_{i}}^{k}\big{)}^{2}\leq
M^{2}\|\boldsymbol{\xi}^{k}(\cdot)\|^{2}_{\ell_{2}}\|\beta_{0}^{k}(\cdot)|_{L_{2}}^{2}$

and

$\mathbb{E}\big{[}\big{(}\widetilde{Y_{i}}^{k}\big{)}^{2}\big{]}=\sigma_{Z}^{2}\sum_{j\in
I_{k}}\xi_{j}^{2}c_{j}^{2}$

Hence by using Hoeffding’s inequality, one concludes that for any $\eta>0$, we
have

$\frac{1}{n}\sum_{i=1}^{n}\big{(}\widetilde{Y_{i}}^{k}\big{)}^{2}-\mathbb{E}\big{[}\big{(}\widetilde{Y_{i}}^{k}\big{)}^{2}\big{]}=\frac{1}{n}\sum_{i=1}^{n}\big{(}\widetilde{Y_{i}}^{k}\big{)}^{2}-\sigma_{Z}^{2}\sum_{j\in
I_{k}}\xi_{j}^{2}c_{j}^{2}\geq-\eta$ (67)

with probability at least
$1-\exp\big{(}-2n\frac{\eta^{2}}{{M}^{4}\|\xi^{l}(\cdot)\|^{2}_{\ell_{2}}\|\beta_{0}^{l}\|_{L_{2}}^{4}}\big{)}$.
Finally, by using the inequalities (64)–(67), together with Parseval’s
equality and some straightforward computations, one gets the desired
inequality (63).

######  Remark 4.

In practice, the deterministic sequence $\boldsymbol{\xi}=(\xi_{j})_{j}$
satisfies a decay condition of the type $\xi_{j}\asymp j^{-s}$ for some $s>0.$
In this case, the combined dyadic decomposition and random pseudo-inverse
scheme estimator $\widehat{\beta}_{n,N}(\cdot)$ is stable in the sense that

$\kappa_{2}(G)=\sum_{k=1}^{K_{N}}\kappa_{2}(G_{k})\lesssim
2^{s}\frac{1.72\cdot\log N}{0.63\cdot\log 2}$ (68)

is a relatively small cumulative condition number.

Next, we are interested in having an $L_{2}$-risk error of a truncated version
of our LFR estimator $\widehat{\beta}_{n,N}(\cdot).$ More precisely, we assume
that the true slope function $\beta_{0}(\cdot)$ is almost every where bounded,
that is there exists $L>0$ such that

$|\beta_{0}(x)|\leq L,\quad a.e.x\in J.$

Then, as for the nonparametric estimator of the previous section, the
truncated $\widehat{\beta}_{n,N}(\cdot)$ estimator, denoted by
$\widetilde{\beta}_{N,L}(\cdot)$ is defined by

$\widetilde{\beta}_{N,L}(x)=\mbox{Sign}\big{(}\widehat{\beta}_{n,N}(x)\big{)}\min\Big{(}L,|\widehat{\beta}_{n,N}(x)|\Big{)}.$

Let $r>0$ and $\eta_{k}>0$ be such that

$\mathbb{P}\big{(}\lambda_{\min}(G_{k})\geq\eta_{k}\big{)}\geq 1-n^{-r},\quad
1\leq k\leq K_{N}.$ (69)

The following theorem provides us with the $L_{2}$-risk of the estimator
$\widetilde{\beta}_{N,M}.$

######  Theorem 7.

Under the previous notations and hypotheses, we have

$\mathbb{E}\Big{[}\|\widetilde{\beta}_{N,M}(\cdot)-\beta_{0}(\cdot)\|_{2}^{2}\Big{]}\leq\frac{\sigma_{Z}^{2}\|(\xi_{i})_{i}\|_{\ell_{2}}^{2}}{n^{2}}\Bigg{(}\sum_{k=1}^{K_{N}}\frac{\sigma_{k}^{2}}{\eta_{k}^{2}}|I_{k}|\Bigg{)}+\frac{4L^{2}}{n^{r}}K_{N}.$
(70)

Proof: Let $\mathbf{\tilde{c}_{n,k}}$ be the expansion coefficients vector of
$\widetilde{\beta}_{N,M}^{k}$, the projection over
$\mbox{Span}\\{\varphi_{j}(\cdot),\,j\in I_{k}\\}$ of
$\widetilde{\beta}_{N,M}.$ For each $1\leq k\leq K_{N},$ let $\Omega_{+,k}$
and $\Omega_{-,k}$ be the set of all possible draw
$(X_{1}(\cdot),\cdots,X_{n}(\cdot))$ for which
$\lambda_{\min}(G_{k})\geq\eta_{k}$ and $0<\lambda_{\min}(G_{k})<\eta_{k}$
respectively. Then , we have

$\mathbb{E}\Big{[}\|\mathbf{\tilde{c}_{n,k}}-\mathbf{c}^{k}\|_{\ell_{2}}^{2}\Big{]}\leq\int_{\Omega_{+,k}}\int_{J}|(\widetilde{\beta}_{N,M}^{k}-\beta_{0}^{k})(x)|^{2}dxd\boldsymbol{\rho}+\int_{\Omega_{-,k}}\int_{J}|(\widetilde{\beta}_{N,M}^{k}-\beta_{0}^{k})(x)|^{2}dxd\boldsymbol{\rho}.$
(71)

where

$\int_{\Omega_{-,k}}\int_{J}|(\widetilde{\beta}_{N,M}^{k}-\beta_{0}^{k})(x)|^{2}dxd\boldsymbol{\rho}\leq\frac{4L^{2}}{n^{r}}.$
(72)

and

$\int_{\Omega_{+,k}}\int_{J}|(\widetilde{\beta}_{N,M}^{k}-\beta_{0}^{k})(x)|^{2}dxd\boldsymbol{\rho}\leq\int_{\Omega_{+,k}}\int_{J}|(\widehat{\beta}_{n,N}^{k}-\beta_{0}^{k})(x)|^{2}dxd\boldsymbol{\rho}.$
(73)

Here, $d\boldsymbol{\rho}$ is the tensor product probability measure, defined
in a similar manner as (46) and associated with the probability law of the
random predictors $X_{i}(\cdot).$ By using Parseval’s equality, we have from
(64)

$\int_{J}|(\widehat{\beta}_{n,N}^{k}-\beta_{0}^{k})(x)|^{2}dx=\|\widehat{\mathbf{c}}_{n,k}-\mathbf{c}^{k}\|_{\ell_{2}}^{2}\leq\frac{1}{n}\|G_{k}^{-1}\|_{2}^{2}\|F_{k}^{\prime}\big{(}\mathbf{\widetilde{Y}}^{k}-\mathbf{Y}^{k}\big{)}\|_{\ell_{2}}^{2}.$

So that on $\Omega_{+,k},$ we obtain

$\int_{J}|(\widehat{\beta}_{n,N}^{k}-\beta_{0}^{k})(x)|^{2}dx\leq\frac{1}{n\eta_{k}^{2}}\|F_{k}^{\prime}\big{(}\mathbf{\widetilde{Y}}^{k}-\mathbf{Y}^{k}\big{)}\|_{\ell_{2}}^{2}.$

But

$F_{k}^{\prime}\big{(}\mathbf{\widetilde{Y}}^{k}-\mathbf{Y}^{k}\big{)}=\frac{1}{\sqrt{n}}\Bigg{[}\displaystyle\sum_{l=1}^{n}-\varepsilon_{l}^{k}\xi_{i}Z_{l,i}\Bigg{]}^{\prime}_{i\in
I_{k}}=\frac{1}{\sqrt{n}}\Big{[}a_{i}\Big{]}^{\prime}_{i\in I_{k}}.$

By using the hypotheses on the noises $\varepsilon_{i},$ it is easy to see
that

$\mathbb{E}\big{[}a_{i}^{2}\big{]}=\mathbb{E}\Bigg{[}\sum_{l,j=1}^{n}\varepsilon^{k}_{l}\varepsilon^{k}_{j}\xi_{l}\xi_{i}Z_{l,i}Z_{j,i}\Bigg{]}=\sigma_{k}^{2}\sigma_{Z}^{2}\sum_{j=1}^{n}\xi_{j}^{2}.$

Consequently, we have

$\int_{J}|(\widehat{\beta}_{n,N}^{k}-\beta_{0}^{k})(x)|^{2}dx\leq\frac{|I_{k}|}{n^{2}\eta_{k}^{2}}\sigma_{k}^{2}\sigma_{Z}^{2}\|(\xi_{i})_{i}\|_{\ell_{2}}^{2}.$

The previous inequality, together with (71), (72) and (73) lead to

$\mathbb{E}\Big{[}\|\mathbf{\tilde{c}_{n,k}}-\mathbf{c}^{k}\|_{\ell_{2}}^{2}\Big{]}\leq\frac{|I_{k}|}{n^{2}\eta_{k}^{2}}\sigma_{k}^{2}\sigma_{Z}^{2}\|(\xi_{i})_{i}\|_{\ell_{2}}^{2}+\frac{4L^{2}}{n^{r}}$

Finally, by using Parseval’s equality and by adding inequalities on subsets
$I_{k}$, one gets the desired $L_{2}-$risk error of the estimator
$\widetilde{\beta}_{N,M}$.

######  Remark 5.

For the special of $K_{N}=1,$ that is our estimator
$\widehat{\beta}_{n,N}(\cdot)$ is constructed without the dyadic partition of
the set $[[1,N]],$ then the $L_{2}-$risk of the previous theorem is simply
given by

$\mathbb{E}\Big{[}\|\widetilde{\beta}_{N,M}(\cdot)-\beta_{0}(\cdot)\|_{2}^{2}\Big{]}\leq\frac{\sigma^{2}\sigma_{Z}^{2}\|\xi\|_{\ell_{2}}^{2}N}{n^{2}\eta^{2}_{N}}+\frac{4L^{2}}{n^{r}}.$

##  5 Numerical Simulations

In this section, we first give some numerical simulations that illustrate the
different results of this paper. Then, we provide real data application of our
nonparametric regression estimator with data corresponding to the Covid-19
daily spread over $15$ months ( starting from March 1st 2020) of four
countries covering different regions of the world.

###  5.1  Numerical examples

Example 1:  In this first example, we illustrate the result of Theorem 3, as
well as an adaptation of this result when the random sampling set is drawn
from a standard normal distribution, as described by (55) and (56). Then, we
illustrate by some numerical test, the inequality (68) which is a consequence
of Theorem 5. For this purpose, we have considered the case of the re-
normalized Jacobi polynomials, so that they are orthonormal over $I_{1}=[0,1]$
with the two special cases of $\alpha=\beta=-\frac{1}{2}$ as well as
$\alpha=\beta=0.$ In this case, the corresponding coefficients
$m_{\alpha,\beta},$ given by Proposition 1 and remark 1, are
$m_{-\frac{1}{2},-\frac{1}{2}}=\frac{2}{\pi},\quad
m_{0,0}=\sqrt{\frac{2+1/\sqrt{2}}{3}}\approx 0.95.$ Then, we have considered
the values of $N=5,10,15,20$ and different values of $n$ the size of the
sampling set $\\{X_{j},j=1,\ldots,n\\}$ following the beta
$B(\alpha,\alpha)-$law. Also, we have computed the average of the condition
number $\kappa_{2}(A_{N})$ over $50$ such realizations of $A_{N}.$ Then, we
have repeated the previous simulations, with the initial random sampling set
$\\{X_{i},\,1\leq i\leq n\\}$ drawn from the standard normal distribution. We
have considered the transformed random sampling ${\tau_{i},\,1\leq i\leq n},$
associated with the $Beta(\alpha+1,\beta+1)$ probability law. We denote by
$\widetilde{A}_{N},$ the counterparts of the random matrix $A_{N}$ when the
transformed sampling points $\tau_{i}$ are used instead of the $X_{i}.$ The
obtained numerical results are given by Table 1.

Next, to illustrate the stability of our proposed LFR estimator
$\widehat{\beta}_{n,N}(\cdot),$ we have computed the cumulative $2-$norm
condition numbers of the positive definite random matrices $G_{k},$ given by
(58), with different values of the integers $N,n$ and the decay rate exponent
$s$ of the deterministic sequence $(\xi_{j})_{j}.$ Also, we assume that the
i.i.d. random variables $Z_{j}$ follow the uniform law
$U(-\sqrt{3},\sqrt{3}).$ In Table 2, we have listed the corresponding values
of the cumulative condition numbers
${\displaystyle\kappa_{2}(G)=\sum_{k=1}^{K_{N}}\kappa_{2}(G_{k}).}$ It is
interesting to note that from these numerical values, the actual values of the
$\kappa_{2}(G)$ are slightly smaller than the theoretical bounds given by the
right-hand side of inequality (68).

$\alpha=\beta$ | $N$ | $n$ | $\kappa_{2}(A_{N})$ | $\kappa_{2}(\widetilde{A}_{N})$ | $\alpha=\beta$ | $N$ | $n$ | $\kappa_{2}(A_{N})$ | $\kappa_{2}(\widetilde{A}_{N})$
---|---|---|---|---|---|---|---|---|---
$-0.5$ | $5$ | $25$ | $7.74$ | $6.48$ | $0.0$ | $5$ | $40$ | $8.36$ | $7.27$
| $10$ | $40$ | $17.57$ | $11.03$ |  | $10$ | $100$ | $19.94$ | $13.36$
| $15$ | $60$ | $24.72$ | $22.03$ |  | $15$ | $125$ | $33.84$ | $22.52$
| $20$ | $100$ | $27.09$ | $12.03$ |  | $20$ | $250$ | $133.07$ | $29.82$
Table 1: Illustration of Theorem 3. $s$ | $N$ | $n$ | $\kappa_{2}(G)$ | $N$ | $n$ | $\kappa_{2}(G)$ | $N$ | $n$ | $\kappa_{2}(G)$ | $N$ | $n$ | $\kappa_{2}(G)$
---|---|---|---|---|---|---|---|---|---|---|---|---
$0.75$ | $20$ | $100$ | $12.05$ | $30$ | $100$ | $15.49$ | $40$ | $100$ | $17.63$ | $50$ | $100$ | $19.66$
|  | $150$ | $11.39$ |  | $150$ | $12.85$ |  | $150$ | $15.66$ |  | $150$ | $18.08$
|  | $200$ | $11.22$ |  | $200$ | $12.26$ |  | $200$ | $15.16$ |  | $200$ | $15.90$
$1.5$ | $20$ | $100$ | $25.59$ | $30$ | $100$ | $29.55$ | $40$ | $100$ | $36.23$ | $50$ | $100$ | $41.01$
|  | $150$ | $23.78$ |  | $150$ | $28.60$ |  | $150$ | $33.32$ |  | $150$ | $37.11$
|  | $200$ | $21.90$ |  | $200$ | $27.97$ |  | $200$ | $31.39$ |  | $200$ | $31.39$
Table 2: Illustration of the cumulative condition number given by (68).

Example 2: A synthetic test function. In this example, for a given real number
$s>0,$ we consider the Weierstrass test function, given by

$W_{s}(x)=\sum_{k\geq 0}\frac{\cos(2^{k}\pi x)}{2^{ks}},\quad-1\leq x\leq 1.$
(74)

Note that $W_{s}\in C^{s-\epsilon}(I),\,\forall\epsilon<s.$ Here, $C^{r}(I)$
is the usual Hölder space defined over the interval $I$ and having a Hölder
smoothness exponent $r>0.$ We have considered the two values of $s=1,2.$ Note
that for $s=1,$ $W_{s}(\cdot)$ is continuous but nowhere differentiable, which
makes this case hard to handle by classical nonparametric regression
estimators. Then, we considered the values of $\alpha=\beta=-\frac{1}{2}$ and
constructed our estimator $\widehat{f}_{n,N}(\cdot),$ based on the first $N+1$
Jacobi polynomials $\widetilde{P}_{k}^{(\alpha,\alpha)}$ with $N=10,20,30$ and
a random sampling set of size $n=100$ and following the beta
$B(\alpha+1,\alpha+1)$ distribution on $I.$ The added centered white noise in
(1) has variance $\sigma^{2},$ with $\sigma=0.1,\,0.05.$ Also, We have
performed $10$ realizations of the numerical approximations of the Weierstrass
function by our estimator.

Moreover, for comparison purpose, we have repeated the above numerical
simulations by applying the Kernel Ridge Regression (KRR) estimator,
associated with the Sinc-kernel, given by ${\displaystyle
K_{c}(x,y)=\frac{\sin(c(x-y))}{\pi(x-y)},\,x,y\in I,}$ where $c>0$ is a tuning
parameter called bandwidth. Note that the popular KRR estimator is briefly
described as follows. Given a reproducing kernel Hilbert space (RKHS),
associated with a positive-definite kernel $K(\cdot,\cdot)$ and given a
convenient regularization parameter $\lambda>0$, the KRR estimator of $f$ is
given by the solution of the following minimization problem,

$\widehat{f}_{\lambda}=\arg\min_{f\in\mathcal{H}}\left\\{\frac{1}{n}\sum_{i=1}^{n}\Big{(}f(X_{i})-Y_{i}\Big{)}^{2}+\lambda\|f\|_{\mathcal{H}}^{2}\right\\},$
(75)

where $\|\cdot\|_{\mathcal{H}}$ is the usual norm of $\mathcal{H}.$ Thanks to
the representer theorem, the solution of the previous minimization problem is
given by

$\widehat{f}_{\lambda}(x)=\sum_{k=1}^{n}c_{k}K(x,X_{k}),\quad\boldsymbol{c}=[c_{1},\ldots,c_{n}]^{\prime}=G_{\lambda}^{-1}\frac{1}{n}[Y_{1},\ldots,Y_{n}]^{\prime},$
(76)

where $G_{\lambda}=\Big{[}\frac{1}{n}K(X_{i},X_{j})\Big{]}_{1\leq i,j\leq
n}+\lambda I_{n}$ is the regularized Gram-matrix, with $I_{n}$ is the $n\times
n$ identity matrix. The convenient value of the regularization parameter is
chosen by the well-known cross validation technique. We denote by
$\widehat{f}_{c},$ the KRR estimator based on the Sinc-kernel with parameter
$c>0.$ The average mean squared estimation errors, associated with the
estimators $\widehat{f}_{n,N}$ and $\widehat{f}_{c}$ for this example and for
different values of $N$ and $c$ are given in Table 3. From these numerical
results, one concludes that both estimators have similar accuracy.
Nonetheless, our estimator is much faster than the KRR based estimator since
it involves the inversion of a relatively small size random matrix. This is
not the case for a KRR based estimator, specially when handling large sampling
dataset. As an example, for a moderate size of the sampling set $n=200,$ and
with the values of $N=c=10,$ the computing time required for the construction
of the KRR estimator $\widehat{f}_{c}$ is approximately 10 times the computing
time for our proposed estimator $\widehat{f}_{n,N}.$

$\sigma$ | $s$ | $N$ | $MSE(\widehat{f}_{n,N})$ | $c$ | $MSE(\widehat{f}_{c})$ | $s$ | $N$ | $MSE(\widehat{f}_{n,N})$ | $c$ | $MSE(\widehat{f}_{c})$
---|---|---|---|---|---|---|---|---|---|---
$0.1$ | $1.0$ | $10$ | $1.53e-1$ | $10$ | $1.60e-1$ | $2.0$ | $10$ | $1.71e-3$ | $10$ | $1.44e-3$
|  | $20$ | $5.88e-3$ | $20$ | $7.53e-3$ |  | $20$ | $2.00e-3$ | $20$ | $1.58e-3$
|  | $30$ | $3.94e-3$ | $30$ | $4.05e-3$ |  | $30$ | $2.93e-3$ | $30$ | $2.12e-3$
$0.05$ | $1.0$ | $10$ | $1.42e-1$ | $10$ | $1.29e-1$ | $2.0$ | $10$ | $8.23e-4$ | $10$ | $9.31e-4$
|  | $20$ | $4.86e-3$ | $20$ | $6.31e-3$ |  | $20$ | $6.35e-4$ | $20$ | $6.98e-4$
|  | $30$ | $1.82e-3$ | $30$ | $1.85e-3$ |  | $30$ | $6.55e-4$ | $30$ | $9.34e-4$
Table 3: Mean squared estimation errors for example 2.

Example 3: In this last example, we illustrate the performance of our LFR
estimator $\widehat{\beta}_{n,N}(\cdot).$ We consider the test LFR problem
that has been proposed in [10] and used later on for comparisons purposes in
[16]. For this test problem, $J=[0,1]$ and the slope function is given by

$\beta_{0}(s)=\sum_{j=1}^{50}c_{j}\varphi_{j}(s)=\sum_{j=1}^{50}4\frac{(-1)^{j+1}}{j^{2}}\varphi_{j}(s),\quad\varphi_{j}(s)=\left\\{\begin{array}[]{ll}1&\mbox{
if }j=1\\\ \sqrt{2}\cos(\pi js)&\mbox{ if }j\geq 1.\end{array}\right.$

The random predictor functional is given by ${\displaystyle
X(\cdot)=\sum_{j=1}^{50}\xi_{j}Z_{j}\varphi_{j}(\cdot),}$ where the
deterministic sequence $(\xi_{j})_{j}$ is given by
${\xi_{j}=\frac{(-1)^{j+1}}{j^{s/2}},\,\,s\geq 0}$ and the $Z_{k}$ are i.i.d.
random sample following the uniform law $U(-\sqrt{3},\sqrt{3}).$ The added
centred white noise $\varepsilon_{i}$ has variance $\sigma^{2}$ with
$\sigma=0.5.$ Then, we have computed our estimator
$\widehat{\beta}_{n,N}(\cdot)$ according our scheme given by (59)-(60), with
values of $N=50,$ $n=100,200,300$ and $s=1.5,2.0,4.0.$ We have computed the
average of the squared prediction and estimation error over $10$ realizations.
These errors are respectively given by

$E_{0}=\|\widehat{\beta}_{0}(\cdot)-\widehat{\beta}_{n,N}(\cdot)\|^{2}_{0}=\sum_{j}\frac{1}{j^{s}}\big{(}c_{j}-\widehat{c}_{j}\big{)}^{2},\quad
E_{2}=\|\widehat{\beta}_{0}-\widehat{\beta}_{n,N}\|^{2}_{2}=\sum_{j}\big{(}c_{j}-\widehat{c}_{j}\big{)}^{2}.$

$s$ | $n$ | $E_{0}$ | $E_{2}$ | $s$ | $n$ | $E_{0}$ | $E_{2}$ | $s$ | $n$ | $E_{0}$ | $E_{2}$
---|---|---|---|---|---|---|---|---|---|---|---
$1.5$ | $100$ | $2.62e-3$ | $7.72e-2$ | $2.0$ | $100$ | $1.52e-3$ | $6.94e-3$ | $4.0$ | $100$ | $2.02e-3$ | $3.32e-1$
| $200$ | $1.28e-3$ | $3.12e-2$ |  | $200$ | $1.11e-3$ | $5.17e-3$ |  | $200$ | $3.92e-4$ | $9.11e-2$
| $300$ | $3.15e-4$ | $1.62e-2$ |  | $300$ | $3.84e-4$ | $1.54e-3$ |  | $300$ | $2.51e-4$ | $1.28e-2$
Table 4: Prediction and estimation errors of the estimator
$\widehat{\beta}_{n,N}(\cdot)$ for $N=50$ and different values of $s$, $n.$

The obtained numerical values of these errors are given by Table 4. These
numerical results are coherent with the theoretical results given by Theorem
5. Moreover, by comparing our results with those obtained in [16] for this
same test problem, we conclude that our estimator
$\widehat{\beta}_{n,N}(\cdot)$ clearly outperforms in terms of accuracy and
convergence speed, the estimators proposed in the previous references.

###  5.2 A Real data Application: Covid-19 world’s countries spread.

In this last section, we use the publicly available daily updated dataset of
world’s countries spread of Covid-19, provided by Our World in Data. This
dataset contains the updated data starting from March 1st 2020 till October
16th and corresponding to the daily confirmed cases as well as the daily
deaths for different world’s countries. The URL of the dataset is
https://github.com/owid/covid-19-data/tree/master/public/data

Figure 1: Plots of the true Covid-19 daily cases (points) versus their Jacobi
based estimates (line) for (a) USA, (b) Algeria, (c) Italy, (d) India.

For the present simulations, we restrict ourselves to the data corresponding
to daily cases during the period March 1st-October 16th and corresponding to
four countries: USA, Algeria, Italy and India. The considered data covers
$m=594$ days and the aim is to build a non-parametric regression estimators
for these countries daily covid-19 cases. To apply our estimator based on
random pseudo-inverse associated with Jacobi polynomials, we have considered
the special values of the parameters $\beta=\alpha=-\frac{1}{2},$ $N=40.$
Then, we have considered a sampling set ${\displaystyle
S=\\{X_{i},\,i=1,\ldots,n\\}}$ with $n=340$ random samples in the interval
$[0,1]$ and following the $B(\alpha+1,\alpha+1)$ law. We found that the
average $2-$condition number of the small size $(N+1)\times(N+1)$ random
matrix $A_{N}=B_{N}^{\prime}B_{N}$ is given by $\kappa_{2}(A_{N})\approx
9.12.$ This ensures the stability of the proposed nonparametric regression
estimator, as predicted by Theorem 3. We have considered $n$ observed outputs
$Y_{i},$ given by $Y_{i}=f([mX_{i}]),$ where $f(k)$ is the observed number of
daily Covid-19 confirmed cases at day $k$ with $1\leq k\leq 594.$ The
coefficients $(\widehat{c}_{k})_{0\leq k\leq N}$ of the Jacobi polynomials
based estimators $\widehat{f}_{n,N}$ are computed in a fast way by using
formula (4). Note that since the random sampling points belong to the interval
$J=[0,1],$ then our estimators are constructed by the use of the re-normalized
Jacobi polynomials given by ${\displaystyle
Q_{k}(x)=\frac{1}{\sqrt{2}}\widetilde{P}_{k}(2x-1),\,x\in J.}$ These last set
of Jacobi polynomials are orthonormal on $J.$ In order to make this estimator
robust (that is not sensitive to data outliers), we have applied the RANSAC
algorithm trick with $10$ iterations and we have chosen the estimator that
best fits the whole data set of $594$ measurements from the $n=340$ the size
of random sampling set used at each iterations. In Figure 1, we have plotted
in points, the graphs of the real data corresponding to the daily confirmed
cases, versus the graph of the associated Jacobi polynomials based estimator
$\widehat{f}_{n,N}(\cdot).$ These numerical simulation results are coherent
with the theoretical results of Theorem 3. In particular, the estimator
$\widehat{f}_{n,N}(\cdot)$ has been made robust by just using few iterations.

## References

  * [1] G. E. Andrews, R. Askey and R. Roy, Special Functions, Cambridge University Press , Cambridge, New York, 1999.
  * [2] D. Benelmadani,· K. Benhenni and · S. Louhichi, The reproducing kernel Hilbert space approach in nonparametric regression problems with correlated observations, Ann. Inst. Stat. Math., 72 (2020), 1479–1500.
  * [3] A. Bonami and A. Karoui, Random Discretization of the Finite Fourier Transform and Related Kernel Random Matrices, J. Fourier Anal. Appl., 26 (2020).
  * [4] B. Bousselmi, J. P. Dupuy and A. Karoui, Random orthogonal projections based schemes for solving nonparametric and linear functional regressions problems, submitted for publication (2020).
  * [5] M. L. Braun, Spectral Properties of the Kernel Matrix and their Relation to Kernel Methods in Machine Learning, 2005.
  * [6] Zi-Cai. Li, Cheng-Sheng Chien and Hung-Tsai. Huang, Effective condition number for finite difference method, Journal of Computational and Applied Mathematics, 198 (1) (2007), 208–235.
  * [7] F. Comte and V. G. Catalot, Regression function estimation as a partly inverse problem, Ann. Inst. Stat. Math., 72 (2020), 1023–1054.
  * [8] B. Funke and C. Palmes, A note on estimating cumulative distribution functions by the use of convolution power kernels, Stat. Probabil. Lett., 121 (2017), 90–98.
  * [9] B. Y. Guo, Jacobi Approximations in Certain Hilbert Spaces and Their Applications to Singular Differential Equations, J. Math. Anal. Appl., 243 (2000), 373-–408.
  * [10] P. Hall and J. L. Horowitz, Methodology and convergence rates for functional linear regression, Ann. Statist., 35 (1), (2007), 70–91.
  * [11] R. A. Horn and C. R. Johnson, Matrix Analysis, second edition, Cambridge University Press, 2013.
  * [12] P. Jaming, A. Karoui and S. Spektor, The Approximation of Almost Time and Band-limited Functions by Their Expansion in Some Orthogonal Polynomials Bases, J. Approx. Theory, 212 (2016), 41–65.
  * [13] J. Jung and L. Sael, Fast and accurate pseudoinverse with sparse matrix reordering and incremental approach, Mach. Learn., 109 (2020), 2333-–2347.
  * [14] J. Prasad and H. Hayashi, On the uniform approximation of smooth functions by Jacobi polynomials, Can. J. Math., 25 (1) (1973), 216–223.
  * [15] Y. Sharon, J. Wright and Y. Ma, Minimum sum of distances estimator: Robustness and stability, Proc. IEEE ACC. St., (2009), 524–530.
  * [16] H. Shin and S. Lee, An RKHS approach to robust functional linear regression, Stat. Sinica, 26 (2016), 255–272.
  * [17] S. Smale and D. X. Zhou, Learning Theory Estimates via Integral Operators and Their Approximations, Constructive Approximation 26 (2) (2007), 153–172.
  * [18] S. Smale and D. X. Zhou, Shannon sampling II: Connections to learning theory, Appl. Comput. Harmon. Anal. 19 (2005), 285–302.
  * [19] J. A. Tropp, An Introduction to Matrix Concentration Inequalities, Foundations and Trends in Machine Learning series, 8 No. 1–2, Now Publishers Inc., 2015.
  * [20] X. Zhao, L. L. Wang and Z. Xie, Sharp error bounds for Jacobi expansions and Gegenbauer-Gauss quadrature of analytic functions, SIAM. J. Numer. Anal., 51 (3) (2013), 1443–1469.
  * [21] G. Zhao and Y. Ma, Robust nonparametric kernel regression estimator, Stat. Probab. Lett., 116 (2016), 72–79.
  * [22] A. Cohen, M.A. Davenport and D.Leviatan, On the stability and accuracy of least square approximations, Found. Comput. Math., 13 (5) (2013), 819–834.