A conjugate-gradient-type rational Krylov subspace method for ill-posed problems

TL;DR

This paper introduces a new regularisation scheme for ill-posed problems by extending conjugate-gradient methods to rational Krylov subspaces, achieving order-optimal regularisation.

Contribution

It presents a novel rational Krylov subspace method based on conjugate gradients for improved regularisation of ill-posed problems.

Findings

Order-optimal regularisation scheme established

Effective for linear inverse problems

Extends conjugate-gradient methods to rational Krylov subspaces

Abstract

Conjugated gradients on the normal equation (CGNE) is a popular method to regularise linear inverse problems. The idea of the method can be summarised as minimising the residuum over a suitable Krylov subspace. It is shown that using the same idea for the shift-and-invert rational Krylov subspace yields an order-optimal regularisation scheme.

Equations

$$Tx=y^{\delta }$$

$$T^{\ast }Tx=T^{\ast }y.$$

$$\underset{\delta \to 0}{\mathrm{limsup}}\{\Vert R_{m(\delta ,y^{\delta })}{y}^{\delta }-x^{+}\Vert |y^{\delta }\in \mathcal{Y},\Vert y^{\delta }-T{x}^{+}\Vert \le \delta \}=0$$

$$m(\delta ,y^{\delta }):=\min \{m\in {\mathbb{N}}_0|\Vert y^{\delta }-T{x}_m^{\delta }\Vert \le \tau \delta \},$$

$$\sup \{\Vert R_{m(\delta ,y^{\delta })}{y}^{\delta }-x^{+}\Vert |x^{+}\in {\mathcal{X}}_{\mu ,\rho },\Vert y^{\delta }-T{x}^{+}\Vert \le \delta \}\le C_{\mu }{\delta }^{\frac{2\mu }{2\mu +1}}{\rho }^{\frac{1}{2\mu +1}},$$

$$x_m^{\delta }=:R_m{y}^{\delta },x_m^{\delta }=\text{mbox}{argmin}_{x\in {\mathcal{K}}_m}\Vert y^{\delta }-Tx\Vert ,m=1,2,\dots ,$$

$${\mathcal{K}}_m={\mathcal{K}}_m(T^{\ast }T,T^{\ast }{y}^{\delta })=\text{mbox}span\{T^{\ast }{y}^{\delta },(T^{\ast }T)T^{\ast }{y}^{\delta },\dots ,(T^{\ast }T{)}^{m-1}{T}^{\ast }{y}^{\delta }\}.$$

$${\mathcal{Q}}_m$$

$$x_m^{\delta }=:R_m{y}^{\delta },x_m^{\delta }=\text{mbox}{argmin}_{x\in {\mathcal{Q}}_m}\Vert y^{\delta }-Tx\Vert ,m=1,2,\dots ,$$

$$(r_1,q_0)=(r_0,q_0)-{\alpha }_0(q_0,q_0)=0,{\alpha }_0=\frac{(r_0,q_0)}{(q_0,q_0)}.$$

$$0=(r_0,q_0)=(y^{\delta },T{T}^{\ast }{r}_0)=(T^{\ast }{r}_0,T^{\ast }{r}_0)$$

$$(q_1,q_0)$$

$$(r_{m+1},q_j)$$

$$0$$

$$(q_{m+1},q_m)$$

$$(q_{m+1},q_j)$$

$$(I+T^{\ast }T/\gamma {)}^{-1}{T}^{\ast }{q}_j$$

$$(q_{m+1},q_j)$$

$$t_0:=(I+T^{\ast }T/\gamma {)}^{-1}{T}^{\ast }{q}_0=c_1{w}_1+c_2{w}_0,c_1,c_2\in \mathbb{R}.$$

$${\mathcal{Q}}_{m+1}=\text{mbox}span\left\{w_0,\dots ,w_m\right\},$$

$$M:=\text{mbox}span\left\{w_0,\dots ,w_m\right\}\subseteq {\mathcal{Q}}_{m+1}.$$

$$0\ne z_m-x_m^{\delta }=\sum _{j=0}^{m-1}{\xi }_j{w}_j,{\xi }_j\in \mathbb{R},$$

$$\Vert y^{\delta }-T{z}_m{\Vert }^2$$

$$0$$

$$T^{\ast }T{x}_{\kappa }^{\delta }=T^{\ast }{y}^{\delta },x_{\kappa }^{\delta }\in {\mathcal{Q}}_{\kappa }\subseteq \mathcal{N}(T{)}^{\perp },$$

$${\mathcal{Q}}_m=\left\{r(T^{\ast }T)T^{\ast }{y}^{\delta }|r\in {\Pi }_{m-1}/(1+\cdot /\gamma {)}^{m-1}\right\}.$$

$$r_m=y^{\delta }-T{x}_m=r_m(T{T}^{\ast })y^{\delta }\text{mbox}or{T}^{\ast }{r}_m=r_m(T^{\ast }T)T^{\ast }{y}^{\delta },$$

$$r_m(\lambda )=\frac{p_m(\lambda )}{(1+\lambda /\gamma {)}^{m-1}}\text{mbox}with{p}_m(\lambda )=\prod _{j=1}^m\left(1-\frac{\lambda }{{\lambda }_{j,m}}\right),m\ge 1,$$

$$0<{\lambda }_{1,m}<{\lambda }_{1,m-1}<{\lambda }_{2,m}<\cdots <{\lambda }_{m-1,m}<{\lambda }_{m,m-1}<{\lambda }_{m,m}\le \Vert T{\Vert }^2.$$

$${\mathcal{Q}}_{\ell }=\text{mbox}span\{w_0,\cdots ,w_{\ell -1}\}=\text{mbox}span\{v_0,\cdots ,v_{\ell -1}\}\text{mbox}for\ell =1,\dots ,m.$$

Full text

• August 2019

1 Introduction

We consider the solution of the linear system

[TABLE]

where the operator $T$ acts continuously between the Hilbert spaces $\mathcal{X}$ and $\mathcal{Y}$. The linear system is assumed to be ill-posed, that is, the range $\mathcal{R}(T)$ is not closed in $\mathcal{Y}$. $y^{\delta}$ is a perturbation of the exact data $y$, such that $\|y^{\delta}-y\|\leq\delta$. $y^{\delta}$ is also called the noisy data and $\delta$ the noise level. For exact data, we assume that $y$ is in the range of $T$, $y\in\mathcal{R}(T)$, which guarantees that there exists a unique $x^{+}\in\mathcal{N}(T)^{\perp}$ such that $Tx^{+}=y$. $\mathcal{N}(T)^{\perp}$ designates the orthogonal complement of the null space $\mathcal{N}(T)$ of $T$. $x^{+}$ can also be characterised as the unique $x^{+}\in\mathcal{N}(T)^{\perp}$ that solves the normal equation

[TABLE]

In fact, the normal equation possesses a unique solution $x^{+}\in\mathcal{N}(T)^{\perp}$ for every $y\in\mathcal{D}(T^{+}):=\mathcal{R}(T)\oplus\mathcal{R}(T)^{\perp}$, where $\mathcal{R}(T)\oplus\mathcal{R}(T)^{\perp}$ designates the direct orthogonal sum of $\mathcal{R}(T)$ and its orthogonal complement $\mathcal{R}(T)^{\perp}=\mathcal{N}(T^{*})$. The linear unbounded map $T^{+}:\mathcal{D}(T^{+})\rightarrow\mathcal{N}(T)^{\perp}$, $y\mapsto x^{+}$, is the Moore–Penrose (generalised) inverse and $x^{+}$ is the minimum-norm solution.

In order to reconstruct the solution $x^{+}$ of the unperturbed problem $Tx^{+}=y$ as good as possible subject to a given noise level $\delta$, special procedures, called regularisation schemes, have to be used. Let $\{R_{m}\}_{m\in\mathbb{N}_{0}}$ be a family of linear or nonlinear operators from $\mathcal{Y}$ to $\mathcal{X}$ with $R_{m}0=0$. If there exists a mapping $m:\mathbb{R}^{+}\times\mathcal{Y}\rightarrow\mathbb{N}_{0}$ such that

[TABLE]

for any $x^{+}\in\mathcal{N}(T)^{\perp}$, then the pair $(R_{m},m(\delta,y^{\delta}))$ is called a (convergent) regularisation scheme for $T$. The mapping $m$ is called parameter choice or stopping rule. We will always use the discrepancy principle as our stopping rule, which is due to Morozov [23]. The discrepancy principle reads: Choose a fixed $\tau>1$ and set:

[TABLE]

where $x_{m}^{\delta}:=R_{m}y^{\delta}$. The discrepancy principle leads to convergent regularisation schemes (cf. [6]). Regularisation schemes might converge arbitrarily slow unless the (unperturbed) data $x^{+}$ satisfies some smoothness assumptions. Convergence rates can be given when $x^{+}$ is in the source set $\mathcal{X}_{\mu,\rho}:=\{x\in\mathcal{X}~{}|~{}x=(T^{*}T)^{\mu}w,~{}\|w\|\leq\rho\}$, $\mu>0$. Regularisation schemes $(R_{m},m(\delta,y^{\delta}))$ that attain the highest possible convergence speed are called of optimal order in $\mathcal{X}_{\mu,\rho}$ if

[TABLE]

where $C_{\mu}$ neither depends on $\delta$ nor on $\rho$.

One of the most popular iterative regularisation schemes is Conjugated Gradients on the Normal Equation (CGNE) that can be stated briefly as

[TABLE]

where $\mathcal{K}_{m}$ is the (polynomial) Krylov subspace

[TABLE]

An efficient algorithm is available to compute these approximations (cf. [15]). CGNE with the discrepancy principle as a stopping rule is an order-optimal regularisation scheme for all $\mu>0$ (cf. Theorem 7.12 in [6],[24]). And, due to its definition, CGNE is the fastest to satisfy the discrepancy principle with respect to all regularisation schemes that compute approximations in the Krylov subspace $\mathcal{K}_{m}$. The analysis of CGNE with respect to its regularisation properties is involved, since the operators $R_{m}$ are nonlinear and not necessarily continuous (cf. Theorem 7.6 in [6], [5]).

In this paper, we will define a method of the same type, but for the shift-and-invert or resolvent Krylov subspace

[TABLE]

where $\gamma>0$ is a fixed real number (e.g. [1, 11, 12, 20, 21, 22, 27, 30]). We define our method by

[TABLE]

combined with the discrepancy principle as its stopping rule. (The minimizer $x_{m}^{\delta}$ is uniquely defined, cf. lemma 2.3 below.) The subspace $\mathcal{Q}_{m}$ belongs to the class of rational Krylov subspaces which have been studied in recent years (e.g. references in [8, 13]). Since this method can be seen as solving the normal equation $(\ref{normalequation})$ approximatively in the shift-and-invert Krylov subspace $\mathcal{Q}_{m}$, the method will be called Shift-and-Invert on the Normal Equation (SINE). Several regularisation schemes have been proposed that compute approximations in the subspace $\mathcal{Q}_{m}$, cf. example 1.1. By definition, our method will be the fastest to stop with respect to the discrepancy principle. SINE is related (but not identical) to CGNE preconditioned by $(I+T^{*}T/\gamma)^{-1}$. Actually, SINE is not a preconditioning technique in the usual sense. Nevertheless, rational Krylov subspaces have been observed of being capable of accelerating the convergence (e.g. [10]). The analysis of SINE with respect to its regularisation properties shares the difficulties of the analysis of CGNE, the family of operators $R_{m}$ is again nonlinear and not continuous in general, which can be seen by generalising the ideas of the proof of theorem 7.6 in [6] and [5].

Example 1.1

Some regularisation schemes with approximations in the subspace $\mathcal{Q}_{m}$.

Iterated Tikhonov-Phillips regularisation (cf. **[3, 7, 17, 18]**). 2. 2.

Applying the implicit Euler method, the implicit midpoint-rule, or the trapezoidal rule with fixed time-step to asymptotic regularisation (Showalter’s regularisation) leads to approximations in $\mathcal{Q}_{m}$ (cf. **[28]**). 3. 3.

The method proposed by Riley in **[29]** applied to the normal equation. 4. 4.

The rational Arnoldi approach proposed in **[2]**.

In section 2, we will show that $x_{m}^{\delta}$ in $(\ref{minprop})$ can be computed efficiently and discuss some basic properties of the method. Convergence for unperturbed data is shown in section 3 before the SINE method is discussed with respect to its regularisation properties in section 4. In section 5, upper bounds on the number of iterations of SINE are discussed by comparing them to the known upper bounds on the number of iterations of CGNE. The findings are illustrated by an experiment in section 6.

Throughout we will use notations identical or closely related to the notations in [6]. In particular, the functional calculus described in section 2.3 of [6] is used without further note. Our proofs will follow closely or sometimes literally the corresponding proofs for CGNE in chapter 7 of [6].

2 Basic properties

We consider algorithm 1, where we choose $x_{0}^{\delta}=0$ without loss of generality. If $x_{0}^{\delta}$ were not zero, the corresponding shift-and-invert Krylov subspace would be spanned with $r_{0}=y^{\delta}-Tx_{0}^{\delta}$ instead of $y^{\delta}$, which allows to use prior information on the solution. First, we aim for the following properties: Algorithm 1 computes $x_{m}^{\delta}$ according to $(\ref{minprop})$, as long as the algorithm does not break down. If Algorithm 1 breaks down in step $\kappa$ with $q_{\kappa}=0$, we have $x_{\kappa}^{\delta}=T^{+}y^{\delta}$ as well as $x_{m}^{\delta}=T^{+}y^{\delta}$ for $m\geq\kappa$ in $(\ref{minprop})$.

Lemma 2.1

As long as $q_{m-1}\not=0$

• (i)

$(r_{m},q_{j})=(T^{*}r_{m},w_{j})=0$, $j=0,\ldots,m-1$

• (ii)

$\alpha_{j}\neq 0$, $j=0,\ldots,m-1$

• (iii)

$(q_{m},q_{j})=0$, $j=0,\ldots,m-1$

Proof. The proof is via induction on $m$. For $m=1$, we have

[TABLE]

Assume $\alpha_{0}=0$ then $(r_{0},q_{0})=0$ (since $(q_{0},q_{0})\neq 0$). It follows

[TABLE]

and therefore $w_{0}=T^{*}r_{0}=0$ and $q_{0}=Tw_{0}=0$. This is a contradiction to our assumption $q_{0}\neq 0$ and therefore $\alpha_{0}\neq 0$ needs to be correct. Further,

[TABLE]

which concludes the proof of our statements for $m=1$. Now we assume that the assertions are satisfied for $m$. Then, we have

[TABLE]

$\alpha_{j}\neq 0$ for $j=0,\ldots,m-1$ follows from the induction hypothesis. We show again by contradiction, that $\alpha_{m}$ can not be zero. Assume $\alpha_{m}=0$, hence $(r_{m},q_{m})=0$. We have

[TABLE]

From this, we immediately have $t_{m}=0$, $\beta_{m-1}=0$, $w_{m}=0$ and finally the contradiction $q_{m}=Tw_{m}=0$. Since this is not true, due to our assumption, $\alpha_{m}\neq 0$ is proved. Now, we target $(iii)$. It follows

[TABLE]

and for $j<m$, we have

[TABLE]

For $1\leq j<m$, it follows from the iteration and $\alpha_{j}\neq 0$ that

[TABLE]

Inserting this into $(\ref{lemmabasicprop1})$ gives

[TABLE]

For $j=0$, a simple but tedious calculation shows that

[TABLE]

With this, $(q_{m+1},q_{0})=0$ follows as above with different coefficients. $\square$

Lemma 2.2

As long as algorithm 1 does not break down with $q_{m}=0$, we have

[TABLE]

where $w_{0},\ldots,w_{m}$ is a basis of $\mathcal{Q}_{m+1}$.

Proof. It follows readily from algorithm 1 that

[TABLE]

Since $q_{0}=Tw_{0},\ldots,q_{m}=Tw_{m}$ is an orthogonal basis of $TM$, due to $(iii)$ of lemma 2.1, $w_{0},\ldots,w_{m}$ are linearly independent. Since the dimension of $\mathcal{Q}_{m+1}$ is smaller or equal to $m+1$, $w_{0},\ldots,w_{m}$ needs to be a basis of $\mathcal{Q}_{m+1}$. $\square$

Lemma 2.3

The iterates $x_{m}^{\delta}$ of algorithm 1 satisfy $(\ref{minprop})$.

Proof. Due to algorithm 1 with $x_{0}^{\delta}=0$, we have $x_{m}^{\delta}\in\mbox{span}\{w_{0},\ldots,w_{m-1}\}=\mathcal{Q}_{m}$. Now, let $z_{m}\in\mathcal{Q}_{m}=\mbox{span}\{w_{0},\ldots,w_{m-1}\}$ such that $z_{m}\neq x_{m}^{\delta}\in\mathcal{Q}_{m}$. Hence, we can write

[TABLE]

and obtain

[TABLE]

by lemma 2.1. The strict inequality, and therefore uniqueness of the minimizer, follows since $\mathcal{Q}_{m}\subseteq\mathcal{N}(T)^{\perp}$ and $\|T\sum_{j=0}^{m-1}\xi_{j}w_{j}\|^{2}>0$ as a consequence. $\square$

Lemma 2.4

If algorithm 1 breaks down in step $\kappa$ with $q_{\kappa}=0$, then $x_{\kappa}^{\delta}=x^{+}=T^{+}y^{\delta}$.

Proof. We first show that $q_{\kappa}=0$ means $T^{*}r_{\kappa}=0$. First assume $\kappa=0$. Then $0=(r_{0},q_{0})=(r_{0},TT^{*}r_{0})=\|T^{*}r_{0}\|^{2}$, hence $T^{*}r_{0}=0$. Now let $\kappa>0$ and assume $(r_{\kappa},q_{\kappa})=0$. Then, with the help of statement $(i)$ of lemma 2.1

[TABLE]

Hence we have $0=(I+T^{*}T/\gamma)t_{\kappa}=T^{*}r_{\kappa}$ in all cases. Since $r_{\kappa}=y^{\delta}-Tx_{k}^{\delta}$, this means

[TABLE]

which characterises the minimum-norm solution, that is $x_{\kappa}^{\delta}=x^{+}$ (cf. Theorem 2.5 and Theorem 2.6 in [6]). $\square$

If the algorithm stops with $q_{\kappa}=0$, it follows from $(\ref{minprop})$ that $x_{m}^{\delta}=x_{\kappa}^{\delta}$, $m\geq\kappa$.

Analogous to the description of the Krylov subspace $\mathcal{K}_{m}$ with the set $\Pi_{m-1}$ of polynomials of degree less than $m$, the shift-and-invert Krylov subspace $\mathcal{Q}_{m}$ can be described with the help of rational functions as

[TABLE]

The functional calculus of section 2.3 in [6] applies to these rational functions. The iterates, residui etc. of algorithm 1 can be identified with the corresponding rational functions (cf. [6, 14]). The following lemma describes some properties of the rational function $r_{m}(\lambda)$ that belongs to the residuum $r_{m}$, i.e., the function $r_{m}(\lambda)$ such that

[TABLE]

respectively.

Lemma 2.5

As long as the stopping index $\kappa$ has not been reached, we have

[TABLE]

where $r_{m}(\lambda)$ is the rational function that describes the residuum $r_{m}$ in algorithm 1. The values $\lambda_{j,m-1}$, $j=1,\ldots,m-1$ of $r_{m-1}(\lambda)$ and the values $\lambda_{j,m}$, $j=1,\ldots,m$ of $r_{m}(\lambda)$ are interlacing, real, and positive, that is

[TABLE]

Proof. Let $v_{0},\cdots,v_{m-1}$ be an orthonormal basis such that

[TABLE]

The basis $v_{j}$ might be obtained from the basis $w_{0},\cdots,w_{m-1}$ by the Gram–Schmidt process (or the Arnoldi process with $T^{*}y^{\delta}$). Then we can represent the iterate $x_{\ell}^{\delta}$ as $x_{\ell}^{\delta}=\sum_{j=0}^{\ell-1}z_{j,\ell}v_{j}$. Since $T^{*}y^{\delta}-T^{*}Tx_{\ell}^{\delta}\perp\mathcal{Q}_{\ell}$ is an equivalent condition to $x_{\ell}^{\delta}$ being the minimizer in $(\ref{minprop})$, we have

[TABLE]

Since $\mathcal{Q}_{\ell}\subset\mathcal{N}(T)^{\perp}=\mathcal{N}(T^{*}T)^{\perp}$, $S_{\ell}$ is invertible, and therefore symmetric positive definite with $\|S_{\ell}\|\leq\|T\|^{2}$ and $z_{\ell}=\beta S_{\ell}^{-1}e_{1}$. Specifically,

[TABLE]

Furthermore, the $(\ell-1,\ell-1)$ submatrix of $S_{\ell}$ is $S_{\ell-1}$ for $\ell=2,\ldots,m$. Inductively, by the interlacing eigenvalue theorem (cf. Theorem 3.6 in [32]), this leads to the result that the eigenvalues of $S_{\ell}$ are separated and interlaced with the eigenvalues of $S_{\ell-1}$. If we designate the eigenvalues of $S_{\ell}$ with $\lambda_{1,\ell}<\cdots<\lambda_{\ell,\ell}$ then we obtain the statement on the interlacing of the numbers in our theorem. We still have to show that the $\lambda_{j,m}$, $j=1,\ldots,m$ are the zeros of our function $r_{m}(\lambda)$. Using the notation of quasi-matrices in the standard representation of the iterate in the rational Krylov subspace (cf. [9, 31]), one obtains

[TABLE]

By comparing $(\ref{rlemma:eq1})$ with $(\ref{rlemma:eq2})$, we obtain

[TABLE]

since $S_{m}$ comes from a Krylov process and the minimal polynomial of $S_{m}$ with respect to $e_{1}$ has therefore degree $m$. Finally, we obtain

[TABLE]

by spectral decomposition. This means that $r_{m}(\lambda)=1-\lambda q_{m-1}(\lambda)/(1+\lambda/\gamma)^{m-1}$ has zeros $\lambda_{j,m}$, $j=1,\ldots,m$. Together with the obvious value $r_{m}(0)=1$, this shows the representation of $r_{m}(\lambda)$ as given in our lemma. $\square$

The inner product introduced in the following lemma will be crucial for the proof of our main theorem. Also, the idea of algorithm 1 can be briefly stated as computing an orthogonal basis $w_{0},\ldots w_{j-1}$ of the rational Krylov subspace $\mathcal{Q}_{j}$ with respect to the inner product $[\cdot,\cdot]$ in $(\ref{rperpQm})$ when the vectors are identified with the corresponding rational functions. In the following, we will often not make a difference between the residuum $r_{m}$ and the rational function $r_{m}(\lambda)$ and denote both by $r_{m}$.

Lemma 2.6

The rational functions $r_{m}$ generated by Algorithm 1 are orthogonal to $\Pi_{m-1}/(1+\cdot/\gamma)^{m-1}$ with respect to the inner product

[TABLE]

where $F_{\lambda}$ designates the spectral family of $TT^{*}$. Among all rational $\varphi\in\Pi_{m}/(1+\cdot/\gamma)^{m-1}$ with $\varphi(0)=1$, $r_{m}$ minimises the functional

[TABLE]

Proof. We have

[TABLE]

which gives the first assertion by lemma 2.1 $(i)$. The second assertion follows by lemma 2.3. $\square$

3 Convergence

The following theorem shows convergence of the iterates $x_{m}$ in algorithm 1 to the minimum-norm solution $x^{+}=T^{+}y$ for data $y\in\mathcal{D}(T^{+})$ and our general choice $x_{0}=0$. For an initial guess $x_{0}\neq 0$, it can be readily shown that the iterates converge to $T^{+}y+P_{\mathcal{N}(T)}x_{0}$, where $P_{\mathcal{N}(T)}$ is the orthogonal projector to the null space of $T$. The superscript $\delta$ has been dropped in this section in order to emphasise that data $y\in\mathcal{D}(T^{+})$ without perturbation is considered.

Theorem 3.1

The sequence of SINE iterates $\{x_{n}\}$ converge to $T^{+}y$ for all $y\in\mathcal{D}(T^{+})$.

Proof. We basically follow the lines of the proof of theorem 7.9 in [6] or [25], respectively. If the iteration terminates after a finite number of steps then the corresponding iterate coincides with $T^{+}y$ according to lemma 2.4. We therefore assume, that the iteration does not terminate. Then we have the ordering

[TABLE]

of the Ritz values according to lemma 2.5. From the representation of the residual rational function $(\ref{rdetail})$ in lemma 2.5, we obtain

[TABLE]

by the same calculations that lead to $(\ref{rmdiffhelp})$ below. Since $r_{m}/(\lambda-\lambda_{1,m})$ is in the space $\Pi_{m-1}/(1+\cdot/\gamma)^{m-1}$, the orthogonality relation $(\ref{rperpQm})$ yields

[TABLE]

which gives

[TABLE]

Since $\lambda/(\lambda-\lambda_{1,m})\geq 1$ for $\lambda\geq\lambda_{1,m}$ we obtain

[TABLE]

And therefore,

[TABLE]

Defining

[TABLE]

we obtain the estimate

[TABLE]

where $E_{\lambda}$ designates the spectral family of $T^{*}T$. For the last inequality, we additionally assumed $y\in\mathcal{R}(T)$, that is, $y=Tx^{+}$. It follows immediately, that $0\leq\varphi_{m}(\lambda)\leq 1$ for $\lambda\in[0,\lambda_{m,1}]$. We further calculate

[TABLE]

By inserting in

[TABLE]

the relation

[TABLE]

we obtain

[TABLE]

and therefore

[TABLE]

Altogether, with $\nu>0$,

[TABLE]

Since $0^{\nu}\varphi^{2}_{m}(0)=0=\lambda_{1,m}^{\nu}\varphi^{2}_{m}(\lambda_{m,1})$, there is at least one $0<\lambda^{*}<\lambda_{m,1}$, such that $(\lambda^{\nu}\varphi^{2}_{m}(\lambda))^{\prime}(\lambda^{*})=0$ and such that the maximum is achieved at this point. Hence the equation

[TABLE]

holds true. We need to distinguish two cases. First case: $\gamma\geq\|T\|^{2}$. Then, we have $0<\lambda^{*}<\lambda_{1,m}\leq\|T\|^{2}\leq\gamma$. Since $\lambda^{*}<\gamma$, we have

[TABLE]

and hence

[TABLE]

and therefore

[TABLE]

Second case: $\gamma<\|T\|^{2}$. We set

[TABLE]

We then have

[TABLE]

We can therefore conclude

[TABLE]

hence we have

[TABLE]

In both cases, we have

[TABLE]

Hence

[TABLE]

We now relax the assumption on $y$ to $y\in\mathcal{D}(T^{+})=\mathcal{R}(T)\oplus\mathcal{R}(T)^{\perp}$. Since $\mathcal{R}(T)^{\perp}=\mathcal{N}(T^{*})$, algorithm 1 produces the same iterates for $y\in\mathcal{R}(T)\oplus\mathcal{R}(T)^{\perp}$ and $P_{\overline{R(T)}}y\in\mathcal{R}(T)$, respectively. Rewriting $P_{\overline{R(T)}}y=Tx^{+}$ with $x^{+}=T^{+}y$, we can apply (19) with $\nu=1$ and obtain

[TABLE]

We now fix $0<\epsilon\leq\lambda_{1,m}$ and obtain the estimates

[TABLE]

We have to consider two cases. If $\lambda_{1,m}\rightarrow 0$ as $m\rightarrow\infty$, then one can choose $\epsilon=\epsilon_{m}=\lambda_{1,m}$. Since $|r_{m}^{\prime}(0)|\geq\lambda_{1,m}^{-1}$, we obtain

[TABLE]

In the second case, if $\lambda_{1,m}\rightarrow\lambda_{1}>0$ as $m\rightarrow\infty$, then we choose $\epsilon_{m}=|r_{m}^{\prime}(0)|^{-\frac{1}{2}}$. Since $|r_{m}^{\prime}(0)|\geq m\|T\|^{-2}$, $\epsilon_{m}\rightarrow 0$ as $m\rightarrow\infty$ and hence, $\epsilon_{m}<\lambda_{1,m}$ for $m$ sufficiently large. In this case we obtain

[TABLE]

Consequently, in any case, we have $\|x^{+}-x_{m}\|\rightarrow 0$ as $m\rightarrow\infty$. $\square$

4 SINE is an order-optimal regularisation method

We assume that

[TABLE]

where the noise level $\delta>0$ is known. Algorithm 1 is stopped with $m=m(\delta,y^{\delta})$ according to the discrepancy principle $(\ref{discprin})$. For the stopping index $m=m(\delta,y^{\delta})\geq 1$,

[TABLE]

is satisfied (with an a-priori chosen $\tau>1$), for $m=0$, only the first inequality holds true. The algorithm always terminates after a finite number of steps, which can be seen as follows. Lemma 4.1 also holds for $\mu=0$ and $\rho=\|x^{+}\|$. Hence

[TABLE]

since $|r_{m}^{\prime}(0)|^{-\frac{1}{2}}\rightarrow 0$, what we already know. The limit of the norm of the residuals exists, since the sequence is non-increasing due to lemma 2.3 or $(\ref{minprop})$, respectively, and bounded from below by zero. Since $\tau\delta>\delta$, the discrepancy principle stops the algorithm after a finite number of steps. If the algorithm has a finite termination index $\kappa$, then $q_{\kappa}=0$. According to lemma 2.4 we have $x_{\kappa}^{\delta}=T^{+}y^{\delta}$, in which case

[TABLE]

and therefore $m(\delta,y^{\delta})\leq\kappa$.

The letter $c$ designates a generic constant in the following lemmata and proofs.

Lemma 4.1

Let $y=Tx^{+}$ with $x^{+}\in\mathcal{X}_{\mu,\rho}$. Then for $0<m\leq\kappa$,

[TABLE]

Proof. The bound $(\ref{resbound})$ proved in theorem 3.1 reads

[TABLE]

As before $\varphi_{m}$ is bounded by $1$ in $[0,\lambda_{1,m}]$ and satisfies the equation $(\ref{basephiest})$ with $\nu=2\mu+1$

[TABLE]

If we insert these estimates and use $y=T(T^{*}T)^{\mu}w$ with $\|w\|\leq\rho$, we obtain

[TABLE]

which gives the assertion. $\square$

The following lemma and its proof are a nearly literal copy of Lemma 7.11 in [6], only that the functions representing the iteration are rational instead of polynomial.

Lemma 4.2

Assume that $y=Tx^{+}$ with $x^{+}\in\mathcal{X}_{\mu,\rho}$. Then for $0\leq m\leq\kappa$,

[TABLE]

where

[TABLE]

Proof. By the interpolation inequality (cf. [6], page 47) and $x_{0}^{\delta}=0$,

[TABLE]

We conclude that the assertion of the lemma is true for $m=0$ by keeping in mind that $r_{0}^{\prime}=0$. Now let $0<m\leq\kappa$. By assumption, we have

[TABLE]

and we choose a positive $\epsilon$ such that

[TABLE]

which in particular implies that $\epsilon$ is smaller than or equal to $\lambda_{1,m}$, cf. $(\ref{rd0})$. Next, we introduce

[TABLE]

where $g_{m}$ is the rational function that represents the $m$-th SINE-iterate in $\mathcal{Q}_{m}$. We obtain

[TABLE]

From here, a literal copy of the proof of Lemma 7.11 in [6] will do. The only difference is that $r_{m}$ is a rational function (cf. $(\ref{rdetail})$) instead of a polynomial. $\square$

Finally, we can prove our main theorem.

Theorem 4.3

If $y\in\mathcal{R}(T)$ and if SINE is stopped according to the discrepancy principle $(\ref{discrep})$ with $m(\delta,y^{\delta})$, then SINE is an order-optimal regularisation method, i.e., if $T^{+}y\in\mathcal{X}_{\mu,\rho}$, then

[TABLE]

Proof. By the definition of the stopping index $m(\delta,y^{\delta})$ by the discrepancy principle, one obtains

[TABLE]

With respect to lemma 4.2, it remains to estimate $|r_{m(\delta,y^{\delta})}^{\prime}(0)|$. For simplicity, we write $m$ instead of $m(\delta,y^{\delta})$ in the following and assume, without loss of generality, that $m\geq 2$. ($m=0$ follows from lemma 4.2 with $r_{0}^{\prime}=0$, $m=1$ refers to the space $\mathcal{Q}_{1}=\mathcal{K}_{1}$ and Theorem 7.12 in [6] applies). By lemma 4.1, we conclude that

[TABLE]

Since $\tau>1$, this implies that

[TABLE]

It remains to estimate

[TABLE]

The rational function

[TABLE]

satisfies

[TABLE]

due to $(\ref{rperpQm})$ in lemma 2.6. Moreover, by definition of $\pi_{m}$ and $(\ref{rd0})$,

[TABLE]

Substituting $u_{m}=\pi_{m}+\lambda\varphi$, then we have $\varphi\in\Pi_{m-2}/(1+\cdot/\gamma)^{m-1}$ and

[TABLE]

We first show that

[TABLE]

Using $(\ref{iterreprat})$, we obtain

[TABLE]

and

[TABLE]

Hence,

[TABLE]

with

[TABLE]

Since

[TABLE]

and

[TABLE]

according to lemma 2.6, we have, again with lemma 2.6,

[TABLE]

By lemma 2.6 and $(\ref{iterreprat})$, we obtain

[TABLE]

and, similarly,

[TABLE]

which finally gives $(\ref{pimab:1})$. Due to

[TABLE]

and by lemma 2.6, we obtain analogously

[TABLE]

and therefore

[TABLE]

Hence, by setting $(\ref{pimab:1})$ and $(\ref{eq9:thm})$ in $(\ref{eq8:thm})$, we obtain

[TABLE]

From here, the proof continues literally as the proof of Theorem 7.12 in [6]. $\square$

5 Upper bounds for the stopping index

The number $m(\delta,y^{\delta})$ of necessary iterations to meet the discrepancy principle reflects the efficiency of the method. Due to construction, SINE will stop faster with respect to the discrepancy principle than any other method that relies on the shift-and-invert Krylov subspace. Here, we will additionally show, that SINE stops earlier than CGNE (or at the same iterate as CGNE, in the worst case) under the same assumptions as in section 4. For this discussion, we designate the stopping index for SINE with parameter $\gamma>0$ as $m^{\gamma}(\delta,y^{\delta})$. When $\gamma$ tends to infinity, SINE turns into CGNE. Therefore we designate the stopping index of CGNE with $m^{\infty}(\delta,y^{\delta})$.

Theorem 5.1

If $y\in\mathcal{R}(T)$ and $\gamma>0$ then $0\leq m^{\gamma}(\delta,y^{\delta})\leq m^{\infty}(\delta,y^{\delta})<\infty\,.$

Proof. For $m\geq 1$, let $p_{m}\in\Pi_{m}$ be an arbitrary polynomial of degree less than or equal to $m$. With $r_{m}(\lambda)=p_{m}(\lambda)/(1+\lambda/\gamma)^{m-1}$, we have

[TABLE]

and therefore

[TABLE]

for all $p_{m}\in\mathcal{P}_{m}$ with $p_{m}(0)=1$, and, finally,

[TABLE]

Hence, by $(\ref{resrep})$, we have

[TABLE]

Due to the definition of CGNE $(\ref{minpropCGNE})$ and SINE $(\ref{minprop})$, it follows immediately that the norm of the residual of the $m$th-SINE iterate $x_{m}^{\mbox{\tiny SINE}}$ is always smaller than or equal to the norm of the residual of the $m$th-CGNE iterate $x_{m}^{\mbox{\tiny CGNE}}$, i.e.,

[TABLE]

holds for all $m\geq 0$, which proves our theorem. (The case $m=0$ is trivial.) $\square$

Theorem 5.1 basically shows that all upper bounds that are known for CGNE also apply to SINE, but SINE might be faster. We explicitly state some corollaries. As a corollary of theorem 7.13 in [6], which is due to [25, 26], and theorem 5.1 we obtain the following statement.

Corollary 5.2

If $y\in\mathcal{R}(T)$, $\gamma\in\mathbb{R}^{+}\cup\{\infty\}$, and $T^{+}y\in\mathcal{X}_{\mu,\rho}$, then

[TABLE]

and this estimate is sharp in the sense that the exponent cannot be replaced by a smaller one and that the bound is supposed to hold true for all possible values of $\gamma$.

Theorem 7.14 and theorem 7.15 in [6] also hold literally for SINE as simple corollaries of theorem 5.1

6 Illustration and discussion

We use the multiplication operator $T:L_{2}(0,1)\rightarrow L_{2}(0,1)$, $Tf(t):=tf(t)$, in order to illustrate the theoretical findings. The range of $T$ is not closed. For example, it can be readily seen that any constant function apart from zero is in the closure $\overline{\mathcal{R}(T)}$ of the range $\mathcal{R}(T)$, but not in the range $\mathcal{R}(T)$ itself. We further have $T^{*}=T$ and $T^{*}Tf(t)=t^{2}f(t)$. For the fractional powers of the operator $T^{*}T$, we obtain $(T^{*}T)^{\mu}f(t)=t^{2\mu}f(t)$, $\mu>0$. We choose the exact solutions $x^{+}_{1}=t$ and $x^{+}_{2}=t^{3}$ with right-hand sides $y_{1}=Tx_{1}^{+}=t^{2}$ and $y_{2}=Tx_{2}^{+}=t^{4}$, respectively. Then $x_{1}^{+}\in\mathcal{X}_{1/2,1}$ and $x_{2}^{+}\in\mathcal{X}_{3/2,1}$. We use the perturbed right-hand sides $y_{i}^{\delta}=y_{i}+\delta$, $i=1,2$, with $\|y_{i}^{\delta}-y_{i}\|=\delta$. The linear systems with the perturbed right-hand sides $y_{i}^{\delta}$, $i=1,2$, do not have a solution and a regularisation is necessary. We now use SINE to compute regularised solutions, where the iteration is stopped according to the discrepancy principle with $\tau=1001/1000$. In figure 1, the $L_{2}$-norm of the error of the computed regularisation is plotted versus $\delta^{-1}$. The red circle-marked line belongs to the error with respect to $x_{1}^{+}=t$ and the green, square-marked line belongs to the error with respect to $x_{2}^{+}=t^{3}$. As predicted by theorem 4.3, the convergence to the exact solution with decreasing perturbation $\delta$ is at least $\delta^{\frac{1}{2}}$ or $\delta^{\frac{3}{2}}$, respectively, which are indicated by the gray lines. The operator is simple enough such that all computations could be conducted exactly by using the computer algebra system Maple.

Due to construction, SINE will stop faster with respect to the discrepancy principle than any other method computing regularisations in the shift-and-invert Krylov subspace $\mathcal{Q}_{m}$. That is, where these methods have been used successfully, SINE should also be a very good choice. That SINE might also be useful with respect to regularisation schemes that do not use the shift-and-invert Krylov subspace $\mathcal{Q}_{m}$, will be illustrated by another simple experiment, where we compare SINE and CGNE. For $\gamma=\frac{1}{1000}$ and $x^{+}=t$, $y^{\delta}=t^{2}+\delta$, $\delta=\frac{1}{1000}$, SINE stops after two steps with

[TABLE]

$c_{1}=-21/5000$, $c_{2}=15070063/15000$, whereas CGNE produces a polynomial of degree $39$ after $19$ steps. Both methods have been stopped according to the discrepancy principle with $\tau=\frac{1001}{1000}$. In figure 2, it can be seen that the SINE regularisation on the left-hand side is qualitatively better than the CGNE regularisation on the right-hand side. The experiment also shows that the stopping index of SINE can be significantly smaller than the stopping index of CGNE. With the designations of section 5, we have

[TABLE]

Altogether, the theory and the experiment suggest that SINE is a valid order-optimal regularisation scheme. It is also hoped, that the given analysis inspires further research on the regularisation properties of rational Krylov subspace methods. For example, it is immediately clear that theorem 5.1 carries over to rational Krylov subspaces with arbitrary negative real poles, when the method is defined analogous to $(\ref{minprop})$. Even choosing negative poles at random can only improve on CGNE with respect to the stopping index. These more general rational Krylov subspace methods might also be seen as accelerations of the nonstationary iterated Tikhonov iteration (e.g. [16]) or of method ${\sl(ii)}$ in example 1.1 with varying step sizes. While there is only one polynomial Krylov subspace, rational Krylov subspaces inspire a wide range of methods that might be adapted to the needs at hand. As a possible application, rational Krylov subspaces have been successfully used to accelerate computations related to seismic imaging (e.g. [4, 19, 33, 34]), which is known to be an ill-posed inverse problem.

Acknowledgments

Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 258734477 – SFB 1173 *