# Extract the information from the big data with randomly distributed noise

Cheng Jin1,2 1\. School of Mathematical Sciences, Fudan University, Shanghai
200433, China.  2\. Shanghai Key Labororary of Contempory Appplied
Mathematics, Shanghai, China.  ,  Zhang Jiantang1 and  Zhong Min3,4,∗ 3\.
School of Mathematics, Southeast University, 210096, Nanjing, Jiangsu
Province.  4\. Nanjing Center for Applied Mathematics, 211135, Nanjing,
Jiangsu Province.  Corresponding author<EMAIL_ADDRESS>
###### Abstract.

In this manuscript, a purely data driven statistical regularization method is
proposed for extracting the information from big data with randomly
distributed noise. Since the variance of the noise maybe large, the method can
be regarded as a general data preprocessing method in ill-posed problems,
which is able to overcome the difficulty that the traditional regularization
method unable to solve, and has superior advantage in computing efficiency.
The unique solvability of the method is proved and a number of conditions are
given to characterize the solution. The regularization parameter strategy is
discussed and the rigorous upper bound estimation of confidence interval of
the error in $L^{2}$ norm is established. Some numerical examples are provided
to illustrate the appropriateness and effectiveness of the method.

##  1\. Introduction

As a rapid development of data collection instruments and information
technology, the generation of big data arises from applications in different
fields as applied mathematics, computer science, geology, biology,
engineering, and even business studies. When deal with these big data, there
are mainly two types of problems. Firstly, there will always be some redundant
data hidden behind, which contains too little effective information or
repeated information. Because of the limitation of computing capacity, too
much attention on these redundant data will seriously affect the storage and
time. Secondly, due to the inevitable measurement errors in the observation,
the big data is generated with randomly distributed noise, the variance of
these random errors can not be very small. Thus, a direct reconstruction based
on these noisy data will lead to unsatisfactory results. To overcome these
problems, it is necessary to propose a method satisfying the following
requirements: on one hand, the approximation should be accurate and stable,
and has strong computing capacity, and low time cost. On the other hand, the
method should also guarantee that the rate of change of the objective function
is not too large, i.e., the accuracy of derivative.

Determining the function $y(x)$ and the first order or higher order
derivatives of $y(x)$ from the random noisy samples of the function values is
called numerical differentiation, which has a widely utilization in many
practical problems. For example, the determination of discontinuous points in
image processing [5], the solution of Abel integral equation [12], and some
inverse problems arising from mathematical and physical equations [13], etc.
There have been many works concerning the convergence analysis of the
numerical algorithms [8, 13, 18], some different methods have been used to get
numerical results [9, 10, 11, 17].

The numerical differentiation is a classical inverse problem in the sense of
unstable dependence of solutions on small perturbation of the data, different
regularization methods for treating such ill-posed problems in one dimension
or higher dimensions were discussed in [2, 3, 6, 7, 14, 15, 16, 21, 22, 23].
However, these methods are not suitable for processing big data with large
variation random noise, this is because they are based on accurate knowledge
of the noise level $\delta$ or extract this information from the nature of the
problem. The prediction or extraction of noiselevel will always over-estimated
or under-estimated, thus the reconstruction accuracy is mainly limited. The
following figure describes the difficulties we encountered, obviously, in the
context of big data, effective information has been hidden, and traditional
regularization methods are unable to handle it.

Figure 1. The noisy samples and exact function for $0.1\times\sin(\pi x)$ with
$\sigma^{2}=0.002$.

There are also several statistics data processing techniques [4, 5, 19, 20],
in which the noise is assumed to be independent and identically distributed,
the reconstruction results converges to sought solution if the degree of
freedom tends to infinity. This greatly increases computational burden when
sample size gets large.

Therefore, it is necessary to consider purely data driven numerical method and
parameter choice rule. In this manuscript, we propose a new simple statistical
method to do the numerical differentiation. One innovation is this method is
can overcome the difficulty that the big data with large variation. What we do
is to separate the big data into $M$ groups, the average values of the $N$
noisy data in each group is calculated and regarded as a processed data. Based
on basic property of distribution, the original variance will be decreased by
$N$ times. The other innovation of our method is superior advantage in
computing, since the dimension and the size of the big data will be greatly
induced after processing. This makes our algorithm run in a relatively small
scale, which greatly saves the storage and improves the efficiency. In
addition, the algorithm can also be utilized as a data preprocessing method
and be applied in different kind of inverse problems.

The manuscript is organize as follows. In section 2, we show the unique
solvability of the method and then present a number of conditions
characterizing the solution. In section 3, we provide the choice strategy for
regularization parameter $\alpha$, and establish the rigorous upper bound
estimation of confidence interval of the error in $L^{2}$ norm, the optimal
choice of $M$ and $N$ will be discussed as well. The algorithm was given in
section 4, and some numerical examples are provided to illustrate the
effectiveness and the computational performance in section 5.

##  2\. Algorithm description and some basic results

###  2.1. Formulation of the problem

Suppose $y=y(x)$ is a function defined on $[0,1]$, and
$0=\hat{x}_{0}<\hat{x}_{1}<\ldots<\hat{x}_{L}=1$ is a uniform grid on $[0,1]$
with meshsize $h=1/L$. Given the noisy observation samples
$\widetilde{y}_{j}=y(\hat{x}_{j})+\eta_{j}$, with the error

$\displaystyle\eta_{j}\sim\mathcal{N}(0,\sigma^{2})\,,$ (2.1)

where the $\sigma^{2}$ denotes the variance of observation noise, we are
interested in reconstructing a function $f(x)$ such that the derivatives of
$f(x)$ are approximations of derivatives of function $y(x)$. Without losing
the generality, we may assume the observations on endpoints are exact, i.e.,

$\displaystyle\widetilde{y}_{0}=y(0)\,,\quad\widetilde{y}_{L}=y(1)\,.$

The main difficulties we encounter in solving such ill-posed numerical
differential problem are the observation size $L$ can be very large, and the
variance $\sigma^{2}$ can be large as well. To overcome these problems, given
a positive integer $M>1$, we separate the big data into $M$ groups. In each
group, there are $N=L/M$ observation points, i.e.,

$\displaystyle G_{1}:=\\{\hat{x}_{1}<\hat{x}_{2}<\ldots<\hat{x}_{N}\\}\,,$
$\displaystyle G_{2}:=\\{\hat{x}_{N+1}<\hat{x}_{N+2}<\ldots<\hat{x}_{2\times
N}\\}\,,$ $\displaystyle\ldots$ $\displaystyle G_{M}:=\\{\hat{x}_{(M-1)\times
N+1}<\hat{x}_{(M-1)\times N+2}<\ldots<\hat{x}_{M\times N}\\}\,.$ (2.2)

Then, we denote the new uniform grid

$\displaystyle\triangle=\\{0=x_{0}<x_{1}<\ldots<x_{M}=1\\}\,$ (2.3)

with $x_{i}=\hat{x}_{iN}(0\leq i\leq M)$ and corresponding mesh size
$h_{M}=1/M$. The new observation samples $\widetilde{Y}_{i}$ is defined as
sample mean of the group $G_{i}$:

$\displaystyle\widetilde{Y}_{i}=\frac{\widetilde{y}_{(i-1)\times
N+1}+\cdots+\widetilde{y}_{i\times N}}{N}\,.$ (2.4)

Let $k>1$ be an integer and $M,N>k$, we denote

$\displaystyle\mathcal{H}:=\\{f|f\in H^{k}(0,1),\ f(0)=y(0),\ f(1)=y(1)\\}\,,$

in which $H^{k}(0,1)$ be the usual Sobolev space consisting of all
$L^{2}(0,1)-$integrable functions whose $k-$order weak derivatives are also
$L^{2}(0,1)$ integrable. Define the cost functional

$\displaystyle\Phi(f):=\frac{1}{M}\sum_{i=1}^{M}\left(\widetilde{Y}_{i}-M_{i}(f)\right)^{2}+\alpha\|f^{(k)}\|_{L^{2}(0,1)}^{2}\,,$
(2.5)

where $\alpha>0$ is regularization parameter, and the

$\displaystyle M_{i}(f)=\frac{1}{h_{M}}\int_{x_{i-1}}^{x_{i}}f(x)dx\,$

is the average value of the function $f$ over the interval $[x_{i-1},x_{i}]$.
It can be shown there exists a unique minimizer $f_{*}\in\mathcal{H}$ for the
functional $\Phi$. The minimizer will be solutions of the numerical
differential problem. In following section 3, a number of conditions are given
to characterize the solution.

######  Remark 2.1.

The similar Tikhonov minimization functional can be found in [14], in which
the author also utilized the average value in the data fitting term, but
$\|f^{\prime}\|_{L^{2}(0,1)}$ in penalty term.

###  2.2. Some results in statistics

In this subsection, we introduce some basic results in statistics.

######  Definition 2.2.

We call elementary probability space a triplet $(\Omega,\mathcal{F},P)$ where
$\Omega$ is a set. $\mathcal{F}$ is an algebra of subsets of $\Omega$ and $P$
is a set of function defined on $\mathbb{R}$ with values in the real interval
$[0,1]$, which satisfies the following relations

  1. (1)

$P(\Omega)=1$,

  2. (2)

if $A_{1},A_{2},\ldots,A_{n}\in\mathcal{F}$ and $A_{i}\bigcap A_{j}=\emptyset$
for $i\neq j$, then

$\displaystyle P(\bigcup_{k=1}^{n}A_{k})=\sum_{k=1}^{n}P(A_{k})\,.$

######  Definition 2.3.

Given the probability space $(\Omega,\mathcal{F},P)$, we call a random
variable a real valued function $\psi:\Omega\rightarrow\mathbb{R}$ such that
$\psi^{-1}(B)\in\mathcal{F}$ for every Borel set $B$.

######  Theorem 2.4.

[Lindeberg-Levy: Central Limit Theorem for iid variables] Let
$X_{1},X_{2},\cdots,X_{n}$ are independent identically distributed (iid)
random variables with finite mean $\mu$ and variance $\sigma^{2}$, then the
random variable

$\displaystyle
Y_{n}=\frac{\sum_{i=1}^{n}X_{i}-n\mu}{\sigma\sqrt{n}}\rightarrow
Z\,,\quad\textrm{as}\ n\rightarrow\infty\,,\ \textrm{in distrbution}\,,$

where the $Z$ denotes the standardized Gaussian random variable whose
probability density function is

$\displaystyle
F_{Z}(x)=\frac{1}{\sqrt{2\pi}}\int_{\infty}^{x}e^{-t^{2}/2}dx\,.$

The theorem can be proving by recalling Levy’s theorem and using character
functions. In different words, the statement of theorem can be expressed by
saying that the variable $\sum_{i=1}^{N}X_{i}$ is asymptotically Gaussian with
mean $n\mu$ and variance $n\sigma^{2}$. This is ,in turn, implies that the
sample mean $\overline{X}:=\frac{1}{n}\sum_{i=1}^{n}X_{i}$ is asymptotically
normal with mean $E(\overline{X})=\mu$ and variance
$D(\overline{X})=\frac{\sigma^{2}}{n}$ and standard deviation
$\frac{\sigma}{\sqrt{n}}$. This fact, we will see, often plays an important
role in cases where a large number of elements is involved.

######  Theorem 2.5.

[Markov’s inequality] If $X$ is any nonnegative random variable and $a>0$,
then

$\displaystyle P(X\geq a)\leq\frac{E(X)}{a}\,.$ (2.6)

More generally, if $\varphi:[0,+\infty)\rightarrow\mathbb{R}$ is a nonnegative
increasing function for $a>0$ and $\varphi(a)>0$, we have

$\displaystyle P(X\geq
a)=P(\varphi(X)\geq\varphi(a))\leq\frac{E(\varphi(X))}{\varphi(a)}\,.$ (2.7)

We also need definition and properties for chi-squared distribution which was
based on the definition of Gamma distribution.

######  Definition 2.6.

Given real numbers $r>0$ and $\lambda>0$, the random variable $X$ is said to
have the gamma probability density function with parameters $r$ and $\lambda$
if

$\displaystyle f_{X}(x)=\frac{\lambda^{r}}{\Gamma(r)}x^{r-1}e^{-\lambda
x}\,,\quad y\geq 0\,,$

in which the $\Gamma(r)$ is gamma function defined by

$\displaystyle\Gamma(r)=\int_{0}^{\infty}x^{r-1}e^{-x}dx$

with properties

  1. (1)

$\Gamma(1)=1$ ,

  2. (2)

$\Gamma(r)=(r-1)\Gamma(r-1)$ ,

  3. (3)

If $r$ is an integer, $\Gamma(r)=(r-1)!$ .

A generalized gamma probability density function plays major role in
statistics, and it can have some different types. Among the most common of all
statistical analyses are procedures known as $\chi^{2}$ (chi-squared)
distribution, which is a special case of the gamma probability density
function with $\lambda=1/2$ and $r=m/2$, and $m$ is a positive integer denotes
the number of degrees of freedom.

######  Proposition 2.7.

Assume $X_{1},\ldots,X_{n}$ be iid random variables,
$X_{i}\sim\mathcal{N}(0,1)$. Then the random variable
$\chi^{2}=\sum_{i=1}^{n}X_{i}^{2}$ has a chi-square distribution with $n$
degrees of freedom denoted by $\chi^{2}\sim\chi^{2}(n)$. The $\chi^{2}-$
distribution have several properties

  1. (1)

The mean and the variance of the $\chi^{2}(n)$distribution are $n$ and $2n$ ,

  2. (2)

If $X\sim\chi^{2}(n_{1})$, $Y\sim\chi^{2}(n_{2})$ and $X,Y$ are independent
variables, then $X+Y\sim\chi^{2}(n_{1}+n_{2})$ ,

  3. (3)

The Cumulative distribution function is

$\displaystyle
F(x,n)=\frac{\gamma(\frac{n}{2},\frac{x}{2})}{\Gamma(\frac{n}{2})}$

in which $\gamma$ is incomplete gamma function defined as

$\displaystyle\gamma\left(\frac{n}{2},\frac{x}{2}\right)=\int_{0}^{x/2}t^{n/2-1}e^{-t}dt\,.$

##  3\. Main theoretical results

For the cost functional $\Phi$ defined in (2.5), we have

######  Theorem 3.1.

There exists one unique function $f_{*}\in\mathcal{H}$ such that

$\displaystyle\Phi(f_{*})\leq\Phi(f)$

for any functions $f\in\mathcal{H}$.

###### Proof.

Step 1 Construction of the minimizer $f_{*}$.

First we assume there exists a minimizer $f_{*}\in\mathcal{H}$. Define a
function $F(\lambda)=\Phi(f_{*}+\lambda g)$ where $g(x)\in H^{k}$ satisfies
$g(0)=g(1)=0$. Due to the minimality of $f_{*}$, we have $F^{\prime}(0)=0$.
Since

$\displaystyle\Phi(f_{*}+\lambda g)$
$\displaystyle=\frac{1}{M}\sum_{i=1}^{M}\left(\widetilde{Y}_{i}-M_{i}(f_{*}+\lambda
g)\right)^{2}+\alpha\|f_{*}^{(k)}+\lambda g^{(k)}\|_{L^{2}(0,1)}^{2}$
$\displaystyle=\frac{1}{M}\sum_{i=1}^{M}\left(\left(\widetilde{Y}_{i}-M_{i}(f_{*})\right)^{2}-2\lambda\left(\widetilde{Y}_{i}-M_{i}(f_{*})\right)M_{i}(g)+\lambda^{2}\left(M_{i}(g)\right)^{2}\right)$
$\displaystyle\quad+\alpha\left(\|f_{*}^{(k)}\|^{2}+2\lambda\int_{0}^{1}f_{*}^{(k)}g^{(k)}dx+\lambda^{2}\|g^{(k)}\|^{2}\right)\,.$

Therefore,

$\displaystyle
F^{\prime}(0)=\frac{1}{M}\sum_{i=1}^{M}\left(-2\left(\widetilde{Y}_{i}-M_{i}(f_{*})\right)M_{i}(g)\right)+2\alpha\int_{0}^{1}f_{*}^{(k)}g^{(k)}dx=0\,.$
(3.1)

Since

$\displaystyle\int_{0}^{1}f_{*}^{(k)}g^{(k)}dx$
$\displaystyle=f_{*}^{(k)}g^{(k-1)}|_{0}^{1}-\int_{0}^{1}f_{*}^{(k+1)}g^{(k-1)}dx$
$\displaystyle=\ldots=\sum_{j=0}^{k-1}(-1)^{j}f_{*}^{(k+j)}g^{(k-1-j)}|_{0}^{1}+(-1)^{k}\int_{0}^{1}f_{*}^{(2k)}gdx\,,$

in addition, $g\in H^{k}$ can be arbitrary, it can be verified that $f_{*}\in
H^{2k-1}(0,1)$ [1] and $f_{*}$ satisfies the following boundary conditions:

$\displaystyle f_{*}^{(k+j)}(0)=f_{*}^{(k+j)}(1)=0\,,\quad
j=0,1,\ldots,k-2\,.$ (3.2)

Therefore,

$\displaystyle F^{\prime}(0)$
$\displaystyle=\frac{1}{M}\sum_{i=1}^{M}\left(-2\left(\widetilde{Y}_{i}-M_{i}(f_{*})\right)M_{i}(g)\right)+2\alpha(-1)^{k}\int_{0}^{1}f_{*}^{(2k)}gdx$
$\displaystyle=-2\sum_{i=1}^{M}\int_{x_{i-1}}^{x_{i}}\left(\widetilde{Y}_{i}-M_{i}(f_{*})\right)gdx+2\alpha(-1)^{k}\int_{0}^{1}f_{*}^{(2k)}gdx=0\,.$

We plug any function $g\in C_{0}^{\infty}(x_{i-1},x_{i})$ into the equation,
it is equivalently

$\displaystyle\int_{x_{i-1}}^{x_{i}}\left[\left(M_{i}(f_{*})-\widetilde{Y}_{i}\right)+(-1)^{k}\alpha
f_{*}^{(2k)}\right]gdx=0\,.$

Consequently,

$\displaystyle
f_{*}^{(2k)}=\frac{(-1)^{k}}{\alpha}\left(\widetilde{Y}_{i}-M_{i}(f_{*})\right)\,,\quad\forall
x\in(x_{i-1},x_{i})$ (3.3)

which is a constant. It can be concluded that $f_{*}(x)|_{(x_{i-1},x_{i})}\in
P_{2k}(x_{i-1},x_{i})$, where $P_{2k}(x_{i-1},x_{i})$ represents the set of
all polynomials of one variable with the order no more than $2k$.

Considering the term $\int_{0}^{1}f_{*}^{(k)}g^{(k)}dx$ by integration by
parts again, for all $g\in C_{0}^{\infty}(0,1)$, since

$\displaystyle\int_{0}^{1}f_{*}^{(k)}g^{(k)}dx$
$\displaystyle=\sum_{i=1}^{M}\int_{x_{i-1}}^{x_{i}}f_{*}^{(k)}g^{(k)}dx$
$\displaystyle=\sum_{i=1}^{M}\left(\sum_{j=0}^{k-1}(-1)^{j}f_{*}^{(k+j)}g^{(k-1-j)}|_{x_{i-1}}^{x_{i}}\right)+(-1)^{k}\sum_{i=1}^{M}\int_{x_{i-1}}^{x_{i}}f_{*}^{(2k)}gdx\,.$

Substituting this into (3.1) and noting equation (3.3), we see

$\displaystyle\sum_{i=1}^{M}\left(\sum_{j=0}^{k-1}(-1)^{j}f_{*}^{(k+j)}g^{(k-1-j)}|_{x_{i-1}}^{x_{i}}\right)=0\,.$

Noting the condition (3.2) and $g(0)=g(1)=0$, it is equivalent that

$\displaystyle\sum_{j=0}^{k-1}(-1)^{j}\left(\sum_{i=1}^{M-1}\left(f_{*}^{(k+j)}(x_{i}-0)-f_{*}^{(k+j)}(x_{i}+0)\right)g^{(k-1-j)}(x_{i})\right)=0\,,$

which implies

$\displaystyle f_{*}^{(k+j)}(x_{i}+0)=f_{*}^{(k+j)}(x_{i}-0)\,,\quad
j=0,\ldots,k-1\,,\quad i=1,\ldots,M-1\,.$ (3.4)

Furthermore, observing that $f_{*}\in H^{2k-1}(0,1)$ which can be continuously
embedded into $C^{2k-2}[0,1]$, we add another $k$ conditions that

$\displaystyle f_{*}^{(j)}(x_{i}+0)=f_{*}^{(j)}(x_{i}-0)\,,\quad
j=0,\ldots,k-1\,,\quad i=1,\ldots,M-1\,.$ (3.5)

Combining the conditions (3.2)-(3.5), the solution $f_{*}(x)$ given as

$\displaystyle f_{*}(x)=c_{1}^{i}+c_{2}^{i}x+\cdots+c_{2k+1}^{i}x^{2k}\,,\quad
x\in(x_{i-1},x_{i}),\quad,i=1,2,\ldots M\,,$ (3.6)

where the $(2k+1)M$ parameters $c_{k}^{i}$, $k=1,2,\ldots 2k+1$,
$i=1,2,\ldots,M$, are determined by solving the following $(2k+1)M$ linear
equations:

$\displaystyle f_{*}^{(j)}(x_{i}+)-f_{*}^{(j)}(x_{i}-)$
$\displaystyle=0\,,\quad j=0,1,\ldots,2k-1,\ i=1,2,\ldots,M-1\,,$
$\displaystyle f_{*}^{(2k)}(x)$
$\displaystyle=\frac{(-1)^{k}}{\alpha}\left(\widetilde{Y}_{i}-M_{i}(f_{*})\right)\,,\quad
x\in(x_{i-1},x_{i})\,,\quad i=1,2,\ldots,M\,,$ $\displaystyle
f_{*}^{(k+j)}(0)$ $\displaystyle=0\,,\quad j=0,1,\ldots,k-2\,,$ $\displaystyle
f_{*}^{(k+j)}(1)$ $\displaystyle=0\,,\quad j=0,1,\ldots,k-2\,,$ $\displaystyle
f(0)$ $\displaystyle=y(0)\,,\quad f(1)=y(1)\,.$ (3.7)

Next, we will prove that these linear equations with respect the unknown
constants are uniquely solvable. We only need to show that the homogenous
equations have only the trivial solution. If we choose $\widetilde{Y}_{i}=0$
for $i=1,2,\ldots,M$, then the linear equations we obtained are just the
homogenous equations. On the other hand, by the definition of functional
$\Phi$, we know that $f_{*}=0$ is the unique minimzer. This means that the
homogenous equation have only the trivial solution.

Step 2: The uniqueness of the minimizer $f_{*}$.

For arbitrary $f\in\mathcal{H}$, we denote $g(x)=f(x)-f_{*}(x)$, it is obvious
that $g(0)=g(1)=0$ and

$\displaystyle\Phi(f)-\Phi(f_{*})$
$\displaystyle=\frac{1}{M}\left[\sum_{i=1}^{M}\left(\widetilde{Y}_{i}-M_{i}(f)\right)^{2}-\left(\widetilde{Y}_{i}-M_{i}(f_{*})\right)^{2}\right]+\alpha\left(\|f^{(k)}\|^{2}-\|f_{*}^{(k)}\|^{2}\right)$
$\displaystyle=\frac{1}{M}\sum_{i=1}^{M}\left(M_{i}(f_{*})-M_{i}(f)\right)\left(2\widetilde{Y}_{i}-M_{i}(f)-M_{i}(f_{*})\right)+\alpha\left(\|f^{(k)}\|^{2}-\|f_{*}^{(k)}\|^{2}\right)\,.$

It is worth to note that

$\displaystyle\|f^{(k)}\|^{2}-\|f_{*}^{(k)}\|^{2}$
$\displaystyle=\|f^{(k)}-f_{*}^{(k)}\|^{2}+2\int_{0}^{1}\left(f^{(k)}-f_{*}^{(k)}\right)f_{*}^{(k)}dx$
$\displaystyle=\|g^{(k)}\|^{2}+2\sum_{j=0}^{k-1}(-1)^{j}f_{*}^{(k+j)}g^{(k-1-j)}|_{0}^{1}+2(-1)^{k}\int_{0}^{1}f_{*}^{(2k)}gdx$
$\displaystyle=\|g^{(k)}\|^{2}+2(-1)^{k}\sum_{i=1}^{M}\int_{x_{i-1}}^{x_{i}}f_{*}^{(2k)}gdx$
$\displaystyle=\|g^{(k)}\|^{2}+2\sum_{i=1}^{M}\int_{x_{i-1}}^{x_{i}}\frac{1}{\alpha}\left(\widetilde{Y}_{i}-M_{i}(f_{*})\right)gdx$
$\displaystyle=\|g^{(k)}\|^{2}+\frac{2}{M}\sum_{i=1}^{M}\frac{1}{\alpha}\left(\widetilde{Y}_{i}-M_{i}(f_{*})\right)M_{i}(g)\,.$

Therefore,

$\displaystyle\Phi(f)-\Phi(f_{*})$
$\displaystyle=\frac{1}{M}\sum_{i=1}^{M}\left(M_{i}(f_{*})-M_{i}(f)\right)\left(2\widetilde{Y}_{i}-M_{i}(f)-M_{i}(f_{*})\right)$
$\displaystyle\quad+\alpha\|g^{(k)}\|^{2}+\frac{2}{M}\sum_{i=1}^{M}\left(\widetilde{Y}_{i}-M_{i}(f_{*})\right)M_{i}(g)$
$\displaystyle=\alpha\|g^{(k)}\|^{2}+\frac{1}{M}\sum_{i=1}^{M}\left(M_{i}(g)\right)^{2}\geq
0\,.$

It means that the function $f_{*}$ is a minimizer of the functional $\Phi(f)$.
If there is another function $f_{1}$ such that $\Phi(f_{*})=\Phi(f_{1})$, then
following the above process, we may have $f_{*}^{(k)}=f_{1}^{(k)}$ and thus
$f_{*}-f_{1}$ is a polynomial of degree $k-1$. In addition since
$M_{i}(f_{*}-f_{1})=0$, implying there exists $\xi\in(x_{i-1},x_{i})$ such
that $f_{*}(\xi)-f_{1}(\xi)=0$, i.e., there exists at least one root in
$(x_{i-1},x_{i})$ and there exist at lease $M$ roots in $(0,1)$. It is a
contradiction provided that $M>k$. Therefore, $f_{1}=f_{*}$ which means that
the minimizer $f_{*}$ is unique.

The proof is complete. ∎

######  Remark 3.2.

For the proof of Theorem 3.1, we know that $f_{*}\in\mathcal{H}$ is a
piecewise polynomial of degree $2k$ which can be determined by $(2k+1)M$
parameters.

In order to obtain the error estimate for the algorithm, we first define an
interpolation operator $Q_{h_{M}}$ from $L^{2}(0,1)$ onto the space of step
functions related to the subdivision $\triangle$ as follows. For a function
$g\in L^{2}(0,1)$, $Q_{h_{M}}g$ is given by

$\displaystyle Q_{h_{M}}g(x)=M_{i}(g),\quad\forall x\in(x_{i-1},x_{i})\,.$
(3.8)

We can obtain the error estimate for the operator $Q_{h_{M}}$ by the usual
scaling argument [c78, hx98], as described in the following result.

######  Lemma 3.3.

For all $g\in H^{1}(0,1)$,

$\displaystyle\|g-Q_{h_{M}}g\|_{L^{2}(0,1)}\leq
h\|g^{\prime}\|_{L^{2}(0,1)}\,.$

######  Lemma 3.4.

Let $y$ be a function in $H^{k}(0,1)$, and let $f_{*}$ be the solution of the
method (LABEL:minimization). Denote $e=f_{*}-y$, and denote the variable

$\displaystyle\triangle_{M}^{2}=\frac{1}{M}\sum_{i=1}^{M}(\widetilde{Y}_{i}-Y_{i})^{2}\,.$
(3.9)

Then, we have the following two estimates

$\displaystyle\|e^{(k)}\|\leq\sqrt{\frac{2\triangle_{M}^{2}}{\alpha}+\frac{2Q^{2}h_{M}^{2}}{\alpha
N^{2}}}+2\|y^{(k)}\|\,,$ (3.10)

and

$\displaystyle\|e\|\leq
h\|e^{\prime}\|+\sqrt{8\triangle_{M}^{2}+\frac{8Q^{2}h_{M}^{2}}{N^{2}}+2\alpha\|y^{(k)}\|^{2}}\,,$
(3.11)

in which $Q$ denote the upper bound of the $\|f^{\prime}\|_{L^{2}(0,1)}$ for
all $f\in H^{k}(0,1)$.

###### Proof.

Putting $y\in\mathcal{H}$ as candidate into the functional $\Phi$ in (2.5), by
the minimality of the $f_{*}$, we have

$\displaystyle\Phi(f_{*})\leq\Phi(y)$
$\displaystyle=\frac{1}{M}\sum_{i=1}^{M}\left(\widetilde{Y}_{i}-M_{i}(y)\right)^{2}+\alpha\|y^{(k)}\|^{2}$
$\displaystyle\leq\frac{2}{M}\sum_{i=1}^{M}\left((\widetilde{Y}_{i}-Y_{i})^{2}+\left(Y_{i}-M_{i}(y)\right)^{2}\right)+\alpha\|y^{(k)}\|^{2}\,.$

Since

$\displaystyle|Y_{i}-M_{i}(y)|$
$\displaystyle=\left|\frac{1}{h_{M}}\int_{x_{i-1}}^{x_{i}}ydx-\frac{y_{(i-1)\times
N+1}+\cdots+y_{i\times N}}{N}\right|$
$\displaystyle=\left|\frac{1}{h_{M}}\int_{x_{i-1}}^{\hat{x}_{(i-1)\times
N+1}}ydx-\frac{y(\hat{x}_{(i-1)\times
N+1})}{N}+\sum_{j=2}^{N}\left(\frac{1}{h_{M}}\int_{\hat{x}_{(i-1)\times
N+j-1}}^{\hat{x}_{(i-1)\times N+j}}ydx-\frac{y(\hat{x}_{(i-1)\times
N+j})}{N}\right)\right|$
$\displaystyle\leq\frac{1}{N}\sum_{j=1}^{N}\left|y(\xi_{j})-y(x_{\hat{(}i-1)\times
N+j})\right|$
$\displaystyle=\frac{1}{N}\sum_{j=1}^{N}\left|y(\hat{x}_{(i-1)\times
N+j})+(\xi_{j}-\hat{x}_{(i-1)\times
N+j})y^{\prime}(\eta_{j})-y(\hat{x}_{(i-1)\times N+j})\right|$
$\displaystyle\leq\frac{Qh_{M}}{N}\,,$

in which $\xi_{1}\in(x_{i-1},\hat{x}_{(i-1)\times N+1}),\
\xi_{j}\in(\hat{x}_{(i-1)\times N+j-1},\hat{x}_{(i-1)\times N+j})\
\textrm{for}\ j\geq 2$ and $\eta_{j}\in(\xi_{j},\hat{x}_{(i-1)\times N+j})$,
this yields

$\displaystyle\frac{2}{M}\sum_{i=1}^{M}|Y_{i}-M_{i}(y)|^{2}\leq\frac{2Q^{2}h_{M}^{2}}{N^{2}}\,.$
(3.12)

Therefore,

$\displaystyle\alpha\|(f_{*})^{(k)}\|^{2}\leq\Phi(f_{*})$
$\displaystyle=\frac{1}{M}\sum_{i=1}^{M}\left(\widetilde{Y}_{i}-M_{i}(f_{*})\right)^{2}+\alpha\|f_{*}^{(k)}\|^{2}$
$\displaystyle\leq
2\triangle_{M}^{2}+\frac{2Q^{2}h_{M}^{2}}{N^{2}}+\alpha\|y^{(k)}\|^{2}\,.$
(3.13)

Form this, it is subsequently that

$\displaystyle\|e^{(k)}\|\leq\|f_{*}^{(k)}\|+\|y^{(k)}\|\leq\sqrt{\frac{2\triangle_{M}^{2}}{\alpha}+\frac{2Q^{2}h_{M}^{2}}{\alpha
N^{2}}}+2\|y^{(k)}\|\,.$

On the other hand, noting that $Q_{h_{M}}$ is a projection operator, we obtain

$\displaystyle\|e\|^{2}=\int_{0}^{1}e^{2}dx$
$\displaystyle=\int_{0}^{1}e(e-Q_{h_{M}}e)dx+\int_{0}^{1}(Q_{h_{M}}e)^{2}dx$
$\displaystyle:=I_{1}+I_{2}\,.$

The $I_{1}$ can be estimated by using the Cauchy-Schwartz inequality and Lemma
3.3,

$\displaystyle\left|\int_{0}^{1}e(e-Q_{h_{M}}e)dx\right|\leq\|e\|\|e-Q_{h_{M}}e\|\leq
h\|e\|\|e^{\prime}\|\,.$

For the $I_{2}$, referring to (3.12) and (3) again, we have that

$\displaystyle\left|\int_{0}^{1}(Q_{h_{M}}e)^{2}dx\right|$
$\displaystyle=\sum_{i=1}^{M}\int_{x_{i-1}}^{x_{i}}(Q_{h_{M}}e)^{2}dx=h_{M}\sum_{i=1}^{M}M_{i}(e)^{2}$
$\displaystyle=h_{M}\sum_{i=1}^{M}\left(M_{i}(y)-M_{i}(f_{*})\right)^{2}$
$\displaystyle\leq\frac{2}{M}\left(\sum_{i=1}^{M}(M_{i}(y)-\widetilde{Y}_{i})^{2}+\sum_{i=1}^{M}(\widetilde{Y}_{i}-M_{i}(f_{*}))^{2}\right)$
$\displaystyle\leq\frac{4}{M}\left(\sum_{i=1}^{M}\left(M_{i}(y)-Y_{i}\right)^{2}\right)+4\triangle_{M}^{2}p+4\triangle_{M}^{2}+\frac{4Q^{2}h_{M}^{2}}{N^{2}}+2\alpha\|y^{(k)}\|^{2}$
$\displaystyle=8\triangle_{M}^{2}+\frac{8Q^{2}h_{M}^{2}}{N^{2}}+2\alpha\|y^{(k)}\|^{2}\,.$

Combining the above two estimates, we have

$\displaystyle\|e\|^{2}\leq
h\|e\|\|e^{\prime}\|+8\triangle_{M}^{2}+\frac{8Q^{2}h_{M}^{2}}{N^{2}}+2\alpha\|y^{(k)}\|^{2}\,,$

hence

$\displaystyle\|e\|\leq
h\|e^{\prime}\|+\sqrt{8\triangle_{M}^{2}+\frac{8Q^{2}h_{M}^{2}}{N^{2}}+2\alpha\|y^{(k)}\|^{2}}\,.$

∎

The above lemma shows that, when the regularization parameter $\alpha$ is
given, the estimates for the $\|e\|_{L^{2}(0,1)}$ and
$\|e^{(k)}\|_{L^{2}(0,1)}$ can be determined and controlled by the statistics
$\triangle_{M}^{2}$, which depend on the random statistics
$\widetilde{Y}_{i}$. Since

$\displaystyle\widetilde{Y}_{i}-Y_{i}\sim\mathcal{N}(0,\frac{\sigma^{2}}{N})\,,$

and due to the independence of the $\widetilde{Y}_{i}$, referring to
Proposition 2.7, we know

$\displaystyle\frac{MN}{\sigma^{2}}\triangle_{M}^{2}\sim\chi^{2}(M)\,.$

For a fixed $p\in(0,1)$, we denote the $(1-p)-$quantile of the $\chi^{2}(M)$
be $Z_{M,p}$, we provide the following lemma:

######  Lemma 3.5.

For $0<p<0.37$, there exists an upper bound for the estimate of
$\frac{Z_{M,p}}{M}$, denoted by $\bar{z}_{M,p}$, which can be determined by
the unique solution of the equation

$\displaystyle xe^{1-x}-p^{\frac{2}{M}}=0$

in $(1,+\infty)$. In addition, the $\bar{z}_{M,p}$ satisfies

$\displaystyle\bar{z}_{M-1,p}>\bar{z}_{M,p}\,,\ \textrm{and}\
\lim_{M\rightarrow\infty}\bar{z}_{M,p}=1\,.$ (3.14)

###### Proof.

Recalling the Markov inequality (2.7), let $\varphi(X)=\exp(tX)$, we have

$\displaystyle P(X\geq
a)=P(\exp(tX)\geq\exp(ta))\leq\frac{\mathbb{E}(\exp(tX))}{\exp(ta)}\,,\quad\forall\
t>0\,,$ (3.15)

in which the numerator of the right hand side of (3.15) is the moment
generating function. If the random variable $X\sim\chi^{2}(M)$, the moment
generating function is defined by

$\displaystyle\mathbb{E}(\exp(tX))=(1-2t)^{-\frac{M}{2}}\,.$

Let $a=\frac{Z_{M,p}}{M}:=z_{M,p}$, where $Z_{M,p}$ is the $(1-p)-$quantile of
the $\chi^{2}(M)$ distribution, the (3.15) becomes

$\displaystyle p=P(X\geq
z_{M,p}M)\leq\frac{(1-2t)^{-\frac{M}{2}}}{\exp(tz_{M,p}M)}\,,\quad\forall
t>0\,.$

Choosing $t=(z_{M,p}-1)/2z_{M,p}$, we have

$\displaystyle p=P(X\geq
z_{M,p}M)\leq\left(z_{M,p}\exp(1-z_{M,p})\right)^{\frac{M}{2}}\,.$ (3.16)

Considering the function $\psi(x)=x\exp(1-x)$, it is obvious that $\psi(x)$ is
strictly decreasing in $(1,+\infty)$ and $\psi(1)=1$,
$\lim_{x\rightarrow+\infty}\psi(x)=0$, therefore, the equation

$\displaystyle\psi(x)=p^{\frac{2}{M}}$

exists a unique solution $\overline{z}_{M,p}\in(1,+\infty)$. Recalling (3.16),
the $\overline{z}_{M,p}$ can be regarded as the upper bound for $z_{M,p}$,
this is because

$\displaystyle
p=\left(\overline{z}_{M,p}\exp(1-\overline{z}_{M,p})\right)^{\frac{M}{2}}\leq\left(z_{M,p}\exp(1-z_{M,p})\right)^{\frac{M}{2}}\,,$

which implies

$\displaystyle z_{M,p}\leq\overline{z}_{M,p}\,.$

In addition, $p^{2/M}$ is increased with respect to $M$, this combining with
the decreasing property of $\psi(x)$, we known $\overline{z}_{M,p}$ is also
decreased with respect to $M$, which yielding (3.14). Considering the
cumulative distribution function in Proposition 2.7, it is necessary demand
$F(z_{M,p}M,M)\geq F(M,M)$ such that $z_{M,p}\geq 1$. Since $F(M,M)$ is
decreased with respect to $M$, and $M$ is a positive integer with $M\geq 2$.
Therefore, we can demand $F(z_{M,p}M,M)\geq F(2,2)\approx 0.63$, yielding

$\displaystyle P(X\leq z_{M,p}M)=F(z_{M,p}M,M)=1-p\geq 0.63\,,$

thus $0<p<0.37$. ∎

######  Remark 3.6.

for fixed $0<p<0.37$, since

$\displaystyle P(X\leq\overline{z}_{M,p}M)\geq P(X\leq z_{M,p}M)=1-p\,,$

for $X=\frac{NM\triangle_{M}^{2}}{\sigma^{2}}\sim\chi^{2}(M)$, the estimate

$\displaystyle\triangle_{M}^{2}\leq\frac{\overline{z}_{M,p}\sigma^{2}}{N}$
(3.17)

is satisfied with probability of at least $1-p$, i.e.,

$\displaystyle
P(\triangle_{M}^{2}\leq\frac{\overline{z}_{M,p}\sigma^{2}}{N})\geq 1-p\,.$

One important thing is how to choose the regularization parameter $\alpha$ in
the functional $\Phi$ so that the minimizer ca be one possible solution of the
numerical differentiation problem. Our consideration is taking
$\alpha=\frac{\bar{c}\sigma^{2}}{N}$ with a constant $\bar{c}$. This is
motivated by the results in previous work in [3]. On one hand, variance
describes the fluctuation level of the random variable, so choosing $\alpha$
be the same order with $\sigma^{2}$ is an intuitive consideration. On the
other hand, based on the central limit theorem, taking the sample mean on a
certain interval as the observation value, the corresponding error variance
will be reduced according to the speed of $1/N$, thus $\alpha$ should also
reduced.

Based on the above discussions, we will establish the convergence results. The
important Sobolev inequality is necessarily be given before the convergence
theorem.

######  Lemma 3.7.

Let $-\infty\leq a<b\leq\infty,1\leq p<\infty$ and $0<\varepsilon_{0}<\infty$,
$f$ is a function that has $m-$th order derivative in $(0,1)$. There exists a
constant $K$ which depends on $\varepsilon_{0}$, $p$, and such that for every
$\varepsilon$, $0<\varepsilon\leq\varepsilon_{0}$, $0\leq j<m$, we have

$\displaystyle\int_{0}^{1}|f^{(j)}|^{p}dt\leq
K\varepsilon\int_{0}^{1}|f^{(m)}|^{p}dt+K\varepsilon^{\frac{-j}{m-j}}\int_{0}^{1}|f|^{p}dt\,.$

######  Theorem 3.8.

Suppose $y=y(x)\in H^{k}(0,1)$ is a function defined on $[0,1]$, given the
noisy observation samples $\widetilde{y}_{j}=y(\hat{x}_{j})+\eta_{j}$ with
$1\leq j\leq L$ of the function $y(x)$ satisfying

$\displaystyle\eta_{i}\sim N(0,\sigma^{2})\,.$ (3.18)

Dividing the samples into $M$ groups as in (2.1) and define the new grid
$\triangle$ as in (2.3). Let $\widetilde{Y}_{i}$ with $1\leq i\leq M$ and
statistics $\triangle_{M}^{2}$ are defined in (2.5) and (3.9) respectively.
For fixed $0<p<0.37$, let $\bar{z}_{M,p}$ be the unique solution for the
function $x\exp(1-x)=p^{2/M}$ on interval $(1,+\infty)$ defined in Lemma 3.5,
then $\triangle_{M}^{2}$ satisfies

$\displaystyle\triangle_{M}^{2}\leq\frac{\overline{z}_{M,p}\sigma^{2}}{N}$

with the probability of at least $1-p$. Let $e=f_{*}-y$, when choosing the
regularization parameter $\alpha=\frac{\bar{c}\sigma^{2}}{N}$, the $L^{2}$
norm of the $e$ and $e^{(k)}$ satisfies the following estimates

$\displaystyle\|e^{(k)}\|\leq\sqrt{\frac{2\overline{z}_{M,p}}{\bar{c}}+\frac{2Q^{2}h_{M}^{2}}{\bar{c}N\sigma^{2}}}+2\|y^{(k)}\|\,,$
(3.19) $\displaystyle\|e\|\leq
h\|e^{\prime}\|+\sqrt{\frac{8\overline{z}_{M,p}\sigma^{2}+2\bar{c}\sigma^{2}\|y^{(k)}\|^{2}}{N}+\frac{8Q^{2}h_{M}^{2}}{N^{2}}}\,,$
(3.20)

with the probability of at least $1-p$. Therefore, for fixed $0\leq j\leq
k-1$, assume $\sigma^{2}>h_{M}^{2}/N$, pthe $L^{2}$ norm of $e^{(j)}$
satisfying the estimate

$\displaystyle\|e^{(j)}\|_{2}\leq
C_{1}h_{M}^{k-j}+C_{2}\left(\frac{\sigma^{2}}{N}\right)^{(k-j)/2k}\,$ (3.21)

with the probability of at least $1-p$. The constants $C_{1}$ and $C_{2}$ are
two constants independent of $h_{M}$, $\sigma$ and $N$.

###### Proof.

The (3.19) and (3.20) is directly obtained from (3.10) and (3.11) by taking
the estimate $\triangle_{M}^{2}$ and $\alpha$. We then prove the (3.21) for
the case $j=1$, i.e.,

$\displaystyle\|e^{\prime}\|_{2}\leq
C_{1}h_{M}^{k-1}+C_{2}\left(\frac{\sigma^{2}}{N}\right)^{\frac{k-1}{2k}}\,,$

Taking $j=1$, $m=k$, $p=2$, $\varepsilon_{0}=1$ and $f=e$ in Lemma 3.7, we can
obtain

$\displaystyle\|e^{\prime}\|^{2}\leq
K\varepsilon\|e^{(k)}\|^{2}+K\varepsilon^{-\frac{1}{k-1}}\|e\|^{2}\,.$

Without losing the generality, we assume that
$\|e\|^{2(k-1)/k}\leq\varepsilon_{0}=1$. Taking
$\varepsilon=\|e\|^{2(k-1)/k}$, it follows that

$\displaystyle\|e^{\prime}\|^{2}\leq K(\|e^{(k)}\|^{2}+1)\|e\|^{2(k-1)/k}\,.$

Referring to the estimates for both $\|e\|$ and $\|e^{(k)}\|$ in (3.19) and
(3.20), we have

$\displaystyle\|e^{\prime}\|^{2}\leq
K\left(h\|e^{\prime}\|+\sqrt{\frac{8\overline{z}_{M,p}\sigma^{2}+2\bar{c}\sigma^{2}\|y^{(k)}\|^{2}}{N}+\frac{8Q^{2}h_{M}^{2}}{N^{2}}}\right)^{2(k-1)/k}\left(1+\frac{2\overline{z}_{M,p}}{\bar{c}}+\frac{2Q^{2}h_{M}^{2}}{\bar{c}N\sigma^{2}}+4\|y^{(k)}\|^{2}\right)\,.$

It is equivalently that

$\displaystyle\|e^{\prime}\|^{k}$ $\displaystyle\leq
K^{k/2}\left(h_{M}\|e^{\prime}\|+\sqrt{\frac{8\overline{z}_{M,p}\sigma^{2}+2\bar{c}\sigma^{2}\|y^{(k)}\|^{2}}{N}+\frac{8Q^{2}h_{M}^{2}}{N^{2}}}\right)^{k-1}\left(1+\frac{2\overline{z}_{M,p}}{\bar{c}}+\frac{2Q^{2}h_{M}^{2}}{\bar{c}N\sigma^{2}}+4\|y^{(k)}\|^{2}\right)^{k/2}$
$\displaystyle\leq
K^{*}\left(1+\frac{2\overline{z}_{M,p}}{\bar{c}}+\frac{2Q^{2}h_{M}^{2}}{\bar{c}N\sigma^{2}}+4\|y^{(k)}\|^{2}\right)^{k/2}$
$\displaystyle\qquad\cdot\left(\|e^{\prime}\|^{k-1}h_{M}^{k-1}+\left(\frac{\left(8\overline{z}_{M,p}+2\bar{c}\|y^{(k)}\|^{2}\right)\sigma^{2}}{N}+\frac{8Q^{2}h_{M}^{2}}{N^{2}}\right)^{(k-1)/2}\right)\,,$

in which the constant $K^{*}:=2^{k-1}K^{k/2}$.

Next we will show that (3.21) can be obtained from the above estimate. We
consider two cases.

Case1: Assume $\|e^{\prime}\|\leq
K^{*}\left(1+\frac{2\overline{z}_{M,p}}{\bar{c}}+\frac{2Q^{2}h_{M}^{2}}{\bar{c}N\sigma^{2}}+4\|y^{(k)}\|^{2}\right)^{k/2}h_{M}^{k-1}$,
since $\overline{z}_{M,p}$ can be bounded, yielding

$\displaystyle\|e^{\prime}\|$ $\displaystyle\leq
K^{*}\left(K_{1}+K_{2}\frac{h_{M}^{2}}{N\sigma^{2}}\right)^{k/2}h_{M}^{k-1}\leq
C_{1}h_{M}^{k-1}\,,$

in which the constants
$K_{1}:=1+\frac{2\overline{z}_{M,p}}{\bar{c}}+4\|y^{(k)}\|^{2}$,
$K_{2}=2Q^{2}/\bar{c}$ and $C_{1}=K^{*}(K_{1}+K_{2})^{k/2}$. The last
inequality is based on the assumption $\sigma^{2}>h_{M}^{2}/N$.

Case2: Assume
$\|e^{\prime}\|>K^{*}\left(1+\frac{2\overline{z}_{M,p}}{\bar{c}}+\frac{2Q^{2}h_{M}^{2}}{\bar{c}pN\sigma^{2}}+4\|y^{(k)}\|^{2}\right)^{k/2}h_{M}^{k-1}$.
Then we can take

$\displaystyle
r=\|e^{\prime}\|-K^{*}\left(1+\frac{2\overline{z}_{M,p}}{\bar{c}}+\frac{2Q^{2}h_{M}^{2}}{\bar{c}N\sigma^{2}}+4\|y^{(k)}\|^{2}\right)^{k/2}h_{M}^{k-1}>0\,,$

then,

$\displaystyle r^{k}$ $\displaystyle\leq\|e^{\prime}\|^{k-1}r$
$\displaystyle\leq\|e^{\prime}\|^{k}-K^{*}\left(1+\frac{2\overline{z}_{M,p}}{\bar{c}}+\frac{2Q^{2}h_{M}^{2}}{\bar{c}N\sigma^{2}}+4\|y^{(k)}\|^{2}\right)^{k/2}h_{M}^{k-1}\|e^{\prime}\|^{k-1}$
$\displaystyle\leq
K^{*}\left(1+\frac{2\overline{z}_{M,p}}{\bar{c}}+\frac{2Q^{2}h_{M}^{2}}{\bar{c}N\sigma^{2}}+4\|y^{(k)}\|^{2}\right)^{k/2}\left(\frac{\left(8\overline{z}_{M,p}+2\|y^{(k)}\|^{2}\right)\sigma^{2}}{N}+\frac{8Q^{2}h_{M}^{2}}{N^{2}}\right)^{(k-1)/2}$
$\displaystyle\leq
C_{1}\left(K_{3}\frac{\sigma^{2}}{N}+K_{4}\frac{h_{M}^{2}}{N^{2}}\right)^{(k-1)/2}$
$\displaystyle\leq C_{2}^{k}\frac{\sigma^{k-1}}{N^{(k-1)/2}}\,,$

where the constants $K_{3}=8\overline{z}_{M,p}+2\|y^{(k)}\|^{2}$,
$K_{4}=8Q^{2}$ and $C_{2}^{k}=C_{1}\max(K_{3},K_{4})^{(k-1)/2}$. Consequently,

$\displaystyle r\leq C_{2}\left(\frac{\sigma^{2}}{N}\right)^{(k-1)/2k}\,.$

The inequality (3.21) for $j=1$ is proved. By a similar method, we can prove
for any $0\leq j<k$, it holds that

$\displaystyle\|e^{(j)}\|_{2}\leq
C_{1}h_{M}^{k-j}+C_{2}\left(\frac{\sigma^{2}}{N}\right)^{(k-j)/2k}\,.$

We will not give the detailed proof here.

∎

######  Remark 3.9.

For particular $k=2$, the estimate (3.21) becomes

$\displaystyle\|e^{\prime}\|_{2}\leq
C_{1}h_{M}+C_{2}\left(\frac{\sigma^{2}}{N}\right)^{1/4}\,.$

When $\sigma$ is fixed, pay attention that $h_{M}=N/J$, this motivate us that
$N=\mathcal{O}(J^{4/5})$ and $M=\mathcal{O}(J^{1/5})$ is the optimal choice.

##  4\. Algorithm

We consider the simple case $k=2$. The other cases can be treated in a similar
way. Since $f_{*}$ is a piece wise polynomial of degree four, we assume that

$\displaystyle
f_{*}(x)=a_{i}+b_{i}(x-x_{i})+c_{i}(x-x_{i})^{2}+d_{i}(x-x_{i})^{3}+e_{i}(x-x_{i})^{4},\quad\textrm{for}\
x\in[x_{i},x_{i+1})\,,$ (4.1)

where there are $5M$ constants $a_{i}$, $b_{i}$, $c_{i}$, $d_{i}$, $e_{i}$ for
$i=0,\ldots,M-1$.

From our reconstruction, we have

$\displaystyle\frac{c_{i-1}-2c_{i}+c_{i+1}}{3h}=3(e_{i-1}+e_{i})h_{M},\quad\textrm{for}\
i=1,\ldots,M-1\,.$

Since $f^{\prime\prime}(0)=f^{\prime\prime}(1)=0$ yields $c_{0}=c_{M}=0$,
denote $c=(c_{1},\ldots,c_{M-1})^{T}$ and
$e=(e_{0},e_{1},\ldots,e_{M-1})^{T}$, we have

$\displaystyle c=6h^{2}(A_{1})^{-1}B_{1}e:=Ce\,,$ (4.2)

with

$\displaystyle A_{1}=\left(\begin{matrix}-2&1&0&\cdots\\\ 1&-2&1&\cdots\\\
\vdots&\vdots&\vdots&\vdots\\\ \cdots&1&-2&1\\\ \cdots&\cdots&1&-2\\\
\end{matrix}\right)_{(M-1)\times(M-1)}B_{1}=\left(\begin{matrix}1&1&0&\cdots\\\
0&1&1&\cdots\\\ \vdots&\vdots&\vdots&\vdots\\\ \cdots&\cdots&1&1\\\
\end{matrix}\right)_{(M-1)\times M}\,.$ (4.3)

Then, since

$\displaystyle
d_{i-1}=\frac{1}{3h_{M}}\left(c_{i}-c_{i-1}-6e_{i-1}h_{M}^{2}\right)\,,$

denote $d=(d_{0},d_{1}\ldots,d_{M-1})^{T}$, it follows that

$\displaystyle d=\frac{1}{3h}Tc-2he:=De\,,\textrm{where}\quad
T=\left(\begin{matrix}1&0&0&\cdots\\\ -1&1&0&\cdots\\\
\vdots&\vdots&\vdots&\vdots\\\ \cdots&\cdots&0&-1\\\
\end{matrix}\right)_{M\times(M-1)}$ (4.4)

Next, since

$\displaystyle
a_{i-1}-2a_{i}+a_{i+1}=(c_{i-1}+c_{i})h_{M}^{2}+(2d_{i-1}+d_{i})h_{M}^{3}+(3e_{i-1}+e_{i})h_{M}^{4},\quad\textrm{for}\
i=1,\ldots,M-1\,,$

this combine with $a_{0}=f(0)$ and $a_{M}=f(1)$, denote
$a=(a_{1},\ldots,a_{M-1})^{T}$, it follows that

$\displaystyle
A_{1}a=h_{M}^{2}Pc+h_{M}^{3}Qd+h_{M}^{3}Re-\frac{1}{h_{M}}v_{1}:=Ae-\bar{v}_{1}\,,$
(4.5)

where

$\displaystyle P=\left(\begin{matrix}1&0&0&\cdots\\\ 1&1&0&\cdots\\\
\vdots&\vdots&\vdots&\vdots\\\ \cdots&\cdots&1&1\\\
\end{matrix}\right)_{(M-1)\times(M-1)}Q=\left(\begin{matrix}2&1&0&\cdots\\\
0&2&1&\cdots\\\ \vdots&\vdots&\vdots&\vdots\\\ \cdots&\cdots&2&1\\\
\end{matrix}\right)_{(M-1)\times M}$ $\displaystyle
R=\left(\begin{matrix}3&1&0&\cdots\\\ 0&3&1&\cdots\\\
\vdots&\vdots&\vdots&\vdots\\\ \cdots&\cdots&3&1\\\
\end{matrix}\right)_{(M-1)\times M}v_{1}=\left(\begin{matrix}f(0)\\\ 0\\\
\vdots\\\ 0\\\ f(1)\\\ \end{matrix}\right)_{(M-1)\times 1}$ (4.6)

In addition, since

$\displaystyle
b_{i-1}=\frac{1}{h_{M}}(a_{i}-a_{i-1})-c_{i-1}h_{M}-d_{i-1}h_{M}^{2}-e_{i-1}h_{M}^{3}\,,\quad\textrm{for}\
i=1\ldots M-1$

Therefore,

$\displaystyle b=\frac{1}{h_{M}}Ta-h_{M}\left(\begin{matrix}0\\\ c\\\
\end{matrix}\right)-h_{M}^{2}d-h_{M}^{3}e-\frac{1}{h_{M}}v_{2}:=Be-\bar{v}_{2}\,,$
(4.7)

where

$\displaystyle R=\left(\begin{matrix}1&0&0&\cdots\\\ -1&1&0&\cdots\\\
\vdots&\vdots&\vdots&\vdots\\\ \cdots&\cdots&0&-1\\\
\end{matrix}\right)_{M\times(M-1)}v_{2}=\left(\begin{matrix}f(0)\\\ 0\\\
\vdots\\\ 0\\\ -f(1)\\\ \end{matrix}\right)_{M\times 1}$ (4.8)

Finally, let $e=(e_{0},e_{1},\ldots,e_{M-1})^{T}$, we have

$\displaystyle 24\alpha e=\widetilde{Y}-\left(\begin{matrix}f(0)\\\ a\\\
\end{matrix}\right)-\frac{h_{M}}{2}b-\frac{h_{M}^{2}}{3}\left(\begin{matrix}0\\\
c\\\ \end{matrix}\right)-\frac{h_{M}^{3}}{4}d-\frac{h_{M}^{4}}{5}e\,.$ (4.9)

This equivalent to

$\displaystyle\left(\left(24\alpha+\frac{h_{M}^{4}}{5}\right)I+\frac{h_{M}}{2}B+\frac{h_{M}^{3}}{4}D\right)e+\left(\begin{matrix}0\\\
\left(A+\frac{h_{M}^{2}}{3}C\right)e\\\
\end{matrix}\right)=\widetilde{Y}+\left(\begin{matrix}-f(0)\\\ \bar{v}_{1}\\\
\end{matrix}\right)+\frac{h_{M}}{2}\bar{v}_{2}\,.$ (4.10)

By solving (4.10) and using the matrix relations, we can get the coefficient
vectors.

##  5\. Numerical example

In this section, some numerical examples are provided to illustrate
computational performance of the method. The regularization parameter is
selected by a posteriori choice L-curve criterion. We divide $[0,1]$ into
$L=1000$ equal subintervals. For a given function
$y=y(x)=x^{3}+2x^{2}-0.5x+1$, we add random noise at each points with normal
distribution $N(0,\sigma^{2})$ with $\sigma^{2}=0.2$ and generate
corresponding noisy observations, see the Figure 2 (left). We divided the
points into $M=10$ groups, so there are $N=100$ points in every group. We
generate the $\\{\widetilde{Y}_{i}\\}_{i=1}^{M}$ be the sample mean of these
$N$ observations, see Figure 2 (right).

Figure 2. The noisy observations $\sigma^{2}=0.2$ (left), the average noisy
observations with $M=10$ (right).

In our numerical example, we fix $k=2$ and do the algorithm in section 4, the
regularization parameter $\alpha$ was suggested be chosen as
$\alpha=\frac{\bar{c}\sigma^{2}}{N}$, now we illustrate the performance of the
heuristic parameter L-curve strategy for determining the constant
$\bar{\alpha}$. We plot the curve
$(\log(\|(f_{*})^{\prime}\|_{L^{2}(0,1)}^{2}),\log(\textrm{residual}))$ with

$\textrm{residual}:=\sqrt{\frac{1}{M}\sum_{i=1}^{M}\left(\widetilde{Y}_{i}-M_{i}(f)\right)^{2}}\,.$

The curve indeed looks like the letter “L”, and we can get the corresponding
constant $\bar{c}=0.0239$.

Figure 3. The loglog figure of the Lcurve.

Then, we present the computational performance of the method and the
regularization parameter is chosen in terms of L-curve method, the
reconstruction for $y(x)$ and $y^{\prime}(x)$ are showed in Figure 4
respectively. Here and in what follows, the blue curve means the exact
function or the exact first derivative, while the red curve means the
reconstructions.

Figure 4. The reconstruction of $y(x)$ (left) and the reconstruction of
$y^{\prime}(x)$ (right) with $M=10$.

Next, we compare the reconstruction results under different choices for $M$,
recalling the Remark 3.9, the optimal choice is $M=\mathcal{O}(J^{1/5})$,
which means $M$ is approximately $5$. we choose $M=5,10,50,100,200$ for
comparison, in all situations, we fix the same constant $\bar{c}=0.0239$ and
thus, due to different $N$, the $\alpha=\bar{c}\sigma^{2}/N$ is also
different. The detailed information was shown in Table 1. It is also necessary
to notice that, the computational complexity was based on the value of $M$,
smaller $M$ represents smaller matrices sizes and thus cheaper computational
costs.

Table 1. The reconstruction comparison for different choices for $M$ $M$ | $N$ | $\|y-f_{*}\|_{L^{2}(0,1)}$ | $\|y-f_{*}\|_{l^{\infty}}$ | $\|y^{\prime}-f_{*}^{\prime}\|_{L^{2}(0,1)}$ | $\|y^{\prime}-f_{*}^{\prime}\|_{l^{\infty}}$
---|---|---|---|---|---
5 | 200 | 0.020805 | 0.036963 | 0.166882 | 0.745254
10 | 100 | 0.027061 | 0.045420 | 0.211428 | 0.815453
50 | 20 | 0.040842 | 0.067272 | 0.249623 | 1.023243
100 | 10 | 0.054420 | 0.085524 | 0.287166 | 1.150439
200 | 5 | 0.079110 | 0.116156 | 0.353859 | 1.333828

Finally, we compare our method with the previous algorithm in [3], in which
they discussed the approximation provided that the noisy observations
satisfying $\|y(x_{j})-y^{\delta}_{j}\|\leq\delta$ for $1\leq j\leq L$, and
the noiselevel $\delta$ should be known a priorily and the regularization
parameter was suggested be chosen as $\alpha=\delta^{2}$. In this example,
since the noise was given as a normal distribution with $\sigma^{2}$ be the
variance, it is appropriate to choose $\alpha=\sigma^{2}$.

Figure 5. The reconstruction of $y(x)$ (left) and the reconstruction of
$y^{\prime}(x)$ (right) using the method in [3].

The error of $\|y-f_{*}\|_{L^{2}(0,1)}$, $\|y-f_{*}\|_{l^{\infty}}$,
$\|y^{\prime}-f_{*}^{\prime}\|_{L^{2}(0,1)}$ and
$\|y^{\prime}-f_{*}^{\prime}\|_{l^{\infty}}$ are $0.059733$, $0.092025$,
$0.299401$ and $1.179538$ respectively. However, the computational complexity
is much higher, since the matrices involved in calculation is $1000\times
1000$.

The advantages of our method will be reflected when the observation data
volume is very large. We provide the second example when the given function is
$y(x)=1+10x^{2}(1-x)^{2}$. Assume $L=10^{6}$, we add random noise at each
points with normal distribution $N(0,\sigma^{2})$ with $\sigma^{2}=0.25$ and
generate corresponding noisy observations, see the Figure 6. This time, we
divided the points into $M=10$ groups, with $N=10^{5}$ points in every group,
and the regularization parameter $\alpha$ was chosen be $\alpha=\sigma^{2}/N$
with constant $\bar{\alpha}=1$ for simplicity. The reconstruction results were
plotted in Figure 6, it is very satisfactory, and the computational complexity
is very low. However, since the matrix size is too large to exceeds the
capacity of the Matlab, the method by [3] was failed to be utilized.

Figure 6. The noisy observations, the average noisy observations with $M=10$,
the approximation of $y(x)$ and the approximation of the $y^{\prime}(x)$.

##  6\. Acknowledgement

J. Cheng is supported by the NSFC (No.11971121), M. Zhong is supported by the
NSFC (No. 11871149) and supported by Zhishan Youth Scholar Program of SEU.

## References

  * [1] Adams, R.A., Sobolev spaces, Pure and Applied Mathematics, Vol. 65 (New York-London: Academic Press).
  * [2] Anderssen, R.S., Hegland, M., 1999, For numerical differentiation, dimensionality can be a blessing! Mathematics of Computing, 68(227), 1121-1141.
  * [3] Cheng J., Jia X.Z. and Wang Y.B., Numerical differentiation and its applications. Inverse Problems in Science and Engineering, 15 (2007) 339-357.
  * [4] Craven P., Wahba G., Smoothing noisy data with spline functions, Numerische Mathematik 31 (1978), 377-403.
  * [5] Deans, S.R., 1983, Radon Transform and its Applications (New York: A Wiley-Interscience Publication, John Wiley $\&$ Sons Inc.).
  * [6] Delhez E., A spline interpolation technique that preserves mass budgets. Appl. Math. Lett. 16 (2003) 16-26.
  * [7] Dinh N.H., Chuong L.H. and Lesnic D. Heuristic regularization methods for numerical differentiation. Computers and Mathematics with Applications, 63 (2012), 816-826.
  * [8] Groetsch C.W., Differentiation of approximately specified functions, Amer. Math. Monthly 98 (1991) 847-850.
  * [9] Groetsch C.W., Optimal order of accuracy in Vasins method for differentiation of noisy functions, J. Optim.Theory Appl. 74 (1992) 373-378.
  * [10] Groetsch C.W., Lanczos generalized derivative, Amer. Math. Monthly 105 (1998) 320-326.
  * [11] Groetsch C.W., Scherzer O., The optimal order of convergence for stable evaluation of differential operators, Electron. J. Differential Equations 4 (1993) 1-10.
  * [12] Gorenflo R., Vessella S., 1991, Abel Integral Equations, Analysis and Applications. Lecture Notes in Mathematics, Vol. 1461 (Berlin: Springer-Verlag).
  * [13] Hanke M., Scherzer O., Inverse problems light: numerical differentiation, Amer. Math. Monthly 108 (2001) 512-521.
  * [14] Huang J.G., Chen Y., A regularization method for the function reconstruction from approximate average fluxes, Inverse Problems 21(2005) 1667-1684.
  * [15] Jia X. Z., Wang Y.B. and Cheng J., The numerical differentiation of scattered data and its error estimate. Mathematics, A Journal of Chinese Universities, 25 (2003), 81-90.
  * [16] Lu S., Wang Y.B., The numerical differentiation of first and second order with Tikhonov regularization. Numerical Mathematics, A Journal of Chinese Universities, 26 (2004) 62-74.
  * [17] D.A. Murio, Automatic numerical differentiation by discrete mollification, Comput. Math. Appl. 13 (1987) 381-386.
  * [18] Ramm A. G., Smirnova A. B., On stable numerical differentiation, Math. Comp. 70 (2001) 1131-1153.
  * [19] Scott L.B., Scott L.R., Efficient methods for data smoothing, SIAM Journal on Numerical Analysis 26 (1989), no. 3, 681-692.
  * [20] Wahba G., Smoothing noisy data with spline functions, Numerische Mathematik 24 (1975), no. 5, 383-393.
  * [21] Wei T., Hon Y.C. and Wang Y.B., Reconstruction of numerical derivatives from scattered noisy data. Inverse Problems, 21 (2005) 657-672.
  * [22] Wang Y.B., Jia X.Z. and Cheng J., A numerical differentiation method and its application to reconstruction of discontinuity. Inverse Problems, 18 (2002) 1461-1476.
  * [23] Wang Y.B., Wei T., Numerical differentiation for two-dimensional scattered data. Journal of Mathematical Analysis and Applications, 312 (2005) 121-137.