Monotonicity of error of regularized solution and its use for parameter choice

Abstract

We consider an ill-posed equation in a Hilbert space with a noisy operator and a noisy right-hand side. The noise level information is given in a general form, as a norm of a certain operator applied to the noise. We derive the monotone error rule (ME-rule) for the choice of the regularization parameter in many methods, giving parameter such that the error is monotonically increasing for larger parameters in the Tikhonov method and for smaller stopping indices in iteration methods. Regularization methods considered include $Y$-scale regularization in (iterated) Tikhonov method and in iteration methods (Landweber method, CG type methods, semi-iterative methods). We also consider modifications of the ME-rule and show in numerical experiments (test problems from Hansen’s Regularization Toolbox, including the sideways heat equation) their advantages over the discrepancy principle.

Keywords:

• ill-posed problems
• regularization methods
• iteration methods
• parameter choice
• monotone error

AMS Subject Classifications:

• 47A52
• 65J20
• 65J10

Introduction

In this paper, we consider linear ill-posed problems(1) $$\begin{array}{c}\hfill A_{0}x=y_0,y_0\in \mathcal{R}(A),\end{array}$$(1) where $A_0$ is a bounded linear operator with non-closed range $\mathcal{R}(A)$ and $X$, $Y$ are infinite-dimensional real Hilbert spaces with inner products $(·,·)$ and norms $\parallel ·\parallel$. We are interested in the minimum-norm solution $x^{†}$ of problem (1(1) $$\begin{array}{c}\hfill A_{0}x=y_0,y_0\in \mathcal{R}(A),\end{array}$$(1) ) and assume that instead of exact data $y_0$ and $A_0$, there are given noisy data $y\in Y$ and $A\in \mathcal{L}(X,Y)$ with(2) $$\begin{array}{c}\hfill \parallel y_0-y\parallel \le \mathit{\delta },\parallel A_0-A\parallel \le h\end{array}$$(2) and known noise levels $\mathit{\delta }$, $h$. Later, we also consider the case of noise level information in more general form (4(4) $$\begin{array}{c}\hfill \parallel D(y_0-y)\parallel \le \mathit{\delta },\parallel D(A_0-A)\parallel \le h.\end{array}$$(4) ).

Ill-posed problems (1(1) $$\begin{array}{c}\hfill A_{0}x=y_0,y_0\in \mathcal{R}(A),\end{array}$$(1) ) arise in a wide variety of problems in applied sciences. For their stable numerical solution, regularization methods are necessary, see [1, 2]. Regularization methods include Tikhonov regularization $x_r={(A^{\ast }A+r^{-1}I)}^{-1}{A}^{\ast }y$ with regularization parameter $r\in {\mathbb{R}}_{+}$ and iterative and projection methods, where the stopping index $r=n\in \mathbb{N}$ is the regularization parameter. Here, $A^{\ast }\in \mathcal{L}(Y,X)$ is the adjoint operator to $A\in \mathcal{L}(X,Y)$. Traditional regularization methods possess the property that in the case of exact data, the error $\parallel x_r^0-x^{†}\parallel$ of regularized solution $x_r^0$, as a function of $r$, is monotonically decreasing for $r\to \infty$. This property is no longer true for the error $\parallel x_r-x^{†}\parallel$ in the case of noisy data. The monotone decrease of the error $\parallel x_r-x^{†}\parallel$ for growing $r$-values can only be guaranteed for small $r$. Typically, $\parallel x_r-x^{†}\parallel$ diverges for $r\to \infty$. Therefore, a rule for the proper choice of the regularization parameter $r$ is necessary.

In the monotone error rule (ME-rule) for choosing a proper regularization parameter, the idea consists in searching for the largest computable regularization parameter $r=r_{\mathrm{ME}}$ for which we can guarantee the ME-property: the error $\parallel x_r-x^{†}\parallel$ is monotonically decreasing for $r\in (0,r_{\mathrm{ME}}]$. For the continuous regularization methods, this means that$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}r}\parallel x_r-x^{†}{\parallel }^2\le 0\forall r\in (0,r_{\mathrm{ME}}],\end{array}$$for the iteration methods and other methods with $r=n\in \mathbb{N}$ this means that(3) $$\begin{array}{c}\hfill \parallel x_n-x^{†}\parallel \le \parallel x_{n-1}-x^{†}\parallel \forall n=1,2,\dots ,n_{\mathrm{ME}}.\end{array}$$(3) From derivation of the ME-rule for concrete regularization methods, one can see for which perturbations of the operator and the right-hand side, the ME-rule gives the optimal regularization parameter.

Consider now the case of the exact operator. In theoretical works on ill-posed problems, often the worst-case error $\sup \{\parallel x_r-x^{†}\parallel :y\in Y,{\mathit{Ax}}^{†}=y_0,\parallel y-y_0\parallel \le \mathit{\delta }\}$ is considered. Here, the exact data $y_0$ are fixed, supremum is over noisy data $y\in Y$. In applications, the exact data $y_0$ are unknown and noisy data $y$ are known. Then, we are interested in finding the parameter $r$ which minimizes the analogue of the worst-case error $\sup \{\parallel x_r-x\parallel :\stackrel{\sim }{y}\in Y,\mathit{Ax}=\stackrel{\sim }{y},\parallel \stackrel{\sim }{y}-y\parallel \le \mathit{\delta }\}$, where given data $y$ are fixed, supremum is over “candidates of exact data” $\stackrel{\sim }{y}\in Y$. This analogue of the worst-case error is minimized by the parameter $r_{\mathrm{ME}}$. This can be seen from derivation of the ME-rule for concrete regularization methods: for $r\in (0,r_{\mathrm{ME}}]$, the error $\parallel x_r-x\parallel$ is decreasing for all “candidates of exact data” $\stackrel{\sim }{y}\in Y$ with $\mathit{Ax}=\stackrel{\sim }{y}$, but for $r>r_{\mathrm{ME}}$ there exists $\stackrel{\sim }{y}\in Y$ with $\mathit{Ax}=\stackrel{\sim }{y}$ where $\parallel x_r-x\parallel >\parallel x_{r_{\mathrm{ME}}}-x\parallel$.

All regularization methods considered in this paper have the property $x_r\in \mathcal{R}(A^{\ast })$ because we need this property for the formulation of the ME rule.

The ME-rule was proposed and studied for continuous methods in [3–8], for iterative methods in [4, 7, 8], analogous rules were proposed in [9, 10]. In this paper, we extend these results in the following directions: we consider the case of noisy operator ($h>0$), the noise level information may be given in general form (4), we derive the ME-rule for a more general class of regularization methods (including the regularization in $Y$-scale,[11] also for iterative and semi-iterative methods). We also give convergence results and error estimates for the ME-rule.

In practical applications, often different regularization parameters are selected and the corresponding regularized solutions are studied online. With our ME-rule, we provide some help for the selection of suitable values of $r$, since we can guarantee that for $r<r_{\mathrm{ME}}$ always $\parallel x_{r_{\mathrm{ME}}}-x^{†}\parallel <\parallel x_r-x^{†}\parallel$ holds. This information may also be used for improving other a posteriori rules R for the choice of the regularization parameter, since the error for the parameter $\max (r_{\mathrm{R}},r_{\text{ME}})$ is less than or equal to the error in rule R. Unfortunately, this observation does not help to improve the discrepancy principle, which gives the parameter $r_{\mathrm{D}}\ge r_{\mathrm{ME}}$, and for a smooth solution $r_{\mathrm{D}}$ is often too large in regularization methods with the finite qualification. In (iterated) Tikhonov method, the ME-rule is in contrast to the discrepancy principle a quasi-optimal parameter choice.

Note that typically the error of the regularized solution decreases monotonically also, somewhat further up to some $r_{\mathrm{opt}}\ge r_{\mathrm{ME}}$. Our numerical experiments suggest to use $r_{\mathrm{MEe}}=c{r}_{\mathrm{ME}}$ (or its integer part) with a certain $c>1$.

The plan of this paper is as follows. In Sections 2 and 3, well-known regularization methods and parameter choice rules are introduced, using noise-level information (4(4) $$\begin{array}{c}\hfill \parallel D(y_0-y)\parallel \le \mathit{\delta },\parallel D(A_0-A)\parallel \le h.\end{array}$$(4) ). In Sections 4 and 5, the ME-rule is derived for the continouous regularization methods and for the iterative methods, respectively. In Section 6 convergence conditions and quasioptimality results of the ME-rule are given. The paper is finished by extensive numerical experiments.

Well-known regularization methods in $Y$-scale

We assume that the noise level information is given in the form(4) $$\begin{array}{c}\hfill \parallel D(y_0-y)\parallel \le \mathit{\delta },\parallel D(A_0-A)\parallel \le h.\end{array}$$(4) where $D$ is a linear injective, possibly unbounded operator in $Y$ with domain $\mathcal{D}(D)$. We assume that $y_0,y\in \mathcal{D}(D)$. The standard case (2(2) $$\begin{array}{c}\hfill \parallel y_0-y\parallel \le \mathit{\delta },\parallel A_0-A\parallel \le h\end{array}$$(2) ) is special case with $D=I$, where $I$ is the identity operator. In [12], the operator $D=L^{-t}$ with $t\in \mathbb{R}$ was used. Here, the operator $L$ is a densely defined, unbounded, self-adjoint and strictly positive operator in the space $Y$, inducing Hilbert scale $Y_s$ with norm $\parallel y{\parallel }_s:=\parallel L^{s}y\parallel$. In the standard case (2(2) $$\begin{array}{c}\hfill \parallel y_0-y\parallel \le \mathit{\delta },\parallel A_0-A\parallel \le h\end{array}$$(2) ), $t=0$. The case $t<0$ corresponds to large noise, the case $t>0$ to small noise. The noise-level information (4(4) $$\begin{array}{c}\hfill \parallel D(y_0-y)\parallel \le \mathit{\delta },\parallel D(A_0-A)\parallel \le h.\end{array}$$(4) ) with different operators $D$ (only the case $h=0$ was considered) was used in works.[13–16]

In many papers, $X_s$-scale is considered with the generating operator $L$ in the space $X$, the equation $\mathit{Ax}=y$ is transformed to the equation $L^{-2s}{A}^{\ast }\mathit{Ax}=L^{-2s}{A}^{\ast }y$ and the regularized solutions are constructed in the form $x_r=g_r(L^{-2s}{A}^{\ast }A)L^{-2s}{A}^{\ast }y$. Here, $r$ is the regularization parameter and the function $g_r(\mathit{\lambda })$ satisfies the conditions(5) $$\begin{array}{cc}\hfill & \underset{0\le \mathit{\lambda }\le \mathrm{\Lambda }}{\sup }|\sqrt{\mathit{\lambda }}{g}_r(\mathit{\lambda })|\le {\mathit{\gamma }}_{\ast }\sqrt{r},r\ge 0,\hfill \\ \hfill & \underset{0\le \mathit{\lambda }\le \mathrm{\Lambda }}{\sup }{\mathit{\lambda }}^p|1-\mathit{\lambda }{g}_r(\mathit{\lambda })|\le {\mathit{\gamma }}_p{r}^{-p},r\ge 0,0\le p\le p_0.\hfill \end{array}$$(5) Here, $p_0$, ${\mathit{\gamma }}_{\ast }$ and ${\mathit{\gamma }}_p$ are positive constants, $\mathrm{\Lambda }$ is at least the norm of the operator in the argument of the function $g_r$, ${\mathit{\gamma }}_0=1$ and the greatest value of $p_0$, for which the inequality (6(6) $$\begin{array}{c}\hfill x_r=g_r(B^{\ast }B)B^{\ast }\overline{y}=B^{\ast }{g}_r({\mathit{BB}}^{\ast })\overline{y}=A^{\ast }{w}_r,w_r=L^{-s}{g}_r({\mathit{BB}}^{\ast })\overline{y}.\end{array}$$(6) ) holds is called the qualification of the method. We can formulate the ME rule for choice of parameter $r$ only for regularized solutions $x_r\in \mathcal{R}(A^{\ast })$. Therefore, we prefer to use the $Y$-scale instead of $X$ scale. The $Y$-scale regularization was proposed by Egger [11, 12]. Here, the equation $\mathit{Ax}=y$ is transformed into the equation$$\begin{array}{c}\hfill \mathit{Bx}=\overline{y},B=L^{-s}A,\overline{y}=L^{-s}y,\end{array}$$assuming $B:X\to Y$ is bounded and $\overline{y}\in Y$. If $s\ge 0$ this is always satisfied; if $s<0$, this means that $\mathcal{R}(A)\in Y_{-s}$ and $y\in Y_{-s}$. The regularized solutions are constructed in the form(6) $$\begin{array}{c}\hfill x_r=g_r(B^{\ast }B)B^{\ast }\overline{y}=B^{\ast }{g}_r({\mathit{BB}}^{\ast })\overline{y}=A^{\ast }{w}_r,w_r=L^{-s}{g}_r({\mathit{BB}}^{\ast })\overline{y}.\end{array}$$(6) Special cases of regularization methods of this form are the following well-known methods (cf. [1, 2, 17]).

1. The Tikhonov method $x_{\mathit{\alpha }}={(\mathit{\alpha }I+B^{\ast }B)}^{-1}{B}^{\ast }\overline{y}$ . Here, $r={\mathit{\alpha }}^{-1}$, $g_r(\mathit{\lambda })={(\mathit{\lambda }+r^{-1})}^{-1}$, $p_0=1$, ${\mathit{\gamma }}_{\ast }=1/2$, ${\mathit{\gamma }}_p=p^p{(1-p)}^{1-p}$.

2. Iterative variants of the Tikhonov method. Let $m\in \mathbb{N}$, $m\ge 1$, $x_{0,\mathit{\alpha }}\in X$ – initial approximation and$$\begin{array}{c}\hfill x_{n,\mathit{\alpha }}={(\mathit{\alpha }I+B^{\ast }B)}^{-1}(\mathit{\alpha }{x}_{n-1,\mathit{\alpha }}+B^{\ast }\overline{y})(n=1,2,\dots ,m).\end{array}$$Here $r={\mathit{\alpha }}^{-1}$, $g_r(\mathit{\lambda })=\frac{1}{\mathit{\lambda }}(1-{\left(\frac{1}{1+r\mathit{\lambda }}\right)}^m)$, $p_0=m$, ${\mathit{\gamma }}_{\ast }=\sqrt{m}$, ${\mathit{\gamma }}_p={\left(\frac{p}{m}\right)}^p{(1-\frac{p}{m})}^{m-p}$.

3. Explicit iteration scheme (Landweber’s method). Let $x_0=0$,$$\begin{array}{c}\hfill x_n=x_{n-1}-\mathit{\mu }(B^{\ast }B{x}_{n-1}-\overline{y}),\mathit{\mu }\in (0,1/\parallel B^{\ast }B\parallel ),n=1,2,\dots \end{array}$$Here, $r=n$, $g_r(\mathit{\lambda })=\frac{1}{\mathit{\lambda }}(1-{(1-\mathit{\mu }\mathit{\lambda })}^r)$, $p_0=\infty$, ${\mathit{\gamma }}_{\ast }=\sqrt{\mathit{\mu }}$, ${\mathit{\gamma }}_p={\left(\frac{p}{\mathit{\mu }e}\right)}^p$.

4. Implicit iteration scheme. Let $\mathit{\alpha }>0$ be a constant, $x_0=0$ and$$\begin{array}{c}\hfill x_n={(\mathit{\alpha }I+B^{\ast }B)}^{-1}(\mathit{\alpha }{x}_{n-1}+B^{\ast }\overline{y}),n=1,2,\dots \end{array}$$Here, $r=n$, $g_r(\mathit{\lambda })=\frac{1}{\mathit{\lambda }}(1-{\left(\frac{\mathit{\alpha }}{\mathit{\alpha }+\mathit{\lambda }}\right)}^r)$, $p_0=\infty$, ${\mathit{\gamma }}_{\ast }=\frac{b_0}{\sqrt{\mathit{\alpha }}}$, where $b_0=\underset{0<\mathit{\lambda }<\infty }{\sup }{\mathit{\lambda }}^{-1/2}(1-e^{-\mathit{\lambda }})\approx 0.6382$ and ${\mathit{\gamma }}_p={\left(\mathit{\alpha }p\right)}^p$.

5. The method of asymptotical regularization: approximation $x_r$ solves the Cauchy problem$$\begin{array}{c}\hfill x^{\prime }(r)+B^{\ast }Bx(r)=B^{\ast }\overline{y},x(0)=0.\end{array}$$Here, $g_r(\mathit{\lambda })=\frac{1}{\mathit{\lambda }}(1-e^{-r\mathit{\lambda }})$, $p_0=\infty$, ${\mathit{\gamma }}_{\ast }=b_0$, ${\mathit{\gamma }}_p={(p/e)}^p$.

In the iterated Tikhonov method, we may use different parameter ${\mathit{\alpha }}_n$ at every iteration step n; then, we get the same approximation as in the extrapolated Tikhonov method [18, 19](7) $$\begin{array}{c}\hfill x_{{\mathit{\alpha }}_1,\dots ,{\mathit{\alpha }}_m}=\sum _{i=1}^m{c}_i{x}_{{\mathit{\alpha }}_i},c_i=\prod _{\frac{j=1}{j\ne i}}^m\frac{{\mathit{\alpha }}_j}{{\mathit{\alpha }}_j-{\mathit{\alpha }}_i},\end{array}$$(7) where $x_{{\mathit{\alpha }}_i}$, $i=1$, ..., $m$ are the approximations of Tikhonov method with the parameters ${\mathit{\alpha }}_i$.

Note that regularized solutions in $X_s$-scales and $Y_s$-scales variants of Tikhonov method have the form(8) $$\begin{array}{c}\hfill x_{\mathit{\alpha }}={(\mathit{\alpha }I+L^{-2s}{A}^{\ast }A)}^{-1}{L}^{-2s}{A}^{\ast }y,x_{\mathit{\alpha }}={(\mathit{\alpha }I+A^{\ast }{L}^{-2s}A)}^{-1}{A}^{\ast }{L}^{-2s}y,\end{array}$$(8) and they minimize the corresponding functionals$$\begin{array}{c}\hfill \parallel \mathit{Ax}-y{\parallel }^2+\mathit{\alpha }\parallel L^{s}x{\parallel }^2,\parallel L^{-s}(\mathit{Ax}-y){\parallel }^2+\mathit{\alpha }\parallel x{\parallel }^2.\end{array}$$If the operator L acts in both spaces X and Y and the operators $A^{\ast }$ and L commute, then these regularized approximations coincide. Much more widespread of these two approximations is the first one (cf.[20]). For the operators L and $L^{-1}$, one may use, for example, the differentiation and integration operators. The fractional powers $L^s$ and $L^{-s}$ can be implemented efficiently, e.g. via FFT or multi-level techniques.

Positive values of $s$ are good for delaying saturation in Tikhonov and iterated Tikhonov methods,[12] negative values of $s$ are good for preconditioning in the iteration methods that drastically decreases the number of iteration steps.[11]

Well-known rules for choice of the regularization parameter

For the choice of the regularization parameter in any regularization method under noise information (4(4) $$\begin{array}{c}\hfill \parallel D(y_0-y)\parallel \le \mathit{\delta },\parallel D(A_0-A)\parallel \le h.\end{array}$$(4) ), some estimate or approximation of value $\parallel D(y-{\mathit{Ax}}^{†})\parallel$ is needed. This value may be estimated as follows:(9) $$\begin{array}{c}\hfill \parallel D(y-{\mathit{Ax}}^{†})\parallel =\parallel D(y-y_0+(A_0-A)x^{†})\parallel \le \mathit{\delta }+h\parallel x^{†}\parallel =:{\mathrm{\Delta }}_1.\end{array}$$(9) Typically, the exact value of $\parallel x^{†}\parallel$ is unknown. Substitution of $\parallel x^{†}\parallel$ by some upper bound $M\ge \parallel x^{†}\parallel$ gives a rough estimate $\parallel D(y-{\mathit{Ax}}^{†})\parallel \le \mathit{\delta }+\mathit{hM}=:{\mathrm{\Delta }}_2.$ If an upper bound $M$ is not available, approximation of $\parallel D(y-{\mathit{Ax}}^{†})\parallel$ by ${\mathrm{\Delta }}_3:=\mathit{\delta }+h\parallel x_{\mathit{\alpha }}\parallel$ (in iteration methods ${\mathrm{\Delta }}_3:=\mathit{\delta }+h\parallel x_n\parallel$) may be used.

In the following, we formulate some a posteriori rules for choosing the regularization parameter which use $\mathrm{\Delta }\in \{{\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2,{\mathrm{\Delta }}_3\}$ and are well known in the case $D=I$, $s=0$. We assume $\mathcal{R}({\mathit{AA}}^{\ast })\subset \mathcal{D}(D)$.

1. Discrepancy principle.[1, 2, 17, 21] In the continuous regularization methods, this principle (D principle) chooses the parameter $\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{D}}$ as the solution of the equation(10) $$\begin{array}{c}\hfill d_{\mathrm{D}}(\mathit{\alpha }):=\parallel D(y-{\mathit{Ax}}_{\mathit{\alpha }})\parallel =C\mathrm{\Delta }\text{with}C\ge 1.\end{array}$$(10) In iteration methods, the discrepancy principle finds $n_{\mathrm{D}}$ as the first index for which $\parallel D(y-{\mathit{Ax}}_n)\parallel \le C\mathrm{\Delta }$.

2. Modified discrepancy principle.[22] In this rule for the methods M1, M2, the parameter $\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{MD}}$ is chosen as the solution of the equation(11) $$\begin{array}{c}\hfill d_{\mathrm{MD}}(\mathit{\alpha }):=\parallel {\mathit{DL}}^s{(I-g_{\mathit{\alpha }}({\mathit{BB}}^{\ast }){\mathit{BB}}^{\ast })}^{1/(2p_0)}{L}^{-s}(y-{\mathit{Ax}}_{\mathit{\alpha }})\parallel =C\mathrm{\Delta }\end{array}$$(11) with $C\ge 1$, assuming $Y_{-s}\subset \mathcal{D}(D)$ (if $D=L^{-t}$ then a sufficient condition is $s\ge t$). Here, $p_0$ is the qualification of the regularization method, see (6(6) $$\begin{array}{c}\hfill x_r=g_r(B^{\ast }B)B^{\ast }\overline{y}=B^{\ast }{g}_r({\mathit{BB}}^{\ast })\overline{y}=A^{\ast }{w}_r,w_r=L^{-s}{g}_r({\mathit{BB}}^{\ast })\overline{y}.\end{array}$$(6) ). One may use also the alternative MD’ rule(12) $$\begin{array}{c}\hfill d_{{\mathrm{MD}}^{\mathrm{\prime }}}(\mathit{\alpha })=(D(y-{\mathit{Ax}}_{\mathit{\alpha }}),D{L}^s{(I-g_{\mathit{\alpha }}({\mathit{BB}}^{\ast }){\mathit{BB}}^{\ast })}^{1/p_0}{L}^{-s}(y-{\mathit{Ax}}_{\mathit{\alpha }}{))}^{1/2}=C\mathrm{\Delta }.\end{array}$$(12)

In case $D=L^{-s}$, the functions $d_{\mathrm{MD}}(\mathit{\alpha })$ and $d_{{\mathrm{MD}}^{\mathrm{\prime }}}(\mathit{\alpha })$ coincide. For regularization methods with $p_0=\infty$, both these forms of the modified discrepancy principle coincide with the discrepancy principle.

ME-rule for the continuous regularization methods

Derivation of the ME-rule

Let us consider the ME-rule for the continuous regularization methods. The prominent example of these methods is Tikhonov method where the regularization parameter is traditionally denoted by $\mathit{\alpha }$. Therefore, in this section we use notation $\mathit{\alpha }=1/r$ instead of r. We assume that the corresponding function $g_{\mathit{\alpha }}(\mathit{\lambda })$ is differentiable with respect to $\mathit{\alpha }$.

Let us consider the regularized approximation(13) $$\begin{array}{c}\hfill x_{\mathit{\alpha }}=A^{\ast }{w}_{\mathit{\alpha }},z_{\mathit{\alpha }}:=\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}{w}_{\mathit{\alpha }}\in Y,\end{array}$$(13) assuming that $z_{\mathit{\alpha }}\in \mathcal{D}({(D^{-1})}^{\ast })$. An example of a regularized approximation of this form is the approximation $x_{\mathit{\alpha }}=A^{\ast }{w}_{\mathit{\alpha }}$ with $w_{\mathit{\alpha }}=L^{-s}{g}_r({\mathit{BB}}^{\ast })\overline{y}$ from (7(7) $$\begin{array}{c}\hfill x_{{\mathit{\alpha }}_1,\dots ,{\mathit{\alpha }}_m}=\sum _{i=1}^m{c}_i{x}_{{\mathit{\alpha }}_i},c_i=\prod _{\frac{j=1}{j\ne i}}^m\frac{{\mathit{\alpha }}_j}{{\mathit{\alpha }}_j-{\mathit{\alpha }}_i},\end{array}$$(7) ). The reformulation of the general idea of ME-rule in terms of $\mathit{\alpha }=1/r$ means: search in the approximation $x_{\mathit{\alpha }}$ for the smallest computable regularization parameter $\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{ME}}$ for which we can guarantee that(14) $$\begin{array}{c}\hfill \frac{\text{d}}{\text{d}\mathit{\alpha }}\parallel x_{\mathit{\alpha }}-x^{†}{\parallel }^2\ge 0\text{for}\text{all}\mathit{\alpha }\in [{\mathit{\alpha }}_{\text{ME}},\infty ).\end{array}$$(14) In order to guarantee this property, we estimate the derivative of the squared error $\parallel x_{\mathit{\alpha }}-x^{†}{\parallel }^2$ with respect to $\mathit{\alpha }$ under condition (10(10) $$\begin{array}{c}\hfill d_{\mathrm{D}}(\mathit{\alpha }):=\parallel D(y-{\mathit{Ax}}_{\mathit{\alpha }})\parallel =C\mathrm{\Delta }\text{with}C\ge 1.\end{array}$$(10) ) as follows:(15) $$\begin{array}{cc}\hfill \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}\parallel x_{\mathit{\alpha }}-x^{†}{\parallel }^2& =(x_{\mathit{\alpha }}-x^{†},A^{\ast }{z}_{\mathit{\alpha }})=(D({\mathit{Ax}}_{\mathit{\alpha }}-y+y-{\mathit{Ax}}^{†}),{(D^{-1})}^{\ast }{z}_{\mathit{\alpha }})\hfill \\ \hfill & \ge \hfill & \hfill \parallel {(D^{-1})}^{\ast }{z}_{\mathit{\alpha }}\parallel (\frac{({\mathit{Ax}}_{\mathit{\alpha }}-y,z_{\mathit{\alpha }})}{\parallel {(D^{-1})}^{\ast }{z}_{\mathit{\alpha }}\parallel }-(\mathit{\delta }+h\parallel x^{†}\parallel )).\end{array}$$(15) This estimate leads us to following ME-rule for the continuous regularization method (14(14) $$\begin{array}{c}\hfill \frac{\text{d}}{\text{d}\mathit{\alpha }}\parallel x_{\mathit{\alpha }}-x^{†}{\parallel }^2\ge 0\text{for}\text{all}\mathit{\alpha }\in [{\mathit{\alpha }}_{\text{ME}},\infty ).\end{array}$$(14) ).

ME rule. Choose $\mathit{\alpha }={\mathit{\alpha }}_{\text{ME}}$ as the largest solution of the equation(16) $$\begin{array}{c}\hfill d_{\text{ME}}(\mathit{\alpha }):=\frac{({\mathit{Ax}}_{\mathit{\alpha }}-y,z_{\mathit{\alpha }})}{\parallel {(D^{-1})}^{\ast }{z}_{\mathit{\alpha }}\parallel }=\mathrm{\Delta }\end{array}$$(16) with $\mathrm{\Delta }\in \{{\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2\}$. If $d_{\text{ME}}(\mathit{\alpha })>\mathrm{\Delta }$ for all $\mathit{\alpha }>0$ take ${\mathit{\alpha }}_{\text{ME}}=0$. If $d_{\text{ME}}(\mathit{\alpha })<\mathrm{\Delta }$ for all $\mathit{\alpha }>0$ take $x_{\mathit{\alpha }}=0$ (it corresponds to $\mathit{\alpha }={\mathit{\alpha }}_{\text{ME}}=\infty$). For regularized approximation (7(7) $$\begin{array}{c}\hfill x_{{\mathit{\alpha }}_1,\dots ,{\mathit{\alpha }}_m}=\sum _{i=1}^m{c}_i{x}_{{\mathit{\alpha }}_i},c_i=\prod _{\frac{j=1}{j\ne i}}^m\frac{{\mathit{\alpha }}_j}{{\mathit{\alpha }}_j-{\mathit{\alpha }}_i},\end{array}$$(7) ), this equation has the form(17) $$\begin{array}{c}\hfill d_{\text{ME}}(\mathit{\alpha }):=\frac{({\mathit{Ax}}_{\mathit{\alpha }}-y,L^{-s}\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}{g}_{\mathit{\alpha }}({\mathit{BB}}^{\ast })L^{-s}y)}{\parallel {(D^{-1})}^{\ast }{L}^{-s}\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}{g}_{\mathit{\alpha }}({\mathit{BB}}^{\ast })L^{-s}y\parallel }=\mathrm{\Delta }.\end{array}$$(17) Here, the assumption $z_{\mathit{\alpha }}\in \mathcal{D}({(D^{-1})}^{\ast })$ is satisfied if $Y_s\subset \mathcal{D}({(D^{-1})}^{\ast })$.

Note that for $\mathrm{\Delta }={\mathrm{\Delta }}_3$, the ME-property is not guaranteed but quasi-optimality retains (see Theorem 6.3). However, $\mathrm{\Delta }={\mathrm{\Delta }}_3$ does not need the norm of the exact solution or a bound of it and in practice may give better results than $\mathrm{\Delta }\in \{{\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2\}$.

From definition of $d_{\text{ME}}(\mathit{\alpha })$ follows the inequality $d_{\text{ME}}(\mathit{\alpha })\le \parallel D(y-{\mathit{Ax}}_{\mathit{\alpha }})\parallel$, therefore ${\mathit{\alpha }}_{\text{D}}\le {\mathit{\alpha }}_{\text{ME}}$ if $C=1$ in the discrepancy principle. The derivation of the ME-rule in (16(16) $$\begin{array}{c}\hfill d_{\text{ME}}(\mathit{\alpha }):=\frac{({\mathit{Ax}}_{\mathit{\alpha }}-y,z_{\mathit{\alpha }})}{\parallel {(D^{-1})}^{\ast }{z}_{\mathit{\alpha }}\parallel }=\mathrm{\Delta }\end{array}$$(16) ) uses only one inequality, which turns to the equality in the case$$\begin{array}{c}\hfill y-A{x}^{†}=-\frac{\mathit{\delta }+h\parallel x^{†}\parallel }{\parallel {(D^{-1})}^{\ast }{z}_{\mathit{\alpha }}\parallel }{(D^{\ast }D)}^{-1}{z}_{\mathit{\alpha }}.\end{array}$$Therefore, in this case the ME-rule gives the optimal parameter, if the equation (18(18) $$\begin{array}{c}\hfill d_{\text{MD}}(\mathit{\alpha })={(D{\mathit{\rho }}_{m,\mathit{\alpha }},D{\mathit{\rho }}_{m+1,\mathit{\alpha }})}^{1/2}\text{and}{d}_{\text{ME}}(\mathit{\alpha })=\frac{({\mathit{\rho }}_{m,\mathit{\alpha }},L^{-2s}{\mathit{\rho }}_{m+1,\mathit{\alpha }})}{\parallel {(D^{-1})}^{\ast }{L}^{-2s}{\mathit{\rho }}_{m+1,\mathit{\alpha }}\parallel }.\end{array}$$(18) ) has a unique solution.

ME-rule and modifications for the (iterated) Tikhonov regularization

Let $x_0=0$. Let ${\mathit{\rho }}_{m,\mathit{\alpha }}$ denote the discrepancy of the (iterated) Tikhonov approximation $x_{\mathit{\alpha }}=x_{m,\mathit{\alpha }}$, i.e. ${\mathit{\rho }}_{m,\mathit{\alpha }}:=y-A{x}_{m,\mathit{\alpha }}$.

Using the identities$$\begin{array}{c}\hfill 1-\mathit{\lambda }{g}_{\mathit{\alpha }}(\mathit{\lambda })={\left(\frac{\mathit{\alpha }}{\mathit{\lambda }+\mathit{\alpha }}\right)}^m\text{and}\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}{g}_{\mathit{\alpha }}(\mathit{\lambda })=-\frac{m}{{\mathit{\alpha }}^2}{\left(\frac{\mathit{\alpha }}{\mathit{\lambda }+\mathit{\alpha }}\right)}^{m+1}\end{array}$$we obtain$$\begin{array}{c}\hfill L^{-s}{\mathit{\rho }}_{m,\mathit{\alpha }}=\overline{y}-{\mathit{Bx}}_{m,\mathit{\alpha }}=[I-{\mathit{BB}}^{\ast }{g}_{\mathit{\alpha }}({\mathit{BB}}^{\ast })]\overline{y}={[\mathit{\alpha }{({\mathit{BB}}^{\ast }+\mathit{\alpha }I)}^{-1}]}^m\overline{y}\end{array}$$and$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}{g}_{\mathit{\alpha }}({\mathit{BB}}^{\ast })\overline{y}=-\frac{m}{{\mathit{\alpha }}^2}{L}^{-s}{\mathit{\rho }}_{m+1,\mathit{\alpha }}.\end{array}$$From these representations and from (18(18) $$\begin{array}{c}\hfill d_{\text{MD}}(\mathit{\alpha })={(D{\mathit{\rho }}_{m,\mathit{\alpha }},D{\mathit{\rho }}_{m+1,\mathit{\alpha }})}^{1/2}\text{and}{d}_{\text{ME}}(\mathit{\alpha })=\frac{({\mathit{\rho }}_{m,\mathit{\alpha }},L^{-2s}{\mathit{\rho }}_{m+1,\mathit{\alpha }})}{\parallel {(D^{-1})}^{\ast }{L}^{-2s}{\mathit{\rho }}_{m+1,\mathit{\alpha }}\parallel }.\end{array}$$(18) ), we conclude that the functions $d_{\text{MD}}(\mathit{\alpha })$ in (13(13) $$\begin{array}{c}\hfill x_{\mathit{\alpha }}=A^{\ast }{w}_{\mathit{\alpha }},z_{\mathit{\alpha }}:=\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}{w}_{\mathit{\alpha }}\in Y,\end{array}$$(13) ) and $d_{\text{ME}}(\mathit{\alpha })$ for the MD’- and ME-rules have the form(18) $$\begin{array}{c}\hfill d_{\text{MD}}(\mathit{\alpha })={(D{\mathit{\rho }}_{m,\mathit{\alpha }},D{\mathit{\rho }}_{m+1,\mathit{\alpha }})}^{1/2}\text{and}{d}_{\text{ME}}(\mathit{\alpha })=\frac{({\mathit{\rho }}_{m,\mathit{\alpha }},L^{-2s}{\mathit{\rho }}_{m+1,\mathit{\alpha }})}{\parallel {(D^{-1})}^{\ast }{L}^{-2s}{\mathit{\rho }}_{m+1,\mathit{\alpha }}\parallel }.\end{array}$$(18) In the case $D=I$, $h=s=0$, convergence and order optimal error estimates for (iterated) Tikhonov method are proved in the case of MD’-rule in [5, 22, 23] and in the case of ME-rule in [3, 5, 6, 8, 24].

One may also use an analog of the ME-rule with the function(19) $$\begin{array}{c}\hfill d_{\text{MEa}}(\mathit{\alpha })=\frac{(D{\mathit{\rho }}_{m,\mathit{\alpha }},D{\mathit{\rho }}_{m+1,\mathit{\alpha }})}{\parallel D{\mathit{\rho }}_{m+1,\mathit{\alpha }}\parallel },\end{array}$$(19) which coincides with ME-rule in the case $D=L^{-s}$.

Consider now some modifications of the ME-rule. We know from (15(15) $$\begin{array}{cc}\hfill \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}\parallel x_{\mathit{\alpha }}-x^{†}{\parallel }^2& =(x_{\mathit{\alpha }}-x^{†},A^{\ast }{z}_{\mathit{\alpha }})=(D({\mathit{Ax}}_{\mathit{\alpha }}-y+y-{\mathit{Ax}}^{†}),{(D^{-1})}^{\ast }{z}_{\mathit{\alpha }})\hfill \\ \hfill & \ge \hfill & \hfill \parallel {(D^{-1})}^{\ast }{z}_{\mathit{\alpha }}\parallel (\frac{({\mathit{Ax}}_{\mathit{\alpha }}-y,z_{\mathit{\alpha }})}{\parallel {(D^{-1})}^{\ast }{z}_{\mathit{\alpha }}\parallel }-(\mathit{\delta }+h\parallel x^{†}\parallel )).\end{array}$$(15) ) that ${\mathit{\alpha }}_{\text{ME}}\ge {\mathit{\alpha }}_{\text{opt}}:=\mathrm{argmin}\{\parallel x_{m,\mathit{\alpha }}-x^{†}\parallel ,\mathit{\alpha }\ge 0\}$. It means that typically somewhat a smaller parameter ${\mathit{\alpha }}_{\text{ME}}/c$ with proper $c>1$ is better parameter than ${\mathit{\alpha }}_{\text{ME}}$. For case $D=I$, $h=s=0$ we optimized the value of $c$ (and other constants below) in extensive numerical experiments and recommend the estimated parameter ${\mathit{\alpha }}_{\text{MEe}}={\mathit{\alpha }}_{\text{ME}}/2.3$. Note that in the case of rough estimate of the noise level, better rules than the ME-rule and the discrepancy principle are rules from the recently derived family of rules,[25] see also [26–29].

In the extrapolated Tikhonov method (8(8) $$\begin{array}{c}\hfill x_{\mathit{\alpha }}={(\mathit{\alpha }I+L^{-2s}{A}^{\ast }A)}^{-1}{L}^{-2s}{A}^{\ast }y,x_{\mathit{\alpha }}={(\mathit{\alpha }I+A^{\ast }{L}^{-2s}A)}^{-1}{A}^{\ast }{L}^{-2s}y,\end{array}$$(8) ) $x_{\mathit{\alpha }}:=x_{{\mathit{\alpha }}_1,\dots ,{\mathit{\alpha }}_m}$, one may use the logarithmically uniform mesh of parameters ${\mathit{\alpha }}_n=\mathit{\alpha }{q}^{n-(m+1)/2}$, $n=1,\dots ,m$; $q<1$ and to choose $\mathit{\alpha }$ by the same rules as in the iterated Tikhonov method.

Note that in [30], an analog of ME-rule in case $D=I$, $h=s=0$ was proposed for the ($m\ge 1$ times iterated) Lavrentiev method, finding the parameter $\mathit{\alpha }={\mathit{\alpha }}_{\text{MEa}}$ as solution of the equation $\parallel {\mathit{\rho }}_{m+1,\mathit{\alpha }}{\parallel }^2/\parallel {\mathit{\rho }}_{m+2,\mathit{\alpha }}\parallel \le C_m\mathit{\delta }$, where $C_1=1.55$, $C_2=1.6$. Numerical experiments in [30, 31] showed that with this parameter choice, the error of regularized solution was in average only 5% larger (in modifications of this rule 3% larger) than with optimal parameter. Note also that several other analogs of ME-rule for the ($m\ge 1$ times iterated) Lavrentiev method were proposed and analysed in [32].

ME-rule for asymptotical regularization

In this regularization method, the regularized solution is given by $x_{\mathit{\alpha }}$ as the solution of the initial value problem$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}x(t)+B^{\ast }Bx(t)=B^{\ast }\overline{y}\text{for}0<t\le r,x(0)=0\end{array}$$with $r=1/\mathit{\alpha }$. In this method, we have $g_{\mathit{\alpha }}(\mathit{\lambda })=(1-e^{-\mathit{\lambda }/\mathit{\alpha }})/\mathit{\lambda }$. Using the identities$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}{g}_{\mathit{\alpha }}(\mathit{\lambda })=-\frac{1}{{\mathit{\alpha }}^2}{e}^{-\mathit{\lambda }/\mathit{\alpha }}=-\frac{1}{{\mathit{\alpha }}^2}(1-\mathit{\lambda }{g}_{\mathit{\alpha }}(\mathit{\lambda }))\end{array}$$we obtain for the method of asymptotical regularization$$\begin{array}{c}\hfill d_{\text{ME}}(\mathit{\alpha })=\frac{({\mathit{Ax}}_{\mathit{\alpha }}-y,L^{-2s}({\mathit{Ax}}_{\mathit{\alpha }}-y))}{\parallel {(D^{-1})}^{\ast }{L}^{-2s}({\mathit{Ax}}_{\mathit{\alpha }}-y)\parallel }.\end{array}$$In the case $s=0$ and $D=I$, the ME-rule coincides with the discrepancy principle for which convergence and order optimal error estimates are well known (cf. [1, 2, 8, 17]).

ME-rule for iterative regularization methods

Derivation of the ME-rule

For solving ill-posed problems (1(1) $$\begin{array}{c}\hfill A_{0}x=y_0,y_0\in \mathcal{R}(A),\end{array}$$(1) ), we now consider iteration methods of the general form(20) $$\begin{array}{c}\hfill x_n=x_{n-1}+A^{\ast }{z}_{n-1},n=1,2,\dots \end{array}$$(20) with $z_n\in Y$. The elements $z_n$ characterize the special iteration method. We assume that $z_n\in \mathcal{D}({(D^{-1})}^{\ast })$. Simple iteration methods are explicit iteration scheme (Landweber method) M3 with $z_n=\mathit{\mu }{L}^{-2s}{\mathit{\rho }}_n$, $\mathit{\mu }\in (0,1/\parallel B^{\ast }B\parallel )$ and implicit iteration scheme M4 with $z_n=L^{-s}{(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-1}{L}^{-s}{\mathit{\rho }}_n={\mathit{\alpha }}^{-1}{L}^{-2s}{\mathit{\rho }}_{n+1}$, $\mathit{\alpha }>0$, where ${\mathit{\rho }}_n=y-{\mathit{Ax}}_n$.

In the monotone error rule (ME-rule), we search for a largest computable iteration number$n_{\text{ME}}$ for which the monotonicity property (3(3) $$\begin{array}{c}\hfill \parallel x_n-x^{†}\parallel \le \parallel x_{n-1}-x^{†}\parallel \forall n=1,2,\dots ,n_{\mathrm{ME}}.\end{array}$$(3) )(21) $$\begin{array}{c}\hfill \parallel x_n-x^{†}\parallel \le \parallel x_{n-1}-x^{†}\parallel \text{for}\text{all}n=1,2,\dots ,n_{\text{ME}}\end{array}$$(21) can be guaranteed. Exploiting (21(21) $$\begin{array}{c}\hfill \parallel x_n-x^{†}\parallel \le \parallel x_{n-1}-x^{†}\parallel \text{for}\text{all}n=1,2,\dots ,n_{\text{ME}}\end{array}$$(21) ) and (10(10) $$\begin{array}{c}\hfill d_{\mathrm{D}}(\mathit{\alpha }):=\parallel D(y-{\mathit{Ax}}_{\mathit{\alpha }})\parallel =C\mathrm{\Delta }\text{with}C\ge 1.\end{array}$$(10) ), we obtain(22) $$\begin{array}{cc}\hfill \parallel x_{n+1}-x^{†}{\parallel }^2-\parallel x_n-x^{†}{\parallel }^2& =\parallel x_n+A^{\ast }{z}_n-x^{†}{\parallel }^2-\parallel x_n-x^{†}{\parallel }^2\hfill \\ \hfill & =2(x_n-x^{†},A^{\ast }{z}_n)+\parallel A^{\ast }{z}_n{\parallel }^2\hfill \\ \hfill & =2(D(A{x}_n-y_0),{(D^{-1})}^{\ast }{z}_n)+\parallel A^{\ast }{z}_n{\parallel }^2\hfill \\ \hfill & =2(D(y-{\mathit{Ax}}^{†}-{\mathit{\rho }}_n),{(D^{-1})}^{\ast }{z}_n)+\parallel A^{\ast }{z}_n{\parallel }^2\hfill \\ \hfill & \le 2\parallel {(D^{-1})}^{\ast }{z}_n\parallel (\mathit{\delta }+h\parallel x^{†}\parallel )-2({\mathit{\rho }}_n,z_n)+\parallel A^{\ast }{z}_n{\parallel }^2\hfill \\ \hfill & =:\mathrm{\Phi }(z_n).\hfill \end{array}$$(22) This estimate leads us to the following ME-rule for iteration methods (21(21) $$\begin{array}{c}\hfill \parallel x_n-x^{†}\parallel \le \parallel x_{n-1}-x^{†}\parallel \text{for}\text{all}n=1,2,\dots ,n_{\text{ME}}\end{array}$$(21) ).

ME rule. Choose $n_{\text{ME}}$ as the first index $n$ satisfying(23) $$\begin{array}{c}\hfill d_{\text{ME}}(n):=\frac{2({\mathit{\rho }}_n,z_n)-\parallel A^{\ast }{z}_n{\parallel }^2}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }\le \mathrm{\Delta }\end{array}$$(23) with $\mathrm{\Delta }\in \{{\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2\}$.

The derivation of the ME-rule in (22(22) $$\begin{array}{cc}\hfill \parallel x_{n+1}-x^{†}{\parallel }^2-\parallel x_n-x^{†}{\parallel }^2& =\parallel x_n+A^{\ast }{z}_n-x^{†}{\parallel }^2-\parallel x_n-x^{†}{\parallel }^2\hfill \\ \hfill & =2(x_n-x^{†},A^{\ast }{z}_n)+\parallel A^{\ast }{z}_n{\parallel }^2\hfill \\ \hfill & =2(D(A{x}_n-y_0),{(D^{-1})}^{\ast }{z}_n)+\parallel A^{\ast }{z}_n{\parallel }^2\hfill \\ \hfill & =2(D(y-{\mathit{Ax}}^{†}-{\mathit{\rho }}_n),{(D^{-1})}^{\ast }{z}_n)+\parallel A^{\ast }{z}_n{\parallel }^2\hfill \\ \hfill & \le 2\parallel {(D^{-1})}^{\ast }{z}_n\parallel (\mathit{\delta }+h\parallel x^{†}\parallel )-2({\mathit{\rho }}_n,z_n)+\parallel A^{\ast }{z}_n{\parallel }^2\hfill \\ \hfill & =:\mathrm{\Phi }(z_n).\hfill \end{array}$$(22) ) uses only one inequality, which turns to the equality in the case$$\begin{array}{c}\hfill y-A{x}^{†}=\frac{\mathit{\delta }+h\parallel x^{†}\parallel }{\parallel {(D^{-1})}^{\ast }{z}_n\parallel }{(D^{\ast }D)}^{-1}{z}_n.\end{array}$$Therefore, in this case the ME-rule gives the optimal parameter.

Note that for $\mathrm{\Delta }={\mathrm{\Delta }}_3$, the ME-property is not guaranteed but for certain methods it is known that quasi-optimality still holds (see Theorem 6.3). However, typically the monotonicity of error also holds up to somewhat larger iteration numbers, hence $\mathrm{\Delta }={\mathrm{\Delta }}_3$ may also give better results than $\mathrm{\Delta }\in \{{\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2\}$.

From definition of $d_{\text{ME}}(n)$ follows the inequality $d_{\text{ME}}(n)\le \parallel D{\mathit{\rho }}_n\parallel =d_{\text{D}}(n)$.

Due to equality ${\mathit{\rho }}_{n+1}={\mathit{\rho }}_n-{\mathit{AA}}^{\ast }{z}_n$, the function $d_{\text{ME}}(n)$ may be presented also in forms(24) $$\begin{array}{c}\hfill d_{\text{ME}}(n):=\frac{({\mathit{\rho }}_n+{\mathit{\rho }}_{n+1},z_n)}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }=\frac{2({\mathit{\rho }}_{n+1},z_n)+\parallel A^{\ast }{z}_n{\parallel }^2}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }.\end{array}$$(24) Inequalities (23(23) $$\begin{array}{c}\hfill d_{\text{ME}}(n):=\frac{2({\mathit{\rho }}_n,z_n)-\parallel A^{\ast }{z}_n{\parallel }^2}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }\le \mathrm{\Delta }\end{array}$$(23) ), (24(24) $$\begin{array}{c}\hfill d_{\text{ME}}(n):=\frac{({\mathit{\rho }}_n+{\mathit{\rho }}_{n+1},z_n)}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }=\frac{2({\mathit{\rho }}_{n+1},z_n)+\parallel A^{\ast }{z}_n{\parallel }^2}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }.\end{array}$$(24) ) may be used:

1. for stopping iterations (21(21) $$\begin{array}{c}\hfill \parallel x_n-x^{†}\parallel \le \parallel x_{n-1}-x^{†}\parallel \text{for}\text{all}n=1,2,\dots ,n_{\text{ME}}\end{array}$$(21) ), if for the first time $d_{\text{ME}}(n)\le \mathrm{\Delta }$;

2. for choosing $z_n$, minimizing the function $\mathrm{\Phi }(z_n)$ in error estimate (23(23) $$\begin{array}{c}\hfill d_{\text{ME}}(n):=\frac{2({\mathit{\rho }}_n,z_n)-\parallel A^{\ast }{z}_n{\parallel }^2}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }\le \mathrm{\Delta }\end{array}$$(23) ) (for example choosing steplength $\mathit{\mu }$ in (26(26) $$\begin{array}{c}\hfill \parallel x_{n+1}-x^{†}{\parallel }^2-\parallel x_n-x^{†}{\parallel }^2\le 2{\mathit{\mu }}_n(\parallel {(D^{\ast })}^{-1}{L}^{-2s}{\mathit{\rho }}_n\parallel \mathrm{\Delta }-\parallel L^{-s}{\mathit{\rho }}_n{\parallel }^2)+{\mathit{\mu }}_n^2\parallel A^{\ast }{L}^{-2s}{\mathit{\rho }}_n{\parallel }^2.\end{array}$$(26) ));

3. for using several iteration methods switching from some faster method (for example CGLS) to a slower method (for example Landweber method): if inequality (23(23) $$\begin{array}{c}\hfill d_{\text{ME}}(n):=\frac{2({\mathit{\rho }}_n,z_n)-\parallel A^{\ast }{z}_n{\parallel }^2}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }\le \mathrm{\Delta }\end{array}$$(23) ) does not guarantee decrease of error in the fast method, the accuracy of this approximation may be further increased by a slower method.

ME-rule in gradient-type methods

Let us consider gradient-type methods of the form (21(21) $$\begin{array}{c}\hfill \parallel x_n-x^{†}\parallel \le \parallel x_{n-1}-x^{†}\parallel \text{for}\text{all}n=1,2,\dots ,n_{\text{ME}}\end{array}$$(21) ) with $z_n={\mathit{\mu }}_n{L}^{-2s}{\mathit{\rho }}_n$, where ${\mathit{\rho }}_n=y-{\mathit{Ax}}_n$ is the discrepancy and ${\mathit{\mu }}_n>0$ is some properly chosen stepsize:(25) $$\begin{array}{c}\hfill x_n=x_{n-1}+{\mathit{\mu }}_{n-1}{A}^{\ast }{L}^{-2s}(y-A{x}_{n-1}),n=1,2,\dots \end{array}$$(25) Special cases of method (26(26) $$\begin{array}{c}\hfill \parallel x_{n+1}-x^{†}{\parallel }^2-\parallel x_n-x^{†}{\parallel }^2\le 2{\mathit{\mu }}_n(\parallel {(D^{\ast })}^{-1}{L}^{-2s}{\mathit{\rho }}_n\parallel \mathrm{\Delta }-\parallel L^{-s}{\mathit{\rho }}_n{\parallel }^2)+{\mathit{\mu }}_n^2\parallel A^{\ast }{L}^{-2s}{\mathit{\rho }}_n{\parallel }^2.\end{array}$$(26) ) are the Landweber method with ${\mathit{\mu }}_n=\mathit{\mu }\in (0,1/\parallel B^{\ast }B\parallel )$, the minimal error method with ${\mathit{\mu }}_n=\frac{\parallel L^{-s}{\mathit{\rho }}_n{\parallel }^2}{\parallel A^{\ast }{L}^{-2s}{\mathit{\rho }}_n{\parallel }^2}$ and the steepest descent method with ${\mathit{\mu }}_n=\frac{\parallel A^{\ast }{L}^{-2s}{\mathit{\rho }}_n{\parallel }^2}{\parallel L^{-s}{\mathit{AA}}^{\ast }{L}^{-2s}{\mathit{\rho }}_n{\parallel }^2}$.

Here, the ME-function (25(25) $$\begin{array}{c}\hfill x_n=x_{n-1}+{\mathit{\mu }}_{n-1}{A}^{\ast }{L}^{-2s}(y-A{x}_{n-1}),n=1,2,\dots \end{array}$$(25) ) has the form$$\begin{array}{c}\hfill d_{\text{ME}}(n)=\frac{({\mathit{\rho }}_n+{\mathit{\rho }}_{n+1},L^{-2s}{\mathit{\rho }}_n)}{2\parallel {(D^{-1})}^{\ast }{L}^{-2s}{\mathit{\rho }}_n\parallel }.\end{array}$$The estimate (23(23) $$\begin{array}{c}\hfill d_{\text{ME}}(n):=\frac{2({\mathit{\rho }}_n,z_n)-\parallel A^{\ast }{z}_n{\parallel }^2}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }\le \mathrm{\Delta }\end{array}$$(23) ) which led us to the ME-rule can also be exploited for finding stepsizes ${\mathit{\mu }}_n>0$ in iteration methods (26(26) $$\begin{array}{c}\hfill \parallel x_{n+1}-x^{†}{\parallel }^2-\parallel x_n-x^{†}{\parallel }^2\le 2{\mathit{\mu }}_n(\parallel {(D^{\ast })}^{-1}{L}^{-2s}{\mathit{\rho }}_n\parallel \mathrm{\Delta }-\parallel L^{-s}{\mathit{\rho }}_n{\parallel }^2)+{\mathit{\mu }}_n^2\parallel A^{\ast }{L}^{-2s}{\mathit{\rho }}_n{\parallel }^2.\end{array}$$(26) ) which guarantee that $x_{n+1}$ is a better approximation for $x^{†}$ than $x_n$. Here, the stepsize ${\mathit{\mu }}_n$ may not only depend on y but also on the noise level $\mathrm{\Delta }$. Exploiting $z_n={\mathit{\mu }}_n{L}^{-2s}{\mathit{\rho }}_n$, the estimate (23(23) $$\begin{array}{c}\hfill d_{\text{ME}}(n):=\frac{2({\mathit{\rho }}_n,z_n)-\parallel A^{\ast }{z}_n{\parallel }^2}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }\le \mathrm{\Delta }\end{array}$$(23) ) obtains form(26) $$\begin{array}{c}\hfill \parallel x_{n+1}-x^{†}{\parallel }^2-\parallel x_n-x^{†}{\parallel }^2\le 2{\mathit{\mu }}_n(\parallel {(D^{\ast })}^{-1}{L}^{-2s}{\mathit{\rho }}_n\parallel \mathrm{\Delta }-\parallel L^{-s}{\mathit{\rho }}_n{\parallel }^2)+{\mathit{\mu }}_n^2\parallel A^{\ast }{L}^{-2s}{\mathit{\rho }}_n{\parallel }^2.\end{array}$$(26) Therefore, $x_{n+1}$ is a better approximation for $x^{†}$ than $x_n$, if$$\begin{array}{c}\hfill \mathrm{\Delta }<\frac{\parallel L^{-s}{\mathit{\rho }}_n{\parallel }^2}{\parallel {(D^{\ast })}^{-1}{L}^{-2s}{\mathit{\rho }}_n\parallel },0<{\mathit{\mu }}_n<\frac{2(\parallel L^{-s}{\mathit{\rho }}_n{\parallel }^2-\parallel {(D^{\ast })}^{-1}{L}^{-2s}{\mathit{\rho }}_n\parallel \mathrm{\Delta })}{\parallel A^{\ast }{L}^{-2s}{\mathit{\rho }}_n{\parallel }^2}.\end{array}$$Minimizing the right-hand side of (27(27) $$\begin{array}{c}\hfill x_n=A^{\ast }{w}_n,n=1,2,\dots ,w_n\in Y\end{array}$$(27) ) with respect to ${\mathit{\mu }}_n$ yields$$\begin{array}{c}\hfill {\mathit{\mu }}_n=\frac{\parallel L^{-s}{\mathit{\rho }}_n{\parallel }^2-\parallel {(D^{\ast })}^{-1}{L}^{-2s}{\mathit{\rho }}_n\parallel \mathrm{\Delta }}{\parallel A^{\ast }{L}^{-2s}{\mathit{\rho }}_n{\parallel }^2}.\end{array}$$Substituting into (27(27) $$\begin{array}{c}\hfill x_n=A^{\ast }{w}_n,n=1,2,\dots ,w_n\in Y\end{array}$$(27) ) shows that for this stepsize, the improvement of the squared error can be estimated by$$\begin{array}{c}\hfill \parallel x_{n+1}-x^{†}{\parallel }^2-\parallel x_n-x^{†}{\parallel }^2\le -\frac{{(\parallel L^{-s}{\mathit{\rho }}_n{\parallel }^2-\parallel (D^{\ast }{)}^{-1}{L}^{-2s}{\mathit{\rho }}_n\parallel \mathrm{\Delta })}^2}{\parallel A^{\ast }{L}^{-2s}{\mathit{\rho }}_n{\parallel }^2}.\end{array}$$In case $D=I$, $h=s=0$ last two relations were given in [8, 10], but convergence results and order optimal error bounds for stopping many gradient methods [2, 9, 10, 17, 33] with the ME-rule were stated in [4, 7–10, 34].

ME-rule in conjugate gradient methods

Some gradient-type methods that do not fit into the class of methods (26(26) $$\begin{array}{c}\hfill \parallel x_{n+1}-x^{†}{\parallel }^2-\parallel x_n-x^{†}{\parallel }^2\le 2{\mathit{\mu }}_n(\parallel {(D^{\ast })}^{-1}{L}^{-2s}{\mathit{\rho }}_n\parallel \mathrm{\Delta }-\parallel L^{-s}{\mathit{\rho }}_n{\parallel }^2)+{\mathit{\mu }}_n^2\parallel A^{\ast }{L}^{-2s}{\mathit{\rho }}_n{\parallel }^2.\end{array}$$(26) ) are conjugate gradient methods (see [9, 10, 33, 36]).

Methods CGLS and CGME are the conjugate gradient method applied to $B^{\ast }\mathit{Bx}=B^{\ast }\overline{y}$ and ${\mathit{BB}}^{\ast }v=\overline{y}$ with $x=B^{\ast }v$, respectively. In both methods, we fix the initial approximation $x_0$, the initial value $u_0=0$ and ${\mathit{\beta }}_0=0$, and compute ${\mathit{\tau }}_0=\overline{y}-{\mathit{Bx}}_0.$ Further, we compute in CGLS $p_0=B^{\ast }{\mathit{\tau }}_0$, and for $n=1,2,\dots$ iteratively$$\begin{array}{cc}\hfill u_n& ={\mathit{\beta }}_{n-1}{u}_{n-1}-{\mathit{\tau }}_{n-1},{\mathit{\gamma }}_n=\frac{\parallel p_{n-1}{\parallel }^2}{\parallel {\mathit{BB}}^{\ast }{u}_n{\parallel }^2},x_n=x_{n-1}-{\mathit{\gamma }}_n{B}^{\ast }{u}_n,\hfill \\ \hfill {\mathit{\tau }}_n& ={\mathit{\tau }}_{n-1}+{\mathit{\gamma }}_n{\mathit{BB}}^{\ast }{u}_n,p_n=B^{\ast }{\mathit{\tau }}_n,{\mathit{\beta }}_n=\frac{\parallel p_n{\parallel }^2}{\parallel p_{n-1}{\parallel }^2}.\hfill \end{array}$$In the method CGME, one computes for $n=1,2,\dots$ iteratively$$\begin{array}{cc}\hfill u_n& ={\mathit{\beta }}_{n-1}{u}_{n-1}-{\mathit{\tau }}_{n-1},{\mathit{\gamma }}_n=\frac{\parallel {\mathit{\tau }}_{n-1}{\parallel }^2}{\parallel B^{\ast }{u}_n{\parallel }^2},\hfill \\ \hfill x_n& =x_{n-1}-{\mathit{\gamma }}_n{B}^{\ast }{u}_n,{\mathit{\tau }}_n={\mathit{\tau }}_{n-1}+{\mathit{\gamma }}_n{\mathit{BB}}^{\ast }{u}_n,{\mathit{\beta }}_n=\frac{\parallel {\mathit{\tau }}_n{\parallel }^2}{\parallel {\mathit{\tau }}_{n-1}{\parallel }^2}\hfill \end{array}$$Both CGLS and CGME method have the form (21(21) $$\begin{array}{c}\hfill \parallel x_n-x^{†}\parallel \le \parallel x_{n-1}-x^{†}\parallel \text{for}\text{all}n=1,2,\dots ,n_{\text{ME}}\end{array}$$(21) ) with $z_{n-1}=-{\mathit{\gamma }}_n{L}^{-s}{u}_n$, hence ME-rule (24(24) $$\begin{array}{c}\hfill d_{\text{ME}}(n):=\frac{({\mathit{\rho }}_n+{\mathit{\rho }}_{n+1},z_n)}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }=\frac{2({\mathit{\rho }}_{n+1},z_n)+\parallel A^{\ast }{z}_n{\parallel }^2}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }.\end{array}$$(24) ) can be used. In case $D=I$, $h=s=0$, the ME-rule for CGLS and CGME methods was proposed and numerically tested in [38].

ME-rule for sequence of approximations

In contrast to iterative methods where approximation on step $n$ is computed using approximation on previous step $n-1$, one may consider an arbitrary sequence of approximations(27) $$\begin{array}{c}\hfill x_n=A^{\ast }{w}_n,n=1,2,\dots ,w_n\in Y\end{array}$$(27) and ask which from these approximations to take as approximate solution.

Theorem 5.1

If in sequence (28(28) $$\begin{array}{c}\hfill d_{\text{ME}}(n)=\frac{({\mathit{\rho }}_n+{\mathit{\rho }}_{n+1},w_{n+1}-w_n)}{2\parallel ({(D^{-1})}^{\ast }(w_{n+1}-w_n)\parallel }\le \mathrm{\Delta }\end{array}$$(28) ) the index $n=n_{\text{ME}}$ is chosen as the first index $n$ satisfying(28) $$\begin{array}{c}\hfill d_{\text{ME}}(n)=\frac{({\mathit{\rho }}_n+{\mathit{\rho }}_{n+1},w_{n+1}-w_n)}{2\parallel ({(D^{-1})}^{\ast }(w_{n+1}-w_n)\parallel }\le \mathrm{\Delta }\end{array}$$(28) with $\mathrm{\Delta }\in \{{\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2\}$, then ME-property (3(3) $$\begin{array}{c}\hfill \parallel x_n-x^{†}\parallel \le \parallel x_{n-1}-x^{†}\parallel \forall n=1,2,\dots ,n_{\mathrm{ME}}.\end{array}$$(3) ) holds.

Proof

This follows from (25(25) $$\begin{array}{c}\hfill x_n=x_{n-1}+{\mathit{\mu }}_{n-1}{A}^{\ast }{L}^{-2s}(y-A{x}_{n-1}),n=1,2,\dots \end{array}$$(25) ) with $z_n=w_{n+1}-w_n$. $\square$

To use the functional $d_{\text{ME}}(n)$, elements $w_n$ are needed. They can be found by first computing $w_n$ and later $x_n=A^{\ast }{w}_n$.

Note that the last theorem also gives an alternative way to derive the ME-rule for continuous regularization methods, using the sequence ${\mathit{\alpha }}_n=q^n$, $q<1$ and taking the limit $q\to 1$.

The last theorem may be applied as an alternative or a complement to any other parameter choice rule for finding a regularized approximation from the sequence in the form $x_n=A^{\ast }{w}_n$, $n=1,2,\dots$ (for example in (iterated) Tikhonov method). Note that computation of the sequence of all approximate solutions for $n=1,2,\dots ,N$ with some $N$ is required, for example, in the balancing principle (see [11, 35, 37, 39]. This theorem guarantees that the stopping index $n_{\text{ME},\text{B}}=\max (n_{\text{ME}},n_{\text{B}})$ is at least as good as the index $n_B$ from arbitrary balancing principle: $\parallel x_{n_{\text{ME},\text{B}}}-x^{†}\parallel \le \parallel x_{n_{\text{B}}}-x^{†}\parallel$.

ME-rule for semi-iterative regularization methods

Let $p_n$ be residual polynoms orthogonal in interval $[0,\parallel B^{\ast }B\parallel ]$ with respect to some measure (see [40]). Then $p_0=1$,$$\begin{array}{c}\hfill p_n(\mathit{\lambda })=p_{n-1}(\mathit{\lambda })+(1-{\mathit{\mu }}_n)(p_{n-2}(\mathit{\lambda })-p_{n-1}(\mathit{\lambda }))-{\mathit{\omega }}_n\mathit{\lambda }{p}_{n-1}(\mathit{\lambda }),n\ge 2\end{array}$$with some constants ${\mathit{\mu }}_n$ and ${\mathit{\omega }}_n$, where ${\mathit{\mu }}_1=1$. Then, semi-iterative approximations are found as$$\begin{array}{cc}\hfill & x_1=x_0+{\mathit{\omega }}_1{A}^{\ast }{L}^{-2s}(y-{\mathit{Ax}}_0),\hfill \\ \hfill & x_n=x_{n-1}+(1-{\mathit{\mu }}_n)(x_{n-2}-x_{n-1})+{\mathit{\omega }}_n{A}^{\ast }{L}^{-2s}(y-{\mathit{Ax}}_{n-1}),n\ge 2.\hfill \end{array}$$In case $x_0=A^{\ast }{w}_0$, we have $x_n=A^{\ast }{w}_n$ with $w_1=w_0+{\mathit{\omega }}_1{L}^{-2s}(y-{\mathit{AA}}^{\ast }{w}_0)$,(29) $$\begin{array}{c}\hfill w_n=w_{n-1}+(1-{\mathit{\mu }}_n)(w_{n-2}-w_{n-1})+{\mathit{\omega }}_n{L}^{-2s}(y-{\mathit{AA}}^{\ast }{w}_{n-1}),n\ge 2.\end{array}$$(29) Thus, the iterations of the semi-iterative methods may be stopped by the ME-rule (29(29) $$\begin{array}{c}\hfill w_n=w_{n-1}+(1-{\mathit{\mu }}_n)(w_{n-2}-w_{n-1})+{\mathit{\omega }}_n{L}^{-2s}(y-{\mathit{AA}}^{\ast }{w}_{n-1}),n\ge 2.\end{array}$$(29) ). We give formula for $w_n$ for the following examples of semi-iterative methods (see [40]).

1. The Chebyshev method of Stiefel$$\begin{array}{c}\hfill w_n=\frac{2n}{n+1}{w}_{n-1}-\frac{n-1}{n+1}{w}_{n-2}+\frac{4n}{n+1}{L}^{-2s}(y-{\mathit{AA}}^{\ast }{w}_{n-1}),n=1,2,\dots \end{array}$$

2. The Chebyshev method of Nemirovskii and Polyak$$\begin{array}{c}\hfill w_1=\frac{4}{3}{L}^{-2s}(y-{\mathit{AA}}^{\ast }{w}_0),\end{array}$$$$\begin{array}{cc}\hfill w_n& =2\frac{2n-1}{2n+1}{w}_{n-1}-\frac{2n-3}{2n+1}{w}_{n-2}+4\frac{2n-1}{2n+1}{L}^{-2s}(y-{\mathit{AA}}^{\ast }{w}_{n-1}),\hfill \\ \hfill & n=2,3,\dots \hfill \end{array}$$

3. The $\mathit{\nu }$-method of Brakhage is method where $w_n$ in (30(30) $$\begin{array}{cc}\hfill x_{n+1}& =x_n+A^{\ast }{L}^{-s}{({\mathit{BB}}^{\ast }+{\mathit{\alpha }}_{n}I)}^{-1}{L}^{-s}(y-{\mathit{Ax}}_n)\hfill \\ \hfill & ={(B^{\ast }B+{\mathit{\alpha }}_{n}I)}^{-1}({\mathit{\alpha }}_{n-1}{x}_n+B^{\ast }{L}^{-s}y),n=0,1,\dots \hfill \end{array}$$(30) ) has for fixed $\mathit{\nu }>0$ the coefficients$$\begin{array}{cc}\hfill {\mathit{\mu }}_1& =1,{\mathit{\mu }}_n=1+\frac{(n-1)(2n-3)(2n+2\mathit{\nu }-1)}{(n+2\mathit{\nu }-1)(2n+4\mathit{\nu }-1)(2n+2\mathit{\nu }-3)},n=2,3,\dots \hfill \\ \hfill {\mathit{\omega }}_1& =\frac{4\mathit{\nu }+2}{4\mathit{\nu }+1},{\mathit{\omega }}_n=\frac{4(2n+2\mathit{\nu }-1)(n+\mathit{\nu }-1)}{(n+2\mathit{\nu }-1)(2n+4\mathit{\nu }-1)},n=2,3,\dots \hfill \end{array}$$

ME-rule in implicit iteration methods

Let us consider implicit iteration methods of the form(30) $$\begin{array}{cc}\hfill x_{n+1}& =x_n+A^{\ast }{L}^{-s}{({\mathit{BB}}^{\ast }+{\mathit{\alpha }}_{n}I)}^{-1}{L}^{-s}(y-{\mathit{Ax}}_n)\hfill \\ \hfill & ={(B^{\ast }B+{\mathit{\alpha }}_{n}I)}^{-1}({\mathit{\alpha }}_{n-1}{x}_n+B^{\ast }{L}^{-s}y),n=0,1,\dots \hfill \end{array}$$(30) This method has form (21(21) $$\begin{array}{c}\hfill \parallel x_n-x^{†}\parallel \le \parallel x_{n-1}-x^{†}\parallel \text{for}\text{all}n=1,2,\dots ,n_{\text{ME}}\end{array}$$(21) ) with(31) $$\begin{array}{c}\hfill z_n=L^{-s}{({\mathit{\alpha }}_{n}I+{\mathit{BB}}^{\ast })}^{-1}{L}^{-s}{\mathit{\rho }}_n={\mathit{\alpha }}_n^{-1}{L}^{-2s}{\mathit{\rho }}_{n+1}.\end{array}$$(31) Therefore, the function $d_{\text{ME}}(n)$ in (25(25) $$\begin{array}{c}\hfill x_n=x_{n-1}+{\mathit{\mu }}_{n-1}{A}^{\ast }{L}^{-2s}(y-A{x}_{n-1}),n=1,2,\dots \end{array}$$(25) ) attains the form$$\begin{array}{c}\hfill d_{\text{ME}}(n)=\frac{({\mathit{\rho }}_n+{\mathit{\rho }}_{n+1},L^{-2s}{\mathit{\rho }}_{n+1})}{\parallel {(D^{-1})}^{\ast }{L}^{-2s}{\mathit{\rho }}_{n+1}\parallel }.\end{array}$$In case $D=I$, $h=s=0$, ME-rule for method (31(31) $$\begin{array}{c}\hfill z_n=L^{-s}{({\mathit{\alpha }}_{n}I+{\mathit{BB}}^{\ast })}^{-1}{L}^{-s}{\mathit{\rho }}_n={\mathit{\alpha }}_n^{-1}{L}^{-2s}{\mathit{\rho }}_{n+1}.\end{array}$$(31) ) was studied in [4, 7, 8, 34].

Let us now consider the question of choosing the parameter ${\mathit{\alpha }}_n$ such that the guaranteed improvement of accuracy of approximation $x_{n+1}$ over $x_n$ would be maximal.

Proposition 5.2

In regularization method (31(31) $$\begin{array}{c}\hfill z_n=L^{-s}{({\mathit{\alpha }}_{n}I+{\mathit{BB}}^{\ast })}^{-1}{L}^{-s}{\mathit{\rho }}_n={\mathit{\alpha }}_n^{-1}{L}^{-2s}{\mathit{\rho }}_{n+1}.\end{array}$$(31) ), the guaranteed improvement $-\mathrm{\Phi }(z_n)$ of accuracy in estimate (23(23) $$\begin{array}{c}\hfill d_{\text{ME}}(n):=\frac{2({\mathit{\rho }}_n,z_n)-\parallel A^{\ast }{z}_n{\parallel }^2}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }\le \mathrm{\Delta }\end{array}$$(23) ) is maximized for $z_n$ in (32(32) $$\begin{array}{c}\hfill \parallel x_r-x^{†}\parallel \le \parallel x_r-x_r^0\parallel +\parallel x_r^0-x^{†}\parallel .\end{array}$$(32) ) with ${\mathit{\alpha }}_n$ as solution $\mathit{\alpha }$ of the equation$$\mathit{\alpha }\parallel {(D^{-1})}^{\ast }{L}^{-s}{(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-1}{L}^{-s}{\mathit{\rho }}_n\parallel \parallel {(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-3/2}{L}^{-s}{\mathit{\rho }}_n{\parallel }^2=({(D^{-1})}^{\ast }{L}^{-s}{(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-1}{L}^{-s}{\mathit{\rho }}_n,{(D^{-1})}^{\ast }{L}^{-s}{(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-2}{L}^{-s}{\mathit{\rho }}_n){\mathrm{\Delta }}_1.$$

Proof

Substitution of $z_n$ in (32(32) $$\begin{array}{c}\hfill \parallel x_r-x^{†}\parallel \le \parallel x_r-x_r^0\parallel +\parallel x_r^0-x^{†}\parallel .\end{array}$$(32) ) to function $\mathrm{\Phi }(z_n)$ in (23(23) $$\begin{array}{c}\hfill d_{\text{ME}}(n):=\frac{2({\mathit{\rho }}_n,z_n)-\parallel A^{\ast }{z}_n{\parallel }^2}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }\le \mathrm{\Delta }\end{array}$$(23) ) shows that minimization of function $\mathrm{\Phi }(z_n)$ is equivalent to minimization of the function$$\begin{array}{cc}\hfill \mathit{\varphi }(\mathit{\alpha })& :=2\parallel {(D^{-1})}^{\ast }{L}^{-s}{(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-1}{\mathit{\tau }}_n\parallel {\mathrm{\Delta }}_1-2({\mathit{\tau }}_n,{(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-1}{\mathit{\tau }}_n)\hfill \\ \hfill & +\parallel B^{\ast }{(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-1}{\mathit{\tau }}_n{\parallel }^2,\hfill \end{array}$$where ${\mathit{\tau }}_n:=L^{-s}{\mathit{\rho }}_n$. Differentiation with use of formula $\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}\parallel ·\parallel =(\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}\parallel ·{\parallel }^2)/(2\parallel ·\parallel )$ gives$$\frac{1}{2}{\mathit{\varphi }}^{\prime }(\mathit{\alpha })=-\frac{({(D^{-1})}^{\ast }{L}^{-s}{(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-1}{\mathit{\tau }}_n,{(D^{-1})}^{\ast }{L}^{-s}{(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-2}{\mathit{\tau }}_n)}{\parallel {(D^{-1})}^{\ast }{L}^{-s}{(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-1}{\mathit{\tau }}_n\parallel }{\mathrm{\Delta }}_1+\parallel {(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-1}{\mathit{\tau }}_n{\parallel }^2-\parallel B^{\ast }{(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-3/2}{\mathit{\tau }}_n{\parallel }^2.$$The equality$$\begin{array}{cc}\hfill \parallel {(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-1}{\mathit{\tau }}_n{\parallel }^2& =((\mathit{\alpha }I+{\mathit{BB}}^{\ast }){\mathit{\tau }}_n,{(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-3}{\mathit{\tau }}_n)\hfill \\ \hfill & =\mathit{\alpha }\parallel {(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-3/2}{\mathit{\tau }}_n{\parallel }^2+\parallel B^{\ast }{(\mathit{\alpha }I+{\mathit{BB}}^{\ast })}^{-3/2}{\mathit{\tau }}_n{\parallel }^2\hfill \end{array}$$shows that ${\mathit{\varphi }}^{\prime }(\mathit{\alpha })=0$ iff $\mathit{\alpha }$ solves the equation (33(33) $$\begin{array}{c}\hfill x_r=B^{\ast }{g}_r({\mathit{BB}}^{\ast })\overline{y}=g_r(B^{\ast }B)B^{\ast }\overline{y}\end{array}$$(33) ). $\square$

Theorem 5.3

In the case $D=L^{-s}$ in regularization method (31(31) $$\begin{array}{c}\hfill z_n=L^{-s}{({\mathit{\alpha }}_{n}I+{\mathit{BB}}^{\ast })}^{-1}{L}^{-s}{\mathit{\rho }}_n={\mathit{\alpha }}_n^{-1}{L}^{-2s}{\mathit{\rho }}_{n+1}.\end{array}$$(31) ), the guaranteed improvement $-\mathrm{\Phi }(z_n)$ of accuracy in estimate (23(23) $$\begin{array}{c}\hfill d_{\text{ME}}(n):=\frac{2({\mathit{\rho }}_n,z_n)-\parallel A^{\ast }{z}_n{\parallel }^2}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }\le \mathrm{\Delta }\end{array}$$(23) ) is maximized, if ${\mathit{\alpha }}_n$ is chosen by the analogue of the discrepancy principle $\parallel L^{-s}{\mathit{\rho }}_{n+1}\parallel ={\mathrm{\Delta }}_1.$

Proof

Indeed, in the case $D=L^{-s}$, the equation (33(33) $$\begin{array}{c}\hfill x_r=B^{\ast }{g}_r({\mathit{BB}}^{\ast })\overline{y}=g_r(B^{\ast }B)B^{\ast }\overline{y}\end{array}$$(33) ) for finding $\mathit{\alpha }=r_n^{-1}$ has the form ${\mathit{\alpha }}_n\parallel {({\mathit{\alpha }}_{n}I+{\mathit{BB}}^{\ast })}^{-1}{L}^{-s}{\mathit{\rho }}_n\parallel ={\mathrm{\Delta }}_1$, which is equivalent to $\parallel L^{-s}{\mathit{\rho }}_{n+1}\parallel ={\mathrm{\Delta }}_1$ due to the equality ${\mathit{\alpha }}_n{({\mathit{\alpha }}_{n}I+{\mathit{BB}}^{\ast })}^{-1}{L}^{-s}{\mathit{\rho }}_n=L^{-s}{\mathit{\rho }}_{n+1}$. $\square$

Convergence and quasi-optimality for the ME-rule

In the ME-rule, considered in previous sections, the regularization parameter $r$ in continuous regularization methods is found as the solution (the largest solution in case of many solutions) of a certain equation $d(r,y,A)=\mathrm{\Delta }$ with $\mathrm{\Delta }\in \{{\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2\}$. In regularization methods, where the regularization parameter is a natural number as in iteration methods, we choose $r_{\text{ME}}$ as the first $r=1$, $2$, ... for which $d(r,y,A)\le \mathrm{\Delta }$. Let us denote by $x_r^0$ the regularized solutions for exact data $y_0$, $A_0$. We give for the ME-rule, the following convergence result.

Theorem 6.1

Let $\mathrm{\Delta }\in \{{\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2\}$. Let $x_r^0$, $x_r$ satisfy the following properties:

1. $\parallel x_r^0-x^{†}\parallel \to 0$ as $r\to \infty$;

2. $\parallel x_{r_1}^0-x_{r_1}\parallel \le \parallel x_{r_2}^0-x_{r_2}\parallel$ (for $r_1\le r_2$);

3. if $r(\mathrm{\Delta })\le \text{const}<\infty$ as $\mathrm{\Delta }\to 0$ , then $\parallel x_{r(\mathrm{\Delta })}-x_{r(\mathrm{\Delta })}^0\parallel \to 0$ as $\mathrm{\Delta }\to 0$;

4. if for some sequence $r_k\le \text{const}<\infty$ ($k\to \infty$) it holds true $d(r_k,y_0,A_0)\to 0$ ($k\to \infty$), then $\parallel x_{r_k}^0-x^{†}\parallel \to 0$ ($k\to \infty$).

Then, $\parallel x_{r_{\text{ME}}}-x^{†}\parallel \to 0$ as $\mathrm{\Delta }\to 0$.

Proof

For every $r$ it holds(32) $$\begin{array}{c}\hfill \parallel x_r-x^{†}\parallel \le \parallel x_r-x_r^0\parallel +\parallel x_r^0-x^{†}\parallel .\end{array}$$(32) First we consider the main case where $r_{\text{ME}}(\mathrm{\Delta })\to \infty$ as $\mathrm{\Delta }\to 0$. Then, $\parallel x_{r_{\text{ME}}}^0-x^{†}\parallel \to 0$ due to the assumption (1). It remains to show that $\parallel x_{r_{\text{ME}}}^0-x_{r_{\text{ME}}}\parallel \to 0$ as $\mathrm{\Delta }\to 0$. Let $r_0(\mathrm{\Delta })$ be a monotonically increasing function giving some a priori regularization parameter guaranteeing convergence $\parallel x_{r_0}-x^{†}\parallel \to 0$ as $\mathrm{\Delta }\to 0$. If $r_{\text{ME}}\le r_0$, then due to assumption (2)$$\begin{array}{c}\hfill \parallel x_{r_{\text{ME}}}^0-x_{r_{\text{ME}}}\parallel \le \parallel x_{r_0}^0-x_{r_0}\parallel .\end{array}$$If $r_{\text{ME}}\ge r_0$, then ME-property says that$$\begin{array}{c}\hfill \parallel x_{r_{\text{ME}}}-x^{†}\parallel \le \parallel x_{r_0}-x^{†}\parallel .\end{array}$$Both estimates converge to 0 as $\mathrm{\Delta }\to 0$.

Consider now the case if $r_{\text{ME}}(\mathrm{\Delta })\le \text{const}$ ($\mathrm{\Delta }\to 0$). Then, for $r=r_{\text{ME}}$ both terms in the right-hand side of (34(34) $$\begin{array}{c}\hfill \parallel x_{r(R)}-x^{†}\parallel \le C\underset{r\ge 0}{\inf }\mathit{\psi }(r)+\mathcal{O}(\mathit{\delta }+h\parallel x^{†}\parallel ),\end{array}$$(34) ) converge as $\mathrm{\Delta }\to 0$ due to assumptions (3) and (4). $\square$

Assumptions (1)–(3) of Theorem 6.1 are general concerning regularization method. Only the last assumption (4) uses concrete form of the function $d(r,y,A)$, for regularization method (7(7) $$\begin{array}{c}\hfill x_{{\mathit{\alpha }}_1,\dots ,{\mathit{\alpha }}_m}=\sum _{i=1}^m{c}_i{x}_{{\mathit{\alpha }}_i},c_i=\prod _{\frac{j=1}{j\ne i}}^m\frac{{\mathit{\alpha }}_j}{{\mathit{\alpha }}_j-{\mathit{\alpha }}_i},\end{array}$$(7) ) it can be proved similarly to Lemma 3.2 in [2, p. 66].

Consider now the case $D=L^{-s}$. The problem $\mathit{Ax}=y$ may be rewritten in form $\mathit{Bx}=\overline{y}$ with operator $B\in L(X,Y)$ and for $x_r$ in (7(7) $$\begin{array}{c}\hfill x_{{\mathit{\alpha }}_1,\dots ,{\mathit{\alpha }}_m}=\sum _{i=1}^m{c}_i{x}_{{\mathit{\alpha }}_i},c_i=\prod _{\frac{j=1}{j\ne i}}^m\frac{{\mathit{\alpha }}_j}{{\mathit{\alpha }}_j-{\mathit{\alpha }}_i},\end{array}$$(7) )(33) $$\begin{array}{c}\hfill x_r=B^{\ast }{g}_r({\mathit{BB}}^{\ast })\overline{y}=g_r(B^{\ast }B)B^{\ast }\overline{y}\end{array}$$(33) the standard regularization theory applies. We get the following results.

Theorem 6.2

Let (4(4) $$\begin{array}{c}\hfill \parallel D(y_0-y)\parallel \le \mathit{\delta },\parallel D(A_0-A)\parallel \le h.\end{array}$$(4) ) hold with $D=L^{-s}$. Let $x_r$ be defined by one of methods M1–M5. Then, choice of $r$ by one of rules D, MD, ME or MEe guarantees convergence $\parallel x_r-x^{†}\parallel \to 0$ as $\mathrm{\Delta }\to 0$ and under additional assumption $x^{†}\in \mathcal{R}({(B^{\ast }B)}^{p/2})$ the order optimal error estimate $\parallel x_r-x^{†}\parallel \le O({\mathrm{\Delta }}^{\frac{p}{p+1}})$ holds with $p\le 2p_0-1$ for rule D, with $p\le 2p_0$ for rules MD, ME and MEe.

proof

Let $D=L^{-s}$. Then, results of this theorem for discrepancy principle follow from the results in [1, 2, 17], for MD rule in [22, 41]. Let $C=1$ in rules D, MD. In methods M1, M2 it holds $d_{\text{MD}}(\mathit{\alpha })\le d_{\text{ME}}(\mathit{\alpha })\le d_{\text{D}}(\mathit{\alpha })$, therefore ${\mathit{\alpha }}_{\text{D}}\le {\mathit{\alpha }}_{\text{ME}}\le {\mathit{\alpha }}_{\text{MD}}$ ([5–8]), in methods M3, M4 $n_{\text{D}}-1\le n_{\text{ME}}\le n_{\text{D}}$ ([4, 7, 8]), in method M5 $d_{\text{ME}}(\mathit{\alpha })=d_{\text{D}}(\mathit{\alpha })$. Therefore, results for ME rule follow from results for D rule and MD rule. The same results hold for parameter $r_{\text{MEe}}=\text{const}{r}_{\text{ME}}$. $\square$

Let us consider the quasi-optimality results for the regularization methods. To characterize the quality of the rule of an a posteriori regularization parameter choice, we use in the following the quasi-optimality property (see [41, 42]). We say that a rule R for a posteriori choice of the regularization parameter $r=r(R)$ is quasi-optimal for a given regularization method $x_r=A^{\ast }{g}_r({\mathit{AA}}^{\ast })y$, if there exists a constant $C$ (not depending on $A_0$, $A$, $x^{†}$, $y$) such that for $\parallel y-y_0\parallel \le \mathit{\delta }$, $\parallel A-A_0\parallel \le h$ the error estimate(34) $$\begin{array}{c}\hfill \parallel x_{r(R)}-x^{†}\parallel \le C\underset{r\ge 0}{\inf }\mathit{\psi }(r)+\mathcal{O}(\mathit{\delta }+h\parallel x^{†}\parallel ),\end{array}$$(34) holds, where the function(35) $$\begin{array}{c}\hfill \mathit{\psi }(r)=\parallel (I-A^{\ast }A{g}_r(A^{\ast }A))x^{†})\parallel +{\mathit{\gamma }}_{\ast }\sqrt{r}(\mathit{\delta }+h\parallel x^{†}\parallel )\end{array}$$(35) is an upper bound of the error $\parallel x_r-x^{†}\parallel$. The following result holds.

Theorem 6.3

Let $D=L^{-s}$. The modified discrepancy principle with $\mathrm{\Delta }\in \{{\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2,{\mathrm{\Delta }}_3\}$ and $C>1$ is quasi-optimal rule for solving the problem $\mathit{Bx}=\overline{y}$ by methods M1–M5: the inequality (36(35) $$\begin{array}{c}\hfill \mathit{\psi }(r)=\parallel (I-A^{\ast }A{g}_r(A^{\ast }A))x^{†})\parallel +{\mathit{\gamma }}_{\ast }\sqrt{r}(\mathit{\delta }+h\parallel x^{†}\parallel )\end{array}$$(35) ) holds where the operators $A$, $A^{\ast }$ are replaced by operators $B$, $B^{\ast }$ in (37). This quasi-optimality is also guaranteed for ME-rule with $\mathrm{\Delta }\in \{{\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2\}$, but also with $\mathrm{\Delta }={\mathrm{\Delta }}_3$ if in (17(17) $$\begin{array}{c}\hfill d_{\text{ME}}(\mathit{\alpha }):=\frac{({\mathit{Ax}}_{\mathit{\alpha }}-y,L^{-s}\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}{g}_{\mathit{\alpha }}({\mathit{BB}}^{\ast })L^{-s}y)}{\parallel {(D^{-1})}^{\ast }{L}^{-s}\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}{g}_{\mathit{\alpha }}({\mathit{BB}}^{\ast })L^{-s}y\parallel }=\mathrm{\Delta }.\end{array}$$(17) ) and (24(24) $$\begin{array}{c}\hfill d_{\text{ME}}(n):=\frac{({\mathit{\rho }}_n+{\mathit{\rho }}_{n+1},z_n)}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }=\frac{2({\mathit{\rho }}_{n+1},z_n)+\parallel A^{\ast }{z}_n{\parallel }^2}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }.\end{array}$$(24) ) $\mathrm{\Delta }$ is replaced by $C\mathrm{\Delta }$ with $C>1$.

Proof

Applying results of [41] (Theorem 3) (see also [34](Theorem 2, Remark 2)) to the approximation (35(35) $$\begin{array}{c}\hfill \mathit{\psi }(r)=\parallel (I-A^{\ast }A{g}_r(A^{\ast }A))x^{†})\parallel +{\mathit{\gamma }}_{\ast }\sqrt{r}(\mathit{\delta }+h\parallel x^{†}\parallel )\end{array}$$(35) ) instead of $x_r=A^{\ast }{g}_r({\mathit{AA}}^{\ast })y$, we get the assertions of theorem about MD rule, but also for ME rule with $\mathrm{\Delta }\in \{{\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2,{\mathrm{\Delta }}_3\}$ if in (17(17) $$\begin{array}{c}\hfill d_{\text{ME}}(\mathit{\alpha }):=\frac{({\mathit{Ax}}_{\mathit{\alpha }}-y,L^{-s}\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}{g}_{\mathit{\alpha }}({\mathit{BB}}^{\ast })L^{-s}y)}{\parallel {(D^{-1})}^{\ast }{L}^{-s}\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}{g}_{\mathit{\alpha }}({\mathit{BB}}^{\ast })L^{-s}y\parallel }=\mathrm{\Delta }.\end{array}$$(17) ) and (24(24) $$\begin{array}{c}\hfill d_{\text{ME}}(n):=\frac{({\mathit{\rho }}_n+{\mathit{\rho }}_{n+1},z_n)}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }=\frac{2({\mathit{\rho }}_{n+1},z_n)+\parallel A^{\ast }{z}_n{\parallel }^2}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }.\end{array}$$(24) ) $C\mathrm{\Delta }$ with $C>1$ is used. In case $\mathrm{\Delta }\in \{{\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2\}$, the increase of $C\ge 1$ in right-hand side $C\mathrm{\Delta }$ in (17(17) $$\begin{array}{c}\hfill d_{\text{ME}}(\mathit{\alpha }):=\frac{({\mathit{Ax}}_{\mathit{\alpha }}-y,L^{-s}\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}{g}_{\mathit{\alpha }}({\mathit{BB}}^{\ast })L^{-s}y)}{\parallel {(D^{-1})}^{\ast }{L}^{-s}\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}{g}_{\mathit{\alpha }}({\mathit{BB}}^{\ast })L^{-s}y\parallel }=\mathrm{\Delta }.\end{array}$$(17) ) and (24(24) $$\begin{array}{c}\hfill d_{\text{ME}}(n):=\frac{({\mathit{\rho }}_n+{\mathit{\rho }}_{n+1},z_n)}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }=\frac{2({\mathit{\rho }}_{n+1},z_n)+\parallel A^{\ast }{z}_n{\parallel }^2}{2\parallel {(D^{-1})}^{\ast }{z}_n\parallel }.\end{array}$$(24) ) leads to decrease of parameter $r\le r_{\text{ME}}$, therefore, to increase of the error $\parallel x_r-x^{†}\parallel$. $\square$

Numerical examples

Our tests are performed on the well-known set of test problems by Hansen [43]: baart, deriv2, foxgood, gravity, heat, ilaplace, phillips, shaw, spikes, wing. In all tests, we used the discretization parameter $n=100$.

For example, the test problem heat (see [43, 44]) represents the heat conduction in one-dimensional semi-infinite rod $[0,\infty )$ with the thermal diffusivity 1. The initial temperature of the rod is 0 and the boundary condition at the point $\mathit{\xi }=0$ is not known. Instead, one measures the temperature at $\mathit{\xi }=1$ for the time $t\in [0,1]$ and tries to determine the temperature at $\mathit{\xi }=0$ from this measurement. (NB! In this paragraph and in the following paragraph, $t$ and $s$ have different meaning from the rest of the paper.)

This results in the Volterra integral equation of the first kind$$\begin{array}{c}\hfill {\int }_0^t\frac{{(t-s)}^{-3/2}}{2\sqrt{\mathit{\pi }}}\exp (-\frac{1}{4(t-s)})x(s)\mathit{ds}=y(t),t\in [0,1],\end{array}$$where $y(t)$ is the measured temperature at $\mathit{\xi }=1$ and $x(t)$ is the temperature at $\mathit{\xi }=0$.

The integral equation is discretized by means of a simple collocation method on the points $t_i=i/n$ and the midpoint rule with $n=100$ points. An exact discrete solution $z$ with components$$\begin{array}{c}\hfill z_i=\{\begin{array}{cc}75t_i^2,\hfill & \text{if}{t}_i\le 0.1\hfill \\ 0.75+(20t_i-2)(3-20t_i),\hfill & \text{if}0.1\le t_i<0.15\hfill \\ 0.75\exp (-2(20t_i-3)),\hfill & \text{if}0.15<t_i\le 0.5\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}\end{array}$$is constructed, and then the right-hand side $b$ is computed as $b=\mathit{Az}$.

Since the performance of methods and rules generally depends on the smoothness $p$ of the exact solution, we complemented the standard solutions $x^{†}$ of the (now discrete) test problems with smoothened solutions ${(A^{\ast }A)}^{p/2}{x}^{†}$ ($p=2$) in all test problems (computing the right-hand side as $A\left({(A^{\ast }A)}^{p/2}{x}^{†}\right)$). After discretization, all problems were scaled (normalized) in such a way that the Euclidian norms of the operator and the right-hand side were 1. We used operator $D=L^{-t}$ with $t\in \{-0.5,0,0.5\}$. We computed $L^{-t}{y}_0$ and added a normally distributed noise such that $\parallel L^{-t}(y-y_0)\parallel$ had the values $0.3$, ${10}^{-1}$, ${10}^{-2}$, ${10}^{-3}$, ${10}^{-4}$, ${10}^{-5}$, ${10}^{-6}$. Here, $\parallel ·\parallel$ means the Euclidean norm.

We generated 10 noise vectors and used these vectors in all of the problems. The problems were regularized by the Tikhonov method in $Y_s$-scales (the second formula in (9(9) $$\begin{array}{c}\hfill \parallel D(y-{\mathit{Ax}}^{†})\parallel =\parallel D(y-y_0+(A_0-A)x^{†})\parallel \le \mathit{\delta }+h\parallel x^{†}\parallel =:{\mathrm{\Delta }}_1.\end{array}$$(9) )), where the regularization parameters were chosen according to the rules that we wanted to compare. We used the $Y_s$-scale with the operator L acting as $\mathit{Lu}=u-u^{\prime \prime }$. In the discrete form, we used the approximating matrix $(l_{\mathit{ij}})$ with the elements $l_{\mathit{ii}}=2n^2+1$, $i=1,\dots ,n$, $l_{i,i+1}=l_{i+1,i}=-n^2$, $i=1,\dots ,n-1$, $l_{\mathit{ij}}=0$ otherwise.

Since in these model equations the exact solution is known, it is possible to find the regularization parameter $\mathit{\alpha }={\mathit{\alpha }}_{\ast }$ which gives the smallest error: $\parallel x_{{\mathit{\alpha }}_{\ast }}-x^{†}\parallel ={\min }_{\mathit{\alpha }>0}\{\parallel x_{\mathit{\alpha }}-x^{†}\parallel \}$. For every rule R, the error ratio $\parallel x_{{\mathit{\alpha }}_{\text{R}}}-x^{†}\parallel /\parallel x_{{\mathit{\alpha }}_{\ast }}-x^{†}\parallel$ describes the performance of the rule R on this particular problem. For better comparison of the cases with different s, in the denominator we always use $s=0$. To compare the rules or to present their properties, the following tables show averages of these error ratios over various parameters of the data set (problems 1–10, noise levels $\mathit{\delta }$, noise vectors).

In Tables 1 and 2 we compare the following rules. Earlier we have introduced the discrepancy principle D (see page 5, formula (11(11) $$\begin{array}{c}\hfill d_{\mathrm{MD}}(\mathit{\alpha }):=\parallel {\mathit{DL}}^s{(I-g_{\mathit{\alpha }}({\mathit{BB}}^{\ast }){\mathit{BB}}^{\ast })}^{1/(2p_0)}{L}^{-s}(y-{\mathit{Ax}}_{\mathit{\alpha }})\parallel =C\mathrm{\Delta }\end{array}$$(11) )), the rule MD’ (page 7, (19(19) $$\begin{array}{c}\hfill d_{\text{MEa}}(\mathit{\alpha })=\frac{(D{\mathit{\rho }}_{m,\mathit{\alpha }},D{\mathit{\rho }}_{m+1,\mathit{\alpha }})}{\parallel D{\mathit{\rho }}_{m+1,\mathit{\alpha }}\parallel },\end{array}$$(19) )), the monotone error rule ME (page 7, (19(19) $$\begin{array}{c}\hfill d_{\text{MEa}}(\mathit{\alpha })=\frac{(D{\mathit{\rho }}_{m,\mathit{\alpha }},D{\mathit{\rho }}_{m+1,\mathit{\alpha }})}{\parallel D{\mathit{\rho }}_{m+1,\mathit{\alpha }}\parallel },\end{array}$$(19) )) and the rule MEa as the analog of the ME-rule (page 7, (20(20) $$\begin{array}{c}\hfill x_n=x_{n-1}+A^{\ast }{z}_{n-1},n=1,2,\dots \end{array}$$(20) )). In rule MEe (recall motivation on page 7), we chose the regularization parameter as ${\mathit{\alpha }}_{\text{MEe}}={\mathit{\alpha }}_{\text{ME}}/c$ with $c=2.3·{10}^{c_1}$, $c_1=(s-t)(1.5+1.3s-0.5t-0.8s^2+2\mathit{st}-t^2)$; this form of the constant was obtained by extensive numerical experiments on different values of $s$ and $t$. In rule MEae, we chose the parameter ${\mathit{\alpha }}_{\text{MEae}}={\mathit{\alpha }}_{\text{MEa}}/c$, where $c=2.3$ for $t=0$ and $c=3$ for $t=-0.5$ or $t=0.5$.

Table 1 Error ratios for different rules.

$t$	$s$	$p$	$\mathit{\delta }$	D	MD’	ME	ME08	ME07	MEe	MEa	MEae
$-0.5$	-0.5	0	l	1.20	1.38	1.28			1.15	1.28	1.17
$-0.5$	-0.5	0	s	1.44	2.41	1.48			1.23	1.48	1.19
$-0.5$	-0.5	2	l	4.10	4.90	4.35			3.73	4.35	3.90
$-0.5$	-0.5	2	s	16.2	18.1	14.9			10.4	14.9	10.2
$-0.5$	0	0	l	3e+2	2e+2	1.97	1.80	1.70	1.26	1e+2	1e+2
$-0.5$	0	0	s	2.84	2.72	2.06	1.86	15.5	1.31	1.89	1.89
$-0.5$	0	2	l	4e+4	2e+4	3.56	2.92	2.58	1.27	2e+4	2e+4
$-0.5$	0	2	s	5e+3	2e+3	2.29	2.03	2.47	1.54	1e+3	1e+3
0	0	0	l	1.15	1.37	1.28			1.13	1.28	1.13
0	0	0	s	1.22	1.51	1.38			1.19	1.38	1.19
0	0	2	l	1.24	1.53	1.29			1.09	1.29	1.09
0	0	2	s	3.66	1.91	1.57			1.11	1.57	1.11
0	0.5	0	l	1.61	1.48	2.13	1.93	1.98	1.32	1.34	1.26
0	0.5	0	s	1.34	1.55	2.00	1.82	3.30	1.36	1.46	1.37
0	0.5	2	l	1.54	1.68	3.20	2.61	2.31	1.19	1.52	1.24
0	0.5	2	s	3.31	3.78	6.97	5.96	5.89	2.95	3.55	3.33
0.5	0.5	0	l	0.92	1.02	0.98			0.90	0.98	0.89
0.5	0.5	0	s	1.18	1.35	1.29			1.15	1.29	1.12
0.5	0.5	2	l	0.61	0.77	0.66			0.60	0.66	0.66
0.5	0.5	2	s	1.18	1.66	1.47			1.13	1.47	1.10
0.5	1	0	l	0.99	1.02	1.56	1.46	1.40	1.14	1.04	0.90
0.5	1	0	s	1.24	1.37	1.77	1.64	2.08	1.22	1.32	1.20
0.5	1	2	l	0.54	0.76	2.37	1.98	1.76	1.31	0.66	0.58
0.5	1	2	s	1.56	1.90	3.22	2.85	3.00	1.70	1.76	1.56

Note that in the case $s>t$, the ME rule (unlike rules D, MD’, MEa) allows mild underestimation of the noise level. In the columns ME08 and ME07, we present the error ratios obtained by using the noise levels $0.8\mathit{\delta }$ and $0.7\mathit{\delta }$ instead of $\mathit{\delta }$. Typically, ME08 and ME07 give better results than ME rule, but sometimes they can give large errors. If instead of the exact noise level approximate noise level is given, we recommend using the MEe-rule with $s>t$. In the case $s=t$, all considered rules failed if the noise level was underestimated.

In Table 1, we present the results separately for large $\mathit{\delta }$’s (l means $\mathit{\delta }\ge {10}^{-2}$) and small $\mathit{\delta }$’s (s means $\mathit{\delta }\le {10}^{-3}$). We present results for $s=t$ and $s=t+0.5$. For a given $t$ in most cases, the best choice of $s$ is $s=t$, in which case all of the considered rules gave similar results. But in the case $s=t=-0.5$, wecould not obtain good results for smooth solutions, especially in case of small noise levels. Then, the error function had a very sharp decrease at its minimum, which was hard to locate. Therefore, in the case $t=-0.5$ we recommend taking $s=0$. Then, the ME and MEe rules gave good results. However, other rules failed in this case in the problems baart, foxgood, and wing.

In most problems, using other values of $s$ and $t$ than $s=t$ and $s=t+0.5$ caused slightly worse results in the rules ME and MEe and essentially worse results in rules D and MD’. However, in the problems deriv2 and phillips in the case of a smooth solution ($p\ge 2$), increasing $s$ increased the accuracy of rules D, MD’, MEe for all $t$ (in problem deriv2 the accuracy of rule D increased for $s$ up to 0.5).

Table 1 shows that the discrepancy principle works well in the case $s=t=0$, $s=t=0.5$ but is not as good for $p>1$. Typically, the MEe-rule gives the best results. We made computations besides MD’ rule (13(13) $$\begin{array}{c}\hfill x_{\mathit{\alpha }}=A^{\ast }{w}_{\mathit{\alpha }},z_{\mathit{\alpha }}:=\frac{\mathrm{d}}{\mathrm{d}\mathit{\alpha }}{w}_{\mathit{\alpha }}\in Y,\end{array}$$(13) ) also with MD rule (12(12) $$\begin{array}{c}\hfill d_{{\mathrm{MD}}^{\mathrm{\prime }}}(\mathit{\alpha })=(D(y-{\mathit{Ax}}_{\mathit{\alpha }}),D{L}^s{(I-g_{\mathit{\alpha }}({\mathit{BB}}^{\ast }){\mathit{BB}}^{\ast })}^{1/p_0}{L}^{-s}(y-{\mathit{Ax}}_{\mathit{\alpha }}{))}^{1/2}=C\mathrm{\Delta }.\end{array}$$(12) ), the results were almost the same.

Table 2 presents the results for the problem heat. In column ‘min’, we present averages of minimal errors over 10 runs for concrete $s$ and $t$. Unlike Table 1, the rules D, MD’, MEa, MEae gave almost as good results as the rule MEe. These rules were not as good as the rule MEe only in the case of smooth solution ($p=2$): rule D for $t=s=0$ and rules MD’, MEa, MEae for $t=-0.5$, $s=0$.

Table 2 Results for the problem heat.

t	s	p	$\mathit{\delta }$	min	D	MD’	ME	ME08	ME07	MEe	MEa	MEae
$-0.5$	-0.5	0	${10}^{-1}$	1.08e–1	1.21	1.58	1.35			1.19	1.35	1.27
$-0.5$	-0.5	0	${10}^{-4}$	9.35e–3	1.05	1.25	1.16			1.05	1.16	1.04
$-0.5$	-0.5	2	${10}^{-1}$	5.56e–2	5.05	5.96	5.12			4.88	5.12	5.37
$-0.5$	-0.5	2	${10}^{-4}$	9.85e–4	24.3	19.4	16.3			13.5	16.3	14.0
$-0.5$	0	0	${10}^{-1}$	9.34e–2	1.14	1.49	1.58	1.34	1.22	1.12	1.33	1.33
$-0.5$	0	0	${10}^{-4}$	9.48e–3	1.12	1.25	1.32	1.26	1.23	1.03	1.18	1.18
$-0.5$	0	2	${10}^{-1}$	1.28e–2	1.04	3.25	2.79	2.07	1.74	1.19	3.10	3.10
$-0.5$	0	2	${10}^{-4}$	7.45e–5	1.29	3.38	1.45	1.27	1.17	1.71	2.77	2.77
0	0	0	${10}^{-1}$	5.97e–1	1.07	1.24	1.15			1.02	1.15	1.02
0	0	0	${10}^{-4}$	3.34e–2	1.00	1.26	1.14			1.01	1.14	1.01
0	0	2	${10}^{-1}$	8.67e–2	1.15	1.22	1.07			1.14	1.07	1.14
0	0	2	${10}^{-4}$	1.15e–3	2.50	1.40	1.18			1.04	1.18	1.04
0	0.5	0	${10}^{-1}$	5.71e–1	1.11	1.24	1.89	1.66	1.51	1.15	1.17	1.04
0	0.5	0	${10}^{-4}$	3.25e–2	1.04	1.27	1.17	1.06	1.17	1.21	1.16	1.01
0	0.5	2	${10}^{-1}$	4.93e–2	0.62	0.81	1.69	1.47	1.35	0.57	0.72	0.60
0	0.5	2	${10}^{-4}$	9.89e–4	0.89	1.19	1.33	1.05	8.52	1.23	1.05	0.88
0.5	0.5	0	${10}^{-1}$	1.18	0.90	0.93	0.91			0.87	0.91	0.87
0.5	0.5	0	${10}^{-4}$	1.64e–1	0.95	1.18	1.08			0.90	1.08	0.88
0.5	0.5	2	${10}^{-1}$	1.41e–1	0.34	0.39	0.36			0.35	0.36	0.37
0.5	0.5	2	${10}^{-4}$	7.21e–3	0.61	0.79	0.71			0.61	0.71	0.63
0.5	1	0	${10}^{-1}$	1.15	0.92	0.94	1.14	1.13	1.11	0.98	0.93	0.89
0.5	1	0	${10}^{-4}$	1.55e–1	1.03	1.20	1.30	1.10	7.04	0.92	1.13	0.92
0.5	1	2	${10}^{-1}$	1.15e–1	0.28	0.36	0.84	0.67	0.59	0.31	0.33	0.29
0.5	1	2	${10}^{-4}$	6.79e–3	0.63	0.80	1.21	1.00	0.86	0.67	0.72	0.58

Both in Tables 1 and 2 in case $t=0.5$, $s=1$, the best results were obtained by using the MEae rule.

Some remarks about numerical realization of algorithms. In our numerical computations, the discretization parameter $n$ was 100. This low dimension allowed us to use the advantages of singular value decomposition of the matrix for fast computations. Note also that instead of the solution of the Equations (11(11) $$\begin{array}{c}\hfill d_{\mathrm{MD}}(\mathit{\alpha }):=\parallel {\mathit{DL}}^s{(I-g_{\mathit{\alpha }}({\mathit{BB}}^{\ast }){\mathit{BB}}^{\ast })}^{1/(2p_0)}{L}^{-s}(y-{\mathit{Ax}}_{\mathit{\alpha }})\parallel =C\mathrm{\Delta }\end{array}$$(11) ), (19(19) $$\begin{array}{c}\hfill d_{\text{MEa}}(\mathit{\alpha })=\frac{(D{\mathit{\rho }}_{m,\mathit{\alpha }},D{\mathit{\rho }}_{m+1,\mathit{\alpha }})}{\parallel D{\mathit{\rho }}_{m+1,\mathit{\alpha }}\parallel },\end{array}$$(19) ) and (20(20) $$\begin{array}{c}\hfill x_n=x_{n-1}+A^{\ast }{z}_{n-1},n=1,2,\dots \end{array}$$(20) ) we made computations on the $\mathit{\alpha }$-sequence ${\mathit{\alpha }}_i=0.9^{-i}$, $i=0,1,\dots$ and used for the corresponding regularization parameter the first ${\mathit{\alpha }}_i$ for which the left hand side of the equation was smaller than the right-hand side.

The computations are made with the exact operator, since the inexact operator is typically accompanied by a significant overestimation of the noise level. In this case, it is preferable to use the rules from the article [25] where the exact operator was considered. We plan to extend these results to the case of a noisy operator in a forthcoming paper.

It is worthwhile to think about the possibility to choose the operator L depending on the operator A.

Note that the ME-rule can also be used in projection methods. It was studied for the exact operator in [45]. Generalizations to noisy operators are to be published in a further paper.

Conclusion

We have considered a linear ill-posed problem $A_{0}x=y_0$ with the operator $A_0\in \mathcal{L}(X,Y)$. A noisy right-hand side $y$ and an operator $A$ with noise levels $\parallel D(y_0-y)\parallel \le \mathit{\delta }$, $\parallel D(A_0-A)\parallel \le h$, where $D$ is certain operator, are given. We derived the monotone error rule (ME rule) for parameter choice in many regularization methods in $Y_{-s}$-scale, generated by powers $s\in \mathbb{R}$ of certain operator $L$. These regularization methods include the (iterated) Tikhonov method and many iterative methods (Landweber method, steepest descent method, conjugate gradient-type methods, semi-iterative methods, implicit iteration method). It was shown for which data the ME rule gives the optimal parameter. For a class of methods, the convergence and quasi-optimality results are given. In extensive numerical experiments, the ME-rule and its modification MEe gave good results, whereby in the case $D=L^{-t}$ with $t<s$ these two rules allowed also moderate underestimation of the noise level (other rules fail in case of underestimated noise level).

Acknowledgments

This work was started while the first author held an appointment as Visiting Professor at the University of Applied Sciences Zittau/Görlitz from June until August 2010 and was supported by DFG (Deutsche Forschungsgemeinschaft) under a cooperative bilateral research project grant. In addition, the work was supported by the Estonian Science Foundation, Grant 9120.