A new iterated Tikhonov regularization method for Fredholm integral equation of first kind

Xiaowei Pang Department of Mathematics, Hebei Normal University, Shijiazhuang, Hebei, China. Hebei Key Laboratory of Computational Mathematics and Applications, Hebei,China [email protected] and Jun Wang Department of Mathematics, Hebei Normal University, Shijiazhuang, Hebei, China. Hebei Research Center of the Basic Discipline Pure Mathematics, Hebei, China. [email protected]
Abstract.

We consider Fredholm integral equation of the first kind, present an efficient new iterated Tikhonov method to solve it. The new Tikhonov iteration method has been proved which can achieve the optimal order under a-priori assumption. In numerical experiments, the new iterated Tikhonov regularization method is compared with the classical iterated Tikhonov method, Landweber iteration method to solve the corresponding discrete problem, which indicates the validity and efficiency of the proposed method.

Key words and phrases:
Fredholm integral equation, Tikhonov regularization, optimal order, discrete ill-posed problem, Landweber iteration.
1991 Mathematics Subject Classification:
 65F10 \cdot 65F22

1. Introduction

In recent years, the research about inverse problems or ill-posed problems draw scientists’ attention. It can be emerged in earth physics, engineering technology and many other fields, such as geophysical problems [19], resistivity inversion problem [17], and computed tomography [18]. Therefore, the investigation of ill-posed problem not only has great scientific innovation significance, but also has certain practical importance.

In general, the inverse problem is much more difficult to solve than the forward problem, owing to its ill-posed feature. In the mid-1960s, the regularization method for dealing with ill-posed problems proposed by Tikhonov, brought the study of ill-posed problems into a new stage. Later, Landweber rewrote the equation (2) into an iteration form. Afterwards, many other technologies applied to regularization method came along successively, including precondition technique [13], adding contraction or penalty [5], multi-parameter regularization methods [6], filter based methods [9], methods coupling of them [10] or other methods [14, 15]. Klann etal. [11] and Hochstenbach etal. [9] discussed measuring the residual error in Tikhonov regularization with a seminorm that uses a fractional power of the Moore-Penrose pseudo inverse of ATAA^{T}A as weighting matrix, which lead to fractional filter methods. The former gave the fractional Landweber method and the latter presented fractional Tikhonov method. In [7], Huckle and Sedlacek also derived seminorms for the Tikhonov–Phillips regularization based on the underlying blur operator, that is using discrete smoothing-norms of the form Lx2\|Lx\|_{2} to substitute the classical 2-norm x2\|x\|_{2} for obtaining regularity, with LL being a discrete approximation to a derivative operator. Using a differential operator in the Tikhonov functional , it will be smoother and get a more accurate reconstruction. In [6], Gazzola and Reichel proposed two multi-parameter regularization methods for linear discrete ill-posed problems, which are based on the projection of a Tikhonov-regularized problem onto Krylov subspaces of increasing dimension. By selecting a proper set of regularization parameters and maximizing a suitable quantity, they can get the approximate solution. Stefano etal. proposed a nested primal–dual method for the efficient solution of regularized convex optimization problem in [1], under a relaxed monotonicity assumption on the scaling matrices and a shrinking condition on the extrapolation parameters, they gave the convergence result for the iteration sequence. Up to now, regularization methods are still powerful tools to settle inverse problems.

In this paper, our goal is to give a new iterated regularization method based on (17), for solving Fredholm integral equation of the first kind. This paper is organized as follows. In Section 2, we recall some basic definitions and preliminaries about the classical Tikhonov and Landweber method, the filter based regularization methods and the optimal order of a regularization method. In Section 3, we present a new iterated Tikhonov method and give the convergence result. Some numerical examples are reported in Section 4. Finally, Section 5 gives the conclusion.

2. Preliminaries

Fredholm integral equation of the first kind will be reviewed in this section firstly. Then, we recall some classical results about Tikhonov method and Landweber iteration.

2.1. Fredholm integral equation of the first kind

Many mathematical physics inverse problems, such as the backwards heat equation problem [2] and the image restoration problem [16], can be reduced to the following Fredholm integral equation of the first kind

(1) abK(s,t)x(s)𝑑s=y(t),s[a,b],\int_{a}^{b}K(s,t)x(s)ds=y(t),~~s\in[a,b],

where a,ba,b is finite or infinite, and x(s)x(s) is unknown, K(s,t)C(a,b)K(s,t)\in C(a,b) is known as a kernel function. If K(s,t)K(s,t) is a continuous kernel function, (1) will be written as the linear operator form

(2) Kx=y.Kx=y.

As we all know, (1) is ill-posed, that is to say, at least one of the existence, uniqueness and stability of the solution is not satisfied. Here it mainly refers to instability, that is, small perturbations in the data on the right hand side will lead to infinite variations in the solution. We consider the following example to illustrate.

Example 2.1.
01(1+ts)etsx(s)𝑑s=y(t)=et,0t1.\int_{0}^{1}(1+ts)e^{ts}x(s)ds=y(t)=e^{t},~~0\leq t\leq 1.

This equation has the unique solution x(t)=1x(t)=1. If we use the Simpson’s rule to approximate the integral, and the step size h=1nh=\frac{1}{n}, then we can get the linear system of equations

j=0nwj(1+titj)etitjx(tj)=y(ih),i=0,1,,n,\sum_{j=0}^{n}w_{j}(1+t_{i}t_{j})e^{t_{i}t_{j}}x(t_{j})=y(ih),~~i=0,1,\cdots,n,

where ww denotes the corresponding weight vector. Table 1 presents the error about x(ih)xix(ih)-x_{i} in different nodal point.

Table 1. The error between numerical solution and true solution at different points
tt n=4n=4 n=8n=8 n=164n=164 n=32n=32
0 0.0774-0.0774 0.1667-0.1667 4.9063-4.9063 1212
14\frac{1}{4} 1.07651.0765 0.4535-0.4535 13.0625-13.0625 32-32
12\frac{1}{2} 0.77300.7730 2.0363-2.0363 16.500016.5000 1313
34\frac{3}{4} 1.07491.0749 0.4393-0.4393 2-2 12-12
11 0.92580.9258 0.83410.8341 0.4063-0.4063 1919

From the above data, we can see that the error is not decreasing as the improvement of the calculation accuracy of the left integral term. It is dangerous to perform numerical calculations at this time. As we stated before, the error of the measure data is used to be inevitable, and we can’t ignore the rounding errors about the computer. So it is difficult to obtain stable numerical solutions for such problems. Based on the above reasons, a stabled method must be adopted—regularization method.

2.2. Tikhonov method and Landweber iteration method

The traditional Tikhonov regularization [8] solves the following minimization problem

(3) minxJαδ(x):=Kxyδ2+αx2,xX.\min_{x}J_{\alpha}^{\delta}(x):=\|Kx-y^{\delta}\|^{2}+\alpha\|x\|^{2},~~x\in X.

If the operator K:XYK:X\rightarrow Y is linear and bounded, the regularization α>0\alpha>0, then the unique minimum xα,δx^{\alpha,\delta} of JαδJ_{\alpha}^{\delta} is also the unique solution of the normal equation

(4) (αI+KK)xα,δ=Ky.\left(\alpha I+K^{*}K\right)x^{\alpha,\delta}=K^{*}y.

Let (μj,xj,yj)(\mu_{j},x_{j},y_{j}) be a singular system for KK, then the solution of Kx=yKx=y is presented by

(5) x=j=11μj(y,yj)xj.x=\sum\limits_{j=1}^{\infty}\frac{1}{\mu_{j}}(y,y_{j})x_{j}.

By the way, Tikhonov method gave a strategy

(6) q(α,μ)=μ2α+μ2q(\alpha,\mu)=\frac{\mu^{2}}{\alpha+\mu^{2}}

to damp the factor 1μj\frac{1}{\mu_{j}} of (5). The function q(α,μ):(0,)×(0,K]q(\alpha,\mu):(0,\infty)\times(0,\|K\|]\rightarrow\mathbb{R} is called as a regularizing filter function. Based on these information, Tikhonov regularization strategy Rα:YX,α>0R_{\alpha}:Y\rightarrow X,\alpha>0 is defined by

(7) Rαy=j=11μjμj2α+μj2(y,yj)xj,yY.R_{\alpha}y=\sum\limits_{j=1}^{\infty}\frac{1}{\mu_{j}}\frac{\mu_{j}^{2}}{\alpha+\mu_{j}^{2}}(y,y_{j})x_{j},~~~y\in Y.

Another methodology to give a regularizing filter function is

(8) q(α,μ)=1(1aμ2)1α,q(\alpha,\mu)=1-(1-a\mu^{2})^{\frac{1}{\alpha}},

if 1α=m\frac{1}{\alpha}=m, and mm represents the iterations, then (8) is the filter function of the Landweber iteration

(9) x0:=0,xm=(IaKK)xm1+aKy.x^{0}:=0,~~x^{m}=(I-aK^{*}K)x^{m-1}+aK^{*}y.

(6) and (8) are regularization filter functions, both of them satisfy the following definition.

Definition 2.2.

[3] Let K:XYK:X\rightarrow Y be compact with singular system (μj,xi,yj)(\mu_{j},x_{i},y_{j}), μ(K)\mu(K) be the closure of j=1{μj}\bigcup\limits_{j=1}^{\infty}\{\mu_{j}\}, and q:μ(K)(0,μ1)q:\mu(K)\subset(0,\mu_{1})\rightarrow\mathbb{R} be a function with the following properties:

(10a) supμj>0|q(α,μj)μj|=c(α)<,\displaystyle\sup_{\mu_{j}>0}|\frac{q(\alpha,\mu_{j})}{\mu_{j}}|=c(\alpha)<\infty,
(10b) |q(α,μj)|c<,c is independent of α,j,\displaystyle|q(\alpha,\mu_{j})|\leq c<\infty,~~\text{c is independent of $\alpha,j$},
(10c) limα0q(α,μj)=1pointwise in μj.\displaystyle\lim_{\alpha\rightarrow 0}q(\alpha,\mu_{j})=1~\text{pointwise in $\mu_{j}$}.

Let Rα:YXR_{\alpha}:Y\rightarrow X be a family operators, α>0\alpha>0, which is defined by

(11) Rαy=j=1q(α,μj)μj(y,yj)xj,yY,R_{\alpha}y=\sum\limits_{j=1}^{\infty}\frac{q(\alpha,\mu_{j})}{\mu_{j}}(y,y_{j})x_{j},~~~y\in Y,

then it is a regularization strategy or a filter based regularization method with Rα=c(α)\|R_{\alpha}\|=c(\alpha), and q(α,μ)q(\alpha,\mu) is called a filter function.

In addition, regularization operators corresponding to (6) and (8) are optimal order under an a-priori assumption [4], or are optimal strategies in the sense of the worst-case error [8]. For the integrity of the article, we give the definition. Next to the definition, there is a sufficient theorem to realize the optimal convergence rate.

Definition 2.3.

For given σ,E>0\sigma,E>0, let

Xσ,E:={xX|zX,zE,x=(KK)σ2z}X.X_{\sigma,E}:=\{x\in X|~\exists~z\in X,~\|z\|\leq E,~x=\left(K^{*}K\right)^{\frac{\sigma}{2}}z\}\subset X.

Define

(δ,σ,Rα):=sup{xxα,δ:xX1,yyδδ},\mathcal{F}(\delta,\sigma,R_{\alpha}):=sup\{\|x-x^{\alpha,\delta}\|:~x\in X_{1},~\|y-y^{\delta}\|\leq\delta\},

for any X1XX_{1}\subset X a subspace, δ>0\delta>0, and for a regularization method RαR_{\alpha}, if

(δ,σ,Rα)cδσσ+1E1σ+1\mathcal{F}(\delta,\sigma,R_{\alpha})\leq c\delta^{\frac{\sigma}{\sigma+1}}E^{\frac{1}{\sigma+1}}

holds, then a regularization method RαR_{\alpha} is called of optimal order under the a-priori assumption xXσ,Ex\in X_{\sigma,E}. If EE is unknown, then redefine a set

Xσ:=σ>0Xσ,E,X_{\sigma}:=\bigcup_{\sigma>0}X_{\sigma,E},

and if

(δ,σ,Rα)cδσσ+1,\mathcal{F}(\delta,\sigma,R_{\alpha})\leq c\delta^{\frac{\sigma}{\sigma+1}},

holds, then we call a regularization method RαR_{\alpha} is of optimal order under the a-priori assumption xXσx\in X_{\sigma}.

Theorem 2.4.

[12] Let K:XYK:X\rightarrow Y be a linear compact operator, Rα:YXR_{\alpha}:Y\rightarrow X is a filter based regularization method, it will be of optimal order under the a-priori assumption xXσ,Ex\in X_{\sigma,E}, σ,E>0\sigma,E>0,

(12a) sup0<μμ1|q(α,μ)μ|cαγ,\displaystyle\sup_{0<\mu\leq\mu_{1}}|\frac{q(\alpha,\mu)}{\mu}|\leq c\alpha^{-\gamma},
(12b) sup0<μμ1|(1q(α,μ))μσ|cσαγσ,\displaystyle\sup_{0<\mu\leq\mu_{1}}|(1-q(\alpha,\mu))\mu^{\sigma}|\leq c_{\sigma}\alpha^{\gamma\sigma},

with the regularization parameter α=c^(δE)1γ(σ+1)\alpha=\hat{c}\left(\frac{\delta}{E}\right)^{\frac{1}{\gamma(\sigma+1)}}, c^=(cσcσ)1γ(σ+1)>0\hat{c}=\left(\frac{c}{\sigma c_{\sigma}}\right)^{\frac{1}{\gamma(\sigma+1)}}>0.

3. Iterated Tikhonov regularization method

In this section, we first look back the standard Tikhonov iteration method, then introduce a new iterated Tikhonov regularization method, it is a generalization of the classical Tikhonov method.

The standard Tikhonov iteration method is

(13) x0,α,δ=0,(αI+KK)xm+1,α,δ=Kyδ+αxm,α,δ.x^{0,\alpha,\delta}=0,~~\left(\alpha I+K^{*}K\right)x^{m+1,\alpha,\delta}=K^{*}y^{\delta}+\alpha x^{m,\alpha,\delta}.

It can be shown that the corresponding regularization filter function is

(14) qm(α,μ)=1(αα+μ2)m,m=1,2,q^{m}(\alpha,\mu)=1-\left(\frac{\alpha}{\alpha+\mu^{2}}\right)^{m},~~m=1,2,...

In [4, 7], the authors introduced a Weighted-II Tikhonov method as the filter based method with the filter function

(15) ql(α,μ)=μ2μ2+α(1(μμ1)2)l,q_{l}(\alpha,\mu)=\frac{\mu^{2}}{\mu^{2}+\alpha\left(1-\left(\frac{\mu}{\mu_{1}}\right)^{2}\right)^{l}},

for α>0\alpha>0 and ll\in\mathbb{N}. Here, we also recall a filter based method—the fractional Tikhonov method [3, 9] with filter function

(16) qr(α,μ)=μ2r(α+μ2)r,q^{r}(\alpha,\mu)=\frac{\mu^{2r}}{\left(\alpha+\mu^{2}\right)^{r}},

for α>0\alpha>0 and r12r\geq\frac{1}{2}. Now, we can introduce a mixed method which combines the filter function of the fractional Tikhonov method and weighted-II Tikhonov method.

Definition 3.1.

Fixing qlr(α,μ)q_{l}^{r}(\alpha,\mu) such that

(17) qlr(α,μ):=μ2r(α(1(μμ1)2)l+μ2)r,q^{r}_{l}(\alpha,\mu):=\frac{\mu^{2r}}{\left(\alpha\left(1-\left(\frac{\mu}{\mu_{1}}\right)^{2}\right)^{l}+\mu^{2}\right)^{r}},

we define the mixed method (fractional weighted Tikhonov method) as the filter based method

(18) Rα,lry:=j=1qlr(α,μ)μj(y,yj)xj.R_{\alpha,l}^{r}y:=\sum\limits_{j=1}^{\infty}\frac{q_{l}^{r}(\alpha,\mu)}{\mu_{j}}(y,y_{j})x_{j}.

It is clear that for l=0l=0 and r=1r=1, it becomes the classical Tikhonov method.

Theorem 3.2.

Let K:XYK:X\rightarrow Y be a linear compact operator with infinite dimensional range and Rα,lrR_{\alpha,l}^{r} be the corresponding family mixed method operator. Then for every given r12r\geq\frac{1}{2}, ll\in\mathbb{N}, Rα,lrR_{\alpha,l}^{r} is a regularization method of optimal order under the a-priori assumption xXσ,Ex\in X_{\sigma,E} with 0<σ20<\sigma\leq 2. Further, if the regularization parameter satisfies α=(δE)2σ+1\alpha=(\frac{\delta}{E})^{\frac{2}{\sigma+1}}, then the best possible rate of convergence with respect to δ\delta is xxα,δ,r=𝒪(δ23)\|x-x^{\alpha,\delta,r}\|=\mathcal{O}(\delta^{\frac{2}{3}}) with σ=2\sigma=2. Moreover, if xxα,δ=𝒪(α)\|x-x^{\alpha,\delta}\|=\mathcal{O}(\alpha), then xX2x\in X_{2}.

Proof.

It has been proved ql(α,μ)q_{l}(\alpha,\mu) is a filter function in [4], so ql(α,μ)q_{l}(\alpha,\mu) satisfies the filter function conditions (10a-10c). By Proposition 12 of [3], for r12r\geq\frac{1}{2}, the function qlr(α,μ)q_{l}^{r}(\alpha,\mu) which meets the condition (10a-10c) can be verified easily. The proof of Qlr(α,μ)Q_{l}^{r}(\alpha,\mu) can meet the requirements (12a-12b) combining the proof of Proposition 12 in [3]. The difference is that Qα1(α,μ)Q_{\alpha}^{1}(\alpha,\mu) is the weighted-II Tikhonov method, and it also has the optimal order 𝒪(δ23)\mathcal{O}(\delta^{\frac{2}{3}}) with γ=12\gamma=\frac{1}{2} in (12b) for every 0<σ20<\sigma\leq 2. ∎

In the following, we will discuss the saturation for the mixed Tikhonov regularization.

Theorem 3.3.

Let K:XYK:X\rightarrow Y be a linear compact operator with infinite dimensional range and let Rα,lrR_{\alpha,l}^{r} be the corresponding family of fractional Tikhonov regularization operators in Definition 3.1 with r12r\geq\frac{1}{2}, ll\in\mathbb{N}. Let α=α(δ,yδ)\alpha=\alpha(\delta,y^{\delta}) be any parameter choice rule, and if

(19) sup{xxα,δ,r:P(yyδ)δ}=o(δ23),\sup\{\|x-x^{\alpha,\delta,r}\|:\|P(y-y^{\delta})\|\leq\delta\}=o\left(\delta^{\frac{2}{3}}\right),

then x=0x=0 with PP is the orthogonal projector onto R(K)¯\overline{R(K)}.

Proof.

For r=1r=1, it is clear that the saturation result follows from weight-II Tikhonov regularization [3]. For r1r\neq 1, we have

(20) xxrα,δ=j=11μj(1qlr(α,μj))(y,yj)xj,x-x_{r}^{\alpha,\delta}=\sum\limits_{j=1}^{\infty}\frac{1}{\mu_{j}}\left(1-q_{l}^{r}(\alpha,\mu_{j})\right)(y,y_{j})x_{j},

and

1qlr(α,μ)=(1+μ2α(1(μμ1)2)l)r(μ2α(1(μμ1)2)l)r(1+μ2α(1(μμ1)2)l)r1(1ql1(α,μ)).\displaystyle\begin{aligned} 1-q_{l}^{r}(\alpha,\mu)&=\frac{\left(1+\frac{\mu^{2}}{\alpha\left(1-\left(\frac{\mu}{\mu_{1}}\right)^{2}\right)^{l}}\right)^{r}-\left(\frac{\mu^{2}}{\alpha\left(1-\left(\frac{\mu}{\mu_{1}}\right)^{2}\right)^{l}}\right)^{r}}{\left(1+\frac{\mu^{2}}{\alpha\left(1-\left(\frac{\mu}{\mu_{1}}\right)^{2}\right)^{l}}\right)^{r-1}}\cdot\left(1-q_{l}^{1}(\alpha,\mu)\right).\end{aligned}

We notice that the above equality will be

1qlr(α,μ)=f(μ2a)(1ql1(α,μ)),1-q_{l}^{r}(\alpha,\mu)=f\left(\frac{\mu^{2}}{a}\right)\cdot\left(1-q_{l}^{1}(\alpha,\mu)\right),

where f(x)=(x+1)rxr(x+1)r1f(x)=\frac{(x+1)^{r}-x^{r}}{(x+1)^{r-1}}, a=α(1(μμ1)2)la=\alpha\left(1-\left(\frac{\mu}{\mu_{1}}\right)^{2}\right)^{l} and ql1(α,μ)=ql(α,μ)q_{l}^{1}(\alpha,\mu)=q_{l}(\alpha,\mu) is the weighted-II Tikhonov. f(x)f(x) is a monotone function, it satisfies f(0)=1f(0)=1 and limxf(x)=r\lim\limits_{x\rightarrow\infty}f(x)=r. Hence, we can get

min{1,r}(1ql1(α,μ))(1qlr(α,μ))max{1,r}(1ql1(α,μ)).\min\{1,r\}\left(1-q_{l}^{1}(\alpha,\mu)\right)\leq\left(1-q_{l}^{r}(\alpha,\mu)\right)\leq\max\{1,r\}\left(1-q_{l}^{1}(\alpha,\mu)\right).

Naturally,

sup{xxα,δ,r:P(yyδ)δ}min{1,r}sup{xxα,δ,1:P(yyδ)δ},\sup\{\|x-x^{\alpha,\delta,r}\|:\|P(y-y^{\delta})\|\leq\delta\}\geq\min\{1,r\}\cdot\sup\{\|x-x^{\alpha,\delta,1}\|:\|P(y-y^{\delta})\|\leq\delta\},

for every yδy^{\delta} satisfies yyδδ\|y-y^{\delta}\|\leq\delta. From Proposition 3.6 in [3], we have

sup{xxα,δ,1:P(yyδ)δ}=o(δ23),\sup\{\|x-x^{\alpha,\delta,1}\|:\|P(y-y^{\delta})\|\leq\delta\}=o(\delta^{\frac{2}{3}}),

Hence, the conclusion follows from the saturation result for Weighted-II Tikhonov (see Corollary 5.3 [3]). ∎

From now on, we propose a new iterated regularization method based on the above mixed Tikhonov. By iterations, we find that a large mm will provide a faster convergence rate (see Theorem 3.5).

Definition 3.4.

Define the iterated fractional weighted Tikhonov method as

(21)

(KK+α(IKKKK)l)rxm,α=(KK)r1Ky+[(KK+α(IKKKK)l)r(KK)r]xm1,α\left(K^{*}K+\alpha\left(I-\frac{K^{*}K}{\|K^{*}K\|}\right)^{l}\right)^{r}x^{m,\alpha}=(K^{*}K)^{r-1}K^{*}y+\left[\left(K^{*}K+\alpha\left(I-\frac{K^{*}K}{\|K^{*}K\|}\right)^{l}\right)^{r}-(K^{*}K)^{r}\right]x^{m-1,\alpha}

with x0,α:=0,r12,α>0x^{0,\alpha}:=0,~r\geq\frac{1}{2},\alpha>0 and ll\in\mathbb{N}. We can define xm,α,δx^{m,\alpha,\delta} as the mm-th iteration of (21) whenever yy is replaced by the noise data yδy^{\delta}.

In whole paper, for convenience, (21) will be called the new iterated Tikhonov method.

Theorem 3.5.

The new iterated Tikhonov in (21) is a filter based regularization method with filter function

(22) Qlm,r(α,μ)=1(1Qlr(α,μ))m,Q^{m,r}_{l}(\alpha,\mu)=1-\left(1-Q^{r}_{l}(\alpha,\mu)\right)^{m},

with Qlr(α,μ)=qlr(α,μ)=(μ2μ2+α(1(μμ1)2)l)rQ^{r}_{l}(\alpha,\mu)=q_{l}^{r}(\alpha,\mu)=\left(\frac{\mu^{2}}{\mu^{2}+\alpha(1-(\frac{\mu}{\mu_{1}})^{2})^{l}}\right)^{r}. Moreover, this method is of optimal order under the a-priori assumption xXσ,Ex\in X_{\sigma,E}, for ll\in\mathbb{N} and 0σ2m0\leq\sigma\leq 2m. Further, regularization parameter α=(δE)2m1+σ\alpha=\left(\frac{\delta}{E}\right)^{\frac{2m}{1+\sigma}}, yields the best convergence rate xxm,α,δ𝒪(δ2m2m+1)\|x-x^{m,\alpha,\delta}\|\leq\mathcal{O}(\delta^{\frac{2m}{2m+1}}) with σ=2m\sigma=2m.

Proof.

Denote C=(KK+α(IKKKK)l)rC=\left(K^{*}K+\alpha\left(I-\frac{K^{*}K}{\|K^{*}K\|}\right)^{l}\right)^{r}, B=C1(KK)r1KyB=C^{-1}(K^{*}K)^{r-1}K^{*}y, and A=C1[C(KK)r]A=C^{-1}\left[C-(K^{*}K)^{r}\right]. By the iteration formulas (21), we have

xm,α,δ=Axm1,α,δ+B=A2xm2,α,δ+(A1+A0)B==k=0m1AkB=k=0m1Ck[C(KK)r]kC1(KK)r1Ky.\displaystyle\begin{aligned} x^{m,\alpha,\delta}&=Ax^{m-1,\alpha,\delta}+B=A^{2}x^{m-2,\alpha,\delta}+(A^{1}+A^{0})B\\ &=\cdot\cdot\cdot\\ &=\sum\limits_{k=0}^{m-1}A^{k}B=\sum\limits_{k=0}^{m-1}C^{-k}\left[C-(K^{*}K)^{r}\right]^{k}C^{-1}(K^{*}K)^{r-1}K^{*}y.\end{aligned}

Let Rαm=k=0m1Ck[C(KK)r]kC1(KK)r1KR_{\alpha}^{m}=\sum\limits_{k=0}^{m-1}C^{-k}\left[C-(K^{*}K)^{r}\right]^{k}C^{-1}(K^{*}K)^{r-1}K^{*}, and the singular system be {μj,xj,yj}\{\mu_{j},x_{j},y_{j}\}, then

Rαmy=j=11μjQlm,r(α,μj)(y,yj)xj,\displaystyle\begin{split}\hskip-14.22636ptR_{\alpha}^{m}y&=\sum\limits_{j=1}^{\infty}\frac{1}{\mu_{j}}Q_{l}^{m,r}(\alpha,\mu_{j})(y,y_{j})x_{j},\end{split}

and

Qlm,r(α,μj)=k=0m1((μj2+α(1(μjμ1)2)l)rμj2r(μj2+α(1(μjμ1)2)l)r)k(μj2(μj2+α(1(μjμ1)2)l))r.\displaystyle\begin{split}\hskip-14.22636ptQ_{l}^{m,r}(\alpha,\mu_{j})=\sum\limits_{k=0}^{m-1}\left(\frac{\left(\mu_{j}^{2}+\alpha\left(1-(\frac{\mu_{j}}{\mu_{1}})^{2}\right)^{l}\right)^{r}-\mu_{j}^{2r}}{\left(\mu_{j}^{2}+\alpha\left(1-(\frac{\mu_{j}}{\mu_{1}})^{2}\right)^{l}\right)^{r}}\right)^{k}\left(\frac{\mu_{j}^{2}}{\left(\mu_{j}^{2}+\alpha\left(1-\left(\frac{\mu_{j}}{\mu_{1}}\right)^{2}\right)^{l}\right)}\right)^{r}.\end{split}

By the definition of ql(α,μ)q_{l}(\alpha,\mu), then

Qlm,r(α,μ)=k=0m1(1(ql(α,μ))r)k(ql(α,μ))r.Q_{l}^{m,r}(\alpha,\mu)=\sum\limits_{k=0}^{m-1}\left(1-\left(q_{l}(\alpha,\mu)\right)^{r}\right)^{k}\left(q_{l}(\alpha,\mu)\right)^{r}.

It is easily to get Qlm,r(α,μ)=1(1(ql(α,μ))r)mQ_{l}^{m,r}(\alpha,\mu)=1-\left(1-(q_{l}(\alpha,\mu))^{r}\right)^{m}, that is the conclusion as we stated.

From the relationship between Qlm,r(α,μ)Q_{l}^{m,r}(\alpha,\mu) and Qlr(α,μ)Q_{l}^{r}(\alpha,\mu), we can deduce

Qlr(α,μ)Qlm,r(α,μ)mQlr(α,μ).Q_{l}^{r}(\alpha,\mu)\leq Q_{l}^{m,r}(\alpha,\mu)\leq mQ_{l}^{r}(\alpha,\mu).

Clearly, ql(α,μ)q_{l}(\alpha,\mu) is weighted-II Tikhonov and it is a regularization filter method. Hence, Qlm,r(α,μ)Q_{l}^{m,r}(\alpha,\mu) satisfies the conditions (10a-10b) and (12a). At the same time, (10c) is easy to check, so it is a filter function naturally. Qlm,r(α,μ)Q_{l}^{m,r}(\alpha,\mu) adapt to the filter based regularization conditions. Finally, we make sure Qlm,r(α,μ)Q_{l}^{m,r}(\alpha,\mu) satisfies condition (12b) for the order optimality.

1Qlm,r(α,μ)=(1Qlr(α,μ))m1Qlr(α,μ)max{1,r}m(1Ql1(α,μ))m=c(1Qlm,1(α,μ)),\displaystyle\begin{aligned} 1-Q_{l}^{m,r}(\alpha,\mu)&=\left(1-Q_{l}^{r}(\alpha,\mu)\right)^{m}\leq 1-Q_{l}^{r}(\alpha,\mu)\\ &\leq\max\{1,r\}^{m}(1-Q_{l}^{1}(\alpha,\mu))^{m}=c\left(1-Q_{l}^{m,1}(\alpha,\mu)\right),\end{aligned}

and notice that Ql1(α,μ)=ql(α,μ)=μ2μ2+α(1(μμ1)2)lQ_{l}^{1}(\alpha,\mu)=q_{l}(\alpha,\mu)=\frac{\mu^{2}}{\mu^{2}+\alpha\left(1-(\frac{\mu}{\mu_{1}})^{2}\right)^{l}} is the weighted-II filter function and Qlm,1(α,μ)=1(1μ2α(1(μμ1)2)l+μ2)mQ_{l}^{m,1}(\alpha,\mu)=1-\left(1-\frac{\mu^{2}}{\alpha\left(1-\left(\frac{\mu}{\mu_{1}}\right)^{2}\right)^{l}+\mu^{2}}\right)^{m} is the filter function of the stationary iterated Tikhonov. So that condition (12b) follows from the properties of stationary Weighted-II iterated Tikhonov, and γ=12\gamma=\frac{1}{2}, 0σ2m0\leq\sigma\leq 2m, therefore, we get the best convergence rate 𝒪(δ2m2m+1)\mathcal{O}(\delta^{\frac{2m}{2m+1}}). ∎

4. Numerical experiments

The purpose of this section is to illustrate the validity from the previous sections with the following example. The classical iterative Tikhonov regularization method, Landweber method and the new iterated Tikhonov method are adopted to get the iterative numerical solutions.

Consider the following integral equation of the first kind:

(23) 0estx(t)𝑑t=h(s),0t<.\int_{0}^{\infty}e^{-st}x(t)dt=h(s),~~0\leq t<\infty.

The kernel operator is given by (Kx)(t)=0estx(s)𝑑s(Kx)(t)=\int_{0}^{\infty}e^{-st}x(s)ds. For numerical computation, we use Gauss-Laguerre quadrature rule with nn points to get the matrix AA corresponding to the kernel. The measure data about the right-hand side function is denoted by yδ=y+δηy^{\delta}=y+\delta\|\eta\|, where η\eta obeys the standard normal distribution, and the perturbation magnitude is δ\delta.

In this example, the right-hand side function h(s)=22s+1h(s)=\frac{2}{2s+1}, hence (23) has the unique solution x(t)=et2x(t)=e^{-\frac{t}{2}}. As mentioned above, we can use regularization method to solve the numerical solution. The operator KK is self-adjoint, so discrete Tikhonov equation takes the following form

(24) (αI+A2)xα,δ=Ayδ.\left(\alpha I+A^{2}\right)x^{\alpha,\delta}=Ay^{\delta}.

4.1. Classical method implementation

First, let the perturbation δ=0\delta=0, that is only the discrete error by quadrature rule will be generated, choose different regularization parameter α=10i,i=1,2,,10\alpha=10^{-i},i=1,2,...,10 by priori, and the quadrature points number n=16,32n=16,32. The numerical discrete errors variation diagram |xxα,δ|l2|x-x^{\alpha,\delta}|_{l^{2}} are showed in Figure 1(a). From Figure 1(a), if α\alpha is small, the error has a big difference between n=16n=16 and n=32n=32 as α<104\alpha<10^{-4} especially.

Refer to caption
(a) Numerical error for different regularization parameter nn in Tikhonov
Refer to caption
(b) Numerical error for different perturbation parameter δ\delta in Tikhonov
Refer to caption
(c) Numerical error for different perturbation parameter δ\delta in Landweber
Refer to caption
(d) The relationship between residual norm and solution norm for different perturbation parameter δ\delta in Tikhonov
Figure 1. Numerical error for nn and δ\delta in Tikhonov and Landweber method

Next, we consider the numerical error for different perturbation δ\delta in Tikhonov method (see Figure 1(b)), and in Landweber method with a=0.5a=0.5 to solve (24) and iteration steps m=1,2,,3000m=1,2,...,3000 (see Figure 1(c)). From the trends of the figures, they show that the numerical error first decrease then increase as α\alpha or mm increase, this coincides with the theory. Besides, we observe that both methods are comparable where precision is concerned.

Figure 1(d) presents the relationship of residual norm and solution norm, when the magnitude of perturbation δ=10j,j=1,2,3,4\delta=10^{-j},j=1,2,3,4 in Tikhonov method. As we can see, the small perturbation will have a small error with the same α\alpha basically. Besides, it looks like a LL curve, that is to say, there is a optimal α\alpha keeping the solution norm and residual norm balance.

4.2. New iterated Tikhonov implementation

Now we use the new iterated Tikhonov regularization method to solve (24). Let α=1e3\alpha=1e-3, a=0.5a=0.5, δ=1e4\delta=1e-4, l=4l=4, and n=32n=32, then we compare the total numeical error by using classical Tikhonov method (13), Landweber iteration method (9) and the new iterated Tikhonov method (21) when iteration steps changes, see Figure 2.

Refer to caption
Figure 2. Iterated Tikhonov error, Landweber iteration error and New iterated Tikhonov error for different mm

From Figure 2, we find that the new iterated Tikhonov method only need less iteration steps to get a smaller error than the other two methods for this problem under these parameters setting, which proves the validity of the proposed method. Finally, let l=2,m=100,r=0.8l=2,m=100,r=0.8, α=1e0,91e1,1e3,1e3\alpha=1e-0,9*1e-1,1e-3,1e-3 for different perturbation δ\delta, the following Table 2 gives the auxiliary specification to prove the efficiency of the new iterated Tikhonov method.

Table 2. The numerical error for different δ\delta by three methods
δ=1e4\delta=1e-4 δ=1e3\delta=1e-3 δ=1e2\delta=1e-2 δ=1e1\delta=1e-1
Iterated Tikhonov 0.2692 0.1035 0.0648 0.0051
Landweber 0.8654 0.1598 0.1336 0.1370
New iterated Tikhonov 0.2418 0.0734 0.0222 0.0046

5. Conclusion

This paper has shown the iterated fractional weight regularization method is an efficient method to solve the Fredholm integral equation of the first kind. The numerical experiments conducted have validated the accuracy of the proposed method and shown that the comparability with the classical Tikhonov method.

References

  • [1] S. Aleotti, S. Bonettini, M. Donatelli, M. Prato and S. Rebegoldi, A nested primal–dual iterated Tikhonov method for regularized convex optimization, Comput. Optim. Appl., (2024), 39pp.
  • [2] K. A. Ames and L. E. Payne, Continuous dependence on modeling for some well-posed perturbations of the backward heat equation, J. Inequal. Appl., 3 (1999), 51–64.
  • [3] D. Bianchi, A. Buccini, M. Donatelli and S. Serra-Capizzano, Iterated fractional Tikhonov regularization, Inverse Problems, 31 (2015), 34pp.
  • [4] D. Bianchi and M. Donatelli, On generalized iterated Tikhonov regularization with operator-dependent seminorms, Electron. Trans. Numer. Anal., 47 (2017), 73–99.
  • [5] Y. H. Chen, Y. Xiang, Z.Y. Shi, J. Lu and Y. J. Wang, Tikhonov regularized penalty matrix construction method based on the magnitude of singular values and its application in near-field acoustic holography, Mech. Syst. Sig. Process., 170 (2022), 108870.
  • [6] S. Gazzola and L. Reiche, A new framework for multi-parameter regularization, BIT, 56 (2016), 919–949.
  • [7] T. K. Huckle and M. Sedlacek, Tikhonov–Phillips regularization with operator dependent seminorms, Numer. Algorithms, 60 (2012), 339–353.
  • [8] A. Kirsch. An introduction to the mathematical theory of inverse problems, Springer, 2011.
  • [9] M. E. Hochstenbach and L. Reichel, Fractional Tikhonov regularization for linear discrete ill-posed problems, BIT, 51 (2011), 197–215.
  • [10] K. Ito, B. Jin and T. Takechi, Multi-parameter Tikhonov regularization - An augmented approach, Chinese Ann. Math. B, 35 (2014), 383–398.
  • [11] E. Klann and R. Ramlau, Regularization by fractional filter methods and data smoothing. Inverse Problems, 24 (2008), 26pp.
  • [12] A. K. Louis, Inverse und Schlecht Gestellte Probleme, B. G. Teubner Stuttgart, 1989.
  • [13] H. B. Li, A preconditioned Krylov subspace method for linear inverse problems with general-form Tikhonov regularization, SIAM J. Sci. Comput., 46 (2024), A2067–A2633.
  • [14] A. Molabahrami, An algorithm based on the regularization and integral mean value methods for the Fredholm integral equations of the first kind, Appl. Math. Model., 37 (2013), 9634–9642.
  • [15] M. Maleki and M. T. Kajani, Numerical approximations for volterra’s population growth model with fractional order via a multi-domain pseudospectral method, Appl. Math. Model. , 39 (2015), 4300–4308.
  • [16] H. Mesgarani and P. Parmour, Application of numerical solution of linear Fredholm integral equation of the first kind for image restoration, Math. Sci., 17 (2023), 371–378.
  • [17] O. Nechaev, V. Glinskikh, I. Mikhaylov and I. Moskaev, Joint inversion of high-frequency induction and lateral logging sounding data in Earth models with tilted principal axes of the electrical resistivity tensor, J. Inverse Ill-Posed Probl., 29 (2021), 295–304.
  • [18] H. S. Park, C. M. Hyun and J. K. Seo, Nonlinear ill-posed problem in low-dose dental cone-beam computed tomography, IMA J. Appl. Math., 89 (2023), 231–253.
  • [19] M. Signorini, S. Zlotnik and P. Díez, Proper generalized decomposition solution of the parameterized Helmholtz problem: application to inverse geophysical problems, Internat. J. Numer. Methods Engrg., 109 (2017), 1085–1102.