Sometimes this projection provides enough regularization to produce a good approximate solution, but often additional regularization is needed. Regularization methods replace an ill-posed problem by a family of well-posed problems, their solution, called regularized solutions, are used as approximations to the desired solution of the inverse problem. These methods always involve some parameter measuring the closeness of the regularized and the original (unregularized) inverse problem, rules (and algorithms) for the choice of these regularization parameters as well as convergence properties of the regularized solutions are central points in the theory of these methods, since they allow to find the right balance between stability and accuracy. The general principles of regularization for ill-posed problems are known. In particular, such principles have been established by A.N. Tikhonov. The literature on various regularization methods based on these general principles is extensive ( Engl, Hanke, and Neubauer [9] , Gilyazov and Gol'dman [11] ).

In statistic, regularization, is associated to penalty based methods or thresholding methods or more generality to “smoothing” techniques. In applications, regularization offers a unifying perspective for many diverse ill-posed inverse problems, a wide range of problems concerned with recovering information from indirect and usually noisy measurements, arising in geophysics, tomography and econometrics. One of the most important, but still insufficiently developed, topics in the theory of ill posed problems is connected with iteration regularization [11] ; i.e, with the utilization of iteration methods of any form for the stable approximate solution of ill-posed problems. Iterative regularization methods tend to be more attractive in terms of numerical cost and implementation, but a number of open questions remain in their theoretical analysis.

In this article we propose an iterative regularized estimator for linear ill-posed problems.

Necessary conditions for convergence are established. These conditions connect the choice of the regularization parameter (i.e., the iteration index) with the projection dimension.

Moreover, we prove that the iterative regularized estimator is optimal in the sense that the estimator achieves the best rate of convergence among all the regularized estimators. A recent work in this direction is developed by Loubes and Luden͂a in [16] , which discusses the problem of estimating inverse nonlinear ill-posed problems with different types of complexity penalties leading either to a model selection estimator or to a regularized estimator.

The choice of the regularization sequence is here crucial, and a lot of work associated with the selection of a good regularization parameter can be found in the literature [10] , [13] . When using iterative methods the problem is finding a good stopping criteria for terminating the iteration procedure. In this article we will use tools developed in the context of model selection via penalization, [2] ,[1] , based on the use of concentration inequalities.

Our article is organized as follows. Section 2 presents basic assumptions and a statement of the discretized inverse problem. In section 3 we discuss regularization methods and we prove consistency of the estimator when the regularization parameter is known. In section 4 we present our main result, prove optimality of an adaptive regularized estimator and give its rate of convergence. Finally, in the last section we introduce regularization by iterative methods for the solution of inverse problems and provide some examples to explain the properties of iterative regularization methods.

2 Preliminaries

2.1 Formulation of the problem and basic assumptions

We assume that the inverse problem is given by $$y=A\tilde{f}+\varepsilon .$$  where $\varepsilon$  is a centered random variable satisfying the moment condition $IE\left(\right|\varepsilon {|}^p/{\sigma }^p)\le p!/2$  and $IE\left({\varepsilon }^2\right)={\sigma }^2$  .

We also need some notations concerning the fixed design settings, $x_i,i=1,\dots ,n$  . Define the empirical measure:

$$P_n=\frac{1}{n}{\sum }_{i=1}^n{\delta }_{x_i}.$$  and the associated empirical norm $$\parallel y{\parallel }_n^2=\parallel y{\parallel }_{P_n}^2=\frac{1}{n}{\sum }_{i=1}^n\left(y\right(x_i){)}^2$$  as well as the empirical scalar product $$<y,\varepsilon {>}_n=\frac{1}{n}{\sum }_{i=1}^n{\varepsilon }_{i}y\left(t_i\right).$$  We assume $A:F\to Y$  , $F,Y$  separable Hilbert spaces. Let $\langle ,\rangle$  stands for the inner product defined over $F$  . We assume that the range of the operator $A$  , $\mathcal{ℛ}\left(A\right),$  is closed in the sense that $\tilde{f}\in \mathcal{N}(A{)}^{\perp }$  , where $\mathcal{N}(A{)}^{\perp }$  is the orthogonal complement of the null space of the operator $A$  .

With this notation let $J\left(f\right)=\parallel y-Af{\parallel }_n^2$  the quadratic risk function. We will denote by $\hat{f}$  the function that minimizes the risk (which may not be unique), defined as

$$\begin{array}{c}\hat{f}=arg{min}_{f\in F}J\left(f\right),\end{array}$$

(2.1)

where the minimum is taken over all functions from $F$  to $Y$  .

The solution of the problem ${min}_{f\in F}J\left(f\right)$  exists if and only if $f$  is a solution of the normal equation

$$\begin{array}{c}{A}^{*}Af=A^{*}y\end{array}$$

(2.2)

where $A^{*}:Y\to F$  is the adjoint operator of $A$  (introduced via the requirement that for all $f\in F$  and $y\in Y$  , $\langle Ax,y{\rangle }_n=\frac{1}{n}\langle x,A^{*}y\rangle$  holds). It is important to remark that the operator $A^{*}$  actually depends on the observation sequence $x_i,i=1,\dots n$  . If $Y$  is generated by $\{{\phi }_j{\}}_{j=1}^m$  and is such that this basis is orthonormal with respect to the $L^2\left(P_n\right)$  norm over $Y$  , and $A$  is the identity then $A^{*}=\frac{1}{n}\left({\phi }_j\right(x_i){)}_{i,j}$  .

It is necessary to mention that the convergence rates can thus be given only over subsets of $F$  , i.e., under a-priori assumptions on the exact solution $\tilde{f}$  . We will formulate such a-priori assumptions, encountered typically in the inverse problem literature, in terms of the exact solutions by considering subsets of $F$  given by some source condition of the form $$A_{\mu ,\rho }=\{f\in F,f=(A^{*}A{)}^{\mu }\omega ,\parallel \omega \parallel \le \rho \}$$  where $0\le \mu \le {\mu }_0$  , ${\mu }_0>0$  and use the notation

$$\begin{array}{c}{A}_{\mu }={\bigcup }_{\rho >0}{A}_{\mu ,\rho }=\mathcal{ℛ}\left(\right(A^{*}A{)}^{\mu })\end{array}$$

(2.3)

These sets are usually called source sets, $f\in A_{\mu ,\rho }$  is said to have a source representation.

The requirement that $f$  be in $A_{\mu ,\rho }$  can be considered as an smoothness condition.

2.2 Projection methods

For numerical calculations, we have to approximate the space $F$  by a finite-dimensional subspace. Estimating over all $F$  is in general not possible. One approach in this direction is regularization by projection, where the regularization is achieved by a finite-dimensional approximation through projection.

Let $M_n$  be a collection of index sets ( $m\in M_n,m=\{j_1\dots ,j_{d_m}\}$  ). We give a sequence $Y_1\subset Y_2\dots \subset Y_m\dots \subset Y$  whose of union is dense in $Y$  . We assume $$dim\left(Y_m\right)=d_m.$$  Let ${\Pi }_{Y_m}^n$  be the orthogonal projector in the empirical norm over the subspace $Y_m$  and let $A_m={\Pi }_{Y_m}^{n}A$  . Define $F_m=A^{*}{Y}_m$  , with $A_m^{*}:Y_m\to F$  the adjoint operator of $A_m$  , and ${\Pi }_{F_m}$  to be the orthogonal projector onto the subspace $F_m$  . Then, by construction $${\Pi }_{F_m}=({\Pi }_{Y_m}^{n}A{)}^{+}{\Pi }_{Y_m}^{n}A.$$  Thus, we shall assume that data are give through an orthogonal design, corresponding to an orthogonal projection ${\Pi }_{Y_m}^n$  as

$$\begin{array}{c}{\Pi }_{Y_m}^{n}y={\Pi }_{Y_m}^{n}A\tilde{f}+{\Pi }_{Y_m}^n\varepsilon .\end{array}$$

(2.4)

With this notation we have that the best-approximate $L^2$  solution has the expression $${\Pi }_{F_m}\tilde{f}=A_m^{†}{y}_m.$$  for $y_m={\Pi }_{Y_m}^{n}y$  in the domain of $A_m^{†}$  . In the following we shall denote ${\tilde{f}}_m={\Pi }_{F_m}\tilde{f}$  .

Our goal is to find the solution of the equation  1.1 in the finite-dimensional subspace $F_m$  of $F$  . We have that for projection without regularization the choice of $F_m$  and of $A_m$  has many advantages. For noisy data and severely ill-posed problems the dimension of the subspace has to be rather low to keep total error estimate small, since for these problems the smallest singular value of $A_m$  decreases rapidly as $m$  increases. To be able to use larger dimensions we have to combine the projection method with additional regularization methods,such as iterative methods [9] ,[11] .

2.3 Singular value decomposition

As often $A_m$  is not of full rank, the singular value decomposition (SVD) of the operator $A_m$  is then a useful tool. Let $({\sigma }_j;{\phi }_j,{\phi }_j{)}_{j\in m}$  be a singular system for a linear operator $A_m$  , that is, $A_m{\phi }_j={\sigma }_j{\phi }_j$  and $A_m^{*}{\phi }_j={\sigma }_j{\phi }_j$  ; where $\{{\sigma }_j^2{\}}_{j\in m}$  are the nonzero eigenvalues of the selfadjoint operator $A_m^{*}{A}_m$  (and also of $A_m{A}_m^{*}$  ), considered in decreasing order. Furthermore, $\{{\phi }_j{\}}_{j\in m}$  and $\{{\phi }_j{\}}_{j\in m}$  are a corresponding complete orthonormal system of eigenvectors of $A_m^{*}{A}_m$  and $A_m{A}_m^{*}$  , respectively. For general linear operators with an SVD decomposition, we can write

$$\begin{array}{c}{A}_{m}f={\sum }_{j\in m}{\sigma }_j\langle f,{\phi }_j\rangle {\phi }_j\end{array}$$

(2.5)

$$\begin{array}{c}{A}_m^{*}{y}_m={\sum }_{j\in m}{\sigma }_j\langle y_m,{\phi }_j\rangle {\phi }_j.\end{array}$$

(2.6)

For $y_m$  in the domain of $A_m^{†}$  , $\mathcal{D}\left(A_m^{†}\right)$  , the best-approximate $L^2$  solution hast the expression $$A_m^{†}{y}_m={\sum }_{j\in m}\frac{\langle y_m,{\phi }_j\rangle }{{\sigma }_j}{\phi }_j={\sum }_{j\in m}\frac{\langle A_m^{*}{y}_m,{\phi }_j\rangle }{{\sigma }_j^2}{\phi }_j.$$  Note that for large $j$  , the term $1/{\sigma }_j$  grows to infinity. Thus, the high frequency errors are strongly amplified. We will asume that ${\sigma }_j=\mathcal{O}\left(j^{-p}\right)$  for some $p>1/2$  , which is clearly related to the ill-posedness of the operator $A$  and the approximation properties of $Y_m$  . For the construction and analysis of regularization methods, we will require some general notation for functions of the operators $A_m^{*}{A}_m$  and $A_m{A}_m^{*}$  .

Let $E_{\lambda }$  be the spectral decomposition of $A_m^{*}{A}_m$  given by $$E_{\lambda }(\cdot )={\sum }_{{\sigma }^2<{\lambda }_j,j\in m}\langle \cdot ,{\phi }_j\rangle {\phi }_j$$  and $H_{\lambda }$  the spectral decomposition of $A_m{A}_m^{*}$  . Then $E_{\lambda }$  is an orthogonal projector and projects onto $$span\left\{{\phi }_j\right|j\in m,{\sigma }^2<\lambda \}.$$  Since $({\sigma }_j^2;{\phi }_j)$  is an eigensystem for the selfadjoint compact operator $A_m^{*}{A}_m$  , $$A_m^{*}{A}_{m}f={\sum }_{j\in m}{\sigma }_j^2\langle f,{\phi }_j\rangle {\phi }_j$$  holds, which will be written (using the definition of the integral below) as $$A_m^{*}{A}_{m}f=\int \lambda d{E}_{\lambda }f$$  for $f\in \mathcal{D}\left(A_m\right)$  . Here the limits of integration could be 0 and $\parallel A_m{\parallel }^2+\epsilon$  for any $\epsilon >0$  . We sometimes omit the limits of integration.

This, motivates the definition

$$\begin{array}{c}G\left(A_m^{*}{A}_m\right):=\int G\left(\lambda \right)d{E}_{\lambda }:={\sum }_{{\sigma }^2={\lambda }_j,j\in m}G\left({\sigma }_j^2\right)\langle \cdot ,{\phi }_j\rangle {\phi }_j\end{array}$$

(2.7)

of a (piecewise) continuous function $G$  of a selfadjoint linear operator on $F_m$  . If $A^{*}A$  is continuously invertible, then $(A^{*}A{)}^{-1}=\int \frac{1}{\lambda }d{E}_{\lambda }$  .

In this case the best-approximate $L^2$  solution, for $y_m$  in the domain of $A_m^{†}$  , can be characterized by the equation

$$\begin{array}{c}{\tilde{f}}_m=A_m^{†}{y}_m=\int \frac{1}{\lambda }d{E}_{\lambda }{A}_m^{*}{y}_m.\end{array}$$

(2.8)

If $G\left(A_m^{*}{A}_m\right)$  is defined via  2.7 , then for $f\in \mathcal{D}\left(G_1\right(A_m^{*}{A}_m\left)\right)$  and $g\in \mathcal{D}\left(G_2\right(A_m^{*}{A}_m\left)\right)$ 

$$\begin{array}{c}\langle G_1\left(A_m^{*}{A}_m\right)f,G_2\left(A_m^{*}{A}_m\right)g\rangle =\int G_1\left(\lambda \right)G_2\left(\lambda \right)d\langle E_{\lambda }f,g\rangle \end{array}$$

(2.9)

and

$$\begin{array}{c}\parallel G\left(A_m^{*}{A}_m\right)f{\parallel }^2=\int G\left(\lambda \right)d\parallel E_{\lambda }f{\parallel }^2.\end{array}$$

(2.10)

The source set, $A_{\mu }$   2.3 , can be characterized via the singular values as follows:

$$(A_m^{*}{A}_m{)}^{\mu }\omega =\int {\lambda }^{\mu }d{E}_{\lambda }\omega ={\sum }_{j\in m}{\sigma }_j^{2\mu }\langle \omega ,{\phi }_j\rangle {\phi }_j$$ 

3 Regularization methods

After the general considerations of the last section, we now explain the construction of a regularization method for the important special case of selfadjoint linear operators. The basic idea for deriving a regularization method is to replace the amplification factors $1/{\lambda }_j$  by a filtered version $Q({\lambda }_j,\alpha )$  , where the filter function is a piecewise continuous, nonnegative and nonincreasing function of $\lambda$  on the segment $[0,\parallel A_m{\parallel }^2]$  for a regularization parameter $\alpha >0$  .

The assumptions over the regularizing coefficients $Q_{\alpha }\left(\lambda \right)$  are technical and are given in order to control fluctuations over set $[0,\parallel A_m{\parallel }^2]$  .

The filter family $\left\{Q\right({\lambda }_j,\alpha ){\}}_{j\in m}$  approximates the function ${\lambda }^{-1}$  for $\alpha \to \infty$  . Intuitively, a regularization on $A_m^{†}$  should then be the replacement of the ill conditioned operator $A_m^{†}$  by a family $\left\{R\right({\lambda }_j,\alpha ){\}}_{j\in m}:Y_m\to F_m$  of continuous operators. Throughout all the article, we shall denote $\left\{Q\right({\lambda }_j,\alpha ){\}}_{j\in m}$  and $\left\{R\right({\lambda }_j,\alpha ){\}}_{j\in m}$  by $Q_{\alpha }$  and $R_{\alpha }$  , respectively. Obviously, for all $\alpha >0,$  $R_{\alpha }$  is bounded.

As the approximation of ${\tilde{f}}_m$  , we then take $$f_{m,\alpha }=Q_{\alpha }\left(A_m^{*}{A}_m\right)A_m^{*}{y}_m=R_{\alpha }{y}_m,$$  where $R_{\alpha }:=\int Q_{\alpha }\left(\lambda \right)d{E}_{\lambda }{A}_m^{*}$  .

Remark 3.1. Note that with the above notation

$$\begin{array}{c}{f}_{m,\alpha }=R_{\alpha }{y}_m=R_{\alpha }{\Pi }_{Y_m}^{n}A\tilde{f}+R_{\alpha }{\Pi }_{Y_m}^n\varepsilon .\end{array}$$

(3.1)

Also that we can write $$R_{\alpha }=Q_{\alpha }\left(A_m^{*}{A}_m\right)A_m^{*}.$$ 

The next theorem gives conditions under which the first term in  3.1 converges to ${\tilde{f}}_m={\Pi }_{F_m}\tilde{f}$  . The proof follows that of [9] , but we include it for the sake of completeness.

Theorem 3.2. Let, for all $\alpha >0$  , $Q_{\alpha }:[0,\parallel A_m{\parallel }^2]\to IR$  be a piecewise continuous and nonincreasing function of $\lambda$  on the segment $[0,\parallel A_m{\parallel }^2]$  . Assume also that there is a $C>0$  such that $$\left|\lambda Q_{\alpha }\right(\lambda \left)\right|\le C,$$  and $${lim}_{\alpha \to \infty }{Q}_{\alpha }\left(\lambda \right)=\frac{1}{\lambda }$$  for all $\lambda \in (0,\parallel A_m{\parallel }^2).$  Then, for all $y\in D\left(A_m^{†}\right)$  , $${lim}_{\alpha \to \infty }{Q}_{\alpha }\left(A_m^{*}{A}_m\right)A_m^{*}{\Pi }_{Y_m}^{n}A\tilde{f}={\tilde{f}}_m$$  holds with ${\tilde{f}}_m=A_m^{†}{y}_m.$ 

Remark 3.3. In order to assume convergence as $\alpha \to \infty$  , it is necessary to choose $Q_{\alpha }$  such that it approximates $1/\lambda$  for all $\lambda \in (0,\parallel A_m{\parallel }^2]$  . Also, note that the condition $\left|\lambda Q_{\alpha }\right(\lambda \left)\right|\le C$  implies that $\parallel A_m{R}_{\alpha }\parallel =\parallel A_m{A}_m^{*}{Q}_{\alpha }\left(A_m^{*}{A}_m\right)\parallel \le C$  , i.e, $\parallel A_m{R}_{\alpha }\parallel$  is uniformly bounded.

• Proof. As in [9] , if ${\tilde{f}}_m$  is defined by  2.8 , then by  2.2 the residual norm has the representation $$\parallel {\tilde{f}}_m-Q_{\alpha }\left(A_m^{*}{A}_m\right)A_m^{*}{A}_m\tilde{f}{\parallel }^2=\parallel (I-Q_{\alpha }(A_m^{*}{A}_m\left)A_m^{*}{A}_m\right){\tilde{f}}_m{\parallel }^2$$  From the formula  2.10 , it follows that $$\parallel {\tilde{f}}_m-Q_{\alpha }\left(A_m^{*}{A}_m\right)A_m^{*}{A}_m\tilde{f}{\parallel }^2={\int }_0^{\parallel A_m{\parallel }^2+}(1-\lambda Q_{\alpha }(\lambda ){)}^{2}d\parallel E_{\lambda }{\tilde{f}}_m{\parallel }^2.$$  Since $(1-\lambda Q_{\alpha }(\lambda ){)}^2$  is bounded by the constant $(1+C{)}^2$  , which is integrable with respect to the measure $d\parallel E_{\lambda }{\tilde{f}}_m{\parallel }^2,$  then by the Dominated Convergence Theorem,

  $$\begin{array}{c}{lim}_{\alpha \to \infty }{\int }_0^{\parallel A_m{\parallel }^2+}(1-\lambda Q_{\alpha }(\lambda ){)}^{2}d\parallel E_{\lambda }{\tilde{f}}_m{\parallel }^2={\int }_0^{\parallel A_m{\parallel }^2+}{lim}_{\alpha \to \infty }(1-\lambda Q_{\alpha }\left(\lambda \right){)}^{2}d\parallel E_{\lambda }{\tilde{f}}_m{\parallel }^2.\end{array}$$

  (3.2)

  Since for $\lambda >0$  , ${lim}_{\alpha \to \infty }(1-\lambda Q_{\alpha }(\lambda \left)\right)=0$  then the integral on the right-hand side of  3.2 equals to $0$  . On the other hand, if $\lambda =0$  , ${lim}_{\alpha \to \infty }(1-\lambda Q_{\alpha }(\lambda \left)\right)=1$  then the equation  3.2 has the form

  $$\begin{array}{c}{lim}_{\alpha \to \infty }{\int }_0^{\parallel A_m{\parallel }^2+}(1-\lambda Q_{\alpha }(\lambda ){)}^{2}d\parallel E_{\lambda }{\tilde{f}}_m{\parallel }^2={lim}_{\lambda \to 0^{+}}\parallel E_{\lambda }{\tilde{f}}_m{\parallel }^2-\parallel E_0{\tilde{f}}_m{\parallel }^2\end{array}$$

  (3.3)

  which is equal the jump of $\parallel E_{\lambda }{\tilde{f}}_m{\parallel }^2$  at $\lambda =0$  . Since ${\tilde{f}}_m\in \mathcal{N}(A_m{)}^{\perp }$  , the term on the right-hand side of  3.3 equals to 0. Thus, $R_{\alpha }{A}_m\tilde{f}$  converges to ${\tilde{f}}_m$  as $\alpha \to \infty$  for $y_m\in D\left(A_m^{†}\right),$  which ends the proof.
  □

Let $Tr\left(B\right)$  the trace of the selfadjoint operator $B^{t}B$  for any square matriz $B$  , which is defined by $$Tr\left(B\right)=\frac{1}{n}{\sum }_{j\in m}{b}_j$$  for $b_j$  eigenvalues of $B^{t}B$  .

We then have the following result,

Theorem 3.4. Let $Q_{\alpha }$  be as in theorem  3.2 . Let $\mu ,\rho >0$  and let ${\omega }_{\mu }:(0,{\alpha }_0)\to R$  be such that for all $\alpha \in (0,{\alpha }_0)$  and $\lambda \in [0,{\sigma }_1^2]$  $${sup}_{0\le \lambda \le {\sigma }_1^2}{\lambda }^{\mu }|1-\lambda Q_{\alpha }(\lambda \left)\right|\le {\omega }_{\mu }\left(\alpha \right)$$  holds. Then for $\tilde{f}\in A_{\mu ,\rho },$  the following inequality holds true

$$\begin{array}{c}IE\parallel {\tilde{f}}_m-f_{\alpha ,m}{\parallel }^2\le 2{\omega }_{\mu }(\alpha {)}^2{\rho }^2+2{\sigma }^{2}Tr(Q_{\alpha }^2\left(A_m^{*}{A}_m\right)A_m{A}_m^{*}).\end{array}$$

(3.4)

• Proof. The proof of this inequality is based on the definition of the estimator $f_{m,\alpha }$  and on the assumptions over this function. We have that the $L^2-$  norm of the difference between the regularized function and the true data function can be bounded by

  $$\begin{array}{c}IE\parallel {\tilde{f}}_m-f_{m,\alpha }{\parallel }^2\le 2IE\parallel {\tilde{f}}_m-R_{\alpha }{A}_m\tilde{f}{\parallel }^2+2IE\parallel R_{\alpha }{A}_m\tilde{f}-f_{m,\alpha }{\parallel }^2\end{array}$$

  (3.5)

  where $f_{m,\alpha }=R_{\alpha }{y}_m$  .
  This is the typical bias-variance decomposition. The first term on the right-hand side is an approximation error, which corresponds to the bias, whereas the second term, variance, is a stability bound on the regularizing operator $R_{\alpha }$  . Note that by the Theorem  3.2 , the first term in  3.5 goes to 0 if $y_m\in \mathcal{D}\left(T_m^{†}\right)$  .

  Let $\omega \in F_m$  with $\parallel \omega \parallel \le \rho$  . Since $\tilde{f}\in F_m$  then ${\Pi }_{F_m}\tilde{f}=(A_m^{*}{A}_m{)}^{\mu }\omega$  . On the other hand, ${\lambda }^{\mu }{sup}_{\lambda }|1-\lambda Q_{\alpha }(\lambda \left)\right|\le {\omega }_{\mu }\left(\alpha \right)$  , then the first term in this equation can be bounded by

  $$\begin{array}{cc}IE\parallel {\tilde{f}}_m-R_{\alpha }{A}_m\tilde{f}{\parallel }^2& =IE\parallel {\tilde{f}}_m-Q_{\alpha }\left(A_m^{*}{A}_m\right)A_m^{*}{A}_m{\tilde{f}}_m{\parallel }^2\end{array}$$

  $$\begin{array}{cc}& =IE\parallel (I-Q_{\alpha }(A_m^{*}{A}_m\left)A_m^{*}{A}_m\right){\tilde{f}}_m{\parallel }^2\end{array}$$

  $$\begin{array}{cc}& \le {\omega }_{\mu }^2{\rho }^2.\end{array}$$

  In order to control the term corresponding to the variance we used that the data perturbation is white noise. Thus,

  $$\begin{array}{cc}IE\parallel R_{\alpha }{A}_m\tilde{f}-f_{m,\alpha }{\parallel }^2& =IE\langle \varepsilon ,\left(Q_{\alpha }\right(A_m^{*}{A}_m\left)A_m^{*}{)}^{*}{Q}_{\alpha }\right(A_m^{*}{A}_m)A_m^{*}\varepsilon \rangle \end{array}$$

  $$\begin{array}{cc}& =IE\langle \varepsilon ,Q_{\alpha }^2\left(A_m^{*}{A}_m\right)A_m{A}_m^{*}\varepsilon \rangle \end{array}$$

  $$\begin{array}{cc}& ={\sigma }^{2}Tr\left(Q_{\alpha }^2\right(A_m^{*}{A}_m\left)A_m{A}_m^{*}\right)\end{array}$$

  which yields the desired result. □

The next result will be useful when studying iterative methods.

Theorem 3.5. Let $Q_{\alpha }$  be as in theorem  3.2 . Assume also that $Q_{\alpha }$  is continuously differentiable and that the function $|1-\lambda Q_{\alpha }(\lambda \left){|}^{\prime }\right|\lambda Q_{\alpha }\left(\lambda \right)-1{|}^{-1}$  doest not decrease. Then the estimates are valid $${sup}_{0\le \lambda \le {\sigma }_1^2}\left|Q_{\alpha }\right(\lambda \left)\right|=Q_{\alpha }\left(0\right),$$  and $${sup}_{0\le \lambda \le {\sigma }_1^2}{\lambda }^{\mu }|1-\lambda Q_{\alpha }(\lambda \left)\right|<{\mu }^{\mu }(\mu +1{)}^{-1}{\omega }_{\mu }(\alpha )$$  where ${\omega }_{\mu }\left(\alpha \right)=Q_{\alpha }(0{)}^{-\mu }.$ 

• Proof. The proof can be carried out by standard techniques. A proof of this result can be found in [11] . □

4 Rates of convergence for the regularized estimator

In any regularization method, the regularization parameter $\alpha$  plays a crucial role. For choosing the parameter, there are general methods of parameter selection. For example, the Discrepancy Principe [20] , Cross-Validation [7] and the L-curve [10] . They differ in the amount of a priori information required as well as in the decision criteria. The appropriate choice of regularization parameter is a difficult problem. We would like too choose $\alpha$  , based on the data in such a way that optimal rates are maintained. This choice should not depend on a priori regularity assumptions.

Our goal is to introduce adaptive methods in the context of statistical inverse problems.

In this section we introduce our adaptive estimator, for a fixed $m=m_0$  . We choose $m_0$  such that $\parallel \tilde{f}-{\Pi }_{F_{m_0}}\tilde{f}{\parallel }^2$  satisfies the optimal rates with high probability since we know $$\parallel \tilde{f}-{\Pi }_{F_{m_0}}\tilde{f}{\parallel }^2<\parallel I-{\Pi }_{Y_{m_0}}{\parallel }^{4\mu }=\mathcal{O}\left(d_{m_0}^{-4\mu p}\right)$$  for a certain $p$  and $0<\mu \le 1/2$  . It is satisfied if the dimension of the set is such that

$$\begin{array}{c}{d}_{m_0}\ge n^{\frac{1}{2p+1}}.\end{array}$$

(4.1)

This leads to the rate $$\parallel \tilde{f}-f_{m,\alpha }{\parallel }^2=\mathcal{O}\left(n^{-\frac{4\mu p}{4\mu p+2p+1}}\right).$$  Analogous results are obtained in the case of Hilbert scales ([12] ,[16] ).

Adaptive model selection is a technique which penalizes the regularization parameter, in such a way that we choose ${\hat{f}}_{m_0,{\alpha }_{\hat{k}}}$  by minimizing $$arg{min}_{k\in \mathcal{K},f\in F_m}\left(\parallel R_{{\alpha }_k}(y_m-A_{m}f){\parallel }^2+pen\left({\alpha }_k\right)\right)$$  where $$\hat{k}=arg{min}_{k\in \mathcal{K}}\left(\parallel R_{{\alpha }_k}(y_m-A_{m}f){\parallel }^2+pen\left({\alpha }_k\right)\right)$$  and $$pen\left({\alpha }_k\right)=r{\sigma }^2(1+L_k)\left[Tr\right(R_{{\alpha }_k}^t{R}_{{\alpha }_k})+{\rho }^2(R_{{\alpha }_k}\left)\right],$$  with $r>2$  and $L_k$  is a sequence which is incorporated in order to control the complexity of the set $\mathcal{K}=\{1,2,\dots ,k_n\}$  , of all possible index up to $k_n$  . Here ${\rho }^2\left(B\right)=\rho \left(B^{t}B\right)$  is the spectral radius of the selfadjoint operator $B^{t}B$  for any square matriz $B$  , which is defined by $${\rho }^2\left(B\right)=\frac{1}{n}{max}_{j\in m}{b}_j$$  for $b_j$  eigenvalues of $B^{t}B$  .

Thus, $\hat{k}$  is selected by minimizing

$$\begin{array}{c}arg{min}_{k\in \mathcal{K}}\left(\parallel R_{{\alpha }_k}(y_m-A_{m}f){\parallel }^2+\frac{r{\sigma }^2(1+L_k)}{n}\left[{\sum }_{j\in m}{Q}^2\left({\lambda }_j\right){\lambda }_j+{max}_{j\in m}{Q}^2\left({\lambda }_j\right){\lambda }_j\right]\right).\end{array}$$

(4.2)

The strategy as proposed in this article automatically provides the optimal order of accuracy.

The regularized estimator has a rate of convergence less or equal than the best rate achieved by the best estimator for a selected model. We have the following result,

Theorem 4.1. For any $f\in F_m$  and any ${\alpha }_k$  the following inequality holds true for $d$  a positive constant that depends on r (as in Lemma  4.4 ),

$$\begin{array}{c}IE\parallel {\tilde{f}}_m-{\hat{f}}_{{\alpha }_{\hat{k}}}{\parallel }^2\le \frac{1}{(1-\nu )}{inf}_{k\in \mathcal{K}}\left[C\right(1+\nu )\parallel {\tilde{f}}_m-f_{{\alpha }_k}{\parallel }^2+2pen({\alpha }_k\left)\right]+\frac{C_1\left(d\right)}{n}\end{array}$$

(4.3)

where $C_1\left(d\right)=4{\sigma }^2{\sum }_k\frac{n{\rho }^2\left(R_{{\alpha }_k}\right)}{d}\left[\sqrt{dr{L}_k\left[\frac{Tr\left(R_{{\alpha }_k}^t{R}_{{\alpha }_k}\right)}{{\rho }^2\left(R_{{\alpha }_k}\right)}+1\right]}+1\right]e^{-\sqrt{dr{L}_k\left[\frac{Tr\left(R_{{\alpha }_k}^t{R}_{{\alpha }_k}\right)}{{\rho }^2\left(R_{{\alpha }_k}\right)}+1\right]}}.$ 

Remark 4.2. An important issue is that equation  4.3 is non asymptotic. The goodness of fit of the estimator is defined by trace, $Tr\left(R_{\alpha }^t{R}_{\alpha }\right)$  , and spectral radius, ${\rho }^2\left(R_{\alpha }\right)$  . Also, the estimator is optimal in the sense that the adaptive estimator achieves the best rate of convergence among all the regularized estimators.

Remark 4.3. Remark that under our assumptions, namely that the basis is orthonormal for the fixed design, both $n{\rho }^2\left(R_k\right)$  and $Tr\left(R_k^t{R}_k\right)/{\rho }^2\left(R_k\right)$  do not depend on n.

• Proof. For any $f_{{\alpha }_k}$  and any $k\in IN$ 

  $$\begin{array}{ccc}& \parallel R_{{\alpha }_{\hat{k}}}(y_m-A_m{\hat{f}}_{{\alpha }_{\hat{k}}}){\parallel }^2+pen\left({\alpha }_{\hat{k}}\right)\le \parallel R_{{\alpha }_k}(y_m-A_m{f}_{{\alpha }_k}){\parallel }^2+pen\left({\alpha }_k\right)& \end{array}$$

  and

  $$\begin{array}{ccc}& \parallel R_{{\alpha }_k}(y_m-A_m{f}_{{\alpha }_k}){\parallel }^2=\parallel R_{{\alpha }_k}{A}_m(\tilde{f}-f_{{\alpha }_k}){\parallel }^2+2\langle R_{{\alpha }_k}{A}_m(\tilde{f}-f_{{\alpha }_k}),R_{{\alpha }_k}{\Pi }_{Y_m}^n\varepsilon \rangle +\parallel R_{{\alpha }_k}{\Pi }_{Y_m}^n\varepsilon {\parallel }^2& \end{array}$$

  Thus, following standard arguments we have

  $$\begin{array}{ccc}& & \parallel R_{{\alpha }_{\hat{k}}}{A}_m(\tilde{f}-{\hat{f}}_{{\alpha }_{\hat{k}}}){\parallel }^2\end{array}$$

  $$\begin{array}{ccc}& \le & \parallel R_{{\alpha }_k}{A}_m(\tilde{f}-f_{{\alpha }_k}){\parallel }^2-2<R_{{\alpha }_{\hat{k}}}{A}_m(\tilde{f}-{\hat{f}}_{{\alpha }_{\hat{k}}}),R_{{\alpha }_{\hat{k}}}{\Pi }_{Y_m}^n\varepsilon >\end{array}$$

  $$\begin{array}{ccc}& & +2<R_{{\alpha }_k}{A}_m(\tilde{f}-f_{{\alpha }_k}),R_{{\alpha }_k}{\Pi }_{Y_m}^n\varepsilon >-\parallel R_{{\alpha }_{\hat{k}}}{\Pi }_{Y_m}^n\varepsilon {\parallel }^2+\parallel R_{{\alpha }_k}{\Pi }_{Y_m}^n\varepsilon {\parallel }^2+pen\left({\alpha }_k\right)+pen\left({\alpha }_{\hat{k}}\right).\end{array}$$

  Let $0<\nu <1$  . Since the algebraic inequality $2ab\le \nu a^2+\frac{1}{\nu }{b}^2$  holds for all $a,b\in IR$  , we find that

  $$\begin{array}{ccc}& & (1-\nu )\parallel R_{{\alpha }_{\hat{k}}}{A}_m(\tilde{f}-{\hat{f}}_{{\alpha }_{\hat{k}}}){\parallel }^2\end{array}$$

  $$\begin{array}{ccc}& \le & (1+\nu )\parallel R_{{\alpha }_k}{A}_m(\tilde{f}-f_{{\alpha }_k}){\parallel }^2+2pen\left({\alpha }_k\right)+2{sup}_{{\alpha }_k}\{\frac{1}{\nu }\parallel R_{{\alpha }_k}{\Pi }_{Y_m}^n\varepsilon {\parallel }^2-pen({\alpha }_k\left)\right\},\end{array}$$

  holds for any $k$  and $f_{{\alpha }_k}\in F_m$  .
  On the other hand, using that is $1\le \parallel R_{{\alpha }_k}A\parallel \le C$  , we have that for any $f_{{\alpha }_k}\in F_{m_0}$  and any $k\in IN$  ,

  $$\begin{array}{ccc}& & (1-\nu )\parallel {\tilde{f}}_m-{\hat{f}}_{{\alpha }_{\hat{k}}}{\parallel }^2\le C(1+\nu )\parallel {\tilde{f}}_m-f_{{\alpha }_k}{\parallel }^2\end{array}$$

  $$\begin{array}{ccc}& +& 2pen\left({\alpha }_k\right)+2C_1{sup}_{{\alpha }_k}\{\parallel R_{{\alpha }_k}{\Pi }_{Y_m}^n\varepsilon {\parallel }^2-pen({\alpha }_k\left)\right\}.\end{array}$$

  The proof then follows directly from the following technical lemma ([3] ,[16] ) which characterizes the supremum of an empirical process by the regularization family.

  Lemma 4.4. Let $\eta \left(A\right)=\sqrt{{\varepsilon }^t{A}^{t}A\varepsilon }=\parallel A\varepsilon \parallel .$  Then, there exists a positive constant $d$  that depends on $r/2$  such that the following inequality holds

  $$\begin{array}{ccc}& & P\left({\eta }^2\right(A)\ge {\sigma }^2[Tr\left(A^{t}A\right)+\rho \left(A^{t}A\right)]r/2(1+L)+{\sigma }^{2}u)\end{array}$$

  (4.4)

  $$\begin{array}{ccc}& & \le exp\{-\sqrt{d(1/\rho (A^{t}A)u+r/2L[Tr\left(A^{t}A\right)/\rho \left(A^{t}A\right)+1\left]\right)}\}.\end{array}$$

  With the above notation, $$\eta \left(R_{{\alpha }_k}\right)=\parallel R_k{\varepsilon }_m\parallel$$  where ${\varepsilon }_m={\Pi }_{Y_m}^n\varepsilon$  .

  Now, with this lemma we have

  $$\begin{array}{ccc}& & P({sup}_{{\alpha }_k}\parallel R_{{\alpha }_k}{\varepsilon }_m{\parallel }^2-pen({\alpha }_k)>{\sigma }^{2}x)\end{array}$$

  $$\begin{array}{ccc}& \le & {\sum }_{k}P\left[{\eta }^2\right(R_{{\alpha }_k})\ge r{\sigma }^2(1+L_k\left)\right[Tr\left(R_{{\alpha }_k}^t{R}_{{\alpha }_k}\right)+{\rho }^2\left(R_{{\alpha }_k}\right)]+{\sigma }^{2}x]\end{array}$$

  $$\begin{array}{ccc}& \le & {\sum }_{k}exp\{-\sqrt{d(1/{\rho }^2(R_{{\alpha }_k})x+r{L}_k[Tr\left(R_{{\alpha }_k}^t{R}_{{\alpha }_k}\right)/{\rho }^2\left(R_{{\alpha }_k}\right)+1\left]\right)}\}\end{array}$$

  Since for $X$  positive $IEX={\int }_0^{\infty }P(X>u)du$  , we then have that

  $$\begin{array}{ccc}& & IE[{sup}_{{\alpha }_k}\parallel R_{{\alpha }_k}{\varepsilon }_m{\parallel }^2-pen({\alpha }_k\left)\right]={\int }_0^{\infty }P[{sup}_{{\alpha }_k}\parallel R_{{\alpha }_k}{\varepsilon }_m{\parallel }^2-pen({\alpha }_k)\ge x]dx\end{array}$$

  $$\begin{array}{ccc}& =& {\sigma }^2{\int }_0^{\infty }P[{sup}_{{\alpha }_k}\parallel R_{{\alpha }_k}{\varepsilon }_m{\parallel }^2-pen({\alpha }_k)\ge {\sigma }^{2}u]du\end{array}$$

  $$\begin{array}{ccc}& \le & {\sigma }^2{\sum }_k{\int }_0^{\infty }exp\{-\sqrt{k_{1}u+k_2}\}du.\end{array}$$

  where $k_1=d/{\rho }^2\left(R_{{\alpha }_k}\right)$  and $k_2=dr{L}_k\left[Tr\right(R_{{\alpha }_k}^T{R}_{{\alpha }_k})/{\rho }^2(R_{{\alpha }_k})+1]$  .

  Let $w=k_{1}u+k_2,$  then

  $$\begin{array}{cc}IE[{sup}_{{\alpha }_k}\parallel R_{{\alpha }_k}{\varepsilon }_m{\parallel }^2-pen({\alpha }_k\left)\right]& \le {\sigma }^2{\sum }_k{\int }_{k_2}^{\infty }\frac{1}{k_1}exp\{-\sqrt{w}\}dw\end{array}$$

  $$\begin{array}{cc}& ={\sigma }^2{\sum }_k\frac{2}{k_1}[-\sqrt{k_2}+1]exp\{-\sqrt{k_2}\}\end{array}$$

  Finally, we have the desired result.

  $$\begin{array}{cc}IE\parallel {\tilde{f}}_m-{\hat{f}}_{{\alpha }_{\hat{k}}}{\parallel }^2\le \frac{1}{(1-\nu )}{inf}_{k\in \mathcal{K}}\left[C\right(1+\nu )\parallel {\tilde{f}}_m-f_{{\alpha }_k}{\parallel }^2+2pen({\alpha }_k\left)\right]+\frac{C_1\left(d\right)}{n}& \end{array}$$

  where $C_1\left(d\right)=4{\sigma }^2{\sum }_k\frac{n{\rho }^2\left(R_{{\alpha }_k}\right)}{d}\left[\sqrt{dr{L}_k\left[\frac{Tr\left(R_{{\alpha }_k}^t{R}_{{\alpha }_k}\right)}{{\rho }^2\left(R_{{\alpha }_k}\right)}+1\right]}+1\right]e^{-\sqrt{dr{L}_k\left[\frac{Tr\left(R_{{\alpha }_k}^t{R}_{{\alpha }_k}\right)}{{\rho }^2\left(R_{{\alpha }_k}\right)}+1\right]}},$  □

5 Regularization by iterative methods

Iterative regularization methods, are very competitive methods for linear inverse problems.

In iterative regularization, one picks an initial guess $f_0$  for the unknown $\tilde{f}$  , and then one iteratively constructs updated approximations via a regularization scheme. The regularization parameter associated with iterative regularization is thus the $“$  stopping point”of the iterative sequence, and an important part of the mathematical theory is the development of stopping criteria for terminating the iteration. In other words, the iteration index plays the role of the regularization parameter $\alpha$  , and the stopping criteria plays of the parameter selection method.

5.1 Descent Methods for Linear Inverse Problems

As an example of iterative regularization, we consider descent methods. Descendent methods have become quite popular in the last years for the solution of linear inverse problems and for nonlinear inverse problems [11] . In this subsection we consider two examples.

As an approximation of ${\tilde{f}}_m$  we will choose $f_{m,\alpha }$  such that

$$\begin{array}{c}{f}_{m,\alpha }=[I-A_m^{*}{A}_m{Q}_{\alpha }(A_m^{*}{A}_m\left)\right]f_0+Q_{\alpha }\left(A_m^{*}{A}_m\right)A_m^{*}{\Pi }_{Y_m}^{n}y\end{array}$$

(5.1)

where $f_0\in F_m$  is an initial approach and this $f_0\in \mathcal{N}(A_m{)}^{\perp }$  [11] .

Most iterative methods for approximating $\tilde{f}$  are based on a transformation of the normal equation into equivalent fixed point equations like $$f=f+A_m^{*}(A_{m}f-y)$$  If $\parallel A_m{\parallel }^2<2$  then the corresponding fixed point operator $I-A_m^{*}{A}_m$  is nonexpansive and one may apply the method of successive approximations. It must be emphasized that $I-A_m^{*}{A}_m$  is no contradiction if our inverse problem is ill-posed, since the spectrum of $A_m^{*}{A}_m$  clusters at the origin.

5.2 Landweber iteration

In this subsection we presented the well-known Landweber iteration, which arises from converting the necessary conditions for minimizing  2.1 into a fixed point iteration. Much development in the last few years has taken place in advancing the theory of Landweber iteration for linear and nonlinear inverse problems.

Using the terminology of the last sections, we introduce the function

$$\begin{array}{c}{Q}_k\left(\lambda \right)={\sum }_{j=0}^{k-1}(1-\lambda {)}^j={\lambda }^{-1}(1-(1-\lambda {)}^k)\end{array}$$

(5.2)

We call $Q_k$  the iteration polynomial of degree $k-1$  . Associated with it is the polynomial $$r_k\left(\lambda \right)=1-\lambda Q_k\left(\lambda \right)=(1-\lambda {)}^k$$  of degree $k$  , which is called the residual polynomial since it determines the residual $y-A_m{f}_{m,k}$  .

Thus, inserting the equation  5.2 in  5.1 we obtain recursively,

$$\begin{array}{c}{f}_{m,k+1}=f_{m,k}-A_m^{*}(A_m{f}_{m,k}-y_m),k=0,1,\dots \end{array}$$

(5.3)

starting from an initial guess $f_0$  . This is a steepest descent method called the linear version of Landweber's iteration. Each step of the iterative process  5.3 is carried out along the direction opposite to the direction of the gradient of the quadratic functional $J\left(f\right)$  in  2.1 . It is known that there is the greatest decrease of the functional along this direction.

If $\parallel A_m\parallel \le 1$  , we considerer $\lambda \in (0,1]$  such that in this interval $\lambda Q_k\left(\lambda \right)$  is uniformly bounded and since $Q_k\left(\lambda \right)$  converge to $1/\lambda$  as $k\to \infty$  then according to Theorem  3.2 the sequences $f_{m,k}$  converge to ${\tilde{f}}_m$  when $y\in \mathcal{D}\left(A_m^{†}\right)$  . If $\parallel A_m\parallel$  is not bounded by one, then we introduce a relaxation parameter $0<\tau <\parallel A_m{\parallel }^{-2}$  in front of $A_m^{*}$  in  5.3 , i.e, we would iterate

$$\begin{array}{c}{f}_{m,k+1}=f_{m,k}-\tau A_m^{*}(A_m{f}_{m,k}-y),k=0,1,\dots \end{array}$$

(5.4)

If $\tau \equiv {\tau }_k$  , one can obtain various variants of the method of steepest descent depending on a choice of the sequence ${\tau }_k$  . The Landweber iteration  5.4 is usually called a method of simple iteration.

In the following we derive a simple estimate for the error propagation in the Landweber iteration. We then have the following result,

Corollary 5.1. Let $\tau =1/(2\parallel A_m{\parallel }^2)<1/{\lambda }_1$  . If $y\in \mathcal{ℛ}\left(A_m\right)$  , then the Landweber iteration is an order optimal regularization method, i.e, $$\parallel {\tilde{f}}_m-f_{k\left(m\right)}{\parallel }^2\le 2c_1{k}^{-2\mu }+2c_2\frac{{\sigma }^2}{n}(\tau k{)}^{(2p+1)/2p},$$  where $c_1={\rho }^2(\frac{\mu }{\tau e}{)}^{\mu }$  and $c_2=\frac{1}{2p+1}(\frac{2p+1}{2p-1}{)}^{(2p+1)/4p}.$ 

• Proof. To apply Theorem  3.4 we have to study the terms of the bias, $IE\parallel {\tilde{f}}_m-R_{\alpha }{A}_m\tilde{f}{\parallel }^2$  , and variance $IE\parallel R_{\alpha }{A}_m\tilde{f}-f_{k\left(m\right)}\parallel$  . By  5.2 we have $${\tilde{f}}_m-R_{\alpha }{A}_m\tilde{f}=(I-A_m^{*}{A}_m{Q}_k(A_m^{*}{A}_m\left)\right){\tilde{f}}_m=(I-A_m^{*}{A}_m{)}^k{\tilde{f}}_m$$  We have to study the residual polynomial $r_k\left(\lambda \right)=(1-\lambda {)}^k$  of the Landweber iteration.
  For $0\le \lambda \le \parallel A_m{\parallel }^2$  the function $${\lambda }^{\mu }|1-\lambda Q_k(\lambda \left)\right|$$  assumes its maximum for $\lambda ={\tau }^{-1}\mu (\mu +k{)}^{-1}$  .

  Thus, we have

  $$\begin{array}{cc}{\lambda }^{\mu }|1-\lambda Q_k(\lambda \left)\right|& \le {max}_{0\le \lambda \le \parallel A_m{\parallel }^2}{\lambda }^{\mu }|1-\lambda Q_k(\lambda \left)\right|\end{array}$$

  $$\begin{array}{cc}& <\frac{{\mu }^{\mu }}{{\tau }^{\mu }(\mu +k{)}^{\mu }}\frac{k^k}{(\mu +k{)}^k}\end{array}$$

  $$\begin{array}{cc}& <{\left(\frac{\mu }{\tau e}\right)}^{\mu }{k}^{-\mu }\end{array}$$

  This leads to numbers ${\omega }_{\mu }\left(k\right)$  as introduced in Theorem  3.4  $${\omega }_{\mu }\left(k\right)={\left(\frac{\mu }{\tau e}\right)}^{\mu }{k}^{-\mu }$$  Thus, the term corresponding to the bias is bounded by $$\parallel {\tilde{f}}_m-R_{\alpha }{A}_m\tilde{f}{\parallel }^2\le {\rho }^2(\frac{\mu }{\tau e}{)}^{2\mu }{k}^{-2\mu }.$$  Next, we establish bounds for the variance term. By assumption, the singular values satisfy ${\lambda }_j\approx j^{-2p}.$  Note that for small values of ${\lambda }_j$  we have

  $Q_k^2\left(\lambda \right)\le (\tau k{)}^2$  $\forall$  $j>m^{\prime }$ 

  and for big values of ${\lambda }_j$  ( ${\lambda }_j\approx {\lambda }_1$  )

  $Q_k^2\left({\lambda }_j\right)\le {\lambda }_j^{-2}$  $\forall$  $j<m^{\prime }$  .

  Consequently,

  $$\begin{array}{cc}nTr\left(Q_k^2\right(A_m{A}_m^{*}\left)A_m{A}_m^{*}\right)={\sum }_{j=1}^m{Q}_k^2\left({\lambda }_j\right){\lambda }_j& \le {\sum }_{j=1}^{m^{\prime }}{\lambda }_j^{-1}+{\sum }_{j>m^{\prime }}(\tau k{)}^2{\lambda }_j\end{array}$$

  $$\begin{array}{cc}& \le {\int }_0^{m^{\prime }}{s}^{2p}ds+(\tau k{)}^2{\int }_0^{m^{\prime }}{s}^{-2p}ds\end{array}$$

  This suggest searching $m^{\prime }\approx c(\tau k{)}^{1/2p}$  for $p>1/2$  , where $c=(\frac{2p+1}{2p+1}{)}^{1/4p}$  . Hence we have,

  $$\begin{array}{cc}IE\parallel R_{\alpha }{A}_m\tilde{f}-f_{k\left(m\right)}{\parallel }^2& =Tr\left(Q_k^2\right(A_m{A}_m^{*}\left)A_m{A}_m^{*}\right)\end{array}$$

  $$\begin{array}{cc}& \le \frac{c^{2p+1}}{2p+1}\frac{(\tau k{)}^{(2p+1)/2p}}{n}.\end{array}$$

  Finally, this implies $$IE\parallel {\tilde{f}}_m-f_{k\left(m\right)}{\parallel }^2\le 2c_1{k}^{-2\mu }+2c_2\frac{{\sigma }^2}{n}(\tau k{)}^{(2p+1)/2p},$$  where $c_1={\rho }^2(\frac{\mu }{\tau e}{)}^{\mu }$  and $c_2=\frac{1}{2p+1}(\frac{2p+1}{2p-1}{)}^{(2p+1)/4p}$  .

  □

Remark 5.2. Note that under the above inequality is satisfied if the dimension of the set is such that $d_{m_0}\approx n^{\frac{1}{4\mu p+2p+1}}$  . Here, the optimal choice of regularization sequence, depending on $p$  and $\mu$  . The optimal rates are of order $IE\parallel \tilde{f}-f_{k\left(m\right)}{\parallel }^2=\mathcal{O}\left(n^{-\frac{4\mu p}{4\mu p+2p+1}}\right)$  . Analogous results are obtained in the ill-posed problem literature, see for example [5] , where typically in a Hilbert scale setting optimal rates are of order $\mathcal{O}\left(n^{-\frac{2s}{2s+2p+1}}\right)$  , with $s=2\mu p$  .

We are ready to state our main result for the Landweber iteration, which bounds the mean squared error of the select estimate ${\hat{f}}_{\hat{k}}$  basically by the smallest mean squared error among the estimates $f_k$  plus a remainder term of order $1/n$  . The result follows from Theorem  4.1 .

Corollary 5.3. Let $\tau =1/(2\parallel A_m{\parallel }^2)<1/{\lambda }_1$  . Next assume $\hat{k}$  as in  4.2 and $d_{m_0}$  as in  4.1 . If $y\in \mathcal{ℛ}\left(A_m\right)$  then for any $f\in F_m$  and any $k$  , the following inequality holds true

$$\begin{array}{cc}IE\parallel {\tilde{f}}_m-{\hat{f}}_{\hat{k}}{\parallel }^2& \le \frac{1}{(1-\nu )}{inf}_{k\in \mathcal{K}}\left[C(1+\nu )\parallel {\tilde{f}}_m-f_k{\parallel }^2+\frac{2r{\sigma }^2(1+L_k)\left(c\right(\tau k{)}^{\frac{2p+1}{2p}}+\tau k)}{n}\right]\end{array}$$

$$\begin{array}{cc}& +\frac{4{\sigma }^2}{n}{\sum }_k\frac{\tau k}{d}\left[\sqrt{dr{L}_k\left[c(\tau k{)}^{1/2p}+1\right]}+1\right]e^{-\sqrt{dr{L}_k\left[c(\tau k{)}^{1/2p}+1\right]}},\end{array}$$

for some $C>0$  and $c=\frac{1}{2p+1}(\frac{2p+1}{2p-1}{)}^{(2p+1)/4p}$  .

• Proof. For fixed ${\lambda }_j$  and $m$  we have that the terms of the trace and spectral radius are bounded by the follows expression

  $$\begin{array}{c}Tr\left(R_k^t{R}_k\right)={\sum }_{j=1}^m{Q}_k^2\left({\lambda }_j\right){\lambda }_j\le {\sum }_{j=1}^{m^{\prime }}{\lambda }_j^{-1}+{\sum }_{j>m^{\prime }}(\tau k{)}^2{\lambda }_j\end{array}$$

  (5.5)

  and

  $$\begin{array}{c}{\rho }^2\left(R_k^t{R}_k\right)={max}_{j\in m}{Q}_k^2\left({\lambda }_j\right){\lambda }_j\le {max}_{j\le m^{\prime }}{\lambda }_j^{-1}+{max}_{j>m^{\prime }}(\tau k{)}^2{\lambda }_j.\end{array}$$

  (5.6)

  Balancing both terms in  5.5 and  5.6 gives the optimal choice of the trace and the spectral radius, respectively. Thus, we have $$Tr\left(R_k^t{R}_k\right)\approx \frac{1}{2p+1}{\left(\frac{2p+1}{2p-1}\right)}^{\frac{(2p+1)}{4p}}\frac{(\tau k{)}^{\frac{2p+1}{2p}}}{n}$$  and $${\rho }^2\left(R_k^t{R}_k\right)\approx \frac{\tau k}{n}$$  Note that the penalization term is roughly proportional to $$\frac{1}{n}\left[\frac{1}{2p+1}{\left(\frac{2p+1}{2p-1}\right)}^{\frac{(2p+1)}{4p}}(\tau k{)}^{\frac{2p+1}{2p}}+\tau k\right]$$  On the other hand $$\frac{Tr\left(R_{\alpha }^t{R}_{\alpha }\right)}{{\rho }^2\left(R_{\alpha }\right)}=\frac{1}{2p+1}{\left(\frac{2p+1}{2p-1}\right)}^{(2p+1)/4p}(\tau k{)}^{1/2p}.$$  The result then follows directly from Theorem  4.1 . □

5.3 Nonlinear multistep iterative process

Many approximate methods widely used in practice are nonlinear. We cite a important example of nonlinear approximate method. We considerer a nonlinear multistep iterative process, which have error residual $$1-\lambda Q_k\left(\lambda \right){=}^k{\prod }_{i=1}(1-{\tau }_{ik}^{-1}\lambda )$$  with ${\tau }_{ik}={\tau }_{ik}(f_0,A,y)>0,0<{\tau }_{1k}\le {\tau }_{2k}\dots \le {\tau }_{kk}\le {\lambda }_1$  . Then for $\lambda >0$  , $Q_k\left(\lambda \right)$  have the following representation $$Q_k\left(\lambda \right)={\lambda }^{-1}\left[1{-}^k{\prod }_{i=1}\right(1-{\tau }_{ik}^{-1}\lambda \left)\right]$$  The following corollary is established.

Corollary 5.4. Let ${\tau }_{ik}={\tau }_{ik}(f_0,A,y)>0,$  with $0<{\tau }_{1k}\le {\tau }_{2k}\dots \le {\tau }_{kk}\le {\lambda }_1$  . If $y\in \mathcal{ℛ}\left(A_m\right)$  , then the nonlinear multistep iterative process is an order optimal regularization method, i.e, $$IE\parallel {\tilde{f}}_m-f_{k\left(m\right)}{\parallel }^2\le 2c_1({\sum }_{i=1}^k{\tau }_{i_k}^{-1}{)}^{-2\mu }+2c_2\frac{{\sigma }^2}{n}({\sum }_{i=1}^k{\tau }_{i_k}^{-1}{)}^{(2p+1)/2p},$$  where $c_1={\rho }^2{\mu }^{\mu }(\mu +1{)}^{-1}$  and $c_2=\frac{1}{2p+1}(\frac{2p+1}{2p-1}{)}^{(2p+1)/4p}.$ 

• Proof. As before, we investigate the behavior of the bias and the variance.
  In the relation $${\lambda }^{\mu }(1-\lambda Q_k(\lambda \left)\right)={{\lambda }^{\mu }}^k{\prod }_{i=1}(1-{\tau }_{ik}^{-1}\lambda )$$  the least upper can not be reached at the points $\lambda =0$  and $\lambda ={\tau }_{1k}$  , since the estimated function is not equal to zero identically.

  On the other hand

  $$\begin{array}{c}[1-\lambda Q_k(\lambda \left){]}^{\prime }\right[\lambda Q_k\left(\lambda \right)-1{]}^{-1}={\sum }_{i=1}^k\frac{1}{{\tau }_{ik}-\lambda }\end{array}$$

  (5.7)

  Since the function in the right-hand of  5.7 does not decrease as a function of $\lambda$  on the half-interval $[0,{\tau }_{1k})$  then, the estimates of the Theorem  3.5 are valid.

  Thus, for $0\le \lambda \le {\tau }_{1k}$  , we have

  $$\begin{array}{cc}{sup}_{0\le \lambda \le {\tau }_{1k}}\left|Q_k\right(\lambda \left)\right|& =Q_k\left(0\right)={\sum }_{i=1}^k{\tau }_{i_k}^{-1}\end{array}$$

  and

  $$\begin{array}{cc}{sup}_{0\le \lambda \le {\tau }_{1k}}{\lambda }^{\mu }|1-\lambda Q_k(\lambda \left)\right|& <{\mu }^{\mu }(\mu +1{)}^{-1}({\sum }_{i=1}^k{\tau }_{i_k}^{-1}{)}^{-\mu }\end{array}$$

  Note that ${\omega }_{\mu }\left(k\right)=({\sum }_{i=1}^k{\tau }_{i_k}^{-1}{)}^{-\mu }.$  Thus, the bias is bounded by $$\parallel {\tilde{f}}_m-R_{\alpha }{A}_m\tilde{f}{\parallel }^2\le c_1({\sum }_{i=1}^k{\tau }_{i_k}^{-1}{)}^{-2\mu },$$  where $c_1={\rho }^2{\mu }^{2\mu }(\mu +1{)}^{-2}$  .

  On the order hand, it is not difficult to see that

  $$\begin{array}{cc}Tr\left(Q_k^2\right(A_m{A}_m^{*}\left)A_m{A}_m^{*}\right)& \le \frac{1}{n}\left[{\sum }_{1\le j\le m^{\prime }}{j}^{2p}+({\sum }_{i=1}^k{\tau }_{ik}^{-1}{)}^2{\sum }_{j>m^{\prime }}{j}^{-2p}\right]\le \frac{1}{n}\left[\frac{{m^{\prime }}^{2p+1}}{2p+1}+({\sum }_{i=1}^k{\tau }_{ik}^{-1}{)}^2\frac{{m^{\prime }}^{-2p+1}}{2p-1}\right]\end{array}$$

  This suggest searching $$m^{\prime }\approx c({\sum }_{i=1}^k{\tau }_{i_k}^{-1}{)}^{1/2p}$$  with $c=(\frac{2p+1}{2p+1}{)}^{1/4p}$  .

  Thus, we have that the term variance is bounded by $$IE\parallel R_{\alpha }{A}_m\tilde{f}-f_{k\left(m\right)}\parallel ={\sigma }^{2}Tr\left(Q_k^2\right(A_m^{*}{A}_m\left)A_m^{*}{A}_m\right)\le c_2\frac{{\sigma }^2}{n}({\sum }_{i=1}^k{\tau }_{i_k}^{-1}{)}^{(2p+1)/2p}$$  where $c_2=\frac{1}{2p+1}(\frac{2p+1}{2p-1}{)}^{(2p+1)/4p}$  Finally we have $$IE\parallel {\tilde{f}}_m-f_{k\left(m\right)}{\parallel }^2\le 2c_1({\sum }_{i=1}^k{\tau }_{i_k}^{-1}{)}^{-2\mu }+2c_2\frac{{\sigma }^2}{n}({\sum }_{i=1}^k{\tau }_{i_k}^{-1}{)}^{(2p+1)/2p}$$  Balancing the bias and variance terms gives the optimal choice $$IE\parallel {\tilde{f}}_m-f_{k\left(m\right)}{\parallel }^2=\mathcal{O}\left(n^{-\frac{4\mu p}{4\mu p+2p+1}}\right).$$  □

We have the following result.

Corollary 5.5. Let ${\tau }_{ik}$  be as in corollary  5.4 . Next assume $\hat{k}$  as in  4.2 and $d_{m_0}$  as in  4.1 . If $y\in \mathcal{ℛ}\left(A_m\right)$  then for any $f\in F_m$  and any $k$  , the following inequality holds true

$$\begin{array}{cc}IE\parallel {\tilde{f}}_m-{\hat{f}}_{\hat{k}}{\parallel }^2& \le \frac{1}{(1-\nu )}{inf}_{k\in \mathcal{K}}\left[C(1+\nu )\parallel {\tilde{f}}_m-f_k{\parallel }^2+\frac{2r{\sigma }^2(1+L_k)\left(c\right({\sum }_{i=1}^k{\tau }_{i_k}^{-1}{)}^{\frac{2p+1}{2p}}+{\sum }_{i=1}^k{\tau }_{i_k}^{-1})}{n}\right]\end{array}$$

$$\begin{array}{cc}& +\frac{4{\sigma }^2}{n}{\sum }_k\frac{{\sum }_{i=1}^k{\tau }_{i_k}^{-1}}{d}\left[\sqrt{dr{L}_k\left[c({\sum }_{i=1}^k{\tau }_{i_k}^{-1}{)}^{1/2p}+1\right]}+1\right]e^{-\sqrt{dr{L}_k\left[c({\sum }_{i=1}^k{\tau }_{i_k}^{-1}{)}^{1/2p}+1\right]}},\end{array}$$

for some $C>0$  and $c=\frac{1}{2p+1}(\frac{2p+1}{2p-1}{)}^{(2p+1)/4p}$  .

• Proof. First observe that $$Tr\left(R_k^t{R}_k\right)\approx \frac{1}{2p+1}{\left(\frac{2p+1}{2p-1}\right)}^{\frac{(2p+1)}{4p}}\frac{({\sum }_{i=1}^k{\tau }_{i_k}^{-1}{)}^{\frac{2p+1}{2p}}}{n}$$  and $${\rho }^2\left(R_k^t{R}_k\right)\approx \frac{{\sum }_{i=1}^k{\tau }_{i_k}^{-1}}{n}$$  Consequently $$\frac{Tr\left(R_k^t{R}_k\right)}{\rho \left(R_k\right)}\approx \frac{1}{2p+1}{\left(\frac{2p+1}{2p-1}\right)}^{\frac{(2p+1)}{4p}}({\sum }_{i=1}^k{\tau }_{i_k}^{-1}{)}^{\frac{1}{2p}}$$  Note that both $n{\rho }^2\left(R_k\right)$  and $Tr\left(R_k^t{R}_k\right)/{\rho }^2\left(R_k\right)$  do not depend on n. The proof then follows directly from theorem  4.1 .
  □

References

1. Barron A. et al., (1999). “Risk Bounds for Model Selection via Penalization ”. Probab. Theory and Related Fields. 113, pp. 467-493.
2. Birg L., and Massart P., (1978). “Minimum Contrast Estimators on Sieves: Exponential Bounds and Rates of Convergence ”. Bernoulli, Vol4.N3, pp 329-395.
3. Bousquet O., (2002). “ Concentration Inequalities for Sub-Additive Functions Using of Entropy Method.”
4. Burger, M., (2001). “A level set method for inverse problems ”. Inverse Problems 17, pp. 13271356.
5. Cavalier L., Golubev G., Picard D., and Tsybakov A., (2002). “Oracle inequalities for inverse problem ”. Ann. Statist, Vol30. N3, pp 843-874.
6. Deuflhard P., Engl H. and Scherzer O., (1998) “A convergence analysis of iterative methods for the solution of nonlinear ill-posed problems under affinely invariant conditions ”. Inverse Problems Vol14, pp. 1081-1106.
7. Dey A.K. et al, (1996) “Cross-Validation for Parameter Selection in inverse estimation problems ”. Scand. J. Statist., Vol23. N4, pp. 609-620.
8. Engl H., (1993) “Regularization methods for the stable solution of inverse problems ”. Surveys on Mathematics for Industry 3, pp. 71-143.
9. Engl H., Hanke M. and Neubauer A., (1996) “Reguralization of Inverse Problems ”. Kluwer Academic Publishers.
10. Engl H. and Grever W., (1994) “Using the L-curve for Determinig Optimal Regularization parameters ”. Numer. MAth. Vol. 69, pp. 25-31.
11. Gilyazov S.F., and Gol'dman N.L., (2000). “Reguralization of Ill-Posed Problems by Iteration Methods”. Kluwer Academic Publishers.
12. Cohen A., Hoffmann M., and Reiss M. “Adaptive Wavelet Galerkin Methods for Linear Inverse Problems”. SIAM J. Numer. Anal. Vol. 42, No. 4, pp. 1479-1501.
13. Kilmer M.E, and O'leary D.P., (2001). “Choosing Reguralization Parameters in iterative Methods for Ill-Posed Problems”. SIAM J. MATRIX ANAL. APPL. Vol. 22, N4, pp. 1204-1221.
14. Lamm P.K., (1999).“Some Recent Developments and Open Problems in Solution Methods for Mathematical Inverse Problems”. Department of Mathematics, Michihan state University, USA.
15. Ledoux, M., and Talagrand, (1996). “Deviation Inequalities for Product Measures ”. ESAIM: Probabilities and Statistics 1, pp. 63-87.
16. Loubes J.-M. and Lude $\tilde{n}$  a C., (2004). “Penalized Estimators for Nonlinear Inverse Problems. ”.
17. Loubes J.-M., and Van De Geer S., (2002). “Adaptive estimation in regression, using soft thresholding type penalties ”. Statistica Neerlandica, 56, pp 453-478.
18. Lude $\tilde{n}$  a C., and Rios,(2003).“Penalized Model Selection for Ill-posed Linear Problems. ”.
19. Mathé P., and Pereverezev S.V., (2003) “Discretization Strategy for Linear Ill-posed problems in variable Hilbert Scales ”. Inverse Problems, Vol.19, N6, pp. 1263-1279.
20. Morozov V.A., (1966). “On the Solution of Functional Equations by the Method of Regularization ”. Soviet Math. Dokl.,7, pp. 414-417.
21. F. O'sullivan., (1986).“A Statistical Perspective on Ill-Posed Inverse Problems ”. Statistical Science. Vol. 1, N4, pp. 502-527.
22. Pereverzev S. and Schock E., (2003). “On the adaptive selection of the parameter in regularization of ill-posed problems ”.
23. Tikhonov, A.,and Arsenin, V., (1977).“Solutions of Ill-Posed Problems ”. Wiley, New York.

Escuela de Matematicas, Facultad de Ciencias, UCV, Av. Los Ilustres, Los Chaguaramos, Codigo Postal 1020-A, Caracas Venezuela. Telf.: (58)212-6051481.
Departamento de Matematicas, IVIC, Carretera Panamericana KM. 11, Aptdo. 21827, Codigo Postal 1020-A, Caracas Venezuela. Telf.: (58)212-5041412.
E-mail address : afermin@euler.ciens.ucv.ve E-mail address : cludena@ivic.ve