Smoothing $\mathcal{L}^2$ gradients in iterative regularization

TL;DR

This paper introduces a new iterative regularization method that smooths the gradient to improve robustness and effectiveness in solving large-scale inverse problems with noisy data, demonstrated through image deblurring and tomography applications.

Contribution

It proposes a novel iterative regularization technique using smoothed descent directions, enhancing stability and performance over traditional gradient-based methods.

Findings

More robust recovery in high-noise scenarios

Improved results compared to Tikhonov, TV, and CGLS methods

Effective in large-scale inverse problems

Abstract

Connected with the rise of interest in inverse problems is the development and analysis of regularization methods, which are a necessity due to the ill-posedness of inverse problems. Tikhonov-type regularization methods are very popular in this regard. However, its direct implementation for large-scale linear or non-linear problems is a non-trivial task. In such scenarios, iterative regularization methods usually serve as a better alternative. In this paper we propose a new iterative regularization method which uses descent directions, different from the usual gradient direction, that enable a more smoother and effective recovery than the later. This is achieved by transforming the original noisy gradient, via a smoothing operator, to a smoother gradient, which is more robust to the noise present in the data. It is also shown that this technique is very beneficial when dealing with data having large noise level. To illustrate the computational efficiency of this method we apply it to numerically solve some classical integral inverse problems, including image deblurring and tomography problems, and compare the results with certain standard regularization methods, such as Tikhonov, TV, CGLS, etc.

Tables

Table 1. Table 1 . Relative errors for Example 6.1 .

Stopping strategy: Morozov’s discrepancy principle, see (5.8)
	$\phi \left(b\right)=0$			$\phi \left(b\right)\ne 0$
Methods	$\Delta x$	$3\Delta x$	$5\Delta x$	$\Delta x$	$3\Delta x$	$5\Delta x$
Tik ($L=I$)	0.1179			0.2748
Tik ($L=|\nabla |$)	0.1822			0.2205
TSVD	0.0695			0.2734
TTLS	0.4544			0.3521
CGLS	0.0814			0.2543
LSQRB	0.0814			0.2543
$\nu$-method	0.1734			0.2651
G (${\nabla }_{{ℒ}^2-{ℒ}^2}$)	0.0803	0.0720	0.0661	0.2574	0.2691	0.2921
G (${\nabla }_{{ℒ}^2-{\mathscr{H}}^1}$)	0.1608	0.1907	0.1861	0.1474	0.2237	0.2218

Table 2. Table 2 . Imaging inverse problems: Example 6.2 and Example 6.3 .

Relative errors in the recoveries using different regularization methods.
	Deblurring		Tomography
Methods	Relative error	Sparsity (# 0s)	Relative error	Sparsity (# 0s)
IRell1	0.3126	0
IRhtv	0.3361	5732	0.5404	225
IRcgls	0.2353	0	0.3436	0
IRfista	0.2349	4058	0.3012	7505
IRmrnsd	0.2257	1	0.2462	1
IRenrich	0.2352	0	0.3424	0
IRirn	0.3199	2035
IRhy._fgmres	0.3126	0
IRhy._gmres	0.2444	0
IRhy._lsqr	0.2793	0	1.025	0
Our (${\nabla }_{{ℒ}^2-{ℒ}^2}$)	0.2357	3797	0.2964	7563

Equations

$$T\phi =g,$$

$$\text{mbox}minimizeF(\psi ):=||T\psi -g|{|}_2^2.$$

$$F(\psi ;\lambda ,L,{\psi }_0,p,q)=||T\psi -g_{\delta }|{|}_p^p+\lambda ||L(\psi -{\psi }_0)|{|}_q^q,$$

$${\psi }_{m+1}^{\delta }={\psi }_m^{\delta }+{\tau }_m{\xi }_m^{\delta },0\le m\le m_0(\delta )<\infty ,$$

$${\psi }_{m+1}^{\delta }=\sum _{i=0}^m{\mu }_{m,i}{\psi }_i^{\delta }+{\tau }_m{\xi }_m^{\delta },$$

$$\sum _{i=0}^m{\mu }_{m,i}=1,{\tau }_m\ne 0,$$

$${\mathcal{K}}_{m+1}(T^{\ast }T,{\xi }_0^{\delta }):=\text{mbox}span\{{\xi }_0^{\delta },(T^{\ast }T){\xi }_0^{\delta },\cdots ,(T^{\ast }T{)}^m{\xi }_0^{\delta }\}.$$

$${\xi }_m^{\delta }=T^{\ast }(I-{\Delta }^{-1})(g_{\delta }-T{\psi }_m^{\delta }),$$

$$-y^{\prime \prime }=f,$$

$$y(a)=0,y(b)=0.$$

$$G(\psi )=||T\psi -g_{\delta }|{|}_{{\mathcal{L}}^2}^2+||u_{\psi }^{\prime }-u_{\delta }^{\prime }|{|}_{{\mathcal{L}}^2}^2,$$

$$u_{\delta }(x):=(-{\Delta }^{-1}{g}_{\delta })(x)={\int }_x^b{\int }_a^{\eta }{g}_{\delta }(\xi )d\xi d\eta -\frac{b-x}{b-a}{\int }_a^b{\int }_a^{\eta }{g}_{\delta }(\xi )d\xi d\eta .$$

$${\left(1+\frac{1}{{\lambda }_1}\right)}^{-1}||u-u_{\psi }|{|}_{{\mathcal{H}}^2}^2\le G(\psi )\le ||u-u_{\psi }|{|}_{{\mathcal{H}}^2}^2,$$

$$G_1(\psi )$$

$$G_2(\psi )$$

$$G_2(\psi )={\int }_a^b({u^{\prime }}^2-{u_{\psi }^{\prime }}^2)-2(T\psi )(u-u_{\psi })dx$$

$$G_2({\psi }_1)-G_2({\psi }_2)={\int }_a^b-2(T({\psi }_1-{\psi }_2))(u-\frac{u_{{\psi }_1}+u_{{\psi }_2}}{2})dx$$

$$G_2^{\prime }(\psi )[h]={\int }_a^b(Th)(-2(u-u_{\psi }))dx$$

$${\nabla }_{{\mathcal{L}}^2}^{\psi }{G}_2=T^{\ast }(-2(u-u_{\psi })),$$

$$G_2^{\prime \prime }(\psi )[h,k]=2(-{\Delta }^{-1}(Th),(Tk){)}_{{\mathcal{L}}^2}$$

$$\underset{ϵ\to 0}{\lim }{u}_{\psi +ϵh}=u_{\psi },$$

$$-u_{\psi }^{\prime \prime }$$

$$-u_{\psi +ϵh}^{\prime \prime }$$

$$-(u_{\psi +ϵh}-u_{\psi }{)}^{\prime \prime }=ϵTh$$

$$((u_{\psi +ϵh}-u_{\psi }{)}^{\prime },(u_{\psi +ϵh}-u_{\psi }{)}^{\prime }{)}_{{\mathcal{L}}^2}=ϵ(Th,u_{\psi +ϵh}-u_{\psi }{)}_{{\mathcal{L}}^2}.$$

$$(1+{\lambda }_1^{-1}{)}^{-1}||u_{\psi +ϵh}-u_{\psi }|{|}_{{\mathcal{H}}^1}^2\le ||(u_{\psi +ϵh}-u_{\psi }{)}^{\prime }|{|}_{{\mathcal{L}}^2}^2=ϵ(Th,u_{\psi +ϵh}-u_{\psi }{)}_{{\mathcal{L}}^2},$$

$$⟹(1+{\lambda }_1^{-1}{)}^{-1}||u_{\psi +ϵh}-u_{\psi }|{|}_{{\mathcal{H}}^1}\le ϵ||Th|{|}_{{\mathcal{L}}^2},$$

$$G_2^{\prime }(\psi )[h]=\underset{ϵ\to 0}{\lim }\frac{G_2(\psi +ϵh)-G_2(\psi )}{ϵ}.$$

$$\frac{G_2(\psi +ϵh)-G_2(\psi )}{ϵ}$$

$$={ϵ}^{-1}{\int }_a^b(u^{\prime }-u_{\psi +ϵh}^{\prime }{)}^2-(u^{\prime }-u_{\psi }^{\prime }{)}^{2}dx$$

$$={ϵ}^{-1}{\int }_a^b(u_{\psi }^{\prime }-u_{\psi +ϵh}^{\prime })(2u^{\prime }-(u_{\psi +ϵh}+u_{\psi }^{\prime }))dx$$

$$={ϵ}^{-1}{\int }_a^b-(u_{\psi }-u_{\psi +ϵh}{)}^{\prime \prime }(2u-(u_{\psi +ϵh}+u_{\psi }))dx$$

$$={ϵ}^{-1}{\int }_a^b-ϵ(Th)(2u-(u_{\psi +ϵh}+u_{\psi }))dx$$

$$=-(Th,2u-(u_{\psi +ϵh}+u_{\psi }){)}_{{\mathcal{L}}^2}.$$

$$G_2^{\prime }(\psi )[h]=(Th,-2(u-u_{\psi }){)}_{{\mathcal{L}}^2}.$$

Full text

1. Introduction

1.1. Inverse Problems and Regularization:

An inverse problem in general is a problem where the output (effext) is known but the input (source) is not, in contrast to a direct problem where we deduce the effect from the source. Mathematically, an inverse problem is often expressed as the problem of finding/estimating a $\varphi$, given $g$, which satisfies the following operator equation:

[TABLE]

where $T$ is also a known operator describing the underlying physical process, the domain and range of $T$ varies depending on the problem. Typically, the solution of (1.1) is approximated by the solution of a least-square problem, i.e.

[TABLE]

However, inverse problems are usually ill-posed, in the sense of violating Hadamard’s third condition: “Continuous dependence of the data”, i.e., even for a slightly perturbed data $g_{\delta}$, such that $||{g-g_{\delta}}||\leq\delta$ (usually small), the inverse recovery becomes unstable, $||{\varphi-\varphi_{\delta}}||>>\delta$ (very large), due to the unboundedness of the (pseudo-) inverse operator $T^{\dagger}$ and the noise present in the data. To counter such instabilities or the ill-posedness of inverse problems, regularization methods have to be employed. In the last few decades, several regularization methods have been established for linear as well as nonlinear inverse problems. Broadly, there exist two kinds of regularization approaches:

1.2. Tikhonov-type regularization:

Tikhonov-type regularization methods are probably the most well known regularization techniques for solving linear as well as nonlinear inverse problems (see [1, 2, 3, 4, 5]), where (instead of minimizing the simple least-square problem (1.2)) one recovers a regularized solution by minimizing a (constrained or) penalized functional

[TABLE]

(for some $p$ and $q$) where $\lambda>0$ is a called the regularization parameter, $||g-g_{\delta}||\leq\delta$ is the error norm, $\psi_{0}$ is an initial guess and $L$ is a regularization operator, with the null spaces of $T$ and $L$ intersecting trivially. The choice of an appropriate parameter value ($\lambda_{0}$) is crucial here, as it balances between the data fitting term $||{T\psi-g_{\delta}}||_{p}^{p}$ and the regularization term $||L(\psi-\psi_{0})||_{q}^{q}$. If $\lambda$ is too small then minimizing (1.3) over-fits the noisy data (thus, leading to a noisy recovery or, statistically speaking, a high-variance recovery), where as, if $\lambda$ is too big then it under-fits the data (thus, leading to an over-smooth recovery, or statistically, a high-bias recovery). Choosing an appropriate $\lambda_{0}$ is not a trivial task, especially when dealing with large-scale problems, and typically one has to compute the solutions for many different $\lambda$ values before choosing an appropriate $\lambda_{0}$.

1.3. (Semi-) Iterative regularization:

Contrary to the above approach, in a (semi-) iterative regularization method one recovers a regularized solution $\varphi^{\delta}$ of (1.1) by simply stopping the minimization process of the least-square problem (1.2) at an appropriate (early) instance. That is, starting from an initial guess $\psi_{0}^{\delta}$, one updates the $m^{th}$ recovered solution ($\psi_{m}^{\delta}$) in a direction ($\xi_{m}^{\delta}$) such that the $(m+1)^{th}$ recovered solution ($\psi_{m+1}^{\delta}$) is better than the previous, meaning $||\psi_{m+1}^{\delta}-\varphi||<||\psi_{m}^{\delta}-\varphi||$, and stops the iteration after certain number of steps (i.e., $m\leq m_{0}$, for some appropriate $m_{0}$, usually depending on the noise level $\delta$). Mathematically,

[TABLE]

where $\tau_{m}$ is the m-th step-length and the descent direction $\xi_{m}^{\delta}$ is usually the negative gradient of the least-square functional $F$, i.e., $\xi_{m}^{\delta}=-T^{*}(T\psi_{m}^{\delta}-g_{\delta})$. When the step-size is fixed ($0<\tau_{m}=\tau<\frac{2}{||T^{*}T||}$) for all $m$, it’s known as Landweber iterations (also known as Richardson iterations) and has been intensively investigated in the literature (see, [6, 7, 1, 2, 8, 9]). The main drawback of Landweber iterations is its slow performance, i.e., it takes a large number of iterations to obtain the optimal convergence rates. To circumvent this drawback many extensions (known as polynomial or accelerated Landweber methods) have been proposed and studied in the framework of regularization (for an overview, see [8, 9, 10, 11]). Such iterative methods are typically known as semi-iterative methods, since when dealing with a noisy data ($g_{\delta}$) one encounters a semi-convergent nature of the recovery process. The advantage of these extended semi-iterative methods over the simple Landweber iteration is that, while Landweber iteration (1.4) uses only the last iterate $\psi_{m}^{\delta}$ to construct the new approximate $\psi_{m+1}^{\delta}$, in a semi-iterative method one make use of the last few iterates (if not all),

[TABLE]

where $\xi_{m}^{\delta}$ is the descent direction (usually $\xi_{m}^{\delta}=-T^{*}(T\psi_{m}^{\delta}-g_{\delta})$). Note that, $\psi_{m+1}^{\delta}-\psi_{0}^{\delta}$ belongs to the Krylov-subspace $\mathcal{K}_{m+1}(T^{*}T,\xi_{0}^{\delta})$, which is defined as

[TABLE]

Hence, these methods are also called as Krylov-subspace methods. Even further acceleration is possible by adapting a conjugate-gradient type method (see [12]), where $\mu_{m,i}$ in (1.5) depend on data $g_{\delta}$ (making it a non-linear method).

1.4. A new iterative regularization method:

First note that, starting from the simple Landweber iterations (1.4), all the generalizations and extensions (1.5) developed in the iterative regularization literature is focused only on improving the speed of the convergence of the descent process. In this paper, we discuss a new iterative regularization method that not only improves the descent rate but also improves the smoothness of the recovery, especially when dealing with data having large noise level, and hence, leads to a much effective and efficient recovery. Note that, in all of the aforementioned iterative methods there is a direct influence of the noisy data $g_{\delta}$ in the descent direction, i.e., $\xi_{m}^{\delta}=T^{*}(g_{\delta}-T\psi_{m}^{\delta})$. Therefore, when the noise level $\delta$ is large then the recovered solution $\varphi^{\delta}$ (though regularized) sometimes still possesses some of the noisy characteristics arising from $g_{\delta}$, especially when $\delta$ is large and $T^{*}$ does not compensate its noisy influence, see Example 6.1. Here, we present a technique to pre-condition the gradient, using a smoothing operator, so that it significantly reduces the noisy influence of $g_{\delta}$, and thus, enhance the smoothness and accuracy of the recovery. First we observe that, a typical real-life data is usually contaminated by additive noise of zero mean, i.e., $g_{\delta}=g+\epsilon_{\delta}$, where $\epsilon_{\delta}$ is a random variable such that $\mathbb{E}(\epsilon_{\delta})=0$ and $||\epsilon_{\delta}||_{2}^{2}=||g_{\delta}-g||_{2}^{2}\leq\delta^{2}$. Also, notice that integrating the noisy data smooths out the noise present in the data, for example, see Figure 1(a) ($g_{\delta}$ vs. $g$, with a relative error $\frac{||g-g_{\delta}||_{2}}{||g||_{2}}\%\approx 10\%$) vs. Figure 1(b) ($\int g_{\delta}$ vs. $\int g$, with $\frac{||\int(g-g_{\delta})||_{2}}{||\int g||_{2}}\%\approx 0.73\%$). This provides a heuristic motivation to incorporate the integrated data (in a regularized manner) during the minimization process to improve the smoothness of the recovery, the theoretical justifications and implementations are demonstrated in the later sections.

We briefly outline our new method here. Let’s assume the operator $T:\mathcal{D}_{T}\rightarrow\mathcal{H}_{2}$ as a linear bounded injective111for simplicity, we work in an unique solution scenario. For a non-injective $T$ one can look for the minimal norm solution or the generalized solution (pseudo-inverse solution). operator, where $\mathcal{D}_{T}\subset\mathcal{H}_{1}$ and $g\in R(T)$ is a known function defined on $[a,b]\subset\mathbb{R}$, a bounded set. Here $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ are Hilbert spaces, that we consider to be subsets of $\mathcal{L}^{2}[a,b]$. Then an approximate solution of problem (1.1) is attained in an iterative fashion using the following descent direction

[TABLE]

where the integral operator $-\Delta^{-1}$ is defined as follows: for any $f\in\mathcal{L}^{2}[a,b]$, $y=-\Delta^{-1}(f)$ is the solution of the following boundary value problem

[TABLE]

Similar to (1.2), we show that the Variational formulation, related to the gradient defined in (1.7), corresponds to the minimization of the following functional

[TABLE]

where $u_{\psi}:=-\Delta^{-1}(T\psi)$ and $u_{\delta}:=-\Delta^{-1}(g_{\delta})$, and the operator $-\Delta^{-1}$ is as defined above in (1.8). Equivalently,

[TABLE]

For the exact $g$, the first term in equation (1.9) is minimized by the solution ($\varphi$) of the inverse problem and in §3 we prove that the second term also has the same minimizer, i.e., the solution of (1.1) is the minimizer of the functional $G$.

Remark 1.1*.*

One may suspect that using an additional integral operator (in the second term of (1.9)) might increase the ill-posedness of the problem. But, as we are not using the noisy data $g_{\delta}$ directly for that term, rather a smoothed version of it (i.e., $u_{\delta}$) such that $||u_{\delta}-u||\approx 0$, the ill-posedness of the problem is not deteriorated. In addition, the first part (involving the original operator $T$ and data $g_{\delta}$) together with a stopping criterion depending on $g_{\delta}$ also influence the recovery of a regularized solution ($\varphi^{\delta}$). In fact, we will show that with proper integration techniques and using both the terms one can attain even lower relative errors during the descent process than using only the first term and can also significantly reduce the semi-convergent nature of the relative errors, see Figures 3(b), 4(g) and 4(h).

Remark 1.2*.*

Note that, since (in practice) we only deal with discrete data sets, to solve an inverse problem in higher dimension ($n\geq 2$) one can convert the n-dimensional array to an one dimensional vector and apply the techniques developed for one dimensional problems. Hence, in this paper we develop the theories only for 1-dimensional problems and apply it solve n-dimensional problems ($n\geq 2$), like image deblurring and tomography, see Examples 6.2 and 6.3.

Here, we follow an improved descent algorithm for the inverse recovery of $\varphi^{\delta}$, where we first start with the normal $\mathcal{L}^{2}$-gradient of $G$ and then upgrade it to various other gradients, including the $\mathcal{H}^{1}$-gradient, to enhance the descent rate and efficiency of the recovery. Using an appropriate gradient is very crucial in the optimization process as it helps to retrieve the features of $\varphi$ more accurately, for example, the $\mathcal{H}^{1}$-gradient not only smooths the noisy $\mathcal{L}^{2}$-gradient but also helps in pre-conditioning certain desired boundary effects, depending on some prior information of the boundary data, (see Example 6.1), details discussed in §4.

In §3 we prove the convexity (and some other properties) of the functional $G$. In §4 we provide a descent algorithm to minimize the functional $G$ by using different gradients (or descent directions) depending on the scenarios, which is crucial for the minimization process. The convergence of the sequence of functions constructed during the descent process to a regularized solution $\varphi^{\delta}$ and its stability are discussed in §5. To validate the numerical viability of the developed theory we perform numerous computational experiments, provided in §6, on some classical inverse problems, like Fredholm type and Volterra type integral equations, and compare our results with the results obtained using certain standard regularization methods, like Tikhonov-type methods, CGLS, LSQR, TV and more, see Tables 1 and 2.

2. Notations and Preliminaries

We adopt the following notations that are used throughout the paper. All functions are real-valued defined on a bounded closed domain $[a,b]\subset\mathbb{R}$. For $1\leq p<\infty$, $\mathcal{L}^{p}[a,b]$ := $(\mathcal{L}^{p},\lvert\lvert.\rvert\rvert_{\mathcal{L}^{p}},[a,b])$ denotes the usual Banach space of $\lvert f\rvert^{p}$ integrable functions on $[a,b]$ and the space $\mathcal{L}^{\infty}[a,b]:=$ $(\mathcal{L}^{\infty},||.||_{\mathcal{L}^{\infty}},[a,b])$ contains the essentially bounded measurable functions. Likewise the Sobolev space $\mathcal{H}^{q}[a,b]:=$ $(\mathcal{H}^{q},||.||_{\mathcal{H}^{q}},[a,b])$ contains all the functions for which $f,\;,f^{\prime},\;,\cdots,\;f^{(q)}\;\in$ $\mathcal{L}^{2}[a,b]$ and the space $\mathcal{H}^{q}_{0}[a,b]:=\{f\in\mathcal{H}^{q}:``f\mbox{ vanishes at the boundary"}\}$. The spaces $\mathcal{L}^{2}[a,b]$ and $\mathcal{H}^{q}$ are Hilbert spaces with inner-products denoted as $(.,.)_{\mathcal{L}^{2}}$ and $(.,.)_{\mathcal{H}^{q}}$, respectively.

Remark 2.1*.*

Note that, integrating the data ($g\in\mathcal{L}^{2}[a,b]$) twice results in $u\in\mathcal{H}^{2}[a,b]\subset\mathcal{L}^{2}[a,b]$. This is particularly very significant in the sense that we are able to upgrade the smoothness of the working data or information space from $\mathcal{L}^{2}[a,b]$ to $\mathcal{H}^{2}[a,b]$ and hence improve the smoothness of the recovery.

Remark 2.2*.*

Also, note that the negative Laplacian operator $-\Delta=-\frac{\partial^{2}}{\partial x^{2}}$ is a positive operator in $\mathcal{L}^{2}[a,b]$ on $\mathcal{D}_{\Delta}=\mathcal{H}^{2}_{0}[a,b]$. And using the positiveness of the operator $-\Delta$ on $\mathcal{H}^{2}_{0}[a,b]$ we can have a bound for $G(\psi)$ as

[TABLE]

where $\lambda_{1}>0$ is the first (smallest) eigenvalue of the positive operator $-\Delta$. Inequality (2.1) also implies, for a sequence of functions $\{\psi_{m}\}\subset\mathcal{L}^{2}[a,b]$, $G(\psi_{m})\rightarrow 0$ if and only if $\lvert\lvert u-u_{\psi_{m}}\rvert\rvert_{\mathcal{H}^{2}}\rightarrow 0$, details in §5.

3. Convexity of the functional G

Now let us define functionals $G_{1}$ and $G_{2}$ corresponding to the first and the second terms of the functional $G$ as follows

[TABLE]

We first analyze the functional $G_{2}$ and then extend it to the functional $G$.

Theorem 3.1**.**

1. (1)

An equivalent form of $G_{2}$, for any $\psi\in\mathcal{L}^{2}[a,b]$, is

[TABLE] 2. (2)

For any two $\psi_{1},\;\psi_{2}\in\mathcal{L}^{2}[a,b]$, we have

[TABLE] 3. (3)

The first Gâteaux differential, at $\psi\in\mathcal{L}^{2}[a,b]$, for $G_{2}$ is given by

[TABLE]

where $h\in\mathcal{L}^{2}[a,b]$. The $\mathcal{L}^{2}$-gradient of $G_{2}$, at $\psi$, is given by

[TABLE]

where the $T^{*}$ is the adjoint of the operator $T$. 4. (4)

The second Gâteaux differential, at any $\psi\in\mathcal{L}^{2}[a,b]$, of $G_{2}$ is given by

[TABLE]

where $h,\;k\in\mathcal{L}^{2}[a,b]$ and $\Delta=\frac{\partial^{2}}{\partial x^{2}}$ on $\mathcal{H}^{2}_{0}[a,b]$. Hence for any $\psi\in\mathcal{L}^{2}[a,b]$, $G_{2}^{\prime\prime}(\psi)$ is a positive definite quadratic form.

Thus $G_{2}$ is a strictly convex222the strict convexity follows if the operator $T$ is injective, otherwise, $G_{2}$ is a convex functional. functional and has a unique global minimum which is attained at $\varphi$.

First we state an ancillary result related to $u_{\psi}$.

Lemma 3.2**.**

For fixed $\psi,\;h\in\mathcal{L}^{2}[a,b]$, we have

[TABLE]

in $\mathcal{H}^{1}[a,b]$.

Proof.

Subtracting the following equations

[TABLE]

we have

[TABLE]

and using $u_{\psi+\epsilon h}-u_{\psi}\in\mathcal{H}^{2}_{0}[a,b]$ we get, via integration by parts,

[TABLE]

Now using the positiveness of the negative laplacian, we have ${\lambda_{1}}\lvert\lvert u_{\psi+\epsilon h}-u_{\psi}\rvert\rvert_{\mathcal{L}^{2}}^{2}\leq\lvert\lvert(u_{\psi+\epsilon h}-u_{\psi})^{\prime}\rvert\rvert_{\mathcal{L}^{2}}^{2}$, where $\lambda_{1}>0$ is the smallest eigenvalue of $-\Delta$, and thus

[TABLE]

where the last line is obtained using the Cauchy-Schwarz inequality. Hence,

$u_{\psi+\epsilon h}\xrightarrow{\epsilon\rightarrow 0}u_{\psi}$ in $\mathcal{H}^{1}[a,b]$, since the operator $T$ is bounded and $\psi,\;h\in\mathcal{L}^{2}[a,b]$ are fixed, which implies the right hand side is of $O(\epsilon)$. ∎

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}^{2}_{0}[a,b]$. In order to prove (3) and (4), we use Lemma 3.2.

• (3)

The Gâteaux derivative of the functional $G_{2}$ at $\psi$ in the direction of $h$, where $h\in\mathcal{L}^{2}[a,b]$, is given by

[TABLE]

Now for a fixed $\epsilon>0$, we have using 3.9,

[TABLE]

Using Lemma 3.2, one obtains the Gâteaux derivative of $G_{2}$ at $\psi\in\mathcal{L}^{2}[a,b]$ in the direction of $h\in\mathcal{L}^{2}[a,b]$ as

[TABLE]

Note that $(Th,-2(u-u_{\psi}))_{\mathcal{L}^{2}}=(h,T^{*}(-2(u-u_{\psi})))_{\mathcal{L}^{2}}$ for all $h\in\mathcal{L}^{2}[a,b]$, where $T^{*}$ is the adjoint of the operator $T$, as $T\in\mathcal{B}(\mathcal{H}_{1},\mathcal{H}_{2})$. Hence, by Riesz representation theorem, the $\mathcal{L}^{2}$-gradient of the functional $G_{2}$ at $\psi$ is given by

[TABLE]

• (4)

Finally, the second Gâteaux derivative for the functional $G_{2}$ at $\psi$ is given by

[TABLE]

Again for a fixed $\epsilon>0$, we have using (3.9)

[TABLE]

Hence from (3.11) and letting $\epsilon\rightarrow 0$ we get

[TABLE]

Here we can see the strict convexity of the functional $G_{2}$, as for any $h\in\mathcal{L}^{2}[a,b]$, we have

[TABLE]

where $-\Delta y=Th$ and $y\in\mathcal{H}^{2}_{0}[a,b]$ (from (3.9)). As $-\Delta$ is a positive operator on $\mathcal{H}^{2}_{0}[a,b]$, $y$ is the trivial solution iff $Th=0$ iff $h\equiv 0$ (T being strictly convex). Thus $G_{2}^{\prime\prime}(\psi)$ is a positive definite form for any $\psi\in\mathcal{L}^{2}[a,b]$. ∎

Now one can similarly prove the above properties for functional $G_{1}$ and extend it for functional $G$ as well, which is as follows:

Theorem 3.3**.**

---

1. (1)

For any two $\psi_{1},\;\psi_{2}\in\mathcal{L}^{2}[a,b]$, we have

[TABLE] 2. (2)

The first Gâteaux differential, at $\psi\in\mathcal{L}^{2}[a,b]$, for G is given by

[TABLE]

where $h\in\mathcal{L}^{2}[a,b]$. The $\mathcal{L}^{2}$-gradient of $G$, at $\psi$, is given by

[TABLE]

where $T^{*}$ is the adjoint of the operator $T$. 3. (3)

The second Gâteaux differential, at any $\psi\in\mathcal{L}^{2}[a,b]$, of G is given by

[TABLE]

where $h,k\in\mathcal{L}^{2}[a,b]$. Hence for any $\psi\in\mathcal{L}^{2}[a,b]$, $G^{\prime\prime}(\psi)$ is a positive definite quadratic form.

4. The Descent Algorithm for G

In this section we discuss the problem of minimizing the functional $G$, via a descent method. One can guess a descent direction by looking at the truncated Taylor expansion of the functional $G$, which is

[TABLE]

for any $\psi$, $h\in\mathcal{L}^{2}[a,b]$, and sufficiently small $\gamma>0$. Thus $G$ is minimized at $\psi$ if the direction $h$ is chosen such that $G^{\prime}(\psi)[h]>0$ for an appropriate $\gamma$.

The followings are some descent directions that can make $G^{\prime}(\psi)[h]>0$.

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

First, notice from Theorem 3.3 that at a given point $\psi\in\mathcal{L}^{2}[a,b]$

[TABLE]

so if we choose the direction $h=\nabla_{\mathcal{L}^{2}}^{\psi}G:=-2T^{*}(u-u_{\psi}+g-T\psi)$ at $\psi$, then $G^{\prime}(\psi)[h]>0$. Though this gradient works well in most situations, there are certain issues associated with it during the descent process. From a theoretical point of view, if $T^{*}(.)(x_{0})=0$ for some $x_{0}\in[a,b]$ at every step of the descent process then $\psi_{0}(x_{0})$ will be invariant during the descent process and if $\psi_{0}(x_{0})\neq\varphi(x_{0})$ then it will lead to severe decay or fluctuation near the point $x_{0}$, for details see [13]. For example, if $T=\int_{a}^{x}(.)dt$ is a Volterra operator then the $\mathcal{L}^{2}$-gradient at any $\psi$ is always zero at the end point $b$, since $T^{*}(.)=\int_{x}^{b}(.)dt$, which implies that at the end point $b$ the recovery will be invariant during the descent process, see Example 6.1.

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

One can circumvent the above problem by opting for the Sobolev or Neuberger gradient, $\nabla_{\mathcal{H}^{1}}^{\psi}G$, instead. It is the solution of the following boundary value problem

[TABLE]

This provides us a gradient $\nabla_{\mathcal{H}^{1}}^{\psi}G:=\phi$, at any $\psi$, with considerably more flexibility at the boundary points $\{a,b\}$. In particular one has

1. (1)

Dirichlet Neuberger gradient : $\phi(a)=0$ and $\phi(b)=0$. 2. (2)

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

Robin or Mixed Neuberger gradient : $\phi(a)=0$ and $\phi^{\prime}(b)=0$ or $\phi^{\prime}(a)=0$ and $\phi(b)=0$.

In addition to the flexibility at the end points, it also enables the new gradient to be a preconditioned (smoothed) version of $\nabla_{\mathcal{L}^{2}}^{\psi}G$, as $\phi=(I-\Delta)^{-1}\;\nabla_{\mathcal{L}^{2}}^{\psi}G$, and hence gives a superior convergence in the steepest descent algorithms when recovering a smooth function, see Example 6.1. One can also exploit the flexibility of the gradient at the end points according to some prior information of $\varphi$ at the end points. For example, if prior knowledge of $\varphi(a)$ and $\varphi(b)$ are known, then one can define $\varphi_{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}$’s during the descent process which leads to a more efficient and faster descent than compared to the normal $\mathcal{L}^{2}$-gradient. Even when $\varphi|_{\{a,b\}}$ is unknown, one can use the Neumann Neuberger gradient that allows free movements at the boundary points rather than gluing to some fixed values, see Example 6.1.

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

Now, one can make use of both the gradients by taking computing the standard Polak-Ribiére conjugate gradient scheme (see [14]) to further boost the descent rate, approximately by a factor of 2. Specifically, the initial search direction at $\psi_{0}$ is $h_{0}=\phi_{0}=\nabla_{\mathcal{H}^{1}}^{\psi_{0}}G$. At $\psi_{m}$ one can use an exact or inexact line search routine to minimize $G$ in the direction of $h_{m}$ resulting in $\psi_{m+1}$. Then $\phi_{m+1}=\nabla_{\mathcal{H}^{1}}^{\psi_{m+1}}G$ and $h_{m+1}=\phi_{m+1}+\gamma_{m}h_{m}$ where

[TABLE]

4.4. The line search method for the functional G

At any given $\psi_{m}\in\mathcal{L}^{2}[a,b]$ the functional $G$ is minimzed in the gradient direction via the single variable function $f_{m}(\gamma)=G(\psi_{m+1}(\gamma))$ and using a line search minimization, where $\psi_{m+1}(\gamma)=\psi_{m}-\gamma\nabla_{\mathcal{H}^{1}}^{\psi_{m}}G$. To bracket the minimum the initial step size $\gamma_{0}$ is chosen by solving the quadratic approximation of the function $f_{m}$ (see [13]), which is given by

[TABLE]

and since the expressions for $G^{\prime\prime}$ and $G^{\prime}$ are known, one can easily compute $\gamma_{0}$.

Hence during the descent process, 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})\geq 0\}$ is strictly decreasing. In the next section, we discuss the convergence of $\{\psi_{m}\}$ to $\varphi$ and the stability of the recovery.

5. Convergence and Regularization

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

5.1. Convergence

Note that, for the sequence $\{\psi_{m}\}$ produced by the steepest descent algorithm (described in Section 4) we have $G(\psi_{m})\rightarrow 0$, as the functional $G$ is non-negative and strictly convex (with the global minimizer $\varphi$, such that G($\varphi$) = 0) and $G(\psi_{m+1})<G(\psi_{m})$. The following theorems describe the convergence of the sequence {$\psi_{m}$} to $\varphi$ in $\mathcal{L}^{2}[a,b]$.

We first prove the convergence result for $G_{2}$ and then extend it for $G$.

Theorem 5.1**.**

Suppose that $\{\psi_{m}\}$ is any sequence of $\mathcal{L}^{2}$-functions such that the sequence $\{G_{2}(\psi_{m})\}$ tends to zero. Then $\{\psi_{m}\}$, and the corresponding sequence $\{g_{m}:=T\psi_{m}\}$, converge in $\mathcal{L}^{2}[a,b]$ to $\varphi$ and $g$, respectively, and the sequence $\{u_{\psi_{m}}\}$ converges to $u$ in $\mathcal{H}^{1}[a,b]$.

Proof.

First note that, using the positivity of $-\Delta$ in $\mathcal{L}^{2}[a,b]$ on $\mathcal{H}^{2}_{0}[a,b]$, we have

[TABLE]

where $\lambda_{1}=\frac{\pi^{2}}{(b-a)^{2}}>0$ is the smallest eigenvalue of $-\Delta$. This immediately gives $u_{\psi_{m}}\xrightarrow{s}u$ in $\mathcal{H}^{1}[a,b]$. Now one can easily see the weak convergence for the sequences $\{\psi_{m}\}$ and $\{g_{m}:=T\psi_{m}\}$ in $\mathcal{L}^{2}[a,b]$ to $\varphi$ and $g=T\varphi$, respectively. Since $u_{\psi_{m}}\xrightarrow{s}u$ in $\mathcal{H}^{1}[a,b]$, this implies $u_{\psi_{m}}\xrightarrow{w}u$ in $\mathcal{L}^{2}[a,b]$. Thus, for any $\psi$ in the range of $T^{*}$ we have

[TABLE]

implying $\psi_{m}\xrightarrow{w}\varphi$ in $R(T^{*})$, and since the range of $T^{*}$ is dense in $\mathcal{L}^{2}[a,b]$ for the linear, bounded and injective $T$, we have $\psi_{m}\xrightarrow{w}\varphi$ in $\mathcal{L}^{2}[a,b]$.

To see the strong convergence, we first express the functional $G_{2}$ in terms of an associated operator, i.e. $G_{2}(\psi)=\lvert\lvert u_{\psi}^{\prime}-u^{\prime}\rvert\rvert_{\mathcal{L}^{2}}^{2}=\lvert\lvert L(\psi)-\psi_{0}\rvert\rvert_{\mathcal{L}^{2}}^{2}$, for some operator $L$ such that $L(\psi)=u_{\psi}^{\prime}$. One can show that, from (1.10), the function $u_{\psi}\in\mathcal{H}^{2}_{0}[a,b]$ can be expressed as, for any $\psi\in\mathcal{L}^{2}[a,b]$ and $x\in[a,b]$,

[TABLE]

and hence, the operator $L(\psi):=u_{\psi}^{\prime}$ can be expressed as

[TABLE]

From (5.4) one can see that the operator $L$ is also linear and bounded in $\mathcal{L}^{2}[a,b]$, and thus, from the convergence theories developed for Landweber iterations or steepest descent methods, see [1, 15], the sequences $\{\psi_{m}\}$ and $\{g_{m}\}$ also converge to $\varphi$ and $g$, respectively, in $\mathcal{L}^{2}[a,b]$. ∎

Similarly, one can extend the convergence results for functional $G$ as

Theorem 5.2**.**

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}\}$, and the corresponding sequence $\{g_{m}:=T\psi_{m}\}$, converge in $\mathcal{L}^{2}[a,b]$ to $\varphi$ and $g$, respectively, and $\{u_{\psi_{m}}\}$ converges to $u$ in $\mathcal{H}^{2}[a,b]$.

5.2. Regularization (through appropriate early stopping)

As explained earlier, when solving an ill-posed problem with noisy data, one can not recover a regularized solution by simply minimizing the associated functional ($F$ in (1.2) or $G_{1}$ in (3.1)) completely. Here, the recovery error follows a semi-convergent nature, i.e., in the initial iterations the sequence $\{\psi_{m}^{\delta}\}$ converges towards the true solution ($\varphi$) but upon further iterations $\psi_{m}^{\delta}$ diverges away from $\varphi$, due to the noise and ill-posedness of the problem. Equivalently, during the initial stage $||\psi_{m}^{\delta}-\varphi||\rightarrow 0$ for $m\leq m_{0}(\delta)<\infty$ (for some appropriate $m_{0}(\delta)$), but upon further iterations $||\psi_{m}^{\delta}-\varphi||\rightarrow\infty$ (or $||\psi_{m}^{\delta}||\rightarrow\infty$), as $m\rightarrow\infty$. Hence such methods are also called as semi-iterative methods and terminating the minimization process at an appropriate early iteration leads to a regularized solution, i.e., the iteration index ($m$) serves as a regularization parameter, for details see [1, 15].

Here we show that the semi-convergent nature is also shown by the functional $G_{2}$ when recovering the true (or original) solution $\varphi$ but using a noisy data $g_{\delta}$. First we see that for an exact $g$ (or equivalently, an exact $u$) and the functional ($G$) constructed based on it, we have the true solution satisfying $G(\varphi)=0$. However, for a given noisy $g_{\delta}$, with $\lvert\lvert g_{\delta}-g\rvert\rvert_{\mathcal{L}^{2}}=\delta>0$, and the functional ($G_{\delta}$) based on it we will have $G_{\delta}(\varphi)>0$, see Theorem 5.3. So if we construct a sequence of functions $\{\psi_{m}^{\delta}\}\subset\mathcal{L}^{2}[a,b]$, using the descent algorithm and based on the noisy data $g_{\delta}$, such that $G_{\delta}(\psi_{m}^{\delta})\rightarrow 0$ then (from Theorem 5.1 or 5.2) we will have $\psi_{m}^{\delta}\rightarrow\varphi_{\delta}$, where $\varphi_{\delta}$ is the recovered noisy solution satisfying $G_{\delta}(\varphi_{\delta})=0$. This implies initially $\psi_{m}^{\delta}\rightarrow\varphi$ and then upon further iterations $\psi_{m}^{\delta}$ diverges away from $\varphi$ and approaches $\varphi_{\delta}$, such that $\lvert\lvert\varphi-\varphi_{\delta}\rvert\rvert_{\mathcal{L}^{2}}>>\delta$. This typical behavior of any ill-posed problem is managed, as stated above, by stopping the descent process at an appropriate iteration $M(\delta)$ such that $G_{\delta}(\psi_{M(\delta)}^{\delta})>0$ but close to zero, this gives us a regularized solution $\varphi^{\delta}:=\psi_{M(\delta)}^{\delta}$ such that $\lvert\lvert\varphi^{\delta}-\varphi\rvert\rvert_{\mathcal{L}^{2}}\leq C(\delta)$.

Following a similar argument as in (2.1) we have a lower bound for $G_{\delta}(\varphi)$ as

Theorem 5.3**.**

Let $G$ and $G_{\delta}$ be the functionals defined corresponding to the given data $g$ and $g_{\delta}$ (or $u,u_{\delta}$) and, $\varphi,\varphi_{\delta}$ denote their respective minimizers (i.e., $G(\varphi)=0$ and $G_{\delta}(\varphi_{\delta})=0$), then we have the following lower bound for $G_{\delta}(\varphi)$

[TABLE]

where $\lambda_{1}>0$ is the smallest eigenvalue of the positive operator $-\Delta$ on $\mathcal{H}^{2}_{0}[a,b]$, and vice-versa for $G(\varphi_{\delta})$.

Therefore, combining (2.1) and Theorem 5.3 we have the following two-sided bound for $G_{\delta}(\varphi)$, for some constants $C_{1}$ and $C_{2}$,

[TABLE]

Thus, when $\delta\rightarrow 0$ we have $G_{\delta}(\varphi)\rightarrow 0$ which implies $\varphi_{\delta}\rightarrow\varphi$ in $\mathcal{L}^{2}[a,b]$.

Though from (5.5) we would like to stop the descent process when

[TABLE]

but as we do not know the exact $g$ (equivalently, the exact $u$) we can not use (5.7) directly as the stopping criteria for the descent process.

5.3. Stopping criterion I

When the error norm $\delta=\lvert\lvert g-g_{\delta}\rvert\rvert_{\mathcal{L}^{2}}$ is known then Morozov’s discrepancy principle, [5], can serve as a stopping criterion for the descent process, i.e., terminate the iteration when

[TABLE]

for an appropriate $\tau>1$333In our experiments, we used the termination condition: $\lvert\lvert T\psi_{m}^{\delta}-g_{\delta}\rvert\rvert_{\mathcal{L}^{2}}<\delta$.. For the convergence and stability of the process see [7].

6. Numerical Results

We follow the following (pseudo-) Algorithm 1 and perform inverse recoveries on some standard integral type equations. A MATLAB code was written to test the numerical viability of the method and the results obtained are then compared with certain standard regularization methods, see Table 1 and Table 2.

6.1. Fredholm Integral Equation of the First Kind

A Fredholm integral equation of the first kind is an integral operator equation that depends on a kernel function $K(s,t)$ and is given by

[TABLE]

Such integral operator equations are classical inverse problems and can be quite ill-posed. The operator equation (6.1) is discretized using either Galerkin or Nystr$\ddot{o}$m methods to yield the following linear discrete ill-posed problem

[TABLE]

where the matrix $A\in\mathbb{R}^{m\times n}$ is the discretized representation of the operator $T$ and the vectors $x\in\mathbb{R}^{n}$, the $b\in\mathbb{R}^{m}$ are the discretized source and effect functions. MATLAB functions in [16, 17, 18] determine the discretizations $A\in\mathbb{R}^{n\times n}$, the scaled discrete approximations of $x\in\mathbb{R}^{n}$ and $b:=Ax\in\mathbb{R}^{m}$. To test the stability of the method, a Gaussian (zero-mean) error vector $\epsilon_{\delta}\in\mathbb{R}^{m}$ is added to $b$ to get a perturbed vector $b_{\delta}\in\mathbb{R}^{m}$, such that $\frac{||\epsilon_{\delta}||_{2}}{||b||_{2}}$ (the relative error norm) is around $10\%$, unless otherwise stated. In particular, when using discrepancy principle to terminate the descent process we assume $\delta=||\epsilon_{\delta}||_{2}$ to be known.

6.2. Numerical Differentiation

In the first example we study the problem of numerical differentiation (of noisy data), which is also governed by an integral equation (also known as a Volterra equation) of the following form

[TABLE]

for $x\in[a,b]$, and comparing (6.3) to (6.1), the kernel function can be considered as $K(s,t)=\chi_{[a,x]}(t)$, where $\chi$ is the standard characteristic function, for details see [13]. Now, for some given discrete noisy data values $\{g_{\delta}(x_{i}))\}_{i=1}^{n}$, we construct a discrete set of values $\{u_{\delta}(x_{i})\}_{i=1}^{n}$ for the function $u_{\delta}$, as defined in (1.10). Note that, when using the formula (1.10) to construct $u_{\delta}$ or $u_{\psi}$, the values of $u_{\delta}(x)$ and $u_{\psi}(x)$ depend on the choice of the value of the integration mesh-size $\Delta x$, because of the discretization of the integral, where for any integrable function $f$,

[TABLE]

Hence, in the following examples we show that increasing the $\Delta x$ values leads to smoother solutions as well as saturating the relative errors in the recovery. Since the kernel function $K(s,t)$ is a simple step-function, which is not smooth, one can notice (see Figure 2) that the recovered solution obtained using the discrepancy principle, though regularized, is also not smooth. However, with increasing $\Delta x$ values the smoothness of the recovery increases, see Figure 3(a). Also, note that here (as explained in §4.1) using $\mathcal{L}^{2}$ or $\mathcal{L}^{2}$-$\mathcal{L}^{2}$ conjugate gradient wouldn’t be a good choice if $\varphi(b)\neq 0$, as $(\nabla_{\mathcal{L}^{2}}G)(b)=0$ (since $T^{*}(b)=0$), which leads to a sharp decay at the end point ‘b’, see Figure 2(b) and Table 1. In such scenarios, either Neumann $\mathcal{H}^{1}$ or $\mathcal{L}^{2}$-$\mathcal{H}^{1}$ conjugate gradient will be very effective.

Example 6.1**.**

We contaminate two functions $g_{1}(x)=\frac{x^{3}}{3}-x$ and $g_{2}(x)=\frac{x^{3}}{3}-\frac{x}{2}$ on $[-1,1]$ with Uniform noise, to produce noisy $g_{1,\delta}$ and $g_{2,\delta}$ with a relative noise level of 10%. We collected evenly spread 200 of those noisy data points, which corresponds to a grid-size of 0.01. Note that, since $g_{2}^{\prime}(x)=x^{2}$ on $[-1,1]$, we have $\varphi_{2}(b)=g_{2}^{\prime}(1)\neq 0$. Now we compute the numerical regularized-derivatives from these noisy $g_{\delta}$’s data points and compare it with the true derivatives $g_{i}^{\prime}$’s. Figure 2 shows the comparison of the recoveries obtained using Tikhonov regularization and our method with the $\mathcal{L}^{2}$-$\mathcal{L}^{2}$ conjugate gradient and the $\mathcal{H}^{1}$-$\mathcal{L}^{2}$ conjugate gradient, with the relative errors in Table 1. Note that, when computing $g_{2,\delta}^{\prime}$ (where $g_{2}^{\prime}(b)\neq 0$) the use of Neumann $\mathcal{L}^{2}$-$\mathcal{H}^{1}$ conjugate gradient leads to a much better and smoother recovery, because of the flexibility at the boundary point $b$ and $\nabla_{\mathcal{H}^{1}}G=(1-\Delta)^{-1}\nabla_{\mathcal{L}^{2}}G$ (leading to a smoother gradient).

In the above computations we used the integration mesh-size $\Delta x=0.01$ for constructing $u_{\delta}$ and $u_{\psi}$. Then we repeat the experiment with the same number of noisy data points but with different $\Delta x$ values for the integration. Note that, here $\Delta x$ represents the mesh-size for the integral to get $u_{\delta}$ or $u_{\psi}$ values and is not connected to the original grid-size of the problem (which determines the number of data points). Figure 3(a) shows the increase in the smoothness level with the increase in the $\Delta x$ values and Figure 3(b) shows the saturating effect shown by increasing $\Delta x$ values. In fact, increasing $\Delta x$ values not only saturates the recovery errors during the descent process but also lowers the minimum error of the recovery (though we cannot obtain that solution directly from the discrepancy principle).

6.3. Imaging Problems

In this subsection we deal with higher dimensional inverse problems. Considering an image as a matrix leads to a 2-dimensional problem. However, we can transform a two dimensional $(n\times n)$-matrix to an one dimensional vector of length $n^{2}$ by simply stacking either the columns or rows over each other and hence, converting the 2-dimensional problem to a 1-dimensional problem. Now, since the integral limits can be arbitrary we assign the limits based on grid or mesh size, i.e., assuming $a=0$ and some value for $\Delta x$ we have $b=n^{2}\Delta x$. Here, since the transformed 1D vector corresponding to an image is highly discontinuous and non-smooth, having a large $\Delta x$ value leads to over-smooth the solution. Hence, we start with the initial choice $\Delta x=10^{-4}$ and show that with conservative increments one can attain a similar saturation of the relative errors as shown in Figure 3(b). We used $\nabla_{\mathcal{L}^{2}}G$ as the descent direction during the recovery process and the dependence of the relative errors on the integration mesh-size ($\Delta x$) is projected in Figure 4.

Example 6.2**.**

**Image deblurring

**In this example we apply our regularization technique to deblurr a noisy blurred image. The mechanism behind blurring an image also results in a Fredholm integral equation of first kind (6.1), with kernel being a Gaussian spread function, i.e., $K(s,t)=\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{\frac{(s-t)^{2}}{2\sigma^{2}}}$, where $\sigma$ determines the spread and the ill-posedness of the problem (ill-posedness increases with the increase in the $\sigma$ value, here we considered $\sigma=3$). A true test image ($\varphi$), see Figure 4(a), is first blurred and then contaminated with Gaussian noise ($\approx$ 10% relative error) to get a blurred noisy image ($g_{\delta}$), see Figure 4(c). Again, the discretized matrix $A\in\mathbb{R}^{n\times n}$ and the vectors $x$, $b$ in $\mathbb{R}^{n}$ (where $n=16384=128\times 128$) are generated using the routine, $blur()$, from [16, 18]. In this problem we performed a constraint minimization444i.e., we project the negative values to zero., $\psi\geq 0$, during the recovery process, since the pixel values are non-negative, and used the $\nabla_{\mathcal{L}^{2}}G$ gradient as the descent direction. Figure 4(e) shows the recovered image (using functional $G$) and Table 2 shows an error comparison between our recovery and the recoveries obtained by using some of the standard iterative regularization methods such as Least Squares minimization with (i) $\mathcal{L}^{1}$-norm (IRell1), (ii) (heuristic) total variation penalization (IRhtv), (iii) Conjugate Gradient (CGLS), (iv) FISTA for constrained LS (IRfista), (v) modified residual norm steepest descent (IRmrnsd), (vi) enriched CGLS (IRenrich), (vii) $\mathcal{L}^{1}$-norm penalization term (IRirn), (viii) hybrid FGMRES for enforcing $\mathcal{L}^{1}$-norm penalization (IRhybrid_fgmres), (ix) hybrid GMRES for square systems (IRhybrid_gmres), (x) hybrid LSQR (IRhybrid_lsqr), (xi) restarted Krylov subspace method (IRrestart), where the routines for the above mentioned methods are taken from [17], with the non-negative constraint (i.e., ‘nonnegativity’,‘on’) and the stopping rule as the discrepancy principle (‘stoprule’,‘DP’).

Example 6.3**.**

**[Tomography test problem]

**In this example we considered the tomography test problem from [18], where the matrix $A\in\mathbb{R}^{n\times n}$ and the vectors $x$, $b\in\mathbb{R}^{n}$ (where $n=16384$) are generated using the routine, $PRtomo()$, defined in [17]. The true image is shown in Figure 4(b) and the noisy image ($\approx$ 10% relative noise) for vector $b=Ax$ is shown in Figure 4(d). Here, we again performed a constrained minimization ($\psi\geq 0$) during the descent process, as in Example 6.2. The recovered image is shown in Figure 4(f) and Table 2 shows an error comparison between our method and some of the standard iterative regularization methods mentioned in Example 6.2.

Remark 6.4*.*

Again, one can notice (from Figure 4(g)) that increasing the $\Delta x$ values saturate the relative errors of the recoveries during the descent process. However, if $\Delta x$ value is increased significantly then it will also decrease the rate of the descent process, as can be seen in Figure 4(h), especially when dealing with highly discontinuous functions (such as images).

7. Conclusion and future research

In this paper we present a new approach for iterative regularization where instead of using the original noisy gradient we opted a smoother gradient, which is more robust to the noise present in the original data ($g_{\delta}$). The numerical results show that the method attains lower minimal relative errors during the recovery process, as well as, saturates the semi-convergent nature of the recovery relative errors. One can also notice that when dealing with high noise level data, where the operator $T^{*}$ fails to smooth the noisy effect of $g_{\delta}$, having smoother gradients is much more effective than the usual gradient direction.

In follow up papers we are trying to generalize this regularization method to solve non-linear inverse problems and also incorporate sparsity in the regularization process.

Acknowledgment

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