On fractional asymptotical regularization of linear ill-posed problems in Hilbert spaces

TL;DR

This paper introduces a fractional-order regularization method called FAR for solving linear ill-posed problems in Hilbert spaces, demonstrating accelerated convergence and numerical efficiency over existing methods.

Contribution

It proposes a novel fractional asymptotical regularization method with theoretical acceleration properties and develops a new iterative scheme based on Adams-Moulton for practical implementation.

Findings

FAR accelerates convergence compared to traditional methods.

Numerical examples confirm the accuracy and efficiency of FAR.

The method is applicable under certain smoothness conditions.

Abstract

In this paper, we study a fractional-order variant of the asymptotical regularization method, called {\it Fractional Asymptotical Regularization (FAR)}, for solving linear ill-posed operator equations in a Hilbert space setting. We assign the method to the general linear regularization schema and prove that under certain smoothness assumptions, FAR with fractional order in the range $(1,2)$ yields an acceleration with respect to comparable order optimal regularization methods. Based on the one-step Adams-Moulton method, a novel iterative regularization scheme is developed for the numerical realization of FAR. Two numerical examples are given to show the accuracy and the acceleration effect of FAR.

Tables

Table 1. Table 1: Comparisons with the Landweber method, Nesterov’s method, Chebyshev method, and CGNE method.

$\delta$	$\theta =0.5$		$\theta =0.8$		Landweber
	$k^{\ast }\left(\delta \right)$	L2Err	$k^{\ast }\left(\delta \right)$	L2Err	$k^{\ast }\left(\delta \right)$	L2Err
Example 1: $y\left(s\right)=s\left(1-s\right)$
0.01	4712	0.2694	1961	0.2684	712	0.2639
0.001	9336	0.1494	7057	0.1484	5997	0.1439
0.0001	$k_{max}$	0.2597	$k_{max}$	0.2073	$k_{max}$	0.1807
Example 2: $y\left(s\right)=s^4{\left(1-s\right)}^3$
0.001	3916	0.3038	2585	0.2918	1347	0.2833
0.0001	8534	0.2003	4947	0.1995	2034	0.1878
0.00001	$k_{max}$	0.3195	$k_{max}$	0.2184	$k_{max}$	0.1509
$\delta$	$\theta =1.2$		$\theta =1.5$		$\theta =1.8$
	$k^{\ast }\left(\delta \right)$	L2Err	$k^{\ast }\left(\delta \right)$	L2Err	$k^{\ast }\left(\delta \right)$	L2Err
Example 1: $y\left(s\right)=s\left(1-s\right)$
0.01	136	0.2476	29	0.2187	11	0.2204
0.001	733	0.1407	97	0.1321	45	0.1325
0.0001	8716	0.0997	1286	0.0673	657	0.0764
Example 2: $y\left(s\right)=s^4{\left(1-s\right)}^3$
0.001	224	0.2520	90	0.2019	58	0.2382
0.0001	883	0.0901	217	0.0767	155	0.0704
0.00001	4582	0.0287	328	0.0067	243	0.0083
$\delta$	Nesterov		Chebyshev		CGNE
	$k^{\ast }\left(\delta \right)$	L2Err	$k^{\ast }\left(\delta \right)$	L2Err	$k^{\ast }\left(\delta \right)$	L2Err
Example 1: $y\left(s\right)=s\left(1-s\right)$
0.01	109	0.2590	127	0.2553	6	0.2213
0.001	704	0.1600	556	0.1496	19	0.1783
0.0001	4160	0.1025	4361	0.0897	38	0.0894
Example 2: $y\left(s\right)=s^4{\left(1-s\right)}^3$
0.001	102	0.2501	75	0.2434	6	0.2919
0.0001	208	0.1176	207	0.0990	13	0.1814
0.00001	1805	0.0280	2226	0.0196	15	0.0847

Equations

$$Ax=y,$$

$$x_{k+1}^{\delta }=x_k^{\delta }+\Delta t{A}^{\ast }(y^{\delta }-A{x}_k^{\delta }),\Delta t\in (0,2/\Vert A{\Vert }^2)(k=0,1,\mathrm{2...})$$

$${\dot{x}}^{\delta }(t)+A^{\ast }A{x}^{\delta }(t)=A^{\ast }{y}^{\delta }$$

$$x^{†}=(A^{\ast }A{)}^{p}v,\Vert v\Vert \le \rho ,$$

$$\Vert x^{\delta }(T_{\ast })-x^{†}\Vert =\mathcal{O}({\delta }^{\frac{2p}{2p+1}})\text{mbox}as\delta \to 0.$$

$$\Vert x^{\delta }(T_{\ast })-x^{†}\Vert =\mathcal{O}({\delta }^{\frac{2p}{2p+1}})\text{mbox}and{T}_{\ast }=\mathcal{O}({\delta }^{-\frac{2}{2p+1}})\text{mbox}as\delta \to 0.$$

$${\ddot{x}}^{\delta }(t)+\eta {\dot{x}}^{\delta }(t)+A^{\ast }A{x}^{\delta }(t)=A^{\ast }{y}^{\delta },\eta >0.$$

$$\left({}^C{D}_{0+}^{\theta }{x}^{\delta }\right)(t)+A^{\ast }A{x}^{\delta }(t)=A^{\ast }{y}^{\delta },D^k{x}^{\delta }(0)=b_k,k=0,...,n-1,$$

$$x_{\alpha }^{\delta }=g_{\alpha }(A^{\ast }A)A^{\ast }{y}^{\delta },$$

$$x_{\alpha }=g_{\alpha }(A^{\ast }A)A^{\ast }A{x}^{†},$$

$$\begin{array}[]{ll}\|x^{\delta}_{\alpha}-x^{\dagger}\|&\leq\|x_{\alpha}-x^{\dagger}\|+\|x^{\delta}_{\alpha}-x_{\alpha}\|\\ &=\|r_{\alpha}(A^{*}A)x^{\dagger}\|+\delta\|g_{\alpha}(A^{*}A)A^{*}\|\\ &\leq\|r_{\alpha}(A^{*}A)x^{\dagger}\|+\gamma_{*}\delta/\sqrt{\alpha}.\end{array}$$

$$\underset{\delta \to 0}{\lim }\alpha =\underset{\delta \to 0}{\lim }{\delta }^2/\alpha =0,$$

$$\Vert r_{\alpha }(A^{\ast }A)x^{†}\Vert \le f(\alpha )(0<\alpha \le \bar{\alpha })$$

$$\Vert x_{{\alpha }_{\ast }}^{\delta }-x^{†}\Vert \le (1+{\gamma }_{\ast })f({\Theta }^{-1}(\delta ))=\mathcal{O}\left(f({\Theta }^{-1}(\delta ))\right)\text{mbox}as\delta \to 0.$$

$$\underset{\lambda \in (0,\Vert A{\Vert }^2]}{\sup }|r_{\alpha }(\lambda )|\phi (\lambda )\le \gamma \phi (\alpha )$$

$$\Vert x_{\alpha }^{\delta }-x^{†}\Vert \le \gamma {\alpha }^p\rho +{\gamma }_{\ast }\delta /\sqrt{\alpha }$$

$$\Vert x_{{\alpha }_{\ast }}^{\delta }-x^{†}\Vert \le \left(\gamma c^p+{\gamma }_{\ast }/\sqrt{c}\right){\delta }^{\frac{2p}{2p+1}}=\mathcal{O}\left({\delta }^{\frac{2p}{2p+1}}\right)\text{mbox}as\delta \to 0.$$

$$\underline{c}\frac{\psi (\alpha )}{\phi (\alpha )}\le \underset{\alpha \le \lambda \le \Vert A{\Vert }^2}{\inf }\frac{\psi (\lambda )}{\phi (\lambda )}(0<\alpha \le \bar{\alpha }).$$

$$x^{†}={\phi }_{\mu }(A^{\ast }A)v,\Vert v\Vert \le \rho ,$$

$$\varphi_{\mu}(\lambda)=\left\{\begin{array}[]{ll}\log^{-\mu}(1/\lambda)\qquad\textrm{~{}for}\quad 0<\lambda\leq e^{-\mu-1},\\ \mbox{arbitrarily extended as index function for}\;\;\lambda>e^{-\mu-1}.\end{array}\right.$$

$$\Vert x_{\alpha }^{\delta }-x^{†}\Vert \le \gamma {\phi }_{\mu }(\alpha )\rho +{\gamma }_{\ast }\delta /\sqrt{\alpha },$$

$$\Vert x_{{\alpha }_{\ast }}^{\delta }-x^{†}\Vert =\mathcal{O}\left({\phi }_{\mu }(\delta )\right)\text{mbox}as\delta \to 0$$

$$⟨\left({}^C{D}_{0+}^{\theta }{x}^{\delta }\right)(t),u_j⟩+{\sigma }_j^2⟨x^{\delta }(t),u_j⟩={\sigma }_j⟨y^{\delta },v_j⟩,j=1,2,....$$

$$\left({}^C{D}_{0+}^{\theta }{\xi }_j\right)(t)+{\sigma }_j^2{\xi }_j(t)={\sigma }_j⟨y^{\delta },v_j⟩,j=1,2,....$$

$$\eta (t)={\int }_0^t(t-\tau {)}^{\theta -1}{E}_{\theta ,\theta }(\lambda (t-\tau {)}^{\theta })\kappa (\tau )d\tau ,$$

$$E_{{\theta }_1,{\theta }_2}(z)=\sum _{k=0}^{\infty }\frac{z^k}{\Gamma ({\theta }_{1}k+{\theta }_2)}.$$

$${\xi }_j(t)={\sigma }_j⟨y^{\delta },v_j⟩{\int }_0^t(t-\tau {)}^{\theta -1}{E}_{\theta ,\theta }(-{\sigma }_j^2(t-\tau {)}^{\theta })d\tau ,$$

$${\int }_0^t{E}_{{\theta }_1,{\theta }_2}(\xi z^{{\theta }_1})z^{{\theta }_2-1}dz=t^{{\theta }_2}{E}_{{\theta }_1,{\theta }_2+1}(\xi t^{{\theta }_1}),{\theta }_2>0,$$

$${\xi }_j(t)={\sigma }_j{t}^{\theta }{E}_{\theta ,\theta +1}(-{\sigma }_j^2{t}^{\theta })⟨y^{\delta },v_j⟩,j=1,2,....$$

$$x^{\delta }(t)=t^{\theta }\sum _{j=1}^{\infty }{E}_{\theta ,\theta +1}(-{\sigma }_j^2{t}^{\theta }){\sigma }_j⟨y^{\delta },v_j⟩u_j=g^{\theta }(t,A^{\ast }A)A^{\ast }{y}^{\delta }$$

Full text

MSC 2010: Primary 47A52; Secondary 26A33, 65J20, 65F22, 65R30,

Key Words and Phrases: Linear ill-posed operator equation, asymptotical regularization, fractional derivatives, stopping rules, source conditions, convergence rates, acceleration

1 Introduction

It belongs to the main areas of competence of the journal “Fractional Calculus & Applied Analysis” to publish papers that handle initial-boundary value problems and abstract initial-value problems in Hilbert spaces for time-fractional diffusion and wave equations, partially also with respect to inverse and ill-posed problems. In this context, we refer for example to the papers [14, 17, 18, 19] by M. Yamamoto, Y. Luchko, and coauthors. The behavior of solutions to such equations plays an important role for finding stable approximations to inverse problems aimed at characterizing parameter functions of the underlying physical processes. As the papers mentioned above show, error estimates and explicit representations of solutions to the fractional equations can also be used to generate new continuous regularization methods. In the present paper, we try to complement this research direction by extending the well-known method of asymptotical regularization (Showalter’s method) in Hilbert spaces to a fractional version using the left-side Caputo fractional derivative. In particular, we are going to derive regularization properties by a stringent analysis and to show at least by means of case studies an acceleration effect based on the fractional derivatives.

Let $A$ be a compact linear operator acting between two infinite dimensional Hilbert spaces $\mathcal{X}$ and $\mathcal{Y}$ such that the range $\mathcal{R}(A)$ of $A$ is an infinite dimensional subspace of $\mathcal{Y}$. It is well-known that then $\mathcal{R}(A)$ is non-closed, i.e. $\overline{\mathcal{R}(A)}^{\mathcal{Y}}\not=\mathcal{R}(A)$, and the linear operator equation

[TABLE]

is ill-posed and requires some kind of regularization. This is typical for operator equations (1.1) which are models for linear inverse problems (cf., e.g., [9] or [10, § 2] and references therein). For simplicity, we denote by $\langle\cdot,\cdot\rangle$ and $\|\cdot\|$ in the sequel the inner products and norms for both Hilbert spaces $\mathcal{X}$ and $\mathcal{Y}$. In this paper, we are interested in stable approaches for solving (1.1) that can be realized by fast iterative algorithms. The corresponding stable approximate solutions should be based on noisy data $y^{\delta}$ of the exact right-hand side $y=Ax^{\dagger}$, where $x^{\dagger}=A^{\dagger}y$ denotes the minimum-norm solution of (1.1). In this context, we consider the noise model $\|y^{\delta}-y\|\leq\delta$ with noise level $\delta>0$.

The simplest iterative regularization approach for solving (1.1) seems to be the Landweber method, which is given by the iteration procedure

[TABLE]

with some starting element $x_{0}\in\mathcal{X}$, where $A^{*}$ denotes the adjoint operator of $A$. The continuous analog to (1.2) as $\Delta t$ tends to zero is known as asymptotic regularization or (see, e.g., [26, 27]). It is written as a first order evolution equation of the form

[TABLE]

with some initial condition, where an artificial scalar time $t$ is introduced. There must be chosen an appropriate finite stopping time $T_{*}=T_{*}(\delta)$ (a priori choice) or $T_{*}=T_{*}(\delta,y^{\delta})$ (a posteriori choice) in order to ensure the regularizing property $x^{\delta}(T_{*})\to x^{\dagger}$ as $\delta\to 0$.

For the asymptotic regularization it is well-known that for all $p>0$ under Hölder-type source conditions

[TABLE]

and for the stopping time $T_{*}$ selected according to the a priori choice $T_{*}=T_{*}(\delta)\sim\delta^{-\frac{2}{2p+1}}$ we have for all $p>0$ (cf., e.g., [26, Theorem 2]) the order optimal convergence rate

[TABLE]

From the analog result of the Landweber iteration in [5, Theorem 6.5] it can be concluded for asymptotical regularization that we have for the stopping time $T_{*}=T_{*}(\delta,y^{\delta})$ chosen according to Morozov’s discrepancy principle the formulas

[TABLE]

Moreover, it has been shown that by using Runge-Kutta integrators, all of the properties of asymptotic regularization (1.3) carry over to its numerical realization [24]. Hence, the continuous model (1.3) is of particular importance for studying the intrinsic properties of a broad class of general linear regularization methods for inverse problems, and can be used for the development of new iterative regularization algorithms by combining some appropriate numerical schemes.

A fatal defect for large-scale problems is the slow performance of Landweber iteration (too many iterations required for optimal stopping) as well as of the asymptotical regularization method, i.e. too excessive stopping times $T_{*}$ are required for obtaining optimal convergence rates (1.5). Therefore, in practice, accelerating strategies are usually used. The well-known methods are semi-iterative methods (e.g. the Brakhage’s $\nu$-method) and the Nesterov acceleration scheme. It has been proven that

• •

for the $\nu$-method, the optimal convergence rates can be obtained with approximately the square root of iterations than those needed for ordinary Landweber iteration [5, § 6.3]. However, in contrast to the Landweber iteration, the $\nu$-methods show a saturation phenomenon; i.e., the optimal convergence rate (1.5) and the asymptotic $k^{*}=\mathcal{O}(\delta^{-\frac{1}{2p+1}})$ hold only for $p\leq\nu$ and $p\leq\nu-0.5$, respectively.

• •

For the Nesterov acceleration scheme, the optimal convergence rates are obtained if $p\leq 1/2$ and if the iteration is terminated according to an a priori stopping rule. If $p>1/2$ or if the iteration is terminated according to the Morozov’s conventional discrepancy principle, only sub-optimal convergence rates can be guaranteed. In both cases, sub-optimal convergence rates can be obtained with the same acceleration speed as by the $\nu$-method [22].

Recently, inspired by these two accelerated regularization methods and the advantage of the continuous model, the authors extended (1.3) in [28] to the damped second order dynamics of the form

[TABLE]

It has been shown that, under condition (1.4), (1.7) exhibits the same convergence rate (1.5) as the asymptotic regularization method for all smoothness index $p>0$. By using the total energy discrepancy principle for choosing the terminating time and the damped symplectic integrators, a new iterative regularization algorithm has been proposed in this regard. Even though numerical experiments demonstrate the acceleration effect of the method (1.7), rigorous mathematical proof is still lacking.

In this paper, we replace the first derivative in the dynamical model (1.3) with appropriate fractional derivatives. Specifically, we consider in Hilbert spaces the following initial value problem

[TABLE]

to an evolution equation of fractional order $\theta$, where $\theta\in(0,2),n=\lfloor\theta\rfloor+1$, and $\lfloor\cdot\rfloor$ is the floor function. $D^{k}$ denotes the usual differential operator of order $k$. The left-side Caputo fractional derivative is defined by $\left({}^{C}D^{\theta}_{0+}x\right)(t):=I^{n-\theta}_{0+}D^{n}x(t)$, where $I^{n-\theta}_{0+}$ is the left-side Riemann-Liouville integral operator, i.e., $(I^{n-\theta}_{0+}x)(t):=\frac{1}{\Gamma(\theta)}\int^{t}_{0}\frac{x(t)}{(t-\tau)^{1-n+\theta}}d\tau$ with the gamma function $\Gamma(\cdot)$. Note that, for $\theta=1$, (1.8) coincides with Showalter’s method (asymptotical regularization). We will call the fractional dynamics (1.8) with an appropriate choice of terminating time Fractional Asymptotical Regularization (FAR). The main goal of this paper is to show that under smoothness assumptions imposed on the exact solution, FAR with $\theta\in(1,2)$ yields an accelerated optimal regularization method.

The paper is structured as follows: in Section 2, we review basic concepts in the general linear regularization theory so that they can be applied to the convergence analysis of FAR. The regularization properties of FAR under Hölder-type and logarithmic source conditions are studied in Sections 3 and 4, respectively. Section 5 is devoted to the numerical realization of FAR. Finally, concluding remarks are given in Section 6.

2 Linear regularization methods revisited

We start with some definitions, taken from [21] and [12], respectively.

Definition 2.1.

A function $\varphi:\mathbb{R}_{+}\to\mathbb{R}_{+}$ is called an index function if it is continuous, strictly increasing and satisfies $\lim_{\lambda\to 0+}\varphi(\lambda)=0$.

Definition 2.2.

A family of functions $g_{\alpha}(\lambda)\;(0<\lambda\leq\|A\|^{2})$, defined for regularization parameters $0<\alpha\leq\bar{\alpha}$, is called a generator function for a linear regularization method to the ill-posed linear operator equation (1.1) if the following three conditions are fulfilled:

• (i)

For the bias function $r_{\alpha}(\lambda)=1-\lambda g_{\alpha}(\lambda)$ we have for any fixed $\lambda\in(0,\|A\|^{2}]$ the limit condition $\lim_{\alpha\to 0}|r_{\alpha}(\lambda)|=0$.

• (ii)

There exists a constant $\gamma_{1}>0$ such that $|r_{\alpha}(\lambda)|\leq\gamma_{1}$ for all $\lambda\in(0,\|A\|^{2}]$ and for all $\alpha\in(0,\bar{\alpha}]$.

• (iii)

There exists a constant $\gamma_{*}>0$ such that $\sqrt{\lambda}|g_{\alpha}(\lambda)|\leq\gamma_{*}/\sqrt{\alpha}$ for all $\alpha\in(0,\bar{\alpha}]$.

Once a generator function family $g_{\alpha}(\lambda)$ is chosen, the approximate solutions to (1.1) based on the noisy data $y^{\delta}$ are calculated by the procedure

[TABLE]

where

[TABLE]

characterizes the noise-free analog to $x_{\alpha}^{\delta}$. By condition (iii) of Definition 2.2 and exploiting the triangle inequality, we obtain the well-known error estimates

[TABLE]

Furthermore, from properties (i) and (ii) of Definition 2.2 we deduce for $\alpha\to 0$ point-wise convergence $r_{\alpha}(A^{*}A)x\to 0$ for any $x\in\mathcal{X}$ (see, e.g., [5, Theorem 4.1]). If the regularization parameters $\alpha=\alpha(\delta)$ or $\alpha=\alpha(\delta,y^{\delta})$ are chosen such that

[TABLE]

one can derive convergence $\|x^{\delta}_{\alpha}-x^{\dagger}\|\to 0$ as $\delta\to 0$ by the estimate (2.3).

If we have an index function $f$ depending on $x^{\dagger}$ as upper bound of the form

[TABLE]

then this function also determines the error profile for the solution $x^{\dagger}$ in the noisy case, and was therefore termed ‘profile function’ in [12]. By using the auxiliary index function $\Theta(\alpha):=\sqrt{\alpha}f(\alpha)$ and choosing the regularization parameter a priori as $\alpha_{*}=\Theta^{-1}(\delta)$, we can derive under additional conditions the convergence rate

[TABLE]

Note that the three requirements in Definition 2.2 are not sufficient to yield profile functions $f$ for classes of solutions $x^{\dagger}$. It is well-known that convergence rates for approximate solutions of form (2.1) are connected with smoothness conditions imposed on the solution $x^{\dagger}$ with respect to the forward operator $A$. In order to measure the sensitivity of a regularization method with respect to possible smoothness assumptions, the concept of qualification in the sense of index functions is introduced:

Definition 2.3.

Let $\varphi$ be an index function. A regularization method (2.1) generated by $g_{\alpha}(\lambda)$ is said to have the qualification $\varphi$ if there are a constant $\gamma>0$, independent of $\alpha$, and a value $\bar{\alpha}>0$ such that for all $\alpha\in(0,\bar{\alpha}]$ the inequality

[TABLE]

is satisfied.

It should be noted that the concept of Definition 2.3 is an alternative to the traditional concept of qualification in the sense of a positive number or infinity used in [27] and [5, Chap. 4]. This alternative concept was originally introduced in [21], but it has been frequently used by other authors recently. Only if the index function $\varphi$ is a qualification in the sense of Definition 2.3 for the regularization method with generator functions $g_{\alpha}$, we can derive the rate result (2.5) from (2.4).

Example 2.1 (Hölder source conditions).

In this example we consider the index functions $\varphi(\lambda)=\lambda^{p}$ associated with the Hölder source conditions (1.4). If such $\varphi$ is for some $p>0$ a qualification for method generated by $g_{\alpha}$, then we obtain from (2.3)

[TABLE]

and hence for $\alpha_{*}=\alpha(\delta)=c\delta^{2p/(2p+1)}\;(c>0)$ the error estimate and convergence rate

[TABLE]

If a benchmark source condition $\psi$ is known to be a qualification for the method generated by $g_{\alpha}$, then other index functions $\varphi$ are also qualifications whenever they are covered by $\psi$, and we refer to [21, Def. 2] and [21, Prop. 3, Remark 5 and Lemma 2] for the following definition and proposition, respectively.

Definition 2.4.

Let $\psi(\lambda)\;(0<\lambda\leq\|A\|^{2})$ be an index function. Then an index function $\varphi(\lambda)\;(0<\lambda\leq\|A\|^{2})$ is said to be covered by $\psi$ if there is $\underline{c}>0$ such that

[TABLE]

Proposition 1.

The index function $\varphi$ is a qualification for the method generated by $g_{\alpha}$ if $\varphi$ is covered by $\psi$ and if $\psi$ is a qualification for that method. If the quotient function $\lambda\mapsto\frac{\psi(\lambda)}{\varphi(\lambda)}$ is increasing for $0<\lambda\leq\bar{\lambda}$ and some $\bar{\lambda}>0$, then $\varphi$ is covered by $\psi$. If, in particular, the index function $\varphi(\lambda)$ is concave for $0<\lambda\leq\bar{\lambda}$, then $\varphi$ is covered by $\psi(\lambda)=\lambda$.

Example 2.2 (Logarithmic source conditions).

The focus of this example is on logarithmic source conditions

[TABLE]

where for exponents $\mu>0$ the function $\varphi_{\mu}$ is defined as

[TABLE]

Note that the condition (2.7) for any $\mu>0$ is significantly weaker than the condition (1.4), even for arbitrarily small $p>0$. As a consequence, the expected logarithmic convergence rates are also significantly lower than in the Hölder case (2.6) and typically occur for severely ill-posed problems (cf., e.g., [13]). If $\varphi_{\mu}$ is for some $\mu>0$ a qualification for method generated by $g_{\alpha}$, then we obtain from (2.3)

[TABLE]

for sufficiently small $\alpha>0$, and hence for $\alpha_{*}=\alpha(\delta)=c\,\delta^{\kappa}\;(c>0,\linebreak 0<\kappa<2)$ the convergence rate

[TABLE]

as a consequence of $\lim_{\delta\to 0}\delta^{\eta}/\varphi_{\mu}(\delta)=0$, which is valid for all $\eta>0$.

3 FAR under Hölder-type source conditions

In this section, we show that the FAR method with fractional order $\theta>0$ introduced in the introduction can be assigned to the general linear regularization schema recalled in Section 2. Therefore, we start with the verification of the generator function $g_{\alpha}$ occurring in formula (2.1) associated with the fractional differential equation (1.8) and the corresponding regularization properties. For simplicity, let $b_{k}=0$ in (1.8) for all $k$ under consideration. The case with non-vanishing initial data can be analyzed similarly, see [28].

Let $\{\sigma_{j};u_{j},v_{j}\}_{j=1}^{\infty}$ be the well-defined singular system for the compact linear operator $A$, i.e. we have $Au_{j}=\sigma_{j}v_{j}$ and $A^{*}v_{j}=\sigma_{j}u_{j}$ with ordered singular values $\|A\|=\sigma_{1}\geq\sigma_{2}\geq\cdot\cdot\cdot\geq\sigma_{j}\geq\sigma_{j+1}\geq\cdot\cdot\cdot\to 0$ as $j\to\infty$. Since the eigenelements $\{u_{j}\}_{j=1}^{\infty}$ and $\{v_{j}\}_{j=1}^{\infty}$ form complete orthonormal systems (with the exception of null-spaces) in $\mathcal{X}$ and $\mathcal{Y}$, respectively, (1.7) is equivalent to

[TABLE]

Using the decomposition $x^{\delta}(t)=\sum\limits_{j=1}^{\infty}\xi_{j}(t)u_{j}$, we obtain

[TABLE]

On the other hand, according to [15, § 4.1.3], the solution to the equation $\left({}^{C}D^{\theta}_{0+}\eta\right)(t)-\lambda\eta(t)=\kappa(t)$ with the vanishing initial data is explicitly given by (see formula (4.1.62) in [15])

[TABLE]

where the two-parametric Mittag-Leffler function $E_{\theta_{1},\theta_{2}}(z)$ is defined as

[TABLE]

Therefore, together with (3.2), we deduce that the solution of (3.1) with the vanishing initial data is given by

[TABLE]

Consequently, by the above result, together with the identity [8, § 4.4]

[TABLE]

we deduce for all $\theta>0$ that

[TABLE]

Thus, by using $x^{\delta}(t)=\sum_{j=1}^{\infty}\xi_{j}(t)u_{j}$, we obtain the explicit formula

[TABLE]

for the solution to (1.7), where the generator function $g^{\theta}(t,\lambda)$ characterizing the FAR method attains the form

[TABLE]

Furthermore, by using the recurrence relations

[TABLE]

we obtain an associated bias function $r^{\theta}(t,\lambda)=1-\lambda g^{\theta}(t,\lambda)$ of the form

[TABLE]

where $E_{\theta}(z)=\sum^{\infty}_{k=0}\frac{z^{k}}{\Gamma(\theta k+1)}$ denotes the classical Mittag-Leffler Function.

Remark 3.1.

For $\theta=0.5$, we have $r(t,\lambda)=e^{-\lambda^{2}t}\left[1+\textrm{erf}(-\lambda\sqrt{t})\right]$, where $\textrm{erf}(z):=\frac{2}{\sqrt{\pi}}\int^{z}_{0}e^{-\tau^{2}}d\tau$ represents the error function. Moreover, $r(t,\lambda)=e^{-\lambda t}$ for $\theta=1$, while $r(t,\lambda)=\cos(\sqrt{\lambda}t)$ when $\theta=2$.

Lemma 3.1.

Let $\theta\in(0,2)$ be a fixed number. Then we have for all $z\in\mathbb{R}_{+}$ the following three inequalities

[TABLE]

where $C_{\theta}>0$ is a constant, depending only on $\theta$.

P r o o f..

Inequalities (3.6) and (3.7) follow from [23, Theorem 1.6] and [8, Corollary 3.7] respectively. For the estimate (3.8), we distinguish three cases: (I) $\theta\in(0,1)$, (II) $\theta=1$, and (III) $\theta\in(1,2)$. In this first case, (3.8) holds according to the uniform estimate [25, Theorem 4]:

[TABLE]

For the second case, we have $E_{1}(-z)=e^{-z}\leq 1$. Now, we consider the last case. To this end, we set $z:=t^{\theta}$ and consider the function $E_{\theta}(-t^{\theta})$. It is well known (see, e.g., [20]) that $E_{\theta}(-t^{\theta})=f_{\theta}(t)+h_{\theta}(t)$, where the completely monotone function $-f_{\theta}$ is defined by

[TABLE]

and the oscillatory part $h_{\theta}$ is given by

[TABLE]

Obviously, for $\theta\in(1,2)$, $f_{\theta}$ is a monotonically increasing function such that $-1<1-\frac{2}{\theta}=\frac{1}{\pi}\int^{\infty}_{0}\frac{r^{\theta-1}\sin(\theta\pi)}{r^{2\theta}+2r^{\theta}\cos(\theta\pi)+1}dr\leq f_{\theta}(t)\leq 0$, while $|h_{\theta}(t)|\leq 2/\theta$ as $\cos(\pi/\theta)<0$. This implies $|E_{\theta}(-t^{\theta})|\leq 3$. $\Box$

We remark that $C_{\theta}\to 1$ as $\theta\searrow 1$, and $C_{\theta}\to\infty$ as $\theta\nearrow 2$.

In order to assign FAR to the general linear regularization schema introduced in Section 2, we exploit the one-to-one correspondence between the artificial time $t>0$ and the conventional regularization parameter $\alpha>0$ by setting $t:=\alpha^{-1/\theta}$.

Theorem 3.1.

For $\theta\in(0,2)$, $g_{\alpha}(\lambda):=g^{\theta}(\alpha^{-1/\theta},\lambda)$ represents a generator function in the sense of Definition 2.2 with $\theta_{*}=C_{\theta}/2$ and $C_{\theta}$ from Lemma 3.1. Consequently, FAR is a linear regularization method for the ill-posed linear operator equation (1.1).

P r o o f..

We have to check the three requirements (i), (ii) and (iii) in Definition 2.2 for all $\theta\in(0,2)$.

Since $r_{\alpha}(\lambda)=r^{\theta}(\alpha^{-1/\theta},\lambda)=E_{\theta}(-\lambda\alpha^{-1})$, it follows from (3.7) that

[TABLE]

which shows that (i) is satisfied. The second condition (ii) is a consequence of (3.8). It remains to find a bound $\gamma_{*}$ in (iii). Note that $g_{\alpha}(\lambda)=\alpha^{-1}E_{\theta,\theta+1}(-\lambda\alpha^{-1})$. By using (3.6), we obtain the inequalities

[TABLE]

and hence $\gamma_{*}=C_{\theta}/2$, which completes the proof. $\Box$

Referring back to the time variable $t$ we obtain a linear regularization method for the ill-posed operator equation (1.1) with the procedure

[TABLE]

The behavior of the bias function $r_{\alpha}(\lambda)$, or equivalently $r_{t^{-\theta}}(\lambda)$, is visualized in Figure 1. It follows from Remark 3.1 that for $\theta=2$, $r_{t^{-2}}(\lambda)=\cos(\sqrt{\lambda t^{2}})$. Obviously, we have for all $\lambda$ under consideration $r_{t^{-2}}(\lambda)\not\to 0$ as $t\to\infty$, which damages the first condition of Definition 2.2. In order to obtain a regularization method in the case $\theta=2$, an additional damping term in the model (1.8) should be introduced, see [28] for details.

Proposition 2.

For all $\theta\in(0,1)\cup(1,2)$ the function $\varphi(\lambda)=\lambda^{p}$ is a qualification of the FAR method if and only if $0<p\leq 1$, i.e. there is a saturation for $p=1$.

P r o o f..

The qualification of the FAR can be deduced from the estimate (3.7) and the following inequalities

[TABLE]

where $\gamma(p,\theta)=C_{\theta}p^{p}(1-p)^{1-p}$ for $0<p<1$, and $\gamma(p,\theta)=C_{\theta}\|A\|^{2(p-1)}$ for $1\leq p<2$. The saturation at $p=1$ follows from the lower bound in (3.9) and the asymptotical behavior, see e.g. [23, § 1.2], $E_{\theta}(-z)z\to 1/\Gamma(1-\theta)$ as $z\to\infty$. $\Box$

Remark 3.2.

By the well-known results of Showalter’s method, see e.g. [28, Remark 2.2], $\varphi(\lambda)=\lambda^{p}$ is a qualification of the FAR method with parameter $\theta=1$ for all $p>0$ without saturation.

Theorem 3.2.

Let $x(t)$ and $x^{\delta}(t)$ be solutions of (1.8) with noise-free and noisy data, respectively. Then, under the assumption $x^{\dagger}\in\mathcal{R}((A^{*}A)^{p})$ (cf. (1.4)), we have for the FAR method in the case $\theta\in(0,1)\cup(1,2)$ and $p\in(0,1]$ the convergence rate

[TABLE]

whenever the terminating time $T_{*}=T_{*}(\delta)$ is chosen according to

[TABLE]

Moreover, the following additional assertions are true if $\theta\in(0,1)$:

• (i)

The convergence rate of method error $\|x(t)-x^{\dagger}\|=\mathcal{O}(t^{-\theta p})$ holds, if and only if $x^{\dagger}$ satisfies a source condition of type (1.4).

• (ii)

The worst-case convergence rate

[TABLE]

is equivalent to the following weaker condition

[TABLE]

• (iii)

If a terminating time $T_{*}$ exists such that

[TABLE]

then $x^{\dagger}=0$.

P r o o f..

The assertions concerning the convergence rate (3.15) follow directly from (2.6) by setting $T_{*}(\delta)=(\alpha_{*})^{-1/\theta}$ in combination with Proposition 2. The remaining assertions (i), (ii) and (iii) of the theorem are based on standard arguments of regularization theory ([5, § 4.2] and [6]) in light of estimates (3.6), (3.9), and (3.11). $\Box$

Remark 3.3.

(a) The first part of Theorem 3.2 means that the optimal convergence rates (3.15) of FAR can be obtained with approximately $1/\theta$ of iterations than needed for ordinary asymptotical regularization method [5]. Therefore, FAR with $\theta\in(1,2)$ yields an accelerated regularization method for all $p>0$. As the estimate (2.6) shows, for fixed $\delta>0$ the error of regularization is proportional to $\gamma$ and $\gamma_{*}$. These constants, however, tend to infinity as $\theta\nearrow 2$ such that values $\theta$ close to $2$ need not be the best choices of this parameter. (b) The source condition (1.4) implies the weaker source condition (3.18). The converse is in general not true, however, (3.18) implies $x^{\dagger}\in\cap_{\nu<p}\mathcal{R}((A^{*}A)^{\nu})$, see [5, Lemma 4.12].

In practice, the a priori stopping rule (3.16) in Theorem 3.2 is not realistic, since a good terminating time $T_{*}$ requires in addition to the knowledge of $\delta$ also the knowledge of $p$ characterizing the specific smoothness of the unknown exact solution $x^{\dagger}$. Therefore, a posteriori parameter choices of the stopping parameter are preferred, and we consider Morozov’s discrepancy principle as the most prominent version exploiting zeros of the discrepancy function

[TABLE]

where we assume $\tau>3$ for the occurring factor of the noise level $\delta$.

Lemma 3.2.

If for $\tau>3$ the condition $\|y^{\delta}\|>\tau\delta$ is satisfied, the function $\chi(T)$ has at least one root.

P r o o f..

By the explicit formula of $x^{\delta}(t)$, c.f. (3.13), we derive together with the estimates (3.8) and (3.14) that

[TABLE]

as $t\to\infty$. The continuity of $\chi(t)$ is obvious as we are dealing with the linear problem. Since $\chi(0)=\|y^{\delta}\|-\tau\delta>0$, the existence of the root of $\chi(t)$ follows from the Bolzano’s theorem. $\Box$

Remark 3.4.

As the proof of Lemma 3.2 indicates, the function $\chi(t)$ is bounded above by a decreasing function $\zeta(t):=\gamma(p,\theta)\|x^{\dagger}\|t^{-\theta/2}-(\tau-3)\delta$. Thus, the trend of $\chi(t)$ is to be a decreasing function, where oscillations may occur, and we refer to Figure 2 for illustration.

Theorem 3.3.

Under the assumption $x^{\dagger}\in\mathcal{R}((A^{*}A)^{p})$, we have for the terminating time $T_{*}=T_{*}(\delta,y^{\delta})$ of FAR chosen as a solution of $\chi(T)=0$ (cf. 3.20) the convergence rates

[TABLE]

for $p\in(0,1]$ if $\theta\in(0,1)\cup(1,2)$ and for all $p>0$ if $\theta=1$.

P r o o f..

Using the interpolation inequality $\|B^{p}u\|\leq\|B^{q}u\|^{p/q}\|u\|^{1-p/q}$ and the source conditions $x^{\dagger}=(A^{*}A)^{p}v$, we deduce that

[TABLE]

Since $T_{*}$ is chosen according to the equation $\chi(T)=0$, we derive that

[TABLE]

Now we combine the estimates (3.25) and (3.29) to obtain, with the source conditions, that

[TABLE]

where $c^{\prime}_{2}:=\left(\tau+3\right)^{2p/(2p+1)}$.

On the other hand, in a similar fashion to (3.29), it is easy to show that

[TABLE]

If we combine the above inequality with the source conditions and the qualification inequality, we obtain

[TABLE]

which yields the estimate for $T_{*}$ in (3.21). Finally, using (2.3), the estimate for $T_{*}$ and (3.33), we conclude that

[TABLE]

This completes the proof. $\Box$

4 FAR under logarithmic source conditions

Now, we turn to the study of FAR method under conditions (2.7).

Proposition 3.

For all $\theta\in(0,2)$ and all $\mu>0$ the index function $\varphi_{\mu}(\lambda)$, defined in (2.8), is a qualification of the FAR method.

P r o o f..

For all $\theta\in(0,2)$ the function $\psi(\lambda)=\lambda$ is a qualification of the FAR method, as shown in Proposition 2. Moreover, for arbitrary $\mu>0$ the index function $\varphi_{\mu}(\lambda)$ is concave for all $0<\lambda\leq e^{-\mu-1}$, hence due to Proposition 1 covered by $\psi$ and consequently also a qualification of the FAR method. This proves the proposition. $\Box$

Based on the above proposition, the following theorem holds as a consequence of the discussions in Example 2.2 (cf. formula (2.9)).

Theorem 4.1.

Let $x^{\delta}(t)$ denote the solution of (1.8) with noisy data $y^{\delta}$ and noise level $\delta>0$. Then, under the assumption $x^{\dagger}\in\mathcal{R}(\varphi_{\mu}(A^{*}A))$ (cf. (2.7)) for arbitrary $\mu>0$, we have for the FAR method and all $0<\theta<2$ the convergence rate

[TABLE]

whenever the terminating time $T_{*}$ of FAR is chosen according to

[TABLE]

for some $0<\kappa<2$,

Now, consider the a posteriori choice of regularization parameter. We start with the following lemma.

Lemma 4.1.

The index function $\sqrt{\lambda}\varphi_{\mu}(\lambda)$ is a qualification of FAR, where $\varphi_{\mu}$ is defined in (2.8). Consequently, a constant $c_{3}$ exists such that

[TABLE]

P r o o f..

The assertion of the lemma follows from Proposition 3 in combination with the fact that $\varphi(\lambda)=\lambda$ is a qualification of FAR for all $\theta\in(0,2)$, and with the following inequality

[TABLE]

$\Box$

Theorem 4.2.

Under the assumptions of Theorem 4.1, but for the terminating time $T_{*}=T_{*}(\delta,y^{\delta})$ of FAR chosen as a solution of $\chi(T)=0$ (cf. 3.20)), we have

[TABLE]

where $W_{a}(z)$ is the unique solution of $\zeta\log^{a}(\zeta)=z$ with respect to $\zeta$, and moreover the convergence rate

[TABLE]

P r o o f..

Without loss of generality we assume that $\delta<e^{-2\mu-1}$. By combining inequality (3.34) and Lemma 4.1 with $\alpha$ replaced by $T^{-\theta}_{*}$, we obtain

[TABLE]

Therefore, for $\delta\in(0,e^{-1})$, and if we set $c_{4}:=c_{3}\theta^{-\mu/2}(\tau-3)^{-1}$, we have

[TABLE]

which yields the estimate (4.2), as well as the inequality

[TABLE]

The estimate (4.3) can be driven according to inequalities (2.3), (4.5), and to the Proposition 3. $\Box$

Remark 4.1.

(a) To the best of our knowledge, under the logarithmic source condition, for obtaining the optimal convergence rate (4.3), existing regularization methods with a posteriori regularization parameter selection methods require the time cost $T_{old}=\mathcal{O}\left(\delta^{-2}\log^{-2-2\mu}(\delta^{-1})\right)$. However, in the case of the FAR method, according to (4.4) we derive for $\theta>1$ that

[TABLE]

for sufficient small $\delta>0$. This means that FAR with $1<\theta<2$ yields a super accelerated optimal regularization method. (b) Similar to Theorem 3.2, the converse result for FAR under logarithmic source conditions can be established by the existing technique, proposed in e.g. [1].

5 Numerical investigation of the FAR method

Roughly speaking, the fractional differential equation (1.8) with appropriate numerical discretization schemes for the artificial time variable and stopping rule of iteration steps yields a concrete iterative regularization method. Rather than performing a rigorous numerical analysis, we provide in this section a numerical example to demonstrate if the regularization property and the acceleration effect of the FAR carry over to its numerical realization. To this end, we consider the one-step Adams-Moulton method [2, 3], namely,

[TABLE]

where $b_{j,k+1}=\frac{\Delta t^{\theta}}{\theta}\left[(k-j+1)^{\theta}-(k-j)^{\theta}\right]$, and $a_{j,k+1}=\frac{\Delta t^{\theta}}{\theta(\theta+1)}\cdot d_{j,k+1}$,

[TABLE]

Remark 5.1.

As the Landweber iterates, the iterates $x^{\delta}_{k+1}$ of (5.1) obviously belong to the Krylov subspace $\textrm{Span}\left\{A^{*}y^{\delta},...,(A^{*}A)^{k-1}A^{*}y^{\delta}\right\}$. Therefore, the solution $x^{\delta}_{k+1}$ of (5.1) can be written as $x^{\delta}_{k}=g_{k}(A^{*}A)A^{*}y^{\delta}$, where $g_{k}$ is a polynomial of degree $k$. For proving the regularization property of (5.1), one has to check that $g_{k}(\lambda)$ fulfills all three conditions in Definition 2.2. Moreover, in order to prove the acceleration effect of 2.2 under source conditions (1.4), one must show that for all $\lambda\in(0,\|A\|^{2}]$: $|r_{k}(\lambda)|(k\lambda)^{p}\to 0$ as $k\to\infty$, where $r_{k}(\lambda)=1-\lambda g_{k}(\lambda)$ represents the residual polynomial. However, all of these issues are not addressed here since it is out of the scope of our current paper. Rather we refer to a similar result in [7] for the second order asymptotical regularization.

Now, by employing the newly developed iterative regularization method (5.1), we present some numerical results for the following integral equation

[TABLE]

In this context, we choose $\mathcal{X}=\mathcal{Y}:=L^{2}[0,1]$ such that the operator $A$ is compact, selfadjoint and injective. Then the operator equation $Ax=y$ can be rewritten as $x=-y^{\prime\prime}$ provided that $y\in H^{2}[0,1]\cap H^{1}_{0}[0,1]$. Moreover, the operator $A$ has the singular system (eigensystem) $\{\sigma_{j};u_{j};u_{j}\}_{j=1}^{\infty}$ with $\sigma_{j}=(j\pi)^{-1}$ and $u_{j}(t)=\sqrt{2}\sin(j\pi t)$. Furthermore, using the interpolation theory (see, e.g., [16]) it is not difficult to show that for $4p-1/2\not\in\mathbb{N}$

[TABLE]

In our simulations, problem (5.2) is numerically solved by the linear finite element method. Let $\mathcal{Y}_{n}$ be the finite element space of piecewise linear functions on a uniform grid with step size $1/(n-1)$. Denote by $P_{n}$ the orthogonal projection operator acting from $\mathcal{Y}$ into $\mathcal{Y}_{n}$. Define $A_{n}:=P_{n}A$ and $\mathcal{X}_{n}:=A^{*}_{n}\mathcal{Y}_{n}$. Let $\{\phi_{j}\}^{n}_{j=1}$ be a basis of $\mathcal{Y}_{n}$, then, instead of (5.2), we solve a system of linear equations $A_{n}x_{n}=y_{n}$ in practice, where $[A_{n}]_{ij}=\int^{1}_{0}\left(\int^{1}_{0}k(s,t)\phi_{i}(s)ds\right)\phi_{j}(t)dt$ and $[y_{n}]_{j}=\int^{1}_{0}y(t)\phi_{j}(t)dt$.

We consider the following two different right-hand sides for (5.2).

Example 5.1.

In this example let $y(s):=s(1-s)$. Then the solution is $x^{\dagger}\equiv 2$ and we have $x^{\dagger}\in\mathcal{R}((A^{*}A)^{p})$ for all $p<1/8$.

Example 5.2.

In this example let $y(s):=s^{4}(1-s)^{3}$. Then the solution attains the form $x^{\dagger}=-6t^{2}(1-t)(2-8t+7t^{2})$, and we have $x^{\dagger}\in\mathcal{R}((A^{*}A)^{p})$ for all $p<5/8$.

Both examples are tested by using the grid size $n=100$. Other parameters are: $\Delta t=19.4850,x_{0}=\dot{x}_{0}=0,\tau=3.1$. Uniformly distributed noises with the magnitude $\delta^{\prime}$ are added to the discretized exact right-hand side:

[TABLE]

where $\textrm{Rand}(x)$ returns a pseudo-random value drawn from a uniform distribution on [0,1]. The noise level of measurement data is calculated by $\delta=\|y^{\delta}_{n}-y_{n}\|_{2}$, where $\|\cdot\|_{2}$ denotes the standard vector norm in $\mathbb{R}^{n}$. The iteration of (5.1) is terminated according to the discrepancy principle ($\tau=3$ for all examples), i.e.

[TABLE]

Finally, to assess the accuracy of the approximate solutions, we define the $L^{2}$-norm relative error for an approximate solution $x^{k^{*}}_{n}$: $\textrm{L2Err}:=\|x^{k^{*}}_{n}-x^{\dagger}\|_{L^{2}[0,1]}/\|x^{\dagger}\|_{L^{2}[0,1]}$, where $x^{\dagger}$ is the exact solution to the corresponding model problem.

The results of the simulation are presented in Table 1, where we can conclude that, in general, the FAR method with $\theta\in(1,2)$ needs fewer iterations and offers more accurate regularized solutions. Concerning the number of iterations, the conjugate gradient method for the normal equation (CGNE, cf., e.g. [11]) performed significantly more effectively than all of other methods. However, the accuracy of the CGNE method is considerably worse than other accelerated regularization methods, since the step size of CGNE is too large to capture the optimal point and the semi-convergence effect disturbs the iteration rather early. Note that we set a maximal iteration number $k_{max}=200,000$ in all of our simulations.

6 Conclusion and outlook

In this paper, we have investigated the Fractional Asymptotical Regularization method (FAR) for solving linear ill-posed operator equations $Ax=y$ with compact forward operators $A$ mapping between infinite dimensional Hilbert spaces. Instead of exact right-hand side $y$, we are given noisy data $y^{\delta}$ obeying the deterministic noise model $\|y^{\delta}-y\|\leq\delta$. We have proven that under both Hölder-type and logarithmic source conditions, FAR with fractional order $\theta\in(0,2)$ exhibits an optimal regularization method. Moreover, if $\theta\in(1,2)$, FAR yields an accelerated method, namely, the optimal convergence rates can be obtained with approximately the $\theta$ root of iterations than needed for ordinary Landweber iteration. Finally, with the help of the one-step Adams-Moulton method, a new iterative regularization algorithm has been introduced.

As has been shown by this manuscript, fractional calculus can play some role for regularization schemes aimed at the stable approximate solution to general linear ill-posed problems. However, to the best of our knowledge, the literature in this direction is quite limited. The initial results in this work might be a bridge between fractional calculus and regularization theory. Of course, there are many remaining interesting problems in this topic. For instance, since the Hilbert scale can be understood as preconditioning for iterative regularization methods [4], it seems to be necessary to derive further assertions on acceleration of the FAR method in Hilbert scales. Certainly, the optimal choice of fractional order $\theta$ and extensions to nonlinear problems are of interest. Yet, in the linear setting, new questions also arise, for example the study of the damped system

[TABLE]

where $\theta>1$, and $\eta>0$ is the damping parameter, which may or may not depend on the artificial time $t$. A remaining natural question is whether the damping term in (6.3) can yield a further acceleration of FAR.

Acknowledgements

We express our gratitude to the anonymous referees whose valuable comments and suggestions allowed us to eliminate weak points of the manuscript and thus to improve the paper.

The work of Y. Zhang is supported by the Alexander von Humboldt Foundation through a postdoctoral researcher fellowship, and the work of B. Hofmann is supported by the German Research Foundation (DFG-grant HO 1454/12-1).