A new regularization approach for numerical differentiation

Abstract

It is well known that the problem of numerical differentiation is an ill-posed problem and one requires regularization methods to approximate the solution. The commonly practiced regularization methods are (external) parameter-based like Tikhonov regularization, which has certain inherent difficulties associated with them. In such scenarios, iterative regularization methods serve as an attractive alternative. In this paper, we propose a novel iterative regularization method where the minimizing functional does not contain the noisy data directly, but rather a smoothed or integrated version of it. The advantage, in addition to circumventing the use of noisy data directly, is that the sequence of functions constructed during the descent process tends to avoid overfitting, and hence, does not corrupt the recovery significantly. To demonstrate the effectiveness of our method we compare the numerical results obtained from our method with the numerical results obtained from certain standard regularization methods such as Tikhonov regularization, Total-variation, etc.

Keywords:

• Numerical differentiation
• Newton-type methods
• Mathematical programming methods
• optimization and variational techniques
• Volterra integral equations
• inverse problems
• Tikhonov regularization
• iterative regularization
• numerical analysis

1991 Mathematics Subject Classifications:

• Primary 65D25
• Secondary 49M15
• 65K05
• 65K10
• 45D05
• 45Q05

1. Introduction

In many applications, we want to calculate the derivative of a function measured experimentally, that is, to differentiate a function obtained from discrete noisy data. The problem consists of calculating stably the derivative of a smooth function g given its noisy data $\tilde{g}$ such that $||\tilde{g}-g||\le \delta$, where the norm $||\cdot ||$ can be either the ${\mathcal{L}}^{\mathrm{\infty }}\mathrm{or}{\mathcal{L}}^2$-norm. This turns out to be an ill-posed problem, since for a small δ such that $||\tilde{g}-g||\le \delta$ we can have $||{\tilde{g}}^{\prime }-g^{\prime }||$ arbitrarily large or even throw $\tilde{g}$ out of the set of differentiable functions. Many methods and techniques have been introduced in the literature regarding this topic (see [1–13] and references therein). They mostly fall into one of three categories: difference methods, interpolation methods and regularization methods. One can use the first two methods to get satisfactory results provided that the function is given precisely, but they fail miserably when they encounter even small amounts of noise, if the step size is not chosen appropriately or in a regularized way. Hence in such scenarios, regularization methods need to be employed to handle the instability arising from noisy data, with Tikhonov's regularization method being very popular in this respect. In a regularization method, the instability is bypassed by reducing the numerical differentiation problem to a family of well-posed problems, depending on a regularization (smoothing) parameter. Once an optimal value for this parameter is found, the corresponding well-posed problem is then solved to obtain an estimate for the derivative. Unfortunately, the computation of the optimal value for this parameter is a nontrivial task, and usually it demands some prior knowledge of the errors involved. Though there are many techniques developed to find an appropriate parameter value, such as the Morzov's discrepancy principle, the L-curve method and others, still sometimes the solution recovered is either over-fitted or over-smoothed; especially when g has sharp edges, discontinuities or the data has extreme noise in it. An attractive alternative is the iterative regularization methods, like the Landweber iterations. In an iterative regularization method (without any explicit dependence on external parameters), one can start the minimization process first and then use the iteration index as an regularization parameter, i.e. one terminates the minimization process at an appropriate time to achieve regularization. More recently iterative methods have been investigated in the frame work of regularization of nonlinear problems, see [14–18].

In this paper, we propose a new smoothing or regularization technique that doesn't involve any external parameters, thus avoiding all the difficulties associated with them. Furthermore, since this method does not involve the noisy data directly, it is very robust to large noise level in the data as well as to errors with non-zero mean. Thus it makes the method computationally feasible when dealing with real-life data, where we do not expect to have any knowledge of the errors involved or when there is extreme noise involved in the data.

Let us briefly outline our method. First, we lay out a brief summary of a typical (external) parameter-dependent regularization method, and then we present our new (external) parameter-independent smoothing technique. For a given differentiable function g the problem of numerical differentiation, finding $\phi =g^{\prime }$, can be expressed via a Volterra integral equation:¹ (1) $T_D\phi =g_1,$(1) where $g_1(x)=g(x)-g(a)$ and (2) $(T_D\phi )(x):={\int }_a^x\phi (\xi )\mathrm{d}\xi ,$(2) for all $x\in [a,b]\subset \mathbb{R}$. Here one can attempt to solve Equation (1(1) $T_D\phi =g_1,$(1) ) by approximating ϕ with the minimizer of the functional $G_1(\psi )=||T_D\psi -g_1|{|}_2^2.$ However, as stated earlier, the recovery becomes unstable for a noisy g. So, to counter the ill-posed nature of this problem regularization techniques are introduced: one instead minimizes the functional² (3) $G_2(\psi ,\beta )=||T_D\psi -g_1|{|}_2^2+\beta ||D\psi ||,$(3) where β is a regularization parameter and D is a differentiation operator of some order (see [1,14]). This converts the ill-posed problem to a family of conditionally well-posed problems depending on β. The first term (fitting term) of the functional $G_2$ ensures that the inverse solution fits well with the given data, when applied by the forward operator $T_D$, and the second term (smoothing term) in (3(3) $G_2(\psi ,\beta )=||T_D\psi -g_1|{|}_2^2+\beta ||D\psi ||,$(3) ) controls the smoothness of the inverse solution. Since the minimizer of the second term is not the solution of the inverse problem, unless it is a trivial solution, one needs to balance between the smoothing and fitting of the inverse recovery by finding an optimal value ${\beta }_0$. Then the exact solution ϕ is approximated by minimizing the corresponding functional $G(\cdot ,{\beta }_0)$.

Given that an important key in any regularization technique is the conversion of an ill-posed problem to a (conditionally) well-posed problem, we make this our main goal. We start with integrating the data twice, which helps in smoothing out the noise. Thus (1(1) $T_D\phi =g_1,$(1) ) can be reformulated as: find a ϕ that satisfies (4) $T_D\phi =-u_1^{″},$(4) where $-u_1^{″}=g_1$ and $g_1^{\prime }=g^{\prime }=\phi$. Equivalently, $u_1(x)=-(-{\int }_x^b{\int }_a^{\eta }{g}_1(\xi )\mathrm{d}\xi \mathrm{d}\eta )$ or (5) $u_1(x)={\int }_x^b{\int }_a^{\eta }g(\xi )\mathrm{d}\xi \mathrm{d}\eta -g(a)[\frac{(b-a{)}^2}{2}-\frac{(x-a{)}^2}{2}].$(5) Now we find the solution of (4(4) $T_D\phi =-u_1^{″},$(4) ) by approximating it with the minimizer of the following functional: (6) $G(\psi )=||u_{\psi }^{\prime }-u_1^{\prime }|{|}_2^2,$(6) where $u_{\psi }$ is the solution, for a given ψ, of the following boundary value problem: (7) $\begin{array}{rl}-u_{\psi }^{″}& =T_D\psi ,\end{array}$(7) (8) $\begin{array}{rl}{u}_{\psi }(a)& =u_1(a),u_{\psi }(b)=u_1(b).\end{array}$(8) Note that unlike the functionals $G_1$ and $G_2$ which are defined on the data g directly, the functional G is defined rather on the transformed data $u_1^{\prime }$. Thus one can see that when the data g has noise in it then working with the smoothed (integrated) data $u_1$ is more effective than working with the noisy data g directly. It will be proved in later sections that this simple change in the working space from g to $u_1$ drastically improves the stability and smoothness of the inverse recovery of ϕ, without adding any further smoothing terms.

Remark 1.1

Typically, real-life data have zero-mean additive noise in it, i.e. (9) $\tilde{g}(x)=g(x)+ϵ(x),$(9) for $x\in [a,b]$, such that ${\mathbb{E}}_{\mu }(ϵ)=0$ and ${\mathbb{E}}_{\mu }({ϵ}^2)\le \delta$ (or equivalently, $||\tilde{g}-g|{|}_{{\mathcal{L}}^2}\le \delta$), where μ is the probability density function of the random variable ε. Therefore, integrating the data smooths out the noise present in the original data ($\tilde{g}$) and hence, the new data ($\widetilde{u_1}$, as defined by (5(5) $u_1(x)={\int }_x^b{\int }_a^{\eta }g(\xi )\mathrm{d}\xi \mathrm{d}\eta -g(a)[\frac{(b-a{)}^2}{2}-\frac{(x-a{)}^2}{2}].$(5) ) for $\tilde{g}$) for the transformed equation (4(4) $T_D\phi =-u_1^{″},$(4) ) has a significantly reduced noise level in it, see the comparison in Figure 5.

We prove, in Section 3, that the functional G is strictly convex and hence has a unique global minimum which satisfies (1(1) $T_D\phi =g_1,$(1) ), that is, $\phi ={\arg \min }_{\psi }G(\psi )$ satisfies $g^{\prime }=\phi$. The global minimum is achieved via an iterative method using an upgraded steepest gradient or conjugate-gradient method, presented in Section 4, where the use of an (Sobolev) ${\mathcal{H}}^1$-gradient for G instead of the commonly used ${\mathcal{L}}^2$-gradient is discussed and is a crucial step in the optimization process. In Section 5, we shed some light on the stability, convergence and well-posedness of this technique. In Section 6, the Volterra operator $T_D$ is improved for better computation and numerical efficiency but keeping it consistent with the theory developed. In Section 7, we provide some numerical results and compare our method of numerical differentiation with some popular regularization methods, namely Tikhonov regularization, total variation, smoothing spline, mollification method and least square polynomial approximation (the comparison is done with results provided in [1,2,19]). Finally, in Section 8 we present an efficient (heuristic) stopping criterion for the descent process when we don't have any prior knowledge of the error norm.

2. Notations and preliminaries

We adopt the following notations that will be used throughout the paper. All functions are real-valued defined on a bounded closed domain $[a,b]\subset \mathbb{R}$. For $1\le p<\mathrm{\infty }$, ${\mathcal{L}}^p[a,b]$ := $({\mathcal{L}}^p,||\cdot |{|}_{{\mathcal{L}}^p},[a,b])$ denotes the usual Banach space of p-integrable functions on $[a,b]$ and the space ${\mathcal{L}}^{\mathrm{\infty }}[a,b]:=$ $({\mathcal{L}}^{\mathrm{\infty }},\Vert \cdot {\Vert }_{{\mathcal{L}}^{\mathrm{\infty }}},[a,b])$ contains the essentially bounded measurable functions. Likewise the Sobolev space ${\mathcal{H}}^q[a,b]:=$ $({\mathcal{H}}^q,\Vert \cdot {\Vert }_{{\mathcal{H}}^q},[a,b])$ contains all the functions for which f, $f^{\prime },\dots ,f^{(q)}\in {\mathcal{L}}^2[a,b]$, with weak differentiation understood, and the space ${\mathcal{H}}_0^q[a,b]:=\{f\in {\mathcal{H}}^q:‘f\text{ vanishes at the boundary}{}^{\prime }\}$. The spaces ${\mathcal{L}}^2[a,b]$ and ${\mathcal{H}}^q$ are Hilbert spaces with inner-products denoted by $(\cdot ,\cdot {)}_{{\mathcal{L}}^2}$ and $(\cdot ,\cdot {)}_{{\mathcal{H}}^q}$, respectively.

Remark 2.1

Note that although the Volterra operator $T_D$ is well defined on the space of integrable functions (${\mathcal{L}}^1[a,b]$), we will restrict it to the space of ${\mathcal{L}}^2[a,b]$ functions, that is, we shall consider ${\mathcal{L}}^2[a,b]$ to be the searching space for the solution of (1(1) $T_D\phi =g_1,$(1) ), since it is a Hilbert space and has a nice inner product $(\cdot ,\cdot {)}_{{\mathcal{L}}^2}$. Hence the domain of the functional G is ${\mathcal{D}}_G={\mathcal{L}}^2[a,b]$. It is not hard to see that the Volterra operator $T_D$ is linear and bounded on ${\mathcal{L}}^2[a,b]$.

Remark 2.2

From Remark 2.1, since $\phi \in {\mathcal{L}}^2[a,b]$ we have $T_D\phi \in {\mathcal{L}}^2[a,b]$. Thus our problem space to be considered is ${\mathcal{H}}^1[a,b]$, that is, $g\in {\mathcal{H}}^1[a,b]$. Note that since we will not be working with g directly but rather with $u_1$, integrating twice makes $u_1\in {\mathcal{H}}^3[a,b]\subset {\mathcal{H}}^1[a,b]$. This is particularly significant in the sense that we are able to upgrade the smoothness of the working data or information space from ${\mathcal{H}}^1[a,b]$ to ${\mathcal{H}}^3[a,b]$ and hence improve the stability and efficiency of numerical computations. This effect can be seen in some of the numerical examples presented in Section 7 where the noise involved in the data is extreme.

Before we proceed to prove the convexity of G and the well-posedness of the method, we will manipulate $u_1$ further to push it to an even nicer space. From Equation (5(5) $u_1(x)={\int }_x^b{\int }_a^{\eta }g(\xi )\mathrm{d}\xi \mathrm{d}\eta -g(a)[\frac{(b-a{)}^2}{2}-\frac{(x-a{)}^2}{2}].$(5) ), we have $u_1(b)=0\text{and}{u}_1(a)={\int }_a^b{\int }_a^{\eta }g(\xi )\mathrm{d}\xi \mathrm{d}\eta -g(a)\frac{(b-a{)}^2}{2}$ Now redefining the working data as $u(x)=u_1(x)-(b-x/b-a)u_1(a)$, we have (10) $\begin{array}{rl}u(x)& ={\int }_x^b{\int }_a^{\eta }g(\xi )\mathrm{d}\xi \mathrm{d}\eta -g(a)[\frac{(b-a{)}^2}{2}-\frac{(x-a{)}^2}{2}]\\ & -\frac{b-x}{b-a}[{\int }_a^b{\int }_a^{\eta }g(\xi )\mathrm{d}\xi \mathrm{d}\eta -g(a)\frac{(b-a{)}^2}{2}],\end{array}$(10) and this gives $u(a)=0$, $u(b)=0$ and $-u^{″}=-u_1^{″}=g_1$. Thus $u\in {\mathcal{H}}_0^3[a,b]$ and the negative Laplace operator $-\mathrm{\Delta }=-{\mathrm{\partial }}^2/\mathrm{\partial }{x}^2$ is a positive operator in ${\mathcal{L}}^2[a,b]$ on the domain ${\mathcal{D}}_{\mathrm{\Delta }}={\mathcal{H}}_0^2[a,b]\supset {\mathcal{H}}_0^3[a,b]$, since for any $v\in {\mathcal{H}}_0^2[a,b]$ (11) $\begin{array}{rl}(-\mathrm{\Delta }v,v{)}_{{\mathcal{L}}^2}& ={\int }_a^b-(\mathrm{\Delta }v)v\mathrm{d}x\\ & =[-v\mathrm{\nabla }v{]}_a^b+{\int }_a^b|\mathrm{\nabla }v{|}^2\mathrm{d}x\\ & =(\mathrm{\nabla }v,\mathrm{\nabla }v{)}_{{\mathcal{L}}^2}\ge 0\end{array}$(11) and equality holds when $\mathrm{\nabla }v=0$, which implies $v\equiv 0$. Here, and until otherwise specified, the notations Δ and ∇, will mean ${\mathrm{\partial }}^2/\mathrm{\partial }{x}^2$, $\mathrm{\partial }/\mathrm{\partial }x$, respectively, throughout the paper. The functional $G(\psi )$ in (6(6) $G(\psi )=||u_{\psi }^{\prime }-u_1^{\prime }|{|}_2^2,$(6) ) is now redefined, with the new u, as (12) $G(\psi )=\Vert u_{\psi }^{\prime }-u^{\prime }{\Vert }_2.$(12) The positivity of the operator $-\mathrm{\Delta }$ in ${\mathcal{L}}^2[a,b]$ on the domain ${\mathcal{D}}_{\mathrm{\Delta }}$ also helps us to obtain an upper and lower bound for $G(\psi )$, for any $\psi \in {\mathcal{L}}^2[a,b]$. First we get the trivial upper bound (13) $G(\psi )=||u^{\prime }-u_{\psi }^{\prime }|{|}_{{\mathcal{L}}^2}^2\le \sum _{j=0}^1{\left||\frac{{\mathrm{\partial }}^j}{\mathrm{\partial }{x}^j}(u-u_{\psi })|\right|}_{{\mathcal{L}}^2}^2=||u-u_{\psi }|{|}_{{\mathcal{H}}^1}^2$(13) To get a lower bound we note from (8(8) $\begin{array}{rl}{u}_{\psi }(a)& =u_1(a),u_{\psi }(b)=u_1(b).\end{array}$(8) ) that $v=u-u_{\psi }\in {\mathcal{H}}_0^2[a,b]$, thus (14) $\begin{array}{rl}G(\psi )& =||u^{\prime }-u_{\psi }^{\prime }|{|}_{{\mathcal{L}}^2}^2\\ & \ge {\lambda }_1||u-u_{\psi }|{|}_{{\mathcal{L}}^2}^2,\end{array}$(14) where ${\lambda }_1>0$ is the smallest eigenvalue of the positive operator $-\mathrm{\Delta }$ on ${\mathcal{D}}_{\mathrm{\Delta }}$. Hence, we get the bounds of $G(\psi )$ in terms of the ${\mathcal{H}}^1[a,b]$ norm of $u-u_{\psi }$ as (15) ${(1+\frac{1}{{\lambda }_1})}^{-1}||u-u_{\psi }|{|}_{{\mathcal{H}}^1}^2\le G(\psi )=||u^{\prime }-u_{\psi }^{\prime }|{|}_{{\mathcal{L}}^2}^2\le ||u-u_{\psi }|{|}_{{\mathcal{H}}^1}^2.$(15) Equation (15(15) ${(1+\frac{1}{{\lambda }_1})}^{-1}||u-u_{\psi }|{|}_{{\mathcal{H}}^1}^2\le G(\psi )=||u^{\prime }-u_{\psi }^{\prime }|{|}_{{\mathcal{L}}^2}^2\le ||u-u_{\psi }|{|}_{{\mathcal{H}}^1}^2.$(15) ) indicates the stability of the method and the convergence of $u_{{\psi }_m}$ to u in ${\mathcal{H}}^1[a,b]$ if $G({\psi }_m)\to 0$, as is explained in detail in Section 5.

3. Convexity of the functional G

In this section, we prove the convexity of the functional G, together with some important properties.

Theorem 3.1

1. An equivalent form of G, for any $\psi \in {\mathcal{L}}^2[a,b],$ is (16) $\begin{array}{rl}G(\psi )& =||u^{\prime }-u_{\psi }^{\prime }|{|}_{{\mathcal{L}}^2}^2\\ & ={\int }_a^b({u^{\prime }}^2-{u_{\psi }^{\prime }}^2)-2(T_D\psi )(u-u_{\psi })\mathrm{d}x.\end{array}$(16)

2. For any ${\psi }_1,{\psi }_2\in {\mathcal{L}}^2[a,b],$ we have (17) $G({\psi }_1)-G({\psi }_2)={\int }_a^b-2(T_D({\psi }_1-{\psi }_2))(u-\frac{u_{{\psi }_1}+u_{{\psi }_2}}{2})\mathrm{d}x.$(17)

3. The first G$\hat{a}$teaux derivative,³ at $\psi \in {\mathcal{L}}^2[a,b],$ for G is given by (18) $G^{\prime }(\psi )[h]={\int }_a^b(T_{D}h)(-2(u-u_{\psi }))\mathrm{d}x,$(18) where $h\in {\mathcal{L}}^2[a,b]$. And the ${\mathcal{L}}^2$-gradient of G, at ψ, is given by (19) ${\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G=T_D^{\ast }(-2(u-u_{\psi })),$(19) where the adjoint $T_D^{\ast }$ of $T_D$ is given by, for any $f\in {\mathcal{L}}^2[a,b],$⁴ (20) $(T_D^{\ast }f)(x)={\int }_x^{b}f(\xi )\mathrm{d}\xi ;\mathrm{\forall }x\in [a,b].$(20)

4. The second Gâteaux derivative,⁵ at any $\psi \in {\mathcal{L}}^2[a,b],$ of G is given by (21) $G^{″}(\psi )[h,k]=2(-{\mathrm{\Delta }}^{-1}(T_{D}h),(T_{D}k){)}_{{\mathcal{L}}^2}$(21) where $h,k\in {\mathcal{L}}^2[a,b]$. Hence for any $\psi \in {\mathcal{L}}^2[a,b],$ $G^{″}(\psi )$ is a positive definite quadratic form.

For the proof of Theorem 3.1, we also need the following ancillary result.

Lemma 3.2

For fixed $\psi ,h\in {\mathcal{L}}^2[a,b]$ we have, in ${\mathcal{H}}^1[a,b],$ (22) $\underset{ϵ\to 0}{lim}{u}_{\psi +ϵh}=u_{\psi }.$(22)

Proof.

Subtracting the following equations: $\begin{array}{rl}-u_{\psi }^{″}& =T_D\psi \\ -u_{\psi +ϵh}^{″}& =T_D(\psi +ϵh),\end{array}$ we have (23) $-\mathrm{\Delta }(u_{\psi +ϵh}-u_{\psi })=ϵT_{D}h$(23) and using $u_{\psi +ϵh}-u_{\psi }\in {\mathcal{H}}_0^2[a,b]$ we get, via integration by parts, $(\mathrm{\nabla }(u_{\psi +ϵh}-u_{\psi }),\mathrm{\nabla }(u_{\psi +ϵh}-u_{\psi }){)}_{{\mathcal{L}}^2}=ϵ(T_{D}h,u_{\psi +ϵh}-u_{\psi }{)}_{{\mathcal{L}}^2}.$ Now from (14(14) $\begin{array}{rl}G(\psi )& =||u^{\prime }-u_{\psi }^{\prime }|{|}_{{\mathcal{L}}^2}^2\\ & \ge {\lambda }_1||u-u_{\psi }|{|}_{{\mathcal{L}}^2}^2,\end{array}$(14) ) we have ${\lambda }_1||u_{\psi +ϵh}-u_{\psi }|{|}_{{\mathcal{L}}^2}^2\le ||\mathrm{\nabla }(u_{\psi +ϵh}-u_{\psi })|{|}_{{\mathcal{L}}^2}^2$, and using the Cauchy–Schwarz inequality, we have $\begin{array}{rl}& (1+{\lambda }_1^{-1}{)}^{-1}||u_{\psi +ϵh}-u_{\psi }|{|}_{{\mathcal{H}}^1}^2\le ||\mathrm{\nabla }(u_{\psi +ϵh}-u_{\psi })|{|}_{{\mathcal{L}}^2}^2=ϵ(T_{D}h,u_{\psi +ϵh}-u_{\psi }{)}_{{\mathcal{L}}^2},\\ & ⟹(1+{\lambda }_1^{-1}{)}^{-1}\Vert u_{\psi +ϵh}-u_{\psi }{\Vert }_{{\mathcal{H}}^1}\le ϵ||T_{D}h|{|}_{{\mathcal{L}}^2},\end{array}$ where ${\lambda }_1>0$. Hence $u_{\psi +ϵh}\stackrel{ϵ\to 0}{\to }{u}_{\psi }$ in ${\mathcal{H}}^1[a,b]$, since the operator $T_D$ is bounded and $\psi ,h\in {\mathcal{L}}^2[a,b]$ are fixed, which implies the right hand side is of $O(ϵ)$.

3.1. Proof of Theorem 3.1

The proof of first two properties (i) and (ii) are straight forward via integration by parts and using the fact that $u-u_{\psi }\in {\mathcal{H}}_0^2[a,b]$.

In order to prove (iii) and (iv), we use Lemma 3.2.

(iii)

The Gâteaux derivative of the functional G at ψ in the direction of $h\in {\mathcal{L}}^2[a,b]$ is given by (24) $G^{\prime }(\psi )[h]=\underset{ϵ\to 0}{lim}\frac{G(\psi +ϵh)-G(\psi )}{ϵ}.$(24) Now for a fixed $ϵ>0$, we have using (23(23) $-\mathrm{\Delta }(u_{\psi +ϵh}-u_{\psi })=ϵT_{D}h$(23) ), $\begin{array}{rl}& \frac{G(\psi +ϵh)-G(\psi )}{ϵ}\\ & ={ϵ}^{-1}{\int }_a^b(u^{\prime }-u_{\psi +ϵh}^{\prime }{)}^2-(u^{\prime }-u_{\psi }^{\prime }{)}^2\mathrm{d}x\\ & ={ϵ}^{-1}{\int }_a^b(u_{\psi }^{\prime }-u_{\psi +ϵh}^{\prime })(2u^{\prime }-(u_{\psi +ϵh}+u_{\psi }^{\prime }))\mathrm{d}x\\ & ={ϵ}^{-1}{\int }_a^b-\mathrm{\Delta }(u_{\psi }-u_{\psi +ϵh})(2u-(u_{\psi +ϵh}+u_{\psi }))\mathrm{d}x\\ & ={ϵ}^{-1}{\int }_a^b-ϵ(T_{D}h)(2u-(u_{\psi +ϵh}+u_{\psi }))\mathrm{d}x\\ & =-(T_{D}h,2u-(u_{\psi +ϵh}+u_{\psi }){)}_{{\mathcal{L}}^2}.\end{array}$Using Lemma 3.2, one obtains the Gâteaux derivative of G at $\psi \in {\mathcal{L}}^2[a,b]$ in the direction of $h\in {\mathcal{L}}^2[a,b]$ as $G^{\prime }(\psi )[h]=(T_{D}h,-2(u-u_{\psi }){)}_{{\mathcal{L}}^2}.$Note that $(T_{D}h,-2(u-u_{\psi }){)}_{{\mathcal{L}}^2}=(h,T_D^{\ast }(-2(u-u_{\psi })){)}_{{\mathcal{L}}^2}$ for all $h\in {\mathcal{L}}^2[a,b]$, where $T_D^{\ast }$ is the adjoint of the operator $T_D$. Hence by Riesz representation theorem, the ${\mathcal{L}}^2$-gradient of the functional G at ψ is given by ${\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G=T_D^{\ast }(-2(u-u_{\psi })).$

(iv)

Finally, the second Gâteaux derivative for the functional G at ψ is given by (25) $G^{″}(\psi )[h,k]=\underset{ϵ\to 0}{lim}\frac{G^{\prime }(\psi +ϵh)[k]-G^{\prime }(\psi )[k]}{ϵ}$(25) Again for a fixed $ϵ>0$, we have using (23(23) $-\mathrm{\Delta }(u_{\psi +ϵh}-u_{\psi })=ϵT_{D}h$(23) ) $\begin{array}{rl}& \frac{G^{\prime }(\psi +ϵh)[k]-G^{\prime }(\psi )[k]}{ϵ}\\ & ={ϵ}^{-1}{\int }_a^b(T_{D}k)(-2(u-u_{\psi +ϵh}))-(T_{D}k)(-2(u-u_{\psi }))\mathrm{d}x\\ & ={ϵ}^{-1}{\int }_a^b-2(T_{D}k)(u_{\psi }-u_{\psi +ϵh})\mathrm{d}x\\ & ={ϵ}^{-1}{\int }_a^b-2(T_{D}k)(ϵ{\mathrm{\Delta }}^{-1}(T_{D}h))\mathrm{d}x\\ & =2{\int }_a^b(T_{D}k)(-{\mathrm{\Delta }}^{-1}(T_{D}h))\mathrm{d}x\\ & =2(-{\mathrm{\Delta }}^{-1}(T_{D}h),T_{D}k{)}_{{\mathcal{L}}^2}.\end{array}$Hence from (25(25) $G^{″}(\psi )[h,k]=\underset{ϵ\to 0}{lim}\frac{G^{\prime }(\psi +ϵh)[k]-G^{\prime }(\psi )[k]}{ϵ}$(25) ) and letting $ϵ\to 0$ we get $G^{″}(\psi )[h,k]=2(-{\mathrm{\Delta }}^{-1}(T_{D}h),T_{D}k{)}_{{\mathcal{L}}^2}.$Here we can see the strict convexity of the functional G, as for any $h\in {\mathcal{L}}^2[a,b]$, we have $\begin{array}{rl}{G}^{″}(\psi )[h,h]& =2(-{\mathrm{\Delta }}^{-1}(T_{D}h),(T_{D}h){)}_{{\mathcal{L}}^2}\\ & =2(y,-\mathrm{\Delta }y{)}_{{\mathcal{L}}^2},\end{array}$where $-\mathrm{\Delta }y=T_{D}h$ and $y\in {\mathcal{H}}_0^2[a,b]$ (from (23(23) $-\mathrm{\Delta }(u_{\psi +ϵh}-u_{\psi })=ϵT_{D}h$(23) )). As $-\mathrm{\Delta }$ is a positive operator on ${\mathcal{H}}_0^2[a,b]$, y is the trivial solution if and only if $T_{D}h=0$. But ${\int }_a^{x}h(\xi )\mathrm{d}\xi =0$for all $x\in [a,b]$ if and only if $h\equiv 0$. Thus $G^{″}(\psi )$ is a positive definite form for any $\psi \in {\mathcal{L}}^2[a,b]$.

In this section, we proved that the functional G is strictly convex and hence has a unique minimizer, which is attained by the solution ϕ of the inverse problem (1(1) $T_D\phi =g_1,$(1) ). We next discuss a descent algorithm that uses the ${\mathcal{L}}^2$-gradient to derive other gradients that provide descent directions, with better and faster descent rates.

4. The descent algorithm

Here we discuss the problem of minimizing the functional G via a descent method. Theorem 3.1 suggests that the minimization of the functional G should be computationally effective in that ϕ is not only the unique global minimum for G but also the unique zero for the gradient, that is, ${\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G\ne 0$ for $\psi \ne \phi$. Now for a given $\psi \in {\mathcal{L}}^2[a,b]$, let $h\in {\mathcal{L}}^2[a,b]$ denote an update direction for ψ. Then Taylor's expansion gives $G(\psi -\alpha h)=G(\psi )-\alpha G^{\prime }(\psi )[h]+O({\alpha }^2)$ or, for sufficiently small $\alpha >0$, we have (26) $G(\psi -\alpha h)-G(\psi )\approx -\alpha G^{\prime }(\psi )[h].$(26) So if we choose the direction h in such a way that $G^{\prime }(\psi )[h]>0$, then we can minimize G along this direction. Thus we can set up a recovery algorithm for ϕ, forming a sequence of values $G({\psi }_{\mathrm{initial}}-{\alpha }_1{h}_1)>\cdots >G({\psi }_{m-1}-{\alpha }_m{h}_m)>G({\psi }_m-{\alpha }_{m+1}{h}_{m+1})\ge 0$ We list a number of different gradient directions that can make $G^{\prime }(\psi )[h]>0$.

1. The ${\mathcal{L}}^2$-Gradient:

   First, notice from Theorem 3.1 that at a given $\psi \in {\mathcal{L}}^2[a,b]$, (27) $G^{\prime }(\psi )[h]=(h,{\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G{)}_{{\mathcal{L}}^2}$(27) so if we choose the direction $h={\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G$ at ψ, then $G^{\prime }(\psi )[h]>0$. However, there are numerical issues associated with ${\mathcal{L}}^2$-gradient of G during the descent process stemming from the fact that it is always zero at b. Consequently, the boundary data at b for the evolving functions ${\psi }_m$ are invariant during the descent and there is no control on the evolving boundary data at b. This can result in severe decay near the boundary point b if ${\psi }_{\mathrm{initial}}(b)\ne \phi (b)$, as the end point for all such ${\psi }_m$ is glued to ${\psi }_{\mathrm{initial}}(b)$, but ${\psi }_m\to \phi$ in ${\mathcal{L}}^2[a,b]$, as is proved in Section 5.

2. The ${\mathcal{H}}^1$-Gradient:

   One can circumvent this problem by opting for the Sobolev gradient ${\mathrm{\nabla }}_{{\mathcal{H}}^1}^{\psi }G$ instead (see [20]), which is also known as the Neuberger gradient. It is defined as follows: for any $h\in {\mathcal{H}}^1[a,b]$ (28) $\begin{array}{rl}{G}^{\prime }(\psi )[h]& =({\mathrm{\nabla }}_{{\mathcal{H}}^1}^{\psi }G,h{)}_{{\mathcal{H}}^1}\\ & =(g^{\prime },h^{\prime }{)}_{{\mathcal{L}}^2}+(g,h{)}_{{\mathcal{L}}^2}\\ & =-(g^{″},h{)}_{{\mathcal{L}}^2}+(g,h{)}_{{\mathcal{L}}^2}+[g^{\prime }h{]}_a^b\\ & =(-g^{″}+g,h{)}_{{\mathcal{L}}^2}+[g^{\prime }h{]}_a^b\end{array}$(28) where $g={\mathrm{\nabla }}_{{\mathcal{H}}^1}^{\psi }G$. Comparing with (27(27) $G^{\prime }(\psi )[h]=(h,{\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G{)}_{{\mathcal{L}}^2}$(27) ) one can obtain the Neuberger gradient g at ψ, by solving the boundary value problem (29) $\begin{array}{rl}-g^{″}+g& ={\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G\\ [g^{\prime }h{]}_a^b& =0.\end{array}$(29) Setting h = g, the boundary condition becomes (30) $[g^{\prime }g{]}_a^b=g^{\prime }(b)g(b)-g^{\prime }(a)g(a)=0.$(30) This provides us a gradient, ${\mathrm{\nabla }}_{{\mathcal{H}}^1}^{\psi }G$, with considerably more flexibility at the boundary points a and b. In particular consider the following cases:

   1. Dirichlet Neuberger gradient : $g(a)=0$ and $g(b)=0$.

   2. Neumann Neuberger gradient : $g^{\prime }(a)=0$ and $g^{\prime }(b)=0$.

   3. Robin or mixed Neuberger gradient : $g(a)=0$ and $g^{\prime }(b)=0$ or $g^{\prime }(a)=0$ and $g(b)=0$.

   This excellent smoothing technique was originally introduced and used by Neuberger. In addition to the flexibility at the end points, it enables the new gradient to be a preconditioned (smoothed) version of ${\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G$, as $g=(I-\mathrm{\Delta }{)}^{-1}{\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G$, and hence gives a superior convergence in the steepest descent algorithms. So now choosing the descent direction h = g at ψ makes $G^{\prime }(\psi )[h]=(g,{\mathrm{\nabla }}_{{\mathcal{H}}^1}^{\psi }G{)}_{{\mathcal{H}}^1}>0$ and hence $G(\psi -\alpha h)-G(\psi )<0$ (from (26(26) $G(\psi -\alpha h)-G(\psi )\approx -\alpha G^{\prime }(\psi )[h].$(26) )). As stated earlier, the greatest advantage of this gradient is the control of boundary data during the descent process, since based on some prior information of the boundary data we can choose any one of the three aforementioned gradients. For example, if some prior knowledge on $\phi (a)$ and $\phi (b)$ are known, then one can define ${\phi }_{\mathrm{initial}}$ as the straight line joining them and use the Dirichlet Neubeger gradient for the descent. Thus the boundary data is preserved in each of the evolving ${\psi }_m$ during the descent process, which leads to a much more efficient, and faster, descent compared to the normal ${\mathcal{L}}^2$-gradient descent. Even when $\phi {|}_{\{a,b\}}$ is unknown, one can use the Neumann Neuberger gradient that allows free movements at the boundary points rather than gluing it to a fixed value. In the latter scenario, one can even take the average of ${\mathrm{\nabla }}_{{\mathcal{H}}^1}^{\psi }G$ and ${\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G$ to make use of both the gradients.

3. The ${\mathcal{L}}^2-{\mathcal{H}}^1$ Conjugate Gradient:

   If one wishes to further boost the descent speed (by, roughly, a factor of 2) and make the best use of both the gradients, then the standard Polak–Ribi$\stackrel{`}{\mathrm{e}}$re conjugate gradient scheme (see [21]) can be implemented. The initial search direction at ${\psi }_0$, is $h_0=g_0={\mathrm{\nabla }}_{{\mathcal{H}}^1}^{{\psi }_0}G$. At ${\psi }_m$ one can use the exact or inexact line search routine to minimize $G(\psi )$ in the direction of $h_m$ resulting in ${\psi }_{m+1}$. Then $g_{m+1}={\mathrm{\nabla }}_{{\mathcal{H}}^1}^{{\psi }_{m+1}}G$ and $h_{m+1}=g_{m+1}+{\gamma }_m{h}_m$, where (31) ${\gamma }_m=\frac{(g_{m+1}-g_m,g_{m+1}{)}_{{\mathcal{H}}^1}}{(g_m,g_m{)}_{{\mathcal{H}}^1}}=\frac{(g_{m+1}-g_m,{\mathrm{\nabla }}_{{\mathcal{L}}^2}^{{\psi }_{m+1}}G{)}_{{\mathcal{L}}^2}}{(g_m,{\mathrm{\nabla }}_{{\mathcal{L}}^2}^{{\psi }_m}G{)}_{{\mathcal{L}}^2}}.$(31)

Remark 4.1

Though the ${\mathcal{L}}^2$ –${\mathcal{H}}^1$ conjugate gradient boost the descent rate, it (sometimes) compromises the accuracy of the recovery, especially when the noise present in the data is extreme, see Table 2 when $\sigma =0.1$. Now one can also construct a ${\mathcal{H}}^1$ –${\mathcal{H}}^1$ conjugate gradient based only on the Sobolev gradient. This is smoother than the ${\mathcal{L}}^2$ –${\mathcal{H}}^1$ conjugate gradient (as Sobolev gradients are smoother, see Equation (28(28) $\begin{array}{rl}{G}^{\prime }(\psi )[h]& =({\mathrm{\nabla }}_{{\mathcal{H}}^1}^{\psi }G,h{)}_{{\mathcal{H}}^1}\\ & =(g^{\prime },h^{\prime }{)}_{{\mathcal{L}}^2}+(g,h{)}_{{\mathcal{L}}^2}\\ & =-(g^{″},h{)}_{{\mathcal{L}}^2}+(g,h{)}_{{\mathcal{L}}^2}+[g^{\prime }h{]}_a^b\\ & =(-g^{″}+g,h{)}_{{\mathcal{L}}^2}+[g^{\prime }h{]}_a^b\end{array}$(28) )) and hence improves the accuracy of the recovered solution (specially for smooth recovery), but it is tad slower than the ${\mathcal{L}}^2$ –${\mathcal{H}}^1$ conjugate gradient. Finally, the simple Sobolev gradient (${\mathrm{\nabla }}_{{\mathcal{H}}^1}^{}G$) provides the most accurate recovery, but at the cost of the descent rate (it the slowest amongst the three), see Tables 2 and 3.

4.1. The line search method

We minimize the single variable function $f_m(\alpha )=G({\psi }_{m+1}(\alpha ))$, where ${\psi }_{m+1}(\alpha )={\psi }_m-\alpha {\mathrm{\nabla }}_{{\mathcal{H}}^1}^{{\psi }_m}G$, via a line search minimization by first bracketing the minimum and then using some well-known optimization techniques like Brent minimization to further approximate it. Note that the function $f_m(\alpha )$ is strictly decreasing in some neighbourhood of $\alpha =0$ as $f_m^{\prime }(0)=-||g_m|{|}_{{\mathcal{H}}^1}^2<0$.

In order to achieve numerical efficiency, we need to carefully choose the initial step size ${\alpha }_0$. For that, we use the quadratic approximation of the function $f_m(\alpha )$ as follows: (32) $f_m(\alpha )\approx G({\psi }_m)-\alpha G^{\prime }({\psi }_m)[g_m]+\frac{1}{2}{\alpha }^2{G}^{″}({\psi }_m)[g_m,g_m],$(32) which gives the minimizing value for α as (33) ${\alpha }_0=\frac{G^{\prime }({\psi }_m)[g_m]}{G^{″}({\psi }_m)[g_m,g_m]}.$(33) Now since ${\alpha }_0$ is derived from the quadratic approximation of the functional G, it is usually very close to the optimal value, thereby reducing the computational time of the descent algorithm significantly. Now if, for $k=0,1,2\dots$, $f_m((k+1){\alpha }_0)>f_m(k{\alpha }_0)$ then we have a bracket, $[max\{k-1,0\}{\alpha }_0,(k+1){\alpha }_0]$, for the minimum and one can use single variable minimization solvers to approximate it.

In this section, we saw a descent algorithm, with various gradients, where starting from an initial guess ${\psi }_0\in {\mathcal{L}}^2[a,b]$, we obtain a sequence of ${\mathcal{L}}^2$-functions ${\psi }_m$ for which the sequence $\{G({\psi }_m)\}$ is strictly decreasing. In the next section, we discuss the convergence of the ${\psi }_m$'s to ϕ and the stability of the recovery.

5. Convergence, stability and conditional well-posedness

Exact data: We first prove that the sequence of functions constructed during the descent process converges to the exact source function in the absence of any error term and then proves the stability of the process in the presence of noise in the data.

5.1. Convergence

First we see that for the sequence $\{{\psi }_m\}$ produced by the steepest descent algorithm, described in Section 4, we have $G({\psi }_m)\to 0$, since the functional G is non-negative and strictly convex (with the global minimizer ϕ, $G(\phi )=0$) and $G({\psi }_{m+1})<G({\psi }_m)$. In this section, we prove that if for any sequence of functions $\{{\psi }_m\}\subset {\mathcal{L}}^2[a,b]$ such that $G({\psi }_m)\to 0$, then ${\psi }_m\stackrel{w}{\to }\phi$ in ${\mathcal{L}}^2[a,b]$, where $\stackrel{w}{\to }$ denotes weak convergence in ${\mathcal{L}}^2[a,b]$.

Theorem 5.1

Suppose that $\{{\psi }_m\}$ is any sequence of ${\mathcal{L}}^2$-functions such that the sequence $\{G({\psi }_m)\}$ tends to zero. Then $\{{\psi }_m\}$ converges weakly to ϕ in ${\mathcal{L}}^2[a,b]$ and $\{u_{{\psi }_m}\}$ converges strongly to u in ${\mathcal{H}}^1[a,b]$. Also, the sequence $\{g_m\}$ converges weakly to $g_1$ in ${\mathcal{L}}^2[a,b],$ where $g_m=-u_{{\psi }_m}^{″}$ and $g_1=-u^{″}$.

Proof.

The proof of $u_{{\psi }_m}\stackrel{s}{\to }u$ in ${\mathcal{H}}^1[a,b]$ is trivial from the bounds of $G(\psi )$ in Equation (15(15) ${(1+\frac{1}{{\lambda }_1})}^{-1}||u-u_{\psi }|{|}_{{\mathcal{H}}^1}^2\le G(\psi )=||u^{\prime }-u_{\psi }^{\prime }|{|}_{{\mathcal{L}}^2}^2\le ||u-u_{\psi }|{|}_{{\mathcal{H}}^1}^2.$(15) ), which gives (34) $||u-u_{{\psi }_m}|{|}_{{\mathcal{H}}^1}^2\le (1+\frac{1}{{\lambda }_1})G({\psi }_m).$(34) To see the weak convergence of $\{{\psi }_m\}$ to ϕ in ${\mathcal{L}}^2[a,b]$, we first prove that the sequence $\{T_D{\psi }_m\}$ converges weakly to $T_D\phi$ in ${\mathcal{L}}^2[a,b]$. Since $\{u_{{\psi }_m}\}$ converges strongly to u in ${\mathcal{H}}^1[a,b]$, this implies $\{u_{{\psi }_m}^{\prime }\}$ and $\{u_{{\psi }_m}\}$ converge weakly to $u^{\prime }$ and u in ${\mathcal{L}}^2[a,b]$, respectively. Now we will use the fact that ${\mathcal{C}}_0^{\mathrm{\infty }}[a,b]$ is dense in ${\mathcal{L}}^2[a,b]$ (in ${\mathcal{L}}^2$-norm), i.e. for any $\psi \in {\mathcal{L}}^2[a,b]$ and $ϵ>0$ there exists a $\tilde{\psi }\in {\mathcal{C}}_0^{\mathrm{\infty }}[a,b]$ such that $||\psi -\tilde{\psi }|{|}_{{\mathcal{L}}^2}\le ϵ$. So for any $\tilde{\psi }\in {\mathcal{C}}_0^{\mathrm{\infty }}[a,b]$ we have, using ${\tilde{\psi }}^{\prime }\in {\mathcal{L}}^2[a,b]$, (35) $\begin{array}{rl}(T_D{\psi }_m-T_D\phi ,\tilde{\psi }{)}_{{\mathcal{L}}^2}& =(-\mathrm{\Delta }(u_{{\psi }_m}-u),\tilde{\psi }{)}_{{\mathcal{L}}^2}\\ & =(\mathrm{\nabla }(u_{{\psi }_m}-u),{\tilde{\psi }}^{\prime }{)}_{{\mathcal{L}}^2},\end{array}$(35) which tends to zero as $m\to \mathrm{\infty }$. Hence by the density of ${\mathcal{C}}_0^{\mathrm{\infty }}[a,b]$ in ${\mathcal{L}}^2[a,b]$, it can be proved that the sequence $\{T_D{\psi }_m\}$ converges weakly to $T_D\phi$ in ${\mathcal{L}}^2[a,b]$.

To prove convergence of $\{{\psi }_m\}$ to $\{\phi \}$ weakly in ${\mathcal{L}}^2[a,b]$, we can use $(T_D({\psi }_m-\phi ),\psi {)}_{{\mathcal{L}}^2}=({\psi }_m-\phi ,T_D^{\ast }\psi {)}_{{\mathcal{L}}^2}$ and (35(35) $\begin{array}{rl}(T_D{\psi }_m-T_D\phi ,\tilde{\psi }{)}_{{\mathcal{L}}^2}& =(-\mathrm{\Delta }(u_{{\psi }_m}-u),\tilde{\psi }{)}_{{\mathcal{L}}^2}\\ & =(\mathrm{\nabla }(u_{{\psi }_m}-u),{\tilde{\psi }}^{\prime }{)}_{{\mathcal{L}}^2},\end{array}$(35) ). Therefore, our proof will be complete if we can show that the range of $T_D^{\ast }$ is dense in ${\mathcal{L}}^2[a,b]$ in ${\mathcal{L}}^2$-norm. Again we start with any $\tilde{\psi }\in {\mathcal{C}}_0^{\mathrm{\infty }}[a,b]$, and using ${\tilde{\psi }}^{\prime }\in {\mathcal{L}}^2[a,b]$, we have $\begin{array}{rl}(T_D^{\ast }(-{\tilde{\psi }}^{\prime }))(x)& =-{\int }_x^b{\tilde{\psi }}^{\prime }(\eta )\mathrm{d}\eta \\ & =\tilde{\psi }(x),\end{array}$ i.e. ${\mathcal{C}}_0^{\mathrm{\infty }}[a,b]\subset Range(T_D^{\ast })$. Hence we have $Range(T_D^{\ast })$ dense in ${\mathcal{L}}^2[a,b]$ in ${\mathcal{L}}^2$-norm. For convergence of $\{g_m\}$ to $g_1$, note that $T_D{\psi }_m=-u_{{\psi }_m}^{″}=g_m$ and $T_D\phi =-u^{″}=g_1$. And since $\{T_D{\psi }_m\}$ converges weakly to $T_D\phi$ in ${\mathcal{L}}^2[a,b]$ implies $\{g_m\}$ converges weakly to $g_1$ in ${\mathcal{L}}^2[a,b]$.

Remark 5.2

It can be further proved that the sequence $\{{\psi }_m\}$ converges strongly to ϕ in ${\mathcal{L}}^2[a,b]$. To prove this, first, one needs to analyse the operator associated with the minimizing functional G as defined in (12(12) $G(\psi )=\Vert u_{\psi }^{\prime }-u^{\prime }{\Vert }_2.$(12) ), i.e. from the definition of $u_{\psi }$ in (7(7) $\begin{array}{rl}-u_{\psi }^{″}& =T_D\psi ,\end{array}$(7) ) together with the criterion $u_{\psi }\in {\mathcal{H}}_0^2[a,b]$ we have (from (10(10) $\begin{array}{rl}u(x)& ={\int }_x^b{\int }_a^{\eta }g(\xi )\mathrm{d}\xi \mathrm{d}\eta -g(a)[\frac{(b-a{)}^2}{2}-\frac{(x-a{)}^2}{2}]\\ & -\frac{b-x}{b-a}[{\int }_a^b{\int }_a^{\eta }g(\xi )\mathrm{d}\xi \mathrm{d}\eta -g(a)\frac{(b-a{)}^2}{2}],\end{array}$(10) ), with $g=T_D\psi$ and $(T_D\psi )(a)=0$, by the definition of $T_D$) (36) $u_{\psi }^{\prime }=-{\int }_a^x(T_D\psi )(\xi )\mathrm{d}\xi +\frac{1}{b-a}{\int }_a^b{\int }_a^{\eta }(T_D\psi )(\xi )\mathrm{d}\xi \mathrm{d}\eta .$(36) From the expression (36(36) $u_{\psi }^{\prime }=-{\int }_a^x(T_D\psi )(\xi )\mathrm{d}\xi +\frac{1}{b-a}{\int }_a^b{\int }_a^{\eta }(T_D\psi )(\xi )\mathrm{d}\xi \mathrm{d}\eta .$(36) ), the operator $L(\psi ):=u_{\psi }^{\prime }$ is both linear and bounded, and hence minimizing the functional G in (12(12) $G(\psi )=\Vert u_{\psi }^{\prime }-u^{\prime }{\Vert }_2.$(12) ) is equivalent to Landweber iterations corresponding to the operator L. Then from the convergence theories developed for Landweber iterations the sequence $\{{\psi }_m\}$ converges to ϕ strongly in ${\mathcal{L}}^2[a,b]$, for details on iterative regularization see [14,16].

Theorem 5.1 proves that for the given function $g_1\in {\mathcal{L}}^2[a,b]$ we are able to construct a sequence of smooth functions, $\{g_m\}\subset {\mathcal{L}}^2[a,b]$, that converges (weakly) to $g_1$ in ${\mathcal{L}}^2[a,b]$. This is critical since, when the data has noise in it one needs to terminate the descent process at an appropriate instance to attain regularization, see Section 5.3, and the (weak) convergence helps us to construct such a stopping criterion in the absence of noise information (δ), see Section 8.

Noisy data: In this section, we consider the data has noise in it and shows that the sequence of functions constructed during the descent process using the noisy data still approximates the exact solution, under some conditions.

5.2. Stability

Here we prove the stability of the process. We will prove this by considering the problem of numerical differentiation as equivalent to finding a unique minimizer of the positive functional G, this makes the problem well-posed. That is, for a given $g\in {\mathcal{H}}^1[a,b]$, and hence a given $u\in {\mathcal{H}}_0^3[a,b]$ and the functional G, the problem of finding (derivative) $\phi \in {\mathcal{L}}^2[a,b]$ such that $g^{\prime }=\phi$ is equivalent to finding the minimizer of the functional G, i.e. a $\psi \in {\mathcal{L}}^2[a,b]$ such that $G(\psi )\le ϵ$, for any small $ϵ>0$, is a conditionally well-posed problem. It is not hard to prove that if two functions g, $\tilde{g}\in {\mathcal{H}}^1[a,b]$ are such that $||g-\tilde{g}|{|}_{{\mathcal{L}}^2}\le \delta$, where $\delta >0$ is small, then the corresponding u, $\tilde{u}\in {\mathcal{H}}_0^3[a,b]$ also satisfy⁶ $||u-\tilde{u}|{|}_{{\mathcal{H}}^1}\le C\delta$, for some constant C.

Theorem 5.3

Suppose the function $\tilde{u}$ in the perturbed version of the function u such that $||u-\tilde{u}|{|}_{{\mathcal{H}}^1}\le \delta ,$ where $\delta >0,$ and let ϕ, $\tilde{\phi }\in {\mathcal{L}}^2[a,b]$ denote their respective recovered functions, such that $T_D\phi =-u^{″}$ and $T_D\tilde{\phi }=-{\tilde{u}}^{″}$. Let the functional G, without loss of generality, be defined based on $\tilde{u},$ that is, $G(\tilde{\phi })=0,$ then we have (37) $0\le G(\phi )\le C{\delta }^2,$(37) where C is some constant.

Proof.

Since $u_{\phi }=u$ and ${\tilde{u}}_{\tilde{\phi }}=\tilde{u}$, the proof follows from the definitions of the corresponding functionals, as $G(\phi )=||u_{\phi }^{\prime }-{\tilde{u}}^{\prime }|{|}_{{\mathcal{L}}^2}^2=||u^{\prime }-{\tilde{u}}^{\prime }|{|}_{{\mathcal{L}}^2}^2\le C{\delta }^2$

In the next theorem, we prove that if a sequence of functions $\{{\psi }_m\}$ converges to $\tilde{\phi }$ in ${\mathcal{L}}^2[a,b]$, then it also approximates ϕ in ${\mathcal{L}}^2[a,b]$, that is, if $G_{\tilde{u}}({\psi }_m)$ is small, then $G_u({\psi }_m)$ is also small where $G_u$ and $G_{\tilde{u}}$ are the functionals formed based on u and $\tilde{u}$, respectively.

Theorem 5.4

Suppose for a sequence of functions $\{{\psi }_m\}\subset {\mathcal{L}}^2[a,b],$ the corresponding sequence $\{G_{\tilde{u}}({\psi }_m)\}\subset \mathbb{R}$ converges to zero, where the functional $G_{\tilde{u}}$ is formed based on $\tilde{u},$ then for the original u such that $||u-\tilde{u}|{|}_{{\mathcal{H}}^1}\le \delta ,$ for small $\delta >0,$ there exists a $M(\delta )\in \mathbb{N}$ such that for all $m\ge M(\delta ),$ $G_u({\psi }_m)\le c{\delta }^2$ for some constant c and $G_u$ is the functional based on u.

Proof.

For u, $\tilde{u}$ and any ${\psi }_m\in {\mathcal{L}}^2[a,b]$, we have $G_u({\psi }_m)=||u^{\prime }-u_{{\psi }_m}^{\prime }|{|}_{{\mathcal{L}}^2}^2$ and $G_{\tilde{u}}({\psi }_m)=||{\tilde{u}}^{\prime }-u_{{\psi }_m}^{\prime }|{|}_{{\mathcal{L}}^2}^2$, then $\begin{array}{rl}{G}_u({\psi }_m)& =||u^{\prime }-u_{{\psi }_m}^{\prime }|{|}_{{\mathcal{L}}^2}^2\\ & \le ||{\tilde{u}}^{\prime }-u^{\prime }|{|}_{{\mathcal{L}}^2}^2+||{\tilde{u}}^{\prime }-u_{{\psi }_m}^{\prime }|{|}_{{\mathcal{L}}^2}^2\end{array}$ and hence from Theorem 5.1, the result follows.

5.3. Conditional well-posedness (iterative-regularization)

As explained earlier, in an external-parameter based regularization method (like Tikhonov-type regularizations) first, one converts the ill-posed problem to a family of well-posed problem (depending on the parameter value β) and then, only after finding an appropriate regularization parameter value (say ${\beta }_0$), one proceeds to the recovery process, i.e. completely minimize the corresponding functional $G(\cdot ,{\beta }_0)$ (as defined in (3(3) $G_2(\psi ,\beta )=||T_D\psi -g_1|{|}_2^2+\beta ||D\psi ||,$(3) )). Whereas in a classical iterative regularization method (not involving any external-parameters like β), one cannot recover a regularized solution by simply minimizing a related functional completely, instead stopping the recovery process at an appropriate instance provides the regularizing effect to the solution. That is, one starts to minimize some related functional (to recover the solution) but then terminates the minimization process at an appropriate iteration before it has reached the minimum (to restrict the influence of the noise), i.e. here the iteration stopping index serves as a regularization parameter, for details see [14,16].

In this section, we further explain the above phenomenon by showing that if one attempts to recover the true (or original) solution ϕ by using a noisy data then it will distort the recovery. First we see that for an exact g (or equivalently, an exact u) and the functional constructed based on it (i.e. $G_u$) we have the true solution (ϕ) satisfying $G_u(\phi )=0$. However, for a given noisy $\tilde{g}$, with $||\tilde{g}-g|{|}_{{\mathcal{L}}^2}=\delta >0$, and the functional based on it (i.e. $G_{\tilde{u}}$) we will have $G_{\tilde{u}}(\phi )>0$, see Theorem 5.5. So if we construct a sequence of functions ${\psi }_m^{\delta }\in {\mathcal{L}}^2[a,b]$, using the descent algorithm and based on the noisy data $\tilde{g}$, such that $G_{\tilde{u}}({\psi }_m^{\delta })\to 0$ then (from Theorem 5.1) we will have ${\psi }_m^{\delta }\to \tilde{\phi }$, where $\tilde{\phi }$ is the recovered noisy solution satisfying $G_{\tilde{u}}(\tilde{\phi })=0$. This implies initially ${\psi }_m^{\delta }\to \phi$ and then upon further iterations ${\psi }_m^{\delta }$ diverges away from ϕ and approaches $\tilde{\phi }$. Hence, the errors in the recoveries $||{\psi }_m^{\delta }-\phi |{|}_{{\mathcal{L}}^2}$ follow a semi-convergence nature, i.e. decreases first and then increases. This is a typical behaviour of any ill-posed problem and is managed, as stated above, by stopping the descent process at an appropriate iteration $M(\delta )$ such that $G_{\tilde{u}}({\psi }_{M(\delta )}^{\delta })>0$ but close to zero (due to the stability Theorems 5.3 and 5.4). Following similar arguments as in (15(15) ${(1+\frac{1}{{\lambda }_1})}^{-1}||u-u_{\psi }|{|}_{{\mathcal{H}}^1}^2\le G(\psi )=||u^{\prime }-u_{\psi }^{\prime }|{|}_{{\mathcal{L}}^2}^2\le ||u-u_{\psi }|{|}_{{\mathcal{H}}^1}^2.$(15) ), we can have a lower bound for $G_{\tilde{u}}(\phi )$.

Theorem 5.5

Given two functions u, $\tilde{u}\in {\mathcal{H}}_0^3[a,b],$ their respective recovery ϕ, $\tilde{\phi }\in {\mathcal{L}}^2[a,b],$ such that $T_D\phi =-u^{″}$ and $T_D\tilde{\phi }=-{\tilde{u}}^{″},$ and let the functional $G_{\tilde{u}}$ be defined based on $\tilde{u},$ that is, $G_{\tilde{u}}(\tilde{\phi })=0,$ then we have (38a) $G_{\tilde{u}}(\phi )=||u^{\prime }-{\tilde{u}}^{\prime }|{|}_{{\mathcal{L}}^2}^2\ge {\lambda }_1||u-\tilde{u}|{|}_{{\mathcal{L}}^2}^2,$(38a) as an ${\mathcal{L}}^2$-lower bound and for a ${\mathcal{H}}^1$-lower bound, we have (38b) $G_{\tilde{u}}(\phi )=||u^{\prime }-{\tilde{u}}^{\prime }|{|}_{{\mathcal{L}}^2}^2\ge {(1+\frac{1}{{\lambda }_1})}^{-1}||u-\tilde{u}|{|}_{{\mathcal{H}}^1}^2$(38b) where ${\lambda }_1={\pi }^2/(b-a{)}^2$ is the smallest eigenvalue of $-\mathrm{\Delta }$ on ${\mathcal{H}}_0^2[a,b]$.

Therefore, combining Theorems 5.3 and 5.5 we have the following two-sided inequality for $G_{\tilde{u}}(\phi )$, for some constants $C_1$ and $C_2$, $0\le C_1||u-\tilde{u}|{|}_{{\mathcal{H}}^1}^2\le G_{\tilde{u}}(\phi )\le C_2||u-\tilde{u}|{|}_{{\mathcal{H}}^1}^2.$ Thus, when $\delta \to 0$ we have $G_{\tilde{u}}(\phi )\to 0$ which implies ${\phi }_{\delta }\to \phi$ in ${\mathcal{L}}^2[a,b]$. Now though we would like to use the bounds in Theorem 5.5 to terminate the descent process, but we do not known the exact g (or equivalently, the exact u), and hence cannot use that as the stopping condition. However, if the error norm $\delta =||g-g_{\delta }|{|}_{{\mathcal{L}}^2}$ is known then one can use Morozov's discrepancy principle, [22], as a stopping criterion for the iteration process, that is, terminate the iteration when (39) $||T{\psi }_m-g_{\delta }|{|}_{{\mathcal{L}}^2}\le \tau \delta$(39) for an appropriate $\tau >1$,⁷ and for unknown δ one usually goes for heuristic approaches to stop the iterations, an example of which is presented in Section 8.

6. Numerical implementation

In this section, we provide an algorithm to compute the derivative numerically. Notice that though one can use the integral operator equation (1(1) $T_D\phi =g_1,$(1) ) to recover ϕ inversely, for computational efficiency we can further improve the operator and the operator equation, but keeping the theory intact. First we see that the adjoint operator $T_D^{\ast }$ can also provide an integral operator equation of interest, (40) $-T_D^{\ast }\phi =g_2,$(40) where $g_2(x)=g(x)-g(b)$ and $(T_D^{\ast }\phi )(x):={\int }_x^b\phi (\xi )\mathrm{d}\xi ,$ for all $x\in [a,b]\subset \mathbb{R}$. Now we can combine both the operator equations (1(1) $T_D\phi =g_1,$(1) ) and (40(40) $-T_D^{\ast }\phi =g_2,$(40) ) to get the following integral operator equation: (41) $T\phi =g_3,$(41) where $g_3(x)=2g(x)-(g(a)+g(b))$ and for all $x\in [a,b]\subset \mathbb{R}$, (42) $\begin{array}{rl}(T\phi )(x)& :=(T_D\phi -T_D^{\ast }\phi )(x)\\ & ={\int }_a^x\phi (\xi )\mathrm{d}\xi -{\int }_x^b\phi (\xi )\mathrm{d}\xi .\end{array}$(42) The advantage of the operator equation (41(41) $T\phi =g_3,$(41) ) over (1(1) $T_D\phi =g_1,$(1) ) or (40(40) $-T_D^{\ast }\phi =g_2,$(40) ) is that it recovers ϕ symmetrically at the end points. For example if we consider the operator equation (40(40) $-T_D^{\ast }\phi =g_2,$(40) ) for the recovery of ϕ, then during the descent process the ${\mathcal{L}}^2$-gradient at ψ will be ${\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G=-T_D(-2(u-u_{\psi }))$ (similar to (19(19) ${\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G=T_D^{\ast }(-2(u-u_{\psi })),$(19) )) which implies ${\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G(a)=0$ for all $\psi \in {\mathcal{L}}^2[a,b]$. Hence the boundary data at a for all the evolving ${\psi }_m$'s are going to be invariant during the descent.⁸ As for (1(1) $T_D\phi =g_1,$(1) ), since ${\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G(b)=0$, the boundary data at b for all the evolving ${\psi }_m$'s are going to be invariant during the descent. Even though one can opt for the Sobolev gradient of G at ψ, ${\mathrm{\nabla }}_{{\mathcal{H}}^1}^{\psi }G$, to counter that problem but, due to the intrinsic decay of the base function ${\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G$ for all ${\psi }_m$'s at a or b, the recovery of ϕ near that respective boundary will not be as good (or symmetric) as at the other end. On the other hand if we use the operator equation (41(41) $T\phi =g_3,$(41) ) for the descent recovery of ϕ then the ${\mathcal{L}}^2$-gradient at ψ is going to be (43) ${\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G=T^{\ast }(-2(u-u_{\psi }))$(43) where (44) $T^{\ast }=T_D^{\ast }-T_D\mathrm{or}{T}^{\ast }=-T.$(44) Thus ${\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G(a)\ne 0$ and ${\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G(b)\ne 0$, and hence the recovery of ϕ at both the end points will be performed symmetrically. Now one can derive other gradients, like the Neuberger or conjugate gradient based on this ${\mathcal{L}}^2$-gradient, for the recovery of ϕ depending on the scenarios, that is, based on the prior knowledge of the boundary information (as explained in Section 4).

Corresponding to the operator equation (41(41) $T\phi =g_3,$(41) ), the smooth or integrated data $u\in {\mathcal{H}}_0^3[a,b]$ will be (45) $\begin{array}{rl}u(x)& =2{\int }_x^b{\int }_a^{\eta }g(\xi )\mathrm{d}\xi \mathrm{d}\eta -(g(a)+g(b))[\frac{(b-a{)}^2}{2}-\frac{(x-a{)}^2}{2}]\\ & -\frac{b-x}{b-a}[2{\int }_a^b{\int }_a^{\eta }g(\xi )\mathrm{d}\xi \mathrm{d}\eta -(g(a)+g(b))\frac{(b-a{)}^2}{2}]\end{array}$(45) Thus our problem set up now is as follows: for a given $g\in {\mathcal{H}}^1[a,b]$ (and hence a given $u\in {\mathcal{H}}_0^3[a,b]$) we want to find a $\phi \in {\mathcal{L}}^2[a,b]$ such that (46) $T\phi =-u^{″}.$(46) Our inverse approach to achieve ϕ will be to minimize the functional G which is defined, for any $\psi \in {\mathcal{L}}^2[a,b]$, by (47) $G(\psi )=||u^{\prime }-u_{\psi }^{\prime }|{|}_{{\mathcal{L}}^2},$(47) where u is as defined in (45(45) $\begin{array}{rl}u(x)& =2{\int }_x^b{\int }_a^{\eta }g(\xi )\mathrm{d}\xi \mathrm{d}\eta -(g(a)+g(b))[\frac{(b-a{)}^2}{2}-\frac{(x-a{)}^2}{2}]\\ & -\frac{b-x}{b-a}[2{\int }_a^b{\int }_a^{\eta }g(\xi )\mathrm{d}\xi \mathrm{d}\eta -(g(a)+g(b))\frac{(b-a{)}^2}{2}]\end{array}$(45) ) and $u_{\psi }$ is the solution of the boundary value problem (48) $\begin{array}{rl}T\psi & =-u_{\psi }^{″},\\ u_{\psi }(a)& =u(a)\mathrm{and}{u}_{\psi }(b)=u(b).\end{array}$(48) Since our new problem set up is almost identical to the old one, the previous theorems and results developed for $T_D$ can be similarly extended to T. Next we provide a pseudo-code, Algorithm 1, for the descent algorithm described earlier.

Remark 6.1

If prior knowledge of $\phi (a)$ and $\phi (b)$ is known then ${\psi }_0$ can be defined as a straight line joining them and we use Dirichlet Neuberger gradient for the descent. If no prior information is known about $\phi (a)$ and $\phi (b)$, then we simply choose ${\psi }_0\equiv 0$ and use the Sobolev gradient or conjugate gradient for the descent. It has been numerically seen that having any information of $\phi (a)$ or $\phi (b)$ and using it, together with appropriate gradient, significantly improves the convergence rate of the descent process and the efficiency of the recovery. In the examples presented here we have not assumed any prior knowledge of $\phi (a)$ or $\phi (b)$ to keep the problem settings as pragmatic as possible.

Remark 6.2

To solve the boundary value problem (29(29) $\begin{array}{rl}-g^{″}+g& ={\mathrm{\nabla }}_{{\mathcal{L}}^2}^{\psi }G\\ [g^{\prime }h{]}_a^b& =0.\end{array}$(29) ) while calculating the Neuberger gradients we used the invariant embedding technique for better numerical results (see [23]). This is very important as this technique enables us to convert the boundary value problem to a system of initial and final value problems and hence one can use the more robust initial value solvers, compared to boundary value solvers, which normally use shooting methods.

Remark 6.3

For all the numerical testings presented in Section 7, we assumed to have prior knowledge on the error norm (δ) and used it as a stopping criteria, as explained in Section 5.3, for the descent process. We also compare the results obtained without any noise information (i.e. using heuristic stopping strategy, see Section 8) with the results obtained using the noise information, see Tables 2 and 3 in Section 8.

7. Results

A MATLAB program was written to test the numerical viability of the method. We take an evenly spaced grid with $h={10}^{-2}$ in all the examples, unless otherwise specified. In all the examples, we used the discrepancy principle (see (39(39) $||T{\psi }_m-g_{\delta }|{|}_{{\mathcal{L}}^2}\le \tau \delta$(39) )) to terminate the iterations when the discrepancy error goes below δ (which is assumed to be known). In Section 8, we discuss the stopping criterion when δ is unknown.

Example 7.1

Comparison with standard regularization methods

In this example, we compare the inverse recovery using our technique with some of the standard regularization methods. We again perturbed the smooth function $g(x)=\cos (x)$, here we consider the regularity of the data as ${\mathcal{H}}^1$ (i.e. $g\in {\mathcal{H}}^{(k+1)}[a,b]$ for k = 0), on $[-.5,.5]$ by random noises to get $\tilde{g}(x)=g(x)+ϵ(x)$, where ε is a normal random variable with mean 0 and standard deviation σ. Like in [1], we generated two data sets, one with $h={10}^{-2}$ (the dense set) and other with $h={10}^{-1}$ (the sparse set). We tested with $\sigma =0.01$ on both the data sets and $\sigma =0.1$ only on the dense set. We compare the relative errors, $||\phi -\tilde{\phi }|{|}_{{\mathcal{L}}^2}/||\phi |{|}_{{\mathcal{L}}^2}$, obtained in our method (using Neumann ${\mathrm{\nabla }}_{{\mathcal{H}}^1}^{}G$-gradient) with the relative errors provided in [1], which are listed in Table 1. Here we can see that our method of numerical differentiation outperforms most of the other methods, in both the dense and sparse situation. Though Tikhonov method performs better for k = 2, that is when g is assumed to be in ${\mathcal{H}}^3[a,b]$, but for the same smoothness consideration, $g\in {\mathcal{H}}^1[a,b]$ or k = 0, it fails miserably. In fact, one can prove that the ill-posed problem of numerical differentiation turns out to be well-posed when $g\in {\mathcal{H}}^{(k+1)}([a,b])$ for k = 1, 2, see [14], which explains the small relative errors in Tikhonov regularization for k = 1, 2. As stated above, the results in Table 1 are obtained using the Sobolev gradient, whereas Tables 2 and 3 show a comparison in the recovery errors using different gradients and different stopping criteria, as well as the descent rates associated with them. To compare with the total variation method, which is very effective in recovering discontinuities in the solution, we perform a test on a sharp-edged function similar to the one presented in [1] and the results are shown in Example 7.4. Hence we can consider this method as an universal approach in every scenarios.

Table 1. Relative errors in derivative approximation.

Methods	$m=100(h=0.01)$, $\sigma =0.01$	$m=100(h=0.01)$, $\sigma =0.1$	$m=10(h=0.1)$, $\sigma =0.01$
Degree-2 polynomial	0.0287	0.3190	0.2786
Tikhonov, k = 0	0.7393	0.8297	0.7062
Tikhonov, k = 1	0.1803	0.3038	0.6420
Tikhonov, k = 2	0.0186	0.0301	0.4432
Cubic Spline	0.1060	1.15	0.3004
Convolution smoothing	0.1059	0.8603	0.2098
Variational method	0.1669	0.7149	0.3419
Our method (k = 0)	0.0607	0.0839	0.1355

Table 2. Relative errors and the descent rates using different gradients.

$m=100(h=0.01)$, extreme error $\sigma =0.1$	Sobolev gradient (${\mathrm{\nabla }}_{{\mathcal{H}}^1}^{}G$)	${\mathcal{L}}^2$–${\mathcal{H}}^1$ conjugate gradient	${\mathcal{H}}^1$–${\mathcal{H}}^1$ conjugate gradient
Recovery errors (using SC-I Section 5.3, known δ)	0.0839	0.4364	0.1522
Recovery errors (using SC-II Section 8, unknown δ)	0.1299	0.1299	0.1299
# iterations (using SC-I Section 5.3, known δ)	39	4	6
# iterations (using SC-II Section 8, unknown δ)	2	2	2

Table 3. Relative errors and the descent rates using different gradients.

$m=100(h=0.01)$, moderate error $\sigma =0.01$	Sobolev gradient (${\mathrm{\nabla }}_{{\mathcal{H}}^1}^{}G$)	${\mathcal{L}}^2$–${\mathcal{H}}^1$ conjugate gradient	${\mathcal{H}}^1$–${\mathcal{H}}^1$ conjugate gradient
Recovery errors (using SC-I Section 5.3, known δ)	0.0607	0.0589	0.0613
Recovery errors (using SC-II Section 8, unknown δ)	0.1129	0.1129	0.1129
# iterations (using SC-I Section 5.3, known δ)	84	5	18
# iterations (using SC-II Section 8, unknown δ)	3	3	3

In the next example, we show that one does not have to assume the normality conditions for the error term, i.e. the assumption that the noise involved should be iid normal random variables is not needed, which is critical for certain other methods, rather it can be a mixture of any random variables.

Example 7.2

Noise as a mixture of random variables

In the previous example, we saw that our method holds out in the presence of large error norm. Here we show that this technique is very effective even in the presence of extreme noise. We perturb the function $g(x)=\sin (x/3)$ on $[0,3\pi ]$ to $\tilde{g}(x)=g(x)+ϵ(x)$,⁹ where ε is the error function obtained from a mixture of uniform($-\delta$, δ) and normal(0, δ) random variables, where $\delta =0.5$. Figure 1 shows the noisy $\tilde{g}$ and the exact g, and Figure 2(a) shows the computed derivative $\tilde{\phi }$ vs. $\phi (x)=\cos (x/3)/3$. The relative error for the recovery of $\tilde{\phi }$ is 0.0071.

Figure 1. Noisy $\tilde{g}$.

Figure 2. Inverse recovery of the derivative $\tilde{\phi }$ and $T\tilde{\phi }$: (a) derivative $\tilde{\phi }$ vs. ϕ and (b) recovery of $\tilde{\phi }$.

In the next example, we further pushed the limits by not having a zero-mean error term, which is again crucial for many other methods.

Example 7.3

Error with non-zero mean

In this example, we will show that this method is impressive even when the noise involved has nonzero mean. We consider the settings of the previous example: $g(x)=\sin (x/3)$ on $[0,3\pi ]$ is perturbed to $\tilde{g}(x)=g(x)+ϵ(x)$ but here the error function ε is a mixture of uniform(-0.8δ, 1.2δ) and normal(0.1, δ), for $\delta =0.1$. Figure 2(b) shows the recovery of the derivative ϕ versus the true derivative. The relative error of the recovery for ϕ is around 0.0719.

In the following two examples, we provide the results of numerical differentiation done on a piece-wise differentiable functions and compare it with the results obtained in [1,2].

Example 7.4

Discontinuous source function

Here we selected a function randomly from the many functions tested in [2]. The selected function has the following definition: $y_2(t)=\{\begin{array}{ll}1-t,& t\in [0,0.5],\\ t,& t\in (0.5,1].\end{array}$ The function $y_2$ is piecewise differentiable except at the point t = 0.5, where it has a sharp edge. The function is then perturbed by a uniform($-\delta$, δ) random variable to get the noisy data $y_{2_{\delta }}$, shown in Figure 4a, where we even increased the error norm in our testing from $\delta =0.001$ in [2] to $\delta =$0.01 in our case. Figure 3 shows the recoveries using the method described here and comparing it with the result obtained in [2], one can see that our recovery outperforms, especially at the boundary points, the recovery in [2]. We also compare it with a similar result obtained in [1]¹⁰ using a total variation regularization method, shown in Figure 4(b).

Figure 3. Recoveries using out method: (a) numerical derivative $\tilde{\phi }$ and (b) smooth approximation $T\tilde{\phi }$.

Figure 4. (a) Noisy data for Example 7.4 and (b) total variation regularization from [1].

Figure 5. Integration smooths out the noise present in the data: (a) noisy $\tilde{g}$ vs exact g and (b) noisy $\tilde{u}$ vs exact u.

8. Stopping criterion II

As explained in Section 5.3, in the presence noise, one has to terminate the descent process at an appropriate iteration to achieve regularization. The discrepancy principle [24–26] provides a stopping condition provided the error norm (δ) is known. However, in many of the practical situations it is very hard to determine an estimate of the error norm. In such cases, heuristic approaches are taken to determine stopping criteria, such as the L-curve method [27,28]. In this section, we present a new heuristic approach to terminate the iterations when the error norm (δ) is unknown. First we notice that the minimizing functional G used here, as defined in (12(12) $G(\psi )=\Vert u_{\psi }^{\prime }-u^{\prime }{\Vert }_2.$(12) ), does not contain the noisy $\tilde{g}$ directly, rather an integrated (smoothed) version of it ($\tilde{u}$), as compared to a minimizing functional (such as $G_2$, defined in (3(3) $G_2(\psi ,\beta )=||T_D\psi -g_1|{|}_2^2+\beta ||D\psi ||,$(3) )) used in any standard regularization method. Hence, in addition to avoiding the noisy data from affecting the recovery, the integration process also helps in constructing a stopping strategy, which is explained below. Figure 5 shows the difference in g and $\tilde{g}$ vs. u and $\tilde{u}$, from Example 7.2. We can see, from Figure 5(b), that the integration smooths out the noise present in $\tilde{g}$ to get $||\tilde{u}-u|{|}_{{\mathcal{L}}^2}/||u|{|}_{{\mathcal{L}}^2}\approx 0.78\mathrm{\%}$, whereas the noise level in $\tilde{g}$ is $||\tilde{g}-g|{|}_{{\mathcal{L}}^2}/||g|{|}_{{\mathcal{L}}^2}\approx 55.44\mathrm{\%}$. Consequently, the sequence $\{{\tilde{g}}_m:=T{\psi }_m{\}}_{m\ge 1}$, constructed during the descent process, converges (weakly) in ${\mathcal{L}}^2[a,b]$ to $\tilde{g}$, rather than strongly to $\tilde{g}$, that is, for any $\varphi \in {\mathcal{L}}^2[a,b]$ the sequence $\{({\tilde{g}}_m-\tilde{g},\varphi {)}_{{\mathcal{L}}^2}\}$ converges to zero. In other words, the integration mitigates the effects of the high oscillations originating from the random variable and also of any outliers (as its support is close to zero measure). Also, since the forward operator T (as defined in (42(42) $\begin{array}{rl}(T\phi )(x)& :=(T_D\phi -T_D^{\ast }\phi )(x)\\ & ={\int }_a^x\phi (\xi )\mathrm{d}\xi -{\int }_x^b\phi (\xi )\mathrm{d}\xi .\end{array}$(42) )) is smooth, the sequence ${\tilde{g}}_m$ first approximate the exact g (as it is also smooth, $g\in {\mathcal{H}}^1[a,b]$), with the corresponding sequence $\{u_{{\psi }_m}\}$ approximating $\tilde{u}\approx u$, and then the sequence $\{{\tilde{g}}_m\}$ attempts to fit the noisy $\tilde{g}$, which leads to a phenomenon known as overfitting. However, when ${\tilde{g}}_m$ tries to overfit the data (i.e. fit $\tilde{g}$) the sequence values $||\tilde{u}-u_{{\psi }_m}|{|}_{{\mathcal{L}}^2}$ increases, since the ovefitting occurs in a smooth fashion (as T is a smooth operator) and, as a result, increases the integral values. This effect can be seen in Figure 6(a), $||u_{{\psi }_m}-u_{\delta }|{|}_{{\mathcal{L}}^2}$ descent for Example 7.1 (when $\sigma =0.1$), and in Figure 6(b), $||T{\psi }_m-g_{\delta }|{|}_{{\mathcal{L}}^2}$ descent for Example 7.2. One can capture the recoveries at these fluctuating¹¹ points (of either $||T{\psi }_m-g_{\delta }|{|}_{{\mathcal{L}}^2}$, $||u_{{\psi }_m}^{\prime }-u_{\delta }^{\prime }|{|}_{{\mathcal{L}}^2}$ or $||u_{{\psi }_m}-u_{\delta }|{|}_{{\mathcal{L}}^2}$) and choose the recovery corresponding to the earliest iteration for which $T{\psi }_m$ fits through $g_{\delta }$. Choosing the early fluctuating iteration is especially important when dealing with data with large error level, such as in Example 7.1 ($\sigma =0.1$) and Example 7.2. For example, from Figure 6(b) if one captures the recovery at iteration 4 then the relative error in the recovery is only 8% (see Figure 7 a). However, even if an appropriate early iteration is not selected, still the recovery errors saturate after certain iterations, rather than blowing up. This is significant when dealing with data having small to moderate error level, such as in Example 7.1 ($\sigma =0.01$), where one can notice (in Figure 7b) that the relative errors of the recoveries attain saturation after recovering the optimal solution, since $||u-u_{\delta }|{|}_{{\mathcal{L}}^2}\approx 0$ for small δ. Tables 2 and 3 show the relative errors of the recoveries obtained using this heuristic stopping criterion.

Figure 6. Fluctuations during the descent process. (a) $||u_{{\psi }_m}-u_{\delta }|{|}_{{\mathcal{L}}^2}$, Example 7.1 ($\sigma =0.1$) and (b) $||\tilde{g}-g_m|{|}_{{\mathcal{L}}^2}$, Example 7.2.

Figure 7. Relative errors in the recoveries during the descent process: (a) Example 7.2 and (b) Example 7.1 ($\sigma =0.01$).

Remark 8.1

Note that this phenomena does not occur in Landweber iterations, since one minimizes the functional containing the noisy data $g_{\delta }$ directly, i.e. $G(\psi )=||T\psi -g_{\delta }|{|}_{{\mathcal{L}}^2}^2$. Figure 8(a) shows the descent of $||T{\psi }_m-g_{\delta }|{|}_{{\mathcal{L}}^2}$, when Landweber iterations are implemented for Example 7.1 ($\sigma =0.01$) and Figure 8(b) shows the corresponding descent of the relative errors of the recovered solutions, where can see no fluctuations in Figure 8(a) but the semi-convergence in Figure 8(b). Therefore, without any prior knowledge of δ it is hard to stop the descent process and avoids the ill-posedness. Whereas, notice the saturation of the relative errors of the recovery (Figure 7(b)) when we implement our method to the same problem.

Figure 8. Landweber iterations on Example 7.1 ($\sigma =0.01$): (a) $||T{\psi }_m-g_{\delta }|{|}_{{\mathcal{L}}^2}$-descent and (b) relative errors descent.

9. Conclusion and future research

This algorithm for numerical differentiation is very effective, even in the presence of extreme noise, as can be seen from the examples presented in Section 7. Furthermore, it serves as a universal method to deal with all scenarios such as when the data set is dense or sparse and when the function g is smooth or not smooth. The key feature in this technique is that we are able to upgrade the working space of the problem from ${\mathcal{H}}^1[a,b]$ to ${\mathcal{H}}_0^3[a,b]$, which is a much smoother space. Additionally, this method also enjoys many advantages of not encountering the involvement of an external regularization parameter, for example one does not have to determine the optimum parameter choice to balance between the fitting and smoothing of the inverse recovery. Even the heuristic approach for the stopping criteria also provides us with a much better recovery, and hence it's very applicable in the absence of the error norm.

In a follow-up paper, we improve this method to calculate derivatives of functions in higher dimensions and for higher order derivatives. Moreover, we can extend this method to encompass any linear inverse problems and thereby generalize the theory , which will be presented in the coming paper where we will apply this method to recover solution of Fredholm Integral Equations, like deconvolution or general Volterra equation.

Acknowledgments

I am very grateful to Prof. Ian Knowles for his support, encouragement and stimulating discussions throughout the preparation of this paper.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Notes

1 The technicality as to how far g can be weakened so that the solution of equation (1(1) $T_D\phi =g_1,$(1) ) is uniquely and stably recovered (finding the ϕ inversely) will be discussed later.

2 Again the domain space for the functionals $G_1$, $G_2$ and G are discussed in details later.

3 It can be further proved that it's also the first Fr$\acute{e}$chet derivative of G at ψ.

4 Hence $T_D^{\ast }$ is also a linear and bounded operator in ${\mathcal{L}}^2[a,b]$.

5 Again, it can be proved that it's the second Fr$\acute{e}$chet derivative of G at ψ.

6 Just for simplicity, we assume $g(a)=\tilde{g}(a)$.

7 In our experiments, we considered $\tau =1$ and the termination condition as $||T{\psi }_m-g_{\delta }|{|}_{{\mathcal{L}}^2}<\delta$.

8 As explained in the ${\mathcal{L}}^2$-gradient version of the descent algorithm for G, in Section 4.

9 To be consistent with ${}^6$, we kept $\tilde{g}(a)=g(a)$ and $\tilde{g}(b)=g(b)$, but can be avoided.

10 Where the test function is $g(x)=|x-0.5|$ on $[0,1]$ and the data set is 100 uniformly distributed points with $\sigma =0.01$.

11 the fluctuating occurs since the values of $||\tilde{g}-g_m|{|}_{{\mathcal{L}}^2}$ tends to decrease first (when approximating the exact g) and then increases (when making a transition from g to $\tilde{g}$) and eventually decreases (when trying to fit the noisy $\tilde{g}$, i.e. overfitting)