# On the convergence analysis of one-shot inversion methods

Marcella Bonazzoli, Houssem Haddar, Tuan Anh Vu

<figure><img src="image_1.png" /></figure>

## To cite this version:

Marcella Bonazzoli, Houssem Haddar, Tuan Anh Vu. On the convergence analysis of one-shot inversion methods. 2024. hal-04151014v2

HAL Id: hal-04151014

https://inria.hal.science/hal-04151014v2

Preprint submitted on 10 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the convergence analysis of one-shot inversion methods

Marcella Bonazzoli, Houssem Haddar, Tuan Anh Vu

<figure><img src="image_1.png" /></figure>

## To cite this version:

Marcella Bonazzoli, Houssem Haddar, Tuan Anh Vu. On the convergence analysis of one-shot inversion methods. SIAM Journal on Applied Mathematics, 2024, 84 (6), pp.2440-2475. 10.1137/23M1585866. hal-04151014v2

HAL Id: hal-04151014

https://inria.hal.science/hal-04151014v2

Submitted on 10 Apr 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the convergence analysis of one-shot inversion methods

Marcella Bonazzoli¹, Houssem Haddar¹, and Tuan-Anh Vu¹,²

¹Inria, UMA, ENSTA Paris, Institut Polytechnique de Paris, 91120 Palaiseau, France

²Institute of Mathematics, Vietnam Academy of Science and Technology, Vietnam

Revised and resubmitted on April 10, 2024

## Abstract

When an inverse problem is solved by a gradient-based optimization algorithm, the corresponding forward and adjoint problems, which are introduced to compute the gradient, can be also solved iteratively. The idea of iterating at the same time on the inverse problem unknown and on the forward and adjoint problem solutions yields the concept of one-shot inversion methods. We are especially interested in the case where the inner iterations for the direct and adjoint problems are incomplete, that is, stopped before achieving a high accuracy on their solutions. Here, we focus on general linear inverse problems and generic fixed-point iterations for the associated forward problem. We analyze variants of the so-called multi-step one-shot methods, in particular semi-implicit schemes with a regularization parameter. We establish sufficient conditions on the descent step for convergence, by studying the eigenvalues of the block matrix of the coupled iterations. Several numerical experiments are provided to illustrate the convergence of these methods in comparison with the classical gradient descent, where the forward and adjoint problems are solved exactly by a direct solver instead. We observe that very few inner iterations are enough to guarantee good convergence of the inversion algorithm, even in the presence of noisy data.

**Keywords.** inverse problems, one-shot methods, convergence analysis, parameter identification

**MSC codes.** 35R30, 65F10, 65N21

## 1 Introduction

For large-scale inverse problems, which often arise in real life applications, the solution of the corresponding forward and adjoint problems is generally computed using an iterative solver, such as preconditioned fixed point or Krylov subspace methods, rather than exactly by a direct solver, such as LU-type solvers (see e.g. [27, 1]). Indeed, the corresponding linear systems could be too large to be handled with direct solvers because of their high memory requirement. In addition, iterative solvers are easier to parallelize on many cores for time speed-up. By coupling the iterative solver with a gradient-based optimization iteration, the idea of *one-step one-shot methods* is to iterate at the same time on the forward problem solution (the state variable), the adjoint problem solution (the adjoint state) and on the inverse problem unknown (the parameter or design variable). If two or more inner iterations are performed on the state and adjoint state before updating the parameter (by starting from the previous iterates as initial guess for the state and adjoint state), we speak of *multi-step one-shot methods*. Our goal is to rigorously analyze the convergence of such inversion methods. In particular, we are interested in those schemes where the inner iterations on the direct and adjoint problems are incomplete, i.e. stopped before achieving convergence. Indeed, solving the forward and adjoint problems exactly by direct solvers or very accurately by iterative solvers could be very time-consuming with little improvement in the accuracy of the inverse problem solution.
# On the convergence analysis of one-shot inversion methods

Marcella Bonazzoli¹, Houssem Haddar¹, and Tuan-Anh Vu¹,²

¹Inria, UMA, ENSTA Paris, Institut Polytechnique de Paris, 91120 Palaiseau, France

²Institute of Mathematics, Vietnam Academy of Science and Technology, Vietnam

Revised and resubmitted on April 10, 2024

## Abstract

When an inverse problem is solved by a gradient-based optimization algorithm, the corresponding forward and adjoint problems, which are introduced to compute the gradient, can be also solved iteratively. The idea of iterating at the same time on the inverse problem unknown and on the forward and adjoint problem solutions yields the concept of one-shot inversion methods. We are especially interested in the case where the inner iterations for the direct and adjoint problems are incomplete, that is, stopped before achieving a high accuracy on their solutions. Here, we focus on general linear inverse problems and generic fixed-point iterations for the associated forward problem. We analyze variants of the so-called multi-step one-shot methods, in particular semi-implicit schemes with a regularization parameter. We establish sufficient conditions on the descent step for convergence, by studying the eigenvalues of the block matrix of the coupled iterations. Several numerical experiments are provided to illustrate the convergence of these methods in comparison with the classical gradient descent, where the forward and adjoint problems are solved exactly by a direct solver instead. We observe that very few inner iterations are enough to guarantee good convergence of the inversion algorithm, even in the presence of noisy data.

**Keywords.** inverse problems, one-shot methods, convergence analysis, parameter identification

**MSC codes.** 35R30, 65F10, 65N21

## 1 Introduction

For large-scale inverse problems, which often arise in real life applications, the solution of the corresponding forward and adjoint problems is generally computed using an iterative solver, such as preconditioned fixed point or Krylov subspace methods, rather than exactly by a direct solver, such as LU-type solvers (see e.g. [27, 1]). Indeed, the corresponding linear systems could be too large to be handled with direct solvers because of their high memory requirement. In addition, iterative solvers are easier to parallelize on many cores for time speed-up. By coupling the iterative solver with a gradient-based optimization iteration, the idea of *one-step one-shot methods* is to iterate at the same time on the forward problem solution (the state variable), the adjoint problem solution (the adjoint state) and on the inverse problem unknown (the parameter or design variable). If two or more inner iterations are performed on the state and adjoint state before updating the parameter (by starting from the previous iterates as initial guess for the state and adjoint state), we speak of *multi-step one-shot methods*. Our goal is to rigorously analyze the convergence of such inversion methods. In particular, we are interested in those schemes where the inner iterations on the direct and adjoint problems are incomplete, i.e. stopped before achieving convergence. Indeed, solving the forward and adjoint problems exactly by direct solvers or very accurately by iterative solvers could be very time-consuming with little improvement in the accuracy of the inverse problem solution.
The concept of one-shot methods was first introduced by Ta'asan [24] for optimal control prob-
lems. Based on this idea, a variety of related methods, such as the all-at-once methods, where
the state equation is included in the misfit functional, were developed for aerodynamic shape op-
timization, see for instance [25, 23, 13, 22, 21] and the literature review in the introduction of
[22]. All-at-once approaches to inverse problems for parameter identification were studied in, e.g.,
[9, 3, 17, 15, 16, 19, 20]. An alternative method, called Wavefield Reconstruction Inversion (WRI),
was introduced for seismic imaging in [28], as an improvement of the classical Full Waveform In-
version (FWI) [26]. WRI is a penalty method which combines the advantages of the all-at-once
approach with those of the reduced approach (where the state equation represents a constraint and
is enforced at each iteration, as in FWI), and was extended to more general inverse problems in [29].

Few convergence proofs, especially for the multi-step one-shot methods, are available in the literature. In particular, for non-linear design optimization problems, Griewank [7] proposed a version of one-step one-shot methods where a Hessian-based preconditioner is used in the design variable iteration. The author proved conditions to ensure that the real eigenvalues of the Jacobian of the coupled iterations are smaller than 1, but these are just necessary and not sufficient conditions to exclude real eigenvalues smaller than −1. In addition, no condition to also bound complex eigenvalues below 1 in modulus was found, and multi-step methods were not investigated. In [10, 11, 5] an exact penalty function of doubly augmented Lagrangian type was introduced to coordinate the coupled iterations, and global convergence of the proposed optimization approach was proved under some assumptions. This particular one-step one-shot approach was later extended to time-dependent problems in [8].

In this work, we consider variants of multi-step one-shot methods where the forward and adjoint
problems are solved using fixed point iterations and the inverse problem is solved using gradient
descent methods. In particular, we analyze semi-implicit schemes with a regularization parameter.
This is a preliminary work where we focus on (discretized) linear inverse problems. The only basic
assumptions we require are the uniqueness of the inverse problem solution and the convergence of the
fixed point iteration for the forward problem. To analyze the convergence of the coupled iterations
we study the real and complex eigenvalues of the block iteration matrix. We prove that if the descent
step is small enough, then the considered one-shot methods converge. Moreover, the upper bounds
for the descent step in this sufficient condition are explicit in the number of inner iterations, in the
norms of the operators involved in the problem, and in the regularization parameter. Note that
previously in our research report [2] we studied (shifted) explicit schemes with no regularization,
while here we include a regularized cost functional and analyze semi-implicit schemes. Moreover, for
the particular scalar case, in [2] we established sufficient and also necessary convergence conditions
on the descent step. In the literature [7, 10, 11, 5] mentioned above, only the one-step one-shot
algorithm has been considered and the analysis for that scheme was far from complete as it was
only known that the real eigenvalues of the iteration system has modulus less than one. In the
current work, we study the eigenvalues in the whole complex plane and consider also the case of
multiple inner iterations on the state variable. In addition, we include the possibility of using some
regularization terms for the inverse variable.

This paper is structured as follows. In Section 2, we introduce the principle of multi-step one-
shot methods and define two variants of these algorithms. Since their analysis is similar, we shall
focus on one of them, namely the semi-implicit scheme. Then, in Section 3, respectively Section 4,
we analyze the convergence of one-step one-shot methods, respectively multi-step one-shot methods:
first, we establish eigenvalue equations for the block matrices of the coupled iterations, then we derive
sufficient convergence conditions on the descent step by studying the location of the eigenvalues in
the complex plane. Finally, in Section 5 we do several numerical tests on the performance of the
different algorithms on a toy 2D Helmholtz inverse problem. These tests are carried out in four cases.
In the first case, the measurements are noise-free and the parameter that we desire to reconstruct
is discretized in a low-dimensional space, while in the second case, the measurements are affected
The concept of one-shot methods was first introduced by Ta'asan [24] for optimal control prob-
lems. Based on this idea, a variety of related methods, such as the all-at-once methods, where
the state equation is included in the misfit functional, were developed for aerodynamic shape op-
timization, see for instance [25, 23, 13, 22, 21] and the literature review in the introduction of
[22]. All-at-once approaches to inverse problems for parameter identification were studied in, e.g.,
[9, 3, 17, 15, 16, 19, 20]. An alternative method, called Wavefield Reconstruction Inversion (WRI),
was introduced for seismic imaging in [28], as an improvement of the classical Full Waveform In-
version (FWI) [26]. WRI is a penalty method which combines the advantages of the all-at-once
approach with those of the reduced approach (where the state equation represents a constraint and
is enforced at each iteration, as in FWI), and was extended to more general inverse problems in [29].

Few convergence proofs, especially for the multi-step one-shot methods, are available in the literature. In particular, for non-linear design optimization problems, Griewank [7] proposed a version of one-step one-shot methods where a Hessian-based preconditioner is used in the design variable iteration. The author proved conditions to ensure that the real eigenvalues of the Jacobian of the coupled iterations are smaller than 1, but these are just necessary and not sufficient conditions to exclude real eigenvalues smaller than −1. In addition, no condition to also bound complex eigenvalues below 1 in modulus was found, and multi-step methods were not investigated. In [10, 11, 5] an exact penalty function of doubly augmented Lagrangian type was introduced to coordinate the coupled iterations, and global convergence of the proposed optimization approach was proved under some assumptions. This particular one-step one-shot approach was later extended to time-dependent problems in [8].

In this work, we consider variants of multi-step one-shot methods where the forward and adjoint
problems are solved using fixed point iterations and the inverse problem is solved using gradient
descent methods. In particular, we analyze semi-implicit schemes with a regularization parameter.
This is a preliminary work where we focus on (discretized) linear inverse problems. The only basic
assumptions we require are the uniqueness of the inverse problem solution and the convergence of the
fixed point iteration for the forward problem. To analyze the convergence of the coupled iterations
we study the real and complex eigenvalues of the block iteration matrix. We prove that if the descent
step is small enough, then the considered one-shot methods converge. Moreover, the upper bounds
for the descent step in this sufficient condition are explicit in the number of inner iterations, in the
norms of the operators involved in the problem, and in the regularization parameter. Note that
previously in our research report [2] we studied (shifted) explicit schemes with no regularization,
while here we include a regularized cost functional and analyze semi-implicit schemes. Moreover, for
the particular scalar case, in [2] we established sufficient and also necessary convergence conditions
on the descent step. In the literature [7, 10, 11, 5] mentioned above, only the one-step one-shot
algorithm has been considered and the analysis for that scheme was far from complete as it was
only known that the real eigenvalues of the iteration system has modulus less than one. In the
current work, we study the eigenvalues in the whole complex plane and consider also the case of
multiple inner iterations on the state variable. In addition, we include the possibility of using some
regularization terms for the inverse variable.

This paper is structured as follows. In Section 2, we introduce the principle of multi-step one-
shot methods and define two variants of these algorithms. Since their analysis is similar, we shall
focus on one of them, namely the semi-implicit scheme. Then, in Section 3, respectively Section 4,
we analyze the convergence of one-step one-shot methods, respectively multi-step one-shot methods:
first, we establish eigenvalue equations for the block matrices of the coupled iterations, then we derive
sufficient convergence conditions on the descent step by studying the location of the eigenvalues in
the complex plane. Finally, in Section 5 we do several numerical tests on the performance of the
different algorithms on a toy 2D Helmholtz inverse problem. These tests are carried out in four cases.
In the first case, the measurements are noise-free and the parameter that we desire to reconstruct
is discretized in a low-dimensional space, while in the second case, the measurements are affected
by noise and the reconstructed parameter is discretized in a higher dimensional space. In the third and fourth cases, we test the robustness with respect to the size of the discretized problem, and the dependence on the norm of the iteration matrix of the forward problem. In particular, we observe that very few inner iterations are enough to guarantee good convergence of the inversion algorithms, even in the presence of noisy data.

Throughout this work, $\langle \cdot, \cdot \rangle$ indicates the usual Hermitian scalar product in $\mathbb{C}^n$, that is $\langle x, y \rangle := \bar{y}^T x, \forall x, y \in \mathbb{C}^n$, and $\| \cdot \|$ the vector/matrix norms induced by $\langle \cdot, \cdot \rangle$. We denote by $A^* = \bar{A}^T$ the adjoint operator of a matrix $A \in \mathbb{C}^{m \times n}$, and likewise by $z^* = \bar{z}$ the conjugate of a complex number $z$. The identity matrix is always denoted by $I$, whose size is understood from context. Finally, for a matrix $T \in \mathbb{C}^{m \times n}$, we denote by $\rho(T)$ the spectral radius of $T$, and when $\rho(T) < 1$, we define

$$s(T) := \sup_{z \in \mathbb{C}, |z| \ge 1} \|(I - T/z)^{-1}\|,$$

which is further studied in Appendix A.

## 2 Multi-step one-shot inversion methods

We focus on (discretized) linear inverse problems, which correspond to a direct (or forward) problem of the form: seek $u \equiv u(\sigma)$ such that

$$u = Bu + M\sigma + F \quad (1)$$

where $u \in \mathbb{R}^{n_u}$, $\sigma \in \mathbb{R}^{n_\sigma}$, $B \in \mathbb{R}^{n_u \times n_u}$, $M \in \mathbb{R}^{n_u \times n_\sigma}$ and $F \in \mathbb{R}^{n_u}$. Here $I - B$ is the invertible matrix associated with the direct problem (e.g. obtained after discretization of a PDE model and applying a preconditioner to the linear system), with parameter $\sigma$. Equation (1) is also referred to as the state equation and $u$ is the state variable. Given $\sigma$, we assume that one solves for $u$ by a fixed point iteration

$$u_{\ell+1} = Bu_{\ell} + M\sigma + F, \quad \ell = 0, 1, \dots \quad (2)$$

We indeed assume $\rho(B) < 1$ so that the fixed point iteration (2) converges for any initial guess $u_0$ (see e.g. [6, Theorem 2.1.1]). Measuring $g = Hu(\sigma)$, where $H \in \mathbb{R}^{n_g \times n_u}$, we consider the linear inverse problem of retrieving $\sigma$ from the knowledge of $g$. Let us set $A := H(I - B)^{-1}M$. The inverse problem can be synthetically written as $A\sigma = g - H(I - B)^{-1}F$, which amounts to inverting the ill-conditioned matrix $A$. We shall assume in the following the uniqueness of the solution for this inverse problem, which is equivalent to the injectivity of $A$. In summary, we set

$$\begin{array}{ll} \text{direct problem:} & u = Bu + M\sigma + F, \\ \text{inverse problem:} & \text{measure } g = Hu(\sigma), \text{ retrieve } \sigma \end{array} \quad (3)$$

with the assumptions:

$$\rho(B) < 1, \quad H(I - B)^{-1}M \text{ is injective.} \quad (4)$$

*Remark 2.1.* Considering real-valued matrices $B$ and $M$ is not a restrictive assumption. Indeed, the case of complex-valued matrices can be rewritten as a system of real-valued equations by doubling the size of the linear system (see [2, Section 5]).

To solve the inverse problem we write its regularized least squares formulation: given some measurements $g$, we seek the regularized solution $\sigma_\alpha^{\text{ex}}$ defined by

$$\sigma_{\alpha}^{\text{ex}} = \underset{\sigma \in \mathbb{R}^{n_{\sigma}}}{\operatorname{argmin}} J(\sigma) \quad \text{where } J(\sigma) := \frac{1}{2} \|Hu(\sigma) - g\|^2 + \frac{\alpha}{2} \|\sigma\|^2, \alpha \ge 0. \quad (5)$$

In case we have exact measurements, i.e. $g = Hu(\sigma^{\text{ex}})$, one can take $\alpha = 0$, so that $\sigma_\alpha^{\text{ex}} = \sigma^{\text{ex}}$. However, in practice, the measurements are often affected by noise, hence in general we cannot
by noise and the reconstructed parameter is discretized in a higher dimensional space. In the third and fourth cases, we test the robustness with respect to the size of the discretized problem, and the dependence on the norm of the iteration matrix of the forward problem. In particular, we observe that very few inner iterations are enough to guarantee good convergence of the inversion algorithms, even in the presence of noisy data.

Throughout this work, $\langle \cdot, \cdot \rangle$ indicates the usual Hermitian scalar product in $\mathbb{C}^n$, that is $\langle x, y \rangle := \bar{y}^T x, \forall x, y \in \mathbb{C}^n$, and $\| \cdot \|$ the vector/matrix norms induced by $\langle \cdot, \cdot \rangle$. We denote by $A^* = \bar{A}^T$ the adjoint operator of a matrix $A \in \mathbb{C}^{m \times n}$, and likewise by $z^* = \bar{z}$ the conjugate of a complex number $z$. The identity matrix is always denoted by $I$, whose size is understood from context. Finally, for a matrix $T \in \mathbb{C}^{m \times n}$, we denote by $\rho(T)$ the spectral radius of $T$, and when $\rho(T) < 1$, we define

$$s(T) := \sup_{z \in \mathbb{C}, |z| \ge 1} \|(I - T/z)^{-1}\|,$$

which is further studied in Appendix A.

## 2 Multi-step one-shot inversion methods

We focus on (discretized) linear inverse problems, which correspond to a direct (or forward) problem of the form: seek $u \equiv u(\sigma)$ such that

$$u = Bu + M\sigma + F \quad (1)$$

where $u \in \mathbb{R}^{n_u}$, $\sigma \in \mathbb{R}^{n_\sigma}$, $B \in \mathbb{R}^{n_u \times n_u}$, $M \in \mathbb{R}^{n_u \times n_\sigma}$ and $F \in \mathbb{R}^{n_u}$. Here $I - B$ is the invertible matrix associated with the direct problem (e.g. obtained after discretization of a PDE model and applying a preconditioner to the linear system), with parameter $\sigma$. Equation (1) is also referred to as the state equation and $u$ is the state variable. Given $\sigma$, we assume that one solves for $u$ by a fixed point iteration

$$u_{\ell+1} = Bu_{\ell} + M\sigma + F, \quad \ell = 0, 1, \dots \quad (2)$$

We indeed assume $\rho(B) < 1$ so that the fixed point iteration (2) converges for any initial guess $u_0$ (see e.g. [6, Theorem 2.1.1]). Measuring $g = Hu(\sigma)$, where $H \in \mathbb{R}^{n_g \times n_u}$, we consider the linear inverse problem of retrieving $\sigma$ from the knowledge of $g$. Let us set $A := H(I - B)^{-1}M$. The inverse problem can be synthetically written as $A\sigma = g - H(I - B)^{-1}F$, which amounts to inverting the ill-conditioned matrix $A$. We shall assume in the following the uniqueness of the solution for this inverse problem, which is equivalent to the injectivity of $A$. In summary, we set

$$\begin{array}{ll} \text{direct problem:} & u = Bu + M\sigma + F, \\ \text{inverse problem:} & \text{measure } g = Hu(\sigma), \text{ retrieve } \sigma \end{array} \quad (3)$$

with the assumptions:

$$\rho(B) < 1, \quad H(I - B)^{-1}M \text{ is injective.} \quad (4)$$

*Remark 2.1.* Considering real-valued matrices $B$ and $M$ is not a restrictive assumption. Indeed, the case of complex-valued matrices can be rewritten as a system of real-valued equations by doubling the size of the linear system (see [2, Section 5]).

To solve the inverse problem we write its regularized least squares formulation: given some measurements $g$, we seek the regularized solution $\sigma_\alpha^{\text{ex}}$ defined by

$$\sigma_{\alpha}^{\text{ex}} = \underset{\sigma \in \mathbb{R}^{n_{\sigma}}}{\operatorname{argmin}} J(\sigma) \quad \text{where } J(\sigma) := \frac{1}{2} \|Hu(\sigma) - g\|^2 + \frac{\alpha}{2} \|\sigma\|^2, \alpha \ge 0. \quad (5)$$

In case we have exact measurements, i.e. $g = Hu(\sigma^{\text{ex}})$, one can take $\alpha = 0$, so that $\sigma_\alpha^{\text{ex}} = \sigma^{\text{ex}}$. However, in practice, the measurements are often affected by noise, hence in general we cannot
reconstruct $\sigma^{\text{ex}}$ but its approximation $\sigma_{\alpha}^{\text{ex}}$. Now we solve the minimization problem in (5). Using the classical Lagrangian technique with real scalar products, we introduce the adjoint state $p \equiv p(\sigma)$, which is the solution of

$$p = B^*p + H^*(Hu - g)$$

and allows us to compute the gradient of the cost functional as

$$\nabla J(\sigma) = M^*p(\sigma) + \alpha\sigma.$$

The classical gradient descent algorithm then reads

$$\text{usual gradient descent:} \quad \begin{cases} \sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^n, \\ u^n = B u^n + M \sigma^n + F, \\ p^n = B^* p^n + H^*(H u^n - g), \end{cases} \quad (6)$$

where $\tau > 0$ is the descent step size, and the state and adjoint state equations are solved exactly at each iteration step for $\sigma$. Notice that when $F = 0$, (6) is equivalent to $\sigma^{n+1} = \sigma^n - \tau A^*(A\sigma^n - g) - \tau\alpha\sigma^n$. One can also consider a slightly different version where an implicit scheme is applied to the regularization term leading to the following semi-implicit gradient scheme

$$\text{semi-implicit gradient descent:} \quad \begin{cases} \sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^{n+1}, \\ u^n = B u^n + M \sigma^n + F, \\ p^n = B^* p^n + H^*(H u^n - g). \end{cases} \quad (7)$$

It can be shown that both algorithms converge for sufficiently small $\tau > 0$: for any initial guess, (6) converges if and only if $\tau < \frac{2}{\rho(A^*A)+\alpha}$ and (7) converges if and only if $(\rho(A^*A) - \alpha)\tau < 2$. These results indicate in particular that we gain more stability with the semi-implicit scheme.

We are interested in methods where the direct and adjoint problems are rather solved iteratively as in (2), and where we iterate at the same time on the forward problem solution and the inverse problem unknown: such methods are called one-shot methods. More precisely, we are interested in two variants of multi-step one-shot methods, defined as follows. Let $n$ be the index of the (outer) iteration on $\sigma$. We update $\sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^n$ as in gradient descent methods (or respectively, $\sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^{n+1}$ as in semi-implicit gradient descent methods), but the state and adjoint state equations are now solved by a fixed point iteration method, using just $k$ inner iterations and as initial guess we naturally choose the information from the previous (outer) step. We then get the two following variants of multi-step one-shot algorithms

$$\text{k-step one-shot:} \quad \begin{cases} \sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^n, \\ u_0^{n+1} = u^n, p_0^{n+1} = p^n, \\ \text{for } \ell = 0, 1, \dots, k-1: \\ \quad \left| \begin{array}{l} u_{\ell+1}^{n+1} = B u_{\ell}^{n+1} + M \sigma^{n+1} + F, \\ p_{\ell+1}^{n+1} = B^* p_{\ell}^{n+1} + H^*(H u_{\ell}^{n+1} - g), \\ u^{n+1} := u_k^{n+1}, p^{n+1} := p_k^{n+1} \end{array} \right. \end{cases} \quad (8)$$

and

$$\text{semi-implicit k-step one-shot:} \quad \begin{cases} \sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^{n+1}, \\ u_0^{n+1} = u^n, p_0^{n+1} = p^n, \\ \text{for } \ell = 0, 1, \dots, k-1: \\ \quad \left| \begin{array}{l} u_{\ell+1}^{n+1} = B u_{\ell}^{n+1} + M \sigma^{n+1} + F, \\ p_{\ell+1}^{n+1} = B^* p_{\ell}^{n+1} + H^*(H u_{\ell}^{n+1} - g), \\ u^{n+1} := u_k^{n+1}, p^{n+1} := p_k^{n+1} \end{array} \right. \end{cases} \quad (9)$$
reconstruct $\sigma^{\text{ex}}$ but its approximation $\sigma_{\alpha}^{\text{ex}}$. Now we solve the minimization problem in (5). Using the classical Lagrangian technique with real scalar products, we introduce the adjoint state $p \equiv p(\sigma)$, which is the solution of

$$p = B^*p + H^*(Hu - g)$$

and allows us to compute the gradient of the cost functional as

$$\nabla J(\sigma) = M^*p(\sigma) + \alpha\sigma.$$

The classical gradient descent algorithm then reads

$$\text{usual gradient descent:} \quad \begin{cases} \sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^n, \\ u^n = B u^n + M \sigma^n + F, \\ p^n = B^* p^n + H^*(H u^n - g), \end{cases} \quad (6)$$

where $\tau > 0$ is the descent step size, and the state and adjoint state equations are solved exactly at each iteration step for $\sigma$. Notice that when $F = 0$, (6) is equivalent to $\sigma^{n+1} = \sigma^n - \tau A^*(A\sigma^n - g) - \tau\alpha\sigma^n$. One can also consider a slightly different version where an implicit scheme is applied to the regularization term leading to the following semi-implicit gradient scheme

$$\text{semi-implicit gradient descent:} \quad \begin{cases} \sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^{n+1}, \\ u^n = B u^n + M \sigma^n + F, \\ p^n = B^* p^n + H^*(H u^n - g). \end{cases} \quad (7)$$

It can be shown that both algorithms converge for sufficiently small $\tau > 0$: for any initial guess, (6) converges if and only if $\tau < \frac{2}{\rho(A^*A)+\alpha}$ and (7) converges if and only if $(\rho(A^*A) - \alpha)\tau < 2$. These results indicate in particular that we gain more stability with the semi-implicit scheme.

We are interested in methods where the direct and adjoint problems are rather solved iteratively as in (2), and where we iterate at the same time on the forward problem solution and the inverse problem unknown: such methods are called one-shot methods. More precisely, we are interested in two variants of multi-step one-shot methods, defined as follows. Let $n$ be the index of the (outer) iteration on $\sigma$. We update $\sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^n$ as in gradient descent methods (or respectively, $\sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^{n+1}$ as in semi-implicit gradient descent methods), but the state and adjoint state equations are now solved by a fixed point iteration method, using just $k$ inner iterations and as initial guess we naturally choose the information from the previous (outer) step. We then get the two following variants of multi-step one-shot algorithms

$$\text{k-step one-shot:} \quad \begin{cases} \sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^n, \\ u_0^{n+1} = u^n, p_0^{n+1} = p^n, \\ \text{for } \ell = 0, 1, \dots, k-1: \\ \quad \left| \begin{array}{l} u_{\ell+1}^{n+1} = B u_{\ell}^{n+1} + M \sigma^{n+1} + F, \\ p_{\ell+1}^{n+1} = B^* p_{\ell}^{n+1} + H^*(H u_{\ell}^{n+1} - g), \\ u^{n+1} := u_k^{n+1}, p^{n+1} := p_k^{n+1} \end{array} \right. \end{cases} \quad (8)$$

and

$$\text{semi-implicit k-step one-shot:} \quad \begin{cases} \sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^{n+1}, \\ u_0^{n+1} = u^n, p_0^{n+1} = p^n, \\ \text{for } \ell = 0, 1, \dots, k-1: \\ \quad \left| \begin{array}{l} u_{\ell+1}^{n+1} = B u_{\ell}^{n+1} + M \sigma^{n+1} + F, \\ p_{\ell+1}^{n+1} = B^* p_{\ell}^{n+1} + H^*(H u_{\ell}^{n+1} - g), \\ u^{n+1} := u_k^{n+1}, p^{n+1} := p_k^{n+1} \end{array} \right. \end{cases} \quad (9)$$
In particular, when $k = 1$, we obtain the two following algorithms

$$ \text{one-step one-shot:} \quad \begin{cases} \sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^n, \\ u^{n+1} = B u^n + M \sigma^{n+1} + F, \\ p^{n+1} = B^* p^n + H^* (H u^n - g) \end{cases} \qquad (10) $$

and

$$ \text{semi-implicit one-step one-shot:} \quad \begin{cases} \sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^{n+1}, \\ u^{n+1} = B u^n + M \sigma^{n+1} + F, \\ p^{n+1} = B^* p^n + H^* (H u^n - g). \end{cases} \qquad (11) $$

Note that when $k \to \infty$, the $k$-step one-shot method (8) formally converges to the usual gradient descent (6), while the semi-implicit $k$-step one-shot method (9) formally converges to the semi-implicit gradient descent (7). Since the analysis of the two schemes (8) and (9) can be done following similar arguments, we choose to concentrate on only one of them, namely the semi-implicit scheme. We refer to [2] for the analysis in the case $\alpha = 0$ (for which the two schemes coincide).

We first analyze the one-step one-shot method ($k = 1$) in Section 3 and then the multi-step one-shot method ($k \ge 1$) in Section 4.

# 3 Convergence of the one-step one-shot method ($k = 1$)

## 3.1 Block iteration matrix and eigenvalue equation

To analyze the convergence of the semi-implicit one-step one-shot method (11), we first express $(\sigma^{n+1}, u^{n+1}, p^{n+1})$ in terms of $(\sigma^n, u^n, p^n)$, by inserting the expression for $\sigma^{n+1}$ into the iteration for $u^{n+1}$ in (11), so that system (11) is rewritten as

$$ \left\{ \begin{aligned} \sigma^{n+1} &= \frac{1}{1+\tau\alpha}\sigma^n - \frac{\tau}{1+\tau\alpha}M^*p^n, \\ u^{n+1} &= B u^n + \frac{1}{1+\tau\alpha}M\sigma^n - \frac{\tau}{1+\tau\alpha}MM^*p^n + F, \\ p^{n+1} &= B^*p^n + H^*Hu^n - H^*g. \end{aligned} \right. \qquad (12) $$

Now, we consider the errors $(\sigma^n - \sigma_\alpha^\text{ex}, u^n - u(\sigma_\alpha^\text{ex}), p^n - p(\sigma_\alpha^\text{ex}))$ with respect to the regularized solution at the $n$-th iteration, and, by abuse of notation, we designate them by $(\sigma^n, u^n, p^n)$. We obtain that these errors satisfy

$$ \left\{ \begin{aligned} \sigma^{n+1} &= \frac{1}{1+\tau\alpha}\sigma^n - \frac{\tau}{1+\tau\alpha}M^*p^n, \\ u^{n+1} &= B u^n + \frac{1}{1+\tau\alpha}M\sigma^n - \frac{\tau}{1+\tau\alpha}MM^*p^n, \\ p^{n+1} &= B^*p^n + H^*Hu^n, \end{aligned} \right. \qquad (13) $$

or equivalently, by putting in evidence the block iteration matrix

$$ \begin{bmatrix} p^{n+1} \\ u^{n+1} \\ \sigma^{n+1} \end{bmatrix} = \begin{bmatrix} B^* & H^*H & 0 \\ -\frac{\tau}{1+\tau\alpha}MM^* & B & \frac{1}{1+\tau\alpha}M \\ -\frac{\tau}{1+\tau\alpha}M^* & 0 & \frac{1}{1+\tau\alpha}I \end{bmatrix} \begin{bmatrix} p^n \\ u^n \\ \sigma^n \end{bmatrix}. \qquad (14) $$

Now recall that a fixed point iteration converges if and only if the spectral radius of its iteration matrix is less than 1. Therefore in the following proposition we establish an eigenvalue equation for the iteration matrix of the semi-implicit one-step one-shot method.

**Proposition 3.1.** Assume that $\lambda \in \mathbb{C}$, $|\lambda| \ge 1$ is an eigenvalue of the iteration matrix in (14). If $\lambda \in \mathbb{C}$, $\lambda \notin \text{Spec}(B)$ then $\exists y \in \mathbb{C}^{n_\sigma}$, $\|y\| = 1$ such that:

$$ (1 + \tau\alpha)\lambda - 1 + \tau\lambda \langle M^*(\lambda I - B^*)^{-1}H^*H(\lambda I - B)^{-1}My, y \rangle = 0. \qquad (15) $$

In particular, $\lambda = 1$ is not an eigenvalue of the iteration matrix.
In particular, when $k = 1$, we obtain the two following algorithms

$$ \text{one-step one-shot:} \quad \begin{cases} \sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^n, \\ u^{n+1} = B u^n + M \sigma^{n+1} + F, \\ p^{n+1} = B^* p^n + H^* (H u^n - g) \end{cases} \qquad (10) $$

and

$$ \text{semi-implicit one-step one-shot:} \quad \begin{cases} \sigma^{n+1} = \sigma^n - \tau M^* p^n - \tau \alpha \sigma^{n+1}, \\ u^{n+1} = B u^n + M \sigma^{n+1} + F, \\ p^{n+1} = B^* p^n + H^* (H u^n - g). \end{cases} \qquad (11) $$

Note that when $k \to \infty$, the k-step one-shot method (8) formally converges to the usual gradient descent (6), while the semi-implicit k-step one-shot method (9) formally converges to the semi-implicit gradient descent (7). Since the analysis of the two schemes (8) and (9) can be done following similar arguments, we choose to concentrate on only one of them, namely the semi-implicit scheme. We refer to [2] for the analysis in the case $\alpha = 0$ (for which the two schemes coincide).

We first analyze the one-step one-shot method ($k = 1$) in Section 3 and then the multi-step one-shot method ($k \ge 1$) in Section 4.

# 3 Convergence of the one-step one-shot method ($k = 1$)

## 3.1 Block iteration matrix and eigenvalue equation

To analyze the convergence of the semi-implicit one-step one-shot method (11), we first express $(\sigma^{n+1}, u^{n+1}, p^{n+1})$ in terms of $(\sigma^n, u^n, p^n)$, by inserting the expression for $\sigma^{n+1}$ into the iteration for $u^{n+1}$ in (11), so that system (11) is rewritten as

$$ \left\{ \begin{aligned} \sigma^{n+1} &= \frac{1}{1+\tau\alpha}\sigma^n - \frac{\tau}{1+\tau\alpha}M^*p^n, \\ u^{n+1} &= Bu^n + \frac{1}{1+\tau\alpha}M\sigma^n - \frac{\tau}{1+\tau\alpha}MM^*p^n + F, \\ p^{n+1} &= B^*p^n + H^*Hu^n - H^*g. \end{aligned} \right. \qquad (12) $$

Now, we consider the errors $(\sigma^n - \sigma_\alpha^{\text{ex}}, u^n - u(\sigma_\alpha^{\text{ex}}), p^n - p(\sigma_\alpha^{\text{ex}}))$ with respect to the regularized solution at the n-th iteration, and, by abuse of notation, we designate them by $(\sigma^n, u^n, p^n)$. We obtain that these errors satisfy

$$ \left\{ \begin{aligned} \sigma^{n+1} &= \frac{1}{1+\tau\alpha}\sigma^n - \frac{\tau}{1+\tau\alpha}M^*p^n, \\ u^{n+1} &= Bu^n + \frac{1}{1+\tau\alpha}M\sigma^n - \frac{\tau}{1+\tau\alpha}MM^*p^n, \\ p^{n+1} &= B^*p^n + H^*Hu^n, \end{aligned} \right. \qquad (13) $$

or equivalently, by putting in evidence the block iteration matrix

$$ \begin{bmatrix} p^{n+1} \\ u^{n+1} \\ \sigma^{n+1} \end{bmatrix} = \begin{bmatrix} B^* & H^*H & 0 \\ -\frac{\tau}{1+\tau\alpha}MM^* & B & \frac{1}{1+\tau\alpha}M \\ -\frac{\tau}{1+\tau\alpha}M^* & 0 & \frac{1}{1+\tau\alpha}I \end{bmatrix} \begin{bmatrix} p^n \\ u^n \\ \sigma^n \end{bmatrix}. \qquad (14) $$

Now recall that a fixed point iteration converges if and only if the spectral radius of its iteration matrix is less than 1. Therefore in the following proposition we establish an eigenvalue equation for the iteration matrix of the semi-implicit one-step one-shot method.

**Proposition 3.1.** Assume that $\lambda \in \mathbb{C}$, $|\lambda| \ge 1$ is an eigenvalue of the iteration matrix in (14). If $\lambda \in \mathbb{C}$, $\lambda \notin \text{Spec}(B)$ then $\exists y \in \mathbb{C}^{n_\sigma}$, $\|y\| = 1$ such that:

$$ (1 + \tau\alpha)\lambda - 1 + \tau\lambda \langle M^*(\lambda I - B^*)^{-1}H^*H(\lambda I - B)^{-1}My, y \rangle = 0. \qquad (15) $$

In particular, $\lambda = 1$ is not an eigenvalue of the iteration matrix.
Proof. Since $\lambda \in \mathbb{C}$ is an eigenvalue of the iteration matrix, there exists a non-zero vector $(\tilde{p}, \tilde{u}, y) \in \mathbb{C}^{n_u+n_u+n_\sigma}$ such that

$$ 
\begin{cases} 
\lambda y = \frac{1}{1+\tau\alpha}y - \frac{\tau}{1+\tau\alpha}M^*\tilde{p}, \\
\lambda \tilde{u} = B\tilde{u} + \frac{1}{1+\tau\alpha}My - \frac{\tau}{1+\tau\alpha}MM^*\tilde{p}, \\
\lambda \tilde{p} = B^*\tilde{p} + H^*H\tilde{u}. 
\end{cases} 
$$

By inserting the first equation into the second equation, we simplify this system of equations as

$$ 
\begin{cases} 
\lambda y = \frac{1}{1+\tau\alpha}y - \frac{\tau}{1+\tau\alpha}M^*\tilde{p}, \\
\lambda \tilde{u} = B\tilde{u} + \lambda My, \\
\lambda \tilde{p} = B^*\tilde{p} + H^*H\tilde{u}. 
\end{cases} 
\quad (16) 
$$

The second equation in (16) gives us directly $\tilde{u}$ in terms of $y$:

$$ \tilde{u} = \lambda(\lambda I - B)^{-1}My. \quad (17) $$

The third equation in (16) at first gives us $\tilde{p}$ in terms of $\tilde{u}$:

$$ \lambda \tilde{p} = (\lambda I - B^*)^{-1}H^*H\tilde{u}, $$

then by combining with equation (17), we obtain $\tilde{p}$ in terms of $y$:

$$ \tilde{p} = \lambda(\lambda I - B)^{-1}H^*H(\lambda I - B)^{-1}My. \quad (18) $$

We also see that $y \neq 0$; indeed if $y = 0$ then equations (17) and (18) yield $u = 0$ and $p = 0$, that is a contradiction. Next, inserting the expression of $\tilde{p}$ in equation (18) back into the first equation in (16), we get

$$ \lambda y = \frac{1}{1 + \tau \alpha} y - \frac{\tau}{1 + \tau \alpha} \lambda M^* (\lambda I - B^*)^{-1} H^* H (\lambda I - B)^{-1} My, $$

which leads to

$$ [(1 + \alpha\tau)\lambda - 1] y + \tau\lambda M^* (\lambda I - B^*)^{-1} H^* H (\lambda I - B)^{-1} My = 0. \quad (19) $$

Finally, by taking scalar product of (19) with $y$, then dividing by $\|y\|^2$, we obtain (15). Now assume that $\lambda = 1$ is an eigenvalue of the iteration matrix, then (15) yields

$$ \alpha + \|H(I - B)^{-1}My\|^2 = 0, $$

which cannot be true due to the injectivity of $H(I - B)^{-1}M$.

□

In the following sections we will show that, for sufficiently small $\tau$, equation (15) cannot hold if $|\lambda| \ge 1$, thus algorithm (11) converges. It is convenient to rewrite (15) as

$$ (1 + \tau\alpha)\lambda^2 - \lambda + \tau\langle M^* (I - B^*/\lambda)^{-1} H^* H (I - B/\lambda)^{-1} My, y \rangle = 0. \quad (20) $$

For the analysis we use auxiliary technical results proved in Appendix A.

## 3.2 Location of the eigenvalues in the complex plane

We now turn our attention to the eigenvalues $\lambda$ for which (20) holds. We would like to derive conditions on the descent step $\tau$ such that all the eigenvalues lie inside the unit circle which would ensure the convergence for the scheme (11). We start with the simple case of real eigenvalues.

**Proposition 3.2 (Real eigenvalues).** Equation (20) admits no solution $\lambda \in \mathbb{R}, \lambda \neq 1, |\lambda| \ge 1$ for all $\tau > 0$.
Proof. Since $\lambda \in \mathbb{C}$ is an eigenvalue of the iteration matrix, there exists a non-zero vector $(\tilde{p}, \tilde{u}, y) \in \mathbb{C}^{n_u+n_u+n_\sigma}$ such that

$$ 
\begin{cases}
\lambda y = \frac{1}{1+\tau\alpha}y - \frac{\tau}{1+\tau\alpha}M^*\tilde{p}, \\
\lambda \tilde{u} = B\tilde{u} + \frac{1}{1+\tau\alpha}My - \frac{\tau}{1+\tau\alpha}MM^*\tilde{p}, \\
\lambda \tilde{p} = B^*\tilde{p} + H^*H\tilde{u}.
\end{cases}
$$

By inserting the first equation into the second equation, we simplify this system of equations as

$$ 
\begin{cases}
\lambda y = \frac{1}{1+\tau\alpha}y - \frac{\tau}{1+\tau\alpha}M^*\tilde{p}, \\
\lambda \tilde{u} = B\tilde{u} + \lambda My, \\
\lambda \tilde{p} = B^*\tilde{p} + H^*H\tilde{u}.
\end{cases}
\quad (16) $$

The second equation in (16) gives us directly $\tilde{u}$ in terms of $y$:

$$ \tilde{u} = \lambda(\lambda I - B)^{-1}My. \quad (17) $$

The third equation in (16) at first gives us $\tilde{p}$ in terms of $\tilde{u}$:

$$ \lambda \tilde{p} = (\lambda I - B^*)^{-1}H^*H\tilde{u}, $$

then by combining with equation (17), we obtain $\tilde{p}$ in terms of $y$:

$$ \tilde{p} = \lambda(\lambda I - B^*)^{-1}H^*H(\lambda I - B)^{-1}My. \quad (18) $$

We also see that $y \neq 0$; indeed if $y = 0$ then equations (17) and (18) yield $u = 0$ and $p = 0$, that is a contradiction. Next, inserting the expression of $\tilde{p}$ in equation (18) back into the first equation in (16), we get

$$ \lambda y = \frac{1}{1 + \tau \alpha} y - \frac{\tau}{1 + \tau \alpha} \lambda M^* (\lambda I - B^*)^{-1} H^* H (\lambda I - B)^{-1} My, $$

which leads to

$$ [(1 + \alpha\tau)\lambda - 1] y + \tau\lambda M^* (\lambda I - B^*)^{-1} H^* H (\lambda I - B)^{-1} My = 0. \quad (19) $$

Finally, by taking scalar product of (19) with $y$, then dividing by $\|y\|^2$, we obtain (15). Now assume that $\lambda = 1$ is an eigenvalue of the iteration matrix, then (15) yields

$$ \alpha + \|H(I - B)^{-1}My\|^2 = 0, $$

which cannot be true due to the injectivity of $H(I - B)^{-1}M$.

□

In the following sections we will show that, for sufficiently small $\tau$, equation (15) cannot hold if $|\lambda| \ge 1$, thus algorithm (11) converges. It is convenient to rewrite (15) as

$$ (1 + \tau\alpha)\lambda^2 - \lambda + \tau \langle M^* (I - B^*/\lambda)^{-1} H^* H (I - B/\lambda)^{-1} My, y \rangle = 0. \quad (20) $$

For the analysis we use auxiliary technical results proved in Appendix A.

## 3.2 Location of the eigenvalues in the complex plane

We now turn our attention to the eigenvalues $\lambda$ for which (20) holds. We would like to derive conditions on the descent step $\tau$ such that all the eigenvalues lie inside the unit circle which would ensure the convergence for the scheme (11). We start with the simple case of real eigenvalues.

**Proposition 3.2 (Real eigenvalues).** Equation (20) admits no solution $\lambda \in \mathbb{R}, \lambda \neq 1, |\lambda| \ge 1$ for all $\tau > 0$.
Proof. When $\lambda \in \mathbb{R}\setminus\{0\}$ equation (20) becomes

$$(1 + \tau\alpha)\lambda^2 - \lambda + \tau\|H(I - B/\lambda)^{-1}My\|^2 = 0.$$

If $\lambda \in \mathbb{R}$, $\lambda \neq 1$, $|\lambda| \ge 1$ then $(1 + \tau\alpha)\lambda^2 - \lambda \ge \lambda^2 - \lambda > 0$, thus the left-hand side of the above equation is positive for any $\tau > 0$. $\square$

For the general case of complex eigenvalues, the study is much more complicated and technical. First, we study separately the very particular and simple case where $B = 0$.

**Proposition 3.3.** When $B = 0$, the eigenvalue equation (20) cannot hold for $\lambda \in \mathbb{C}$, $|\lambda| \ge 1$ if $\tau > 0$ and

$$(\|H\|^2\|M\|^2 - \alpha)\tau < 1.$$

*Proof.* When $B = 0$, $HM$ is injective by (4), and equation (20) becomes $(1+\tau\alpha)\lambda^2 - \lambda + \tau\|HMy\|^2 = 0$. One can prove the proposition by explicitly computing the roots of the second order polynomial with respect to $\lambda$. One can also apply the following elementary Lemma 3.4 below, which can be deduced from Marden's work [18] or Appendix C of [2]. Indeed, the previous polynomial can be written as $P(\lambda) = a_0 + a_1\lambda + \lambda^2$ with $a_0 = \frac{\tau\|HMy\|^2}{1+\tau\alpha}$ and $a_1 = \frac{-1}{1+\tau\alpha}$. We have $a_0 - a_1 + 1 = \frac{\tau(\|HMy\|^2+\alpha)+2}{1+\tau\alpha} > 0$ and $a_0 + a_1 + 1 = \frac{\tau(\|HMy\|^2+\alpha)}{1+\tau\alpha} > 0$ since $\|HMy\|^2 > 0$ thanks to the injectivity of $HM$. Therefore, according to Lemma 3.4, the roots of $P$ stay strictly inside the unit circle of the complex plane if and only if $|a_0| = \frac{\tau\|HMy\|^2}{1+\tau\alpha} < 1$, which is equivalent to $(\|HMy\|^2 - \alpha)\tau < 1$. Since $\|y\| = 1$, the later inequality is verified if $(\|H\|^2\|M\|^2 - \alpha)\tau < 1$, which is the statement of the proposition. $\square$

**Lemma 3.4.** Let $a_0, a_1 \in \mathbb{R}$, then all roots of $\mathcal{P}(z) = a_0 + a_1z + z^2$ stay strictly inside the unit circle of the complex plane if and only if

$$|a_0| < 1 \quad \text{and} \quad (a_0 - a_1 + 1)(a_0 + a_1 + 1) > 0.$$

Now we consider the more complicated case where $0 \neq \rho(B) < 1$. The following proposition summarizes the results we obtained.

**Proposition 3.5 (Complex eigenvalues).** If $0 \neq \rho(B) < 1$, there exists $\tau > 0$ sufficiently small such that equation (20) admits no solution $\lambda \in \mathbb{C}\setminus\mathbb{R}$, $|\lambda| \ge 1$. In particular, if $\|B\| < 1$, given any $\delta_0 > 0$ and $0 < \theta_0 \le \frac{\pi}{4}$, one can choose

$$\tau < \min_{1 \le i \le 3} \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|^2)^4} \varphi_i(\|B\|) + C_i \alpha \right)^{-1},$$

where the functions $\varphi_i$, $1 \le i \le 3$, are defined by

$$\varphi_1(b) := 4b^2, \quad \varphi_2(b) := \frac{1}{2 \sin \frac{\theta_0}{2}} (1+b)^2 (1-b)^2 \quad \text{and} \quad \varphi_3(b) := \frac{2c}{\delta_0} b^2,$$

and the pure constants $C_i$, $1 \le i \le 3$, are defined by

$$C_1 := \sqrt{2} - 1, \quad C_2 := \sqrt{2} + \frac{1}{2 \sin \frac{\theta_0}{2}} - 1 \quad \text{and} \quad C_3 := \frac{\sqrt{c}}{\delta_0} - 1$$

with

$$c := \frac{1 + 2\delta_0 \sin \frac{3\theta_0}{2} + \delta_0^2}{\cos^2 \frac{3\theta_0}{2}}.$$
Proof. When $\lambda \in \mathbb{R}\setminus\{0\}$ equation (20) becomes

$$(1 + \tau\alpha)\lambda^2 - \lambda + \tau\|H(I - B/\lambda)^{-1}My\|^2 = 0.$$

If $\lambda \in \mathbb{R}$, $\lambda \neq 1$, $|\lambda| \ge 1$ then $(1 + \tau\alpha)\lambda^2 - \lambda \ge \lambda^2 - \lambda > 0$, thus the left-hand side of the above equation is positive for any $\tau > 0$. $\square$

For the general case of complex eigenvalues, the study is much more complicated and technical. First, we study separately the very particular and simple case where $B = 0$.

**Proposition 3.3.** When $B = 0$, the eigenvalue equation (20) cannot hold for $\lambda \in \mathbb{C}$, $|\lambda| \ge 1$ if $\tau > 0$ and

$$(\|H\|^2\|M\|^2 - \alpha)\tau < 1.$$

*Proof.* When $B = 0$, $HM$ is injective by (4), and equation (20) becomes $(1+\tau\alpha)\lambda^2 - \lambda + \tau\|HMy\|^2 = 0$. One can prove the proposition by explicitly computing the roots of the second order polynomial with respect to $\lambda$. One can also apply the following elementary Lemma 3.4 below, which can be deduced from Marden's work [18] or Appendix C of [2]. Indeed, the previous polynomial can be written as $P(\lambda) = a_0 + a_1\lambda + \lambda^2$ with $a_0 = \frac{\tau\|HMy\|^2}{1+\tau\alpha}$ and $a_1 = \frac{-1}{1+\tau\alpha}$. We have $a_0 - a_1 + 1 = \frac{\tau(\|HMy\|^2+\alpha)+2}{1+\tau\alpha} > 0$ and $a_0 + a_1 + 1 = \frac{\tau(\|HMy\|^2+\alpha)}{1+\tau\alpha} > 0$ since $\|HMy\|^2 > 0$ thanks to the injectivity of $HM$. Therefore, according to Lemma 3.4, the roots of $P$ stay strictly inside the unit circle of the complex plane if and only if $|a_0| = \frac{\tau\|HMy\|^2}{1+\tau\alpha} < 1$, which is equivalent to $(\|HMy\|^2 - \alpha)\tau < 1$. Since $\|y\| = 1$, the later inequality is verified if $(\|H\|^2\|M\|^2 - \alpha)\tau < 1$, which is the statement of the proposition. $\square$

**Lemma 3.4.** Let $a_0, a_1 \in \mathbb{R}$, then all roots of $\mathcal{P}(z) = a_0 + a_1z + z^2$ stay strictly inside the unit circle of the complex plane if and only if

$$|a_0| < 1 \quad \text{and} \quad (a_0 - a_1 + 1)(a_0 + a_1 + 1) > 0.$$

Now we consider the more complicated case where $0 \neq \rho(B) < 1$. The following proposition summarizes the results we obtained.

**Proposition 3.5 (Complex eigenvalues).** If $0 \neq \rho(B) < 1$, there exists $\tau > 0$ sufficiently small such that equation (20) admits no solution $\lambda \in \mathbb{C}\setminus\mathbb{R}$, $|\lambda| \ge 1$. In particular, if $\|B\| < 1$, given any $\delta_0 > 0$ and $0 < \theta_0 \le \frac{\pi}{4}$, one can choose

$$\tau < \min_{1 \le i \le 3} \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|^2)^4} \varphi_i(\|B\|) + C_i \alpha \right)^{-1},$$

where the functions $\varphi_i$, $1 \le i \le 3$, are defined by

$$\varphi_1(b) := 4b^2, \quad \varphi_2(b) := \frac{1}{2 \sin \frac{\theta_0}{2}} (1+b)^2 (1-b)^2 \quad \text{and} \quad \varphi_3(b) := \frac{2c}{\delta_0} b^2,$$

and the pure constants $C_i$, $1 \le i \le 3$, are defined by

$$C_1 := \sqrt{2} - 1, \quad C_2 := \sqrt{2} + \frac{1}{2 \sin \frac{\theta_0}{2}} - 1 \quad \text{and} \quad C_3 := \frac{\sqrt{c}}{\delta_0} - 1$$

with

$$c := \frac{1 + 2\delta_0 \sin \frac{3\theta_0}{2} + \delta_0^2}{\cos^2 \frac{3\theta_0}{2}}.$$
### Proof. Step 1. Rewrite equation (20) by separating real and imaginary parts.

Let $\lambda = R(\cos \theta + i \sin \theta)$ in polar form where $R = |\lambda| \ge 1$ and $\theta \in (-\pi, \pi)$, $\theta \ne 0$. Write $1/\lambda = r(\cos \phi + i \sin \phi)$ in polar form where $r = 1/|\lambda| = 1/R \le 1$ and $\phi = -\theta \in (-\pi, \pi)$. By Lemma A.3, we have

$$ \left(I - \frac{B}{\lambda}\right)^{-1} = P(\lambda) + iQ(\lambda), \quad \left(I - \frac{B^*}{\lambda}\right)^{-1} = P(\lambda)^* + iQ(\lambda)^* $$

where $P(\lambda)$ and $Q(\lambda)$ are $\mathbb{C}^{n_u \times n_u}$ matrices that satisfy the following bounds in the case $\|B\| < 1$ for all $|\lambda| \ge 1$:

$$ \|P(\lambda)\| \le p := \frac{1}{1 - \|B\|}, \qquad (21) $$

$$ \|Q(\lambda)\| \le q_1 := \frac{\|B\|}{1 - \|B\|} \quad \text{and} \quad \|Q(\lambda)\| \le |\sin \theta|q_2 \quad \text{with} \quad q_2 := \frac{\|B\|}{(1 - \|B\|)^2}. \qquad (22) $$

These bounds still hold in the case $0 \ne \rho(B) < 1$ with

$$ p = (1 + \|B\|)s(B)^2, \quad q_1 = \|B\|s(B)^2 \quad \text{and} \quad q_2 = \|B\|s(B)^2. \qquad (23) $$

To simplify the notation, we will not explicitly write the dependence of $P$ and $Q$ on $\lambda$. Now we rewrite (20) as

$$ (1 + \tau\alpha)(\lambda^2 - \lambda) + \tau\alpha(\lambda - 1) + \tau\alpha + \tau G(P^* + iQ^*, P + iQ) = 0 \qquad (24) $$

where

$$ G(X, Y) := \langle M^* X H^* H Y M y, y \rangle \in \mathbb{C}, \quad X, Y \in \mathbb{C}^{n_u \times n_u}. $$

Notice that $G$ is a bilinear form and $G(X, Y) = G(Y^*, X^*)^*$ so that $G(X, Y) + G(X^*, Y^*)$ is real. With these properties of $G$, we expand (24) and take its real and imaginary parts, which yields

$$ (1 + \tau\alpha) \operatorname{Re}(\lambda^2 - \lambda) + \tau\alpha \operatorname{Re}(\lambda - 1) + \tau\alpha + \tau[G(P^*, P) - G(Q^*, Q)] = 0, \qquad (25) $$

and

$$ (1 + \tau\alpha) \operatorname{Im}(\lambda^2 - \lambda) + \tau\alpha \operatorname{Im}(\lambda - 1) + \tau[G(P^*, Q) + G(Q^*, P)] = 0. \qquad (26) $$

### Step 2. Use a suitable combination of equations (25) and (26).

Let $\gamma \in \mathbb{R}$. Multiplying equation (26) with $\gamma$ then summing it with equation (25), we obtain:

$$ (1 + \tau\alpha)[\operatorname{Re}(\lambda^2 - \lambda) + \gamma \operatorname{Im}(\lambda^2 - \lambda)] + \tau\alpha[\operatorname{Re}(\lambda - 1) + \gamma \operatorname{Im}(\lambda - 1)] \\ + \tau\alpha + \tau[G(P^*, P) - G(Q^*, Q) + \gamma G(P^*, Q) + \gamma G(Q^*, P)] = 0, $$

or equivalently,

$$ (1 + \tau\alpha)[\operatorname{Re}(\lambda^2 - \lambda) + \gamma \operatorname{Im}(\lambda^2 - \lambda)] + \tau\alpha[\operatorname{Re}(\lambda - 1) + \gamma \operatorname{Im}(\lambda - 1)] \\ + \tau\alpha + \tau G(P^* + \gamma Q^*, P + \gamma Q) - (1 + \gamma^2)\tau G(Q^*, Q) = 0. \qquad (27) $$

Now we consider four cases for $\lambda$ as in Lemma A.4 (see Figure 1):

* Case 1. $\operatorname{Re}(\lambda^2 - \lambda) \ge 0$;

* Case 2. $\operatorname{Re}(\lambda^2 - \lambda) < 0$ and $\theta \in [\theta_0, \pi - \theta_0] \cup [-\pi + \theta_0, -\theta_0]$ for fixed $0 < \theta_0 \le \frac{\pi}{4}$;

* Case 3. $\operatorname{Re}(\lambda^2 - \lambda) < 0$ and $\theta \in (-\theta_0, \theta_0)$ for fixed $0 < \theta_0 \le \frac{\pi}{4}$;

* Case 4. $\operatorname{Re}(\lambda^2 - \lambda) < 0$ and $\theta \in (\pi - \theta_0, \pi) \cup (-\pi, -\pi + \theta_0)$ for fixed $0 < \theta_0 \le \frac{\pi}{4}$.
### Proof. Step 1. Rewrite equation (20) by separating real and imaginary parts.

Let $\lambda = R(\cos \theta + i \sin \theta)$ in polar form where $R = |\lambda| \ge 1$ and $\theta \in (-\pi, \pi)$, $\theta \ne 0$. Write $1/\lambda = r(\cos \phi + i \sin \phi)$ in polar form where $r = 1/|\lambda| = 1/R \le 1$ and $\phi = -\theta \in (-\pi, \pi)$. By Lemma A.3, we have

$$ \left(I - \frac{B}{\lambda}\right)^{-1} = P(\lambda) + iQ(\lambda), \quad \left(I - \frac{B^*}{\lambda}\right)^{-1} = P(\lambda)^* + iQ(\lambda)^* $$

where $P(\lambda)$ and $Q(\lambda)$ are $\mathbb{C}^{n_u \times n_u}$ matrices that satisfy the following bounds in the case $\|B\| < 1$ for all $|\lambda| \ge 1$:

$$ \|P(\lambda)\| \le p := \frac{1}{1 - \|B\|}, \qquad (21) $$

$$ \|Q(\lambda)\| \le q_1 := \frac{\|B\|}{1 - \|B\|} \quad \text{and} \quad \|Q(\lambda)\| \le |\sin \theta|q_2 \quad \text{with} \quad q_2 := \frac{\|B\|}{(1 - \|B\|)^2}. \qquad (22) $$

These bounds still hold in the case $0 \ne \rho(B) < 1$ with

$$ p = (1 + \|B\|)s(B)^2, \quad q_1 = \|B\|s(B)^2 \quad \text{and} \quad q_2 = \|B\|s(B)^2. \qquad (23) $$

To simplify the notation, we will not explicitly write the dependence of $P$ and $Q$ on $\lambda$. Now we rewrite (20) as

$$ (1 + \tau\alpha)(\lambda^2 - \lambda) + \tau\alpha(\lambda - 1) + \tau\alpha + \tau G(P^* + iQ^*, P + iQ) = 0 \qquad (24) $$

where

$$ G(X, Y) := \langle M^*XH^*HYM_y, y \rangle \in \mathbb{C}, \quad X, Y \in \mathbb{C}^{n_u \times n_u}. $$

Notice that $G$ is a bilinear form and $G(X, Y) = G(Y^*, X^*)^*$ so that $G(X, Y) + G(X^*, Y^*)$ is real. With these properties of $G$, we expand (24) and take its real and imaginary parts, which yields

$$ (1 + \tau\alpha) \operatorname{Re}(\lambda^2 - \lambda) + \tau\alpha \operatorname{Re}(\lambda - 1) + \tau\alpha + \tau[G(P^*, P) - G(Q^*, Q)] = 0, \qquad (25) $$

and

$$ (1 + \tau\alpha) \operatorname{Im}(\lambda^2 - \lambda) + \tau\alpha \operatorname{Im}(\lambda - 1) + \tau[G(P^*, Q) + G(Q^*, P)] = 0. \qquad (26) $$

### Step 2. Use a suitable combination of equations (25) and (26).

Let $\gamma \in \mathbb{R}$. Multiplying equation (26) with $\gamma$ then summing it with equation (25), we obtain:

$$ (1 + \tau\alpha)[\operatorname{Re}(\lambda^2 - \lambda) + \gamma \operatorname{Im}(\lambda^2 - \lambda)] + \tau\alpha[\operatorname{Re}(\lambda - 1) + \gamma \operatorname{Im}(\lambda - 1)] \\ + \tau\alpha + \tau[G(P^*, P) - G(Q^*, Q) + \gamma G(P^*, Q) + \gamma G(Q^*, P)] = 0, $$

or equivalently,

$$ (1 + \tau\alpha)[\operatorname{Re}(\lambda^2 - \lambda) + \gamma \operatorname{Im}(\lambda^2 - \lambda)] + \tau\alpha[\operatorname{Re}(\lambda - 1) + \gamma \operatorname{Im}(\lambda - 1)] \\ + \tau\alpha + \tau G(P^* + \gamma Q^*, P + \gamma Q) - (1 + \gamma^2)\tau G(Q^*, Q) = 0. \qquad (27) $$

Now we consider four cases for $\lambda$ as in Lemma A.4 (see Figure 1):

* Case 1. $\operatorname{Re}(\lambda^2 - \lambda) \ge 0$;

* Case 2. $\operatorname{Re}(\lambda^2 - \lambda) < 0$ and $\theta \in [\theta_0, \pi - \theta_0] \cup [-\pi + \theta_0, -\theta_0]$ for fixed $0 < \theta_0 \le \frac{\pi}{4}$;

* Case 3. $\operatorname{Re}(\lambda^2 - \lambda) < 0$ and $\theta \in (-\theta_0, \theta_0)$ for fixed $0 < \theta_0 \le \frac{\pi}{4}$;

* Case 4. $\operatorname{Re}(\lambda^2 - \lambda) < 0$ and $\theta \in (\pi - \theta_0, \pi) \cup (-\pi, -\pi + \theta_0)$ for fixed $0 < \theta_0 \le \frac{\pi}{4}$.
<figure><img src="image_3.png" /><figcaption>Figure 1: Illustration of the regions in the complex plane associated with the cases 1, 2 and 3 indicated in the proof of Proposition 3.5 for $\theta_0 = \frac{\pi}{4}$. The center circle is the unit circle.</figcaption></figure>

Cases 1, 2 and 3 are respectively treated in Lemmas 3.6, 3.7 and 3.8 below. Notice that case 4 corresponds to an empty set, according to Lemma A.4 (iv). The statement of the proposition easily follows from the combination of Lemmas 3.6, 3.7 and 3.8. In particular, the parameter $\theta_0$ is mainly used for case 2 to obtain a lower bound for $|\lambda - 1|$. This is why we require $\theta_0 > 0$. We need $\theta_0 \le \frac{\pi}{4}$ to design a lower bound for a suitable combination of $\mathrm{Re}(\lambda^2 - \lambda)$ and $\mathrm{Im}(\lambda^2 - \lambda)$ in case 3. $\square$

In the lemmas below we shall make use of the following obvious property where $\|y\| = 1$:

$$0 \le G(X^*, X) = \|HXMy\|^2 \le (\|H\|\|M\|\|X\|)^2, \quad \forall X \in \mathbb{C}^{n_u \times n_u}. \quad (28)$$

**Lemma 3.6 (Case 1).** Let $0 \ne \rho(B) < 1$ and let $|\lambda| \ge 1$ with $\mathrm{Re}(\lambda^2 - \lambda) \ge 0$. Then equation (20) cannot hold if one chooses

$$\tau < \left( \|H\|^2 \|M\|^2 s(B)^4 \cdot 4 \|B\|^2 + (\sqrt{2} - 1)\alpha \right)^{-1}.$$

Moreover, if $\|B\| < 1$, the result is also true if

$$\tau < \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|^2)^4} \cdot 4 \|B\|^2 + (\sqrt{2} - 1)\alpha \right)^{-1}.$$

*Proof.* Define

$$\gamma_1 = \gamma_1(\lambda) := \begin{cases} 1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) \ge 0, \\ -1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) < 0 \end{cases}$$

as in Lemma A.4. Writing (27) for $\gamma = \gamma_1$ and using $\gamma_1^2 = 1$ we obtain

$$(1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha[\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)] \\ + \tau\alpha + \tau G(P^* + \gamma_1 Q^*, P + \gamma_1 Q) - 2\tau G(Q^*, Q) = 0. \quad (29)$$

Since $G(P^* + \gamma_1 Q^*, P + \gamma_1 Q) \ge 0$ by (28) and $\tau\alpha \ge 0$, then the left-hand side of (29) is positive if $\tau$ satisfies

$$(1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)] - \tau\alpha|\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)| - 2\tau G(Q^*, Q) > 0,$$

or equivalently,

$$1 + \tau\alpha - \tau\alpha \frac{|\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} - 2\tau \frac{G(Q^*, Q)}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} > 0. \quad (30)$$
<figure><img src="image_3.png" /><figcaption>Figure 1: Illustration of the regions in the complex plane associated with the cases 1, 2 and 3 indicated in the proof of Proposition 3.5 for $\theta_0 = \frac{\pi}{4}$. The center circle is the unit circle.</figcaption></figure>

Cases 1, 2 and 3 are respectively treated in Lemmas 3.6, 3.7 and 3.8 below. Notice that case 4 corresponds to an empty set, according to Lemma A.4 (iv). The statement of the proposition easily follows from the combination of Lemmas 3.6, 3.7 and 3.8. In particular, the parameter $\theta_0$ is mainly used for case 2 to obtain a lower bound for $|\lambda - 1|$. This is why we require $\theta_0 > 0$. We need $\theta_0 \le \frac{\pi}{4}$ to design a lower bound for a suitable combination of $\mathrm{Re}(\lambda^2 - \lambda)$ and $\mathrm{Im}(\lambda^2 - \lambda)$ in case 3. $\square$

In the lemmas below we shall make use of the following obvious property where $\|y\| = 1$:

$$0 \le G(X^*, X) = \|HXMy\|^2 \le (\|H\|\|M\|\|X\|)^2, \quad \forall X \in \mathbb{C}^{n_u \times n_u}. \quad (28)$$

**Lemma 3.6 (Case 1).** Let $0 \ne \rho(B) < 1$ and let $|\lambda| \ge 1$ with $\mathrm{Re}(\lambda^2 - \lambda) \ge 0$. Then equation (20) cannot hold if one chooses

$$\tau < (\|H\|^2 \|M\|^2 s(B)^4 \cdot 4 \|B\|^2 + (\sqrt{2} - 1)\alpha)^{-1}.$$

Moreover, if $\|B\| < 1$, the result is also true if

$$\tau < \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|^2)^4} \cdot 4 \|B\|^2 + (\sqrt{2} - 1)\alpha \right)^{-1}.$$

*Proof.* Define

$$\gamma_1 = \gamma_1(\lambda) := \begin{cases} 1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) \ge 0, \\ -1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) < 0 \end{cases}$$

as in Lemma A.4. Writing (27) for $\gamma = \gamma_1$ and using $\gamma_1^2 = 1$ we obtain

$$(1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha[\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)] \\ + \tau\alpha + \tau G(P^* + \gamma_1 Q^*, P + \gamma_1 Q) - 2\tau G(Q^*, Q) = 0. \quad (29)$$

Since $G(P^* + \gamma_1 Q^*, P + \gamma_1 Q) \ge 0$ by (28) and $\tau\alpha \ge 0$, then the left-hand side of (29) is positive if $\tau$ satisfies

$$(1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)] - \tau\alpha |\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)| - 2\tau G(Q^*, Q) > 0,$$

or equivalently,

$$1 + \tau\alpha - \tau\alpha \frac{|\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} - 2\tau \frac{G(Q^*, Q)}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} > 0. \quad (30)$$
Notice that the choice of $\gamma_1$ ensures $\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda) > 0$. By Lemma A.4 we have

$$ \frac{|\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\sqrt{1 + \gamma_1^2 |\lambda - 1|}}{|\lambda(\lambda - 1)|} = \frac{\sqrt{2}}{|\lambda|} \le \sqrt{2}. $$

Using again Lemma A.4 and (28), we have

$$ \frac{G(Q^*, Q)}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\|H\|\|M\| \sin\theta |q_2|^2}{2|\sin(\theta/2)|} = 2 \left|\sin\frac{\theta}{2}\right| \cos^2\frac{\theta}{2} \|H\|^2 \|M\|^2 q_2^2 \\ \le 2 \|H\|^2 \|M\|^2 q_2^2. $$

Inserting the two previous inequalities in (30) gives the desired results using definitions (22) and (23) of $q_2$. $\square$

**Lemma 3.7 (Case 2).** Let $0 \ne \rho(B) < 1$ and let $|\lambda| \ge 1$, $\mathrm{Re}(\lambda^2 - \lambda) < 0$, $\theta \in [\theta_0, \pi - \theta_0] \cup [-\pi + \theta_0, -\theta_0]$ for given $0 < \theta_0 \le \frac{\pi}{4}$. Then equation (20) cannot hold if one chooses

$$ \tau < \left( \|H\|^2 \|M\|^2 s(B)^4 \cdot \frac{(1 + 2\|B\|)^2}{2 \sin \frac{\theta_0}{2}} + \left( \sqrt{2} + \frac{1}{2 \sin \frac{\theta_0}{2}} - 1 \right) \alpha \right)^{-1}. $$

Moreover, if $\|B\| < 1$, the result is also true if

$$ \tau < \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|)^2} \cdot \frac{(1 + \|B\|)^2}{2 \sin \frac{\theta_0}{2}} + \left( \sqrt{2} + \frac{1}{2 \sin \frac{\theta_0}{2}} - 1 \right) \alpha \right)^{-1}. $$

*Proof.* Define

$$ \gamma_2 = \gamma_2(\lambda) := \begin{cases} -1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) \ge 0, \\ 1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) < 0 \end{cases} $$

as in Lemma A.4 (ii). Writing (27) for $\gamma = \gamma_2$ and using $\gamma_2^2 = 1$ we obtain

$$ (1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha[\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)] \\ + \tau\alpha + \tau G(P^* + \gamma_2 Q^*, P + \gamma_2 Q) - 2\tau G(Q^*, Q) = 0. \quad (31) $$

Since $G(Q^*, Q) \ge 0$ by (28), the left-hand side of (31) is negative if $\tau$ satisfies

$$ (1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha|\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)| \\ + \tau\alpha + \tau G(P^* + \gamma_2 Q^*, P + \gamma_2 Q) < 0, $$

or equivalently,

$$ -1 - \tau\alpha + \tau\alpha \frac{|\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)|}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} + \frac{\tau\alpha}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \\ + \tau \frac{G(P^* + \gamma_2 Q^*, P + \gamma_2 Q)}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} < 0. \quad (32) $$

Notice that the choice of $\gamma_2$ ensures $\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda) < 0$. In the following we derive upper bounds independent of $\lambda$ for the terms appearing with the negative sign in (32). By Lemma A.4 (ii) we have

$$ \frac{|\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)|}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \le \frac{\sqrt{1 + \gamma_2^2 |\lambda - 1|}}{|\lambda(\lambda - 1)|} = \frac{\sqrt{2}}{|\lambda|} \le \sqrt{2}. $$
Notice that the choice of $\gamma_1$ ensures $\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda) > 0$. By Lemma A.4 we have

$$ \frac{|\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\sqrt{1 + \gamma_1^2 |\lambda - 1|}}{|\lambda(\lambda - 1)|} = \frac{\sqrt{2}}{|\lambda|} \le \sqrt{2}. $$

Using again Lemma A.4 and (28), we have

$$ \frac{G(Q^*, Q)}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\|(H)\|\|M\| \sin\theta |q_2|^2}{2|\sin(\theta/2)|} = 2 \left|\sin\frac{\theta}{2}\right| \cos^2\frac{\theta}{2} \|H\|^2 \|M\|^2 q_2^2 \\ \le 2 \|H\|^2 \|M\|^2 q_2^2. $$

Inserting the two previous inequalities in (30) gives the desired results using definitions (22) and (23) of $q_2$. □

**Lemma 3.7 (Case 2).** Let $0 \ne \rho(B) < 1$ and let $|\lambda| \ge 1$, $\mathrm{Re}(\lambda^2 - \lambda) < 0$, $\theta \in [\theta_0, \pi - \theta_0] \cup [-\pi + \theta_0, -\theta_0]$ for given $0 < \theta_0 \le \frac{\pi}{4}$. Then equation (20) cannot hold if one chooses

$$ \tau < \left( \|H\|^2 \|M\|^2 s(B)^4 \cdot \frac{(1 + 2\|B\|)^2}{2 \sin \frac{\theta_0}{2}} + \left( \sqrt{2} + \frac{1}{2 \sin \frac{\theta_0}{2}} - 1 \right) \alpha \right)^{-1}. $$

Moreover, if $\|B\| < 1$, the result is also true if

$$ \tau < \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|)^2} \cdot \frac{(1 + \|B\|)^2}{2 \sin \frac{\theta_0}{2}} + \left( \sqrt{2} + \frac{1}{2 \sin \frac{\theta_0}{2}} - 1 \right) \alpha \right)^{-1}. $$

*Proof.* Define

$$ \gamma_2 = \gamma_2(\lambda) := \begin{cases} -1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) \ge 0, \\ 1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) < 0 \end{cases} $$

as in Lemma A.4 (ii). Writing (27) for $\gamma = \gamma_2$ and using $\gamma_2^2 = 1$ we obtain

$$ (1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha[\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)] \\ + \tau\alpha + \tau G(P^* + \gamma_2 Q^*, P + \gamma_2 Q) - 2\tau G(Q^*, Q) = 0. \quad (31) $$

Since $G(Q^*, Q) \ge 0$ by (28), the left-hand side of (31) is negative if $\tau$ satisfies

$$ (1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha|\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)| \\ + \tau\alpha + \tau G(P^* + \gamma_2 Q^*, P + \gamma_2 Q) < 0, $$

or equivalently,

$$ -1 - \tau\alpha + \tau\alpha \frac{|\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)|}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} + \frac{\tau\alpha}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \\ + \tau \frac{G(P^* + \gamma_2 Q^*, P + \gamma_2 Q)}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} < 0. \quad (32) $$

Notice that the choice of $\gamma_2$ ensures $\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda) < 0$. In the following we derive upper bounds independent of $\lambda$ for the terms appearing with the negative sign in (32). By Lemma A.4 (ii) we have

$$ \frac{|\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)|}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \le \frac{\sqrt{1 + \gamma_2^2 |\lambda - 1|}}{|\lambda(\lambda - 1)|} = \frac{\sqrt{2}}{|\lambda|} \le \sqrt{2}. $$
and

$$ \frac{1}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \le \frac{1}{2 \sin \frac{\theta_0}{2}}. $$

Using again Lemma A.4 (ii) and (28), we have

$$ \frac{G(P^* + \gamma_2 Q^*, P + \gamma_2 Q)}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \le \frac{\|H\|^2 \|M\|^2 (p + q_1)^2}{2 \sin \frac{\theta_0}{2}}. $$

Inserting these previous inequalities in (32) gives the desired results using definitions (21), (22) and (23) of p and $q_1$. □

**Lemma 3.8 (Case 3).** Let $0 \ne \rho(B) < 1$ and let $|\lambda| \ge 1$, $\mathrm{Re}(\lambda^2 - \lambda) < 0$, $\theta \in (-\theta_0, \theta_0)$ for given $0 < \theta_0 \le \frac{\pi}{4}$. For any $\delta_0 > 0$, equation (20) cannot hold if one chooses

$$ \tau < \left( \|H\|^2 \|M\|^2 s(B)^4 \cdot \frac{2c}{\delta_0} \|B\|^2 + \left( \frac{\sqrt{c}}{\delta_0} - 1 \right) \alpha \right)^{-1} $$

where $c = c(\theta_0, \delta_0) := (1 + 2\delta_0 \sin \frac{3\theta_0}{2} + \delta_0^2) / \cos^2 \frac{3\theta_0}{2}$. Moreover, if $0 < \|B\| < 1$, the result is also true if

$$ \tau < \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|)^{-4}} \cdot \frac{2c}{\delta_0} \|B\|^2 + \left( \frac{\sqrt{c}}{\delta_0} - 1 \right) \alpha \right)^{-1}. $$

*Proof.* Define

$$ \gamma_3 = \gamma_3(\mathrm{sign}(\theta)) := \begin{cases} \left(\delta_0 + \sin \frac{3\theta_0}{2}\right) / \cos \frac{3\theta_0}{2} & \text{if } \theta > 0, \\ -\left(\delta_0 + \sin \frac{3\theta_0}{2}\right) / \cos \frac{3\theta_0}{2} & \text{if } \theta < 0 \end{cases} $$

as in Lemma A.4 (iii). Writing (27) for $\gamma = \gamma_3$ we obtain

$$ (1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha[\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)] \\ + \tau\alpha + \tau G(P^* + \gamma_3 Q^*, P + \gamma_3 Q) - \tau(1 + \gamma_3^2)G(Q^*, Q) = 0. \quad (33) $$

Since $G(P^* + \gamma_3 Q^*, P + \gamma_3 Q) \ge 0$ by (28) and $\tau\alpha \ge 0$, the left-hand side of (33) is positive if $\tau$ satisfies

$$ (1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)] - \tau\alpha|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)| - \tau(1 + \gamma_3^2)G(Q^*, Q) > 0, $$

or equivalently,

$$ 1 + \tau\alpha - \tau\alpha \frac{|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} - \tau(1 + \gamma_3^2) \frac{G(Q^*, Q)}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} > 0. \quad (34) $$

Notice that the choice of $\gamma_3$ ensures $\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda) > 0$, also

$$ 1 + \gamma_3^2 = 1 + \frac{\left(\delta_0 + \sin \frac{3\theta_0}{2}\right)^2}{\cos^2 \frac{3\theta_0}{2}} = \frac{1 + 2\delta_0 \sin \frac{3\theta_0}{2} + \delta_0^2}{\cos^2 \frac{3\theta_0}{2}} =: c $$

is a constant greater than $\delta_0^2$. By Lemma A.4 (iii) we have

$$ \frac{|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\sqrt{1 + \gamma_3^2}}{\delta_0} = \frac{\sqrt{c}}{\delta_0}. $$
and

$$ \frac{1}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \le \frac{1}{2 \sin \frac{\theta_0}{2}}. $$

Using again Lemma A.4 (ii) and (28), we have

$$ \frac{G(P^* + \gamma_2 Q^*, P + \gamma_2 Q)}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \le \frac{\|H\|^2 \|M\|^2 (p + q_1)^2}{2 \sin \frac{\theta_0}{2}}. $$

Inserting these previous inequalities in (32) gives the desired results using definitions (21), (22) and (23) of p and $q_1$. □

**Lemma 3.8 (Case 3).** Let $0 \ne \rho(B) < 1$ and let $|\lambda| \ge 1$, $\mathrm{Re}(\lambda^2 - \lambda) < 0$, $\theta \in (-\theta_0, \theta_0)$ for given $0 < \theta_0 \le \frac{\pi}{4}$. For any $\delta_0 > 0$, equation (20) cannot hold if one chooses

$$ \tau < \left( \|H\|^2 \|M\|^2 s(B)^4 \cdot \frac{2c}{\delta_0} \|B\|^2 + \left( \frac{\sqrt{c}}{\delta_0} - 1 \right) \alpha \right)^{-1} $$

where $c = c(\theta_0, \delta_0) := (1 + 2\delta_0 \sin \frac{3\theta_0}{2} + \delta_0^2) / \cos^2 \frac{3\theta_0}{2}$. Moreover, if $0 < \|B\| < 1$, the result is also true if

$$ \tau < \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|)^{-4}} \cdot \frac{2c}{\delta_0} \|B\|^2 + \left( \frac{\sqrt{c}}{\delta_0} - 1 \right) \alpha \right)^{-1}. $$

*Proof.* Define

$$ \gamma_3 = \gamma_3(\mathrm{sign}(\theta)) := \begin{cases} \left(\delta_0 + \sin \frac{3\theta_0}{2}\right) / \cos \frac{3\theta_0}{2} & \text{if } \theta > 0, \\ -\left(\delta_0 + \sin \frac{3\theta_0}{2}\right) / \cos \frac{3\theta_0}{2} & \text{if } \theta < 0 \end{cases} $$

as in Lemma A.4 (iii). Writing (27) for $\gamma = \gamma_3$ we obtain

$$ (1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha[\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)] \\ + \tau\alpha + \tau G(P^* + \gamma_3 Q^*, P + \gamma_3 Q) - \tau(1 + \gamma_3^2)G(Q^*, Q) = 0. \quad (33) $$

Since $G(P^* + \gamma_3 Q^*, P + \gamma_3 Q) \ge 0$ by (28) and $\tau\alpha \ge 0$, the left-hand side of (33) is positive if $\tau$ satisfies

$$ (1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)] - \tau\alpha|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)| - \tau(1 + \gamma_3^2)G(Q^*, Q) > 0, $$

or equivalently,

$$ 1 + \tau\alpha - \tau\alpha \frac{|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} - \tau(1 + \gamma_3^2) \frac{G(Q^*, Q)}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} > 0. \quad (34) $$

Notice that the choice of $\gamma_3$ ensures $\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda) > 0$, also

$$ 1 + \gamma_3^2 = 1 + \frac{\left(\delta_0 + \sin \frac{3\theta_0}{2}\right)^2}{\cos^2 \frac{3\theta_0}{2}} = \frac{1 + 2\delta_0 \sin \frac{3\theta_0}{2} + \delta_0^2}{\cos^2 \frac{3\theta_0}{2}} =: c $$

is a constant greater than $\delta_0^2$. By Lemma A.4 (iii) we have

$$ \frac{|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\sqrt{1 + \gamma_3^2}}{\delta_0} = \frac{\sqrt{c}}{\delta_0}. $$
Using again Lemma A.4 (iii) and (28), we have

$$ \frac{G(Q^*, Q)}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{(\|H\|\|M\|\sin\theta|q_2|^2}{2\delta_0|\sin(\theta/2)|} = \frac{2}{\delta_0} \left|\sin\frac{\theta}{2}\right| \cos^2\frac{\theta}{2} \|H\|^2 \|M\|^2 q_2^2 \\ \le \frac{2}{\delta_0} \|H\|^2 \|M\|^2 q_2^2. $$

Inserting the two previous inequalities in (34) gives the desired results using definitions (22) and (23) of $q_2$. □

### 3.3 Final result ($k = 1$)

Considering Proposition 3.3 (for $B = 0$), Proposition 3.2 (for real eigenvalues and $B \ne 0$) and taking the bound in Proposition 3.5 (for complex eigenvalues and $B \ne 0$), we obtain a sufficient condition on the descent step $\tau$ to ensure convergence of the semi-implicit one-step one-shot method.

**Theorem 3.9 (Convergence of semi-implicit one-step one-shot).** Under assumption (4), the one-step one-shot method (11) converges for sufficiently small $\tau$. In particular, for $\|B\| < 1$, there exists an explicit piecewise (at most a 4th order) polynomial function $P_1$ and a pure constant $C > 0$ such that $P_1 > 0$ on $[0, 1)$ and it is enough to take

$$ \tau < \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|^4)^4} P_1(\|B\|) + C\alpha \right)^{-1}. $$

We emphasize that the bound on $\tau$ in Theorem 3.9 depends only on $\|B\|$, $\|M\|$, $\|H\|$ and the regularization parameter $\alpha$. This bound does not depend on the dimensions of $\sigma$, $u$ and $g$.

## 4 Convergence of the multi-step one-shot method ($k \ge 1$)

We now tackle the general case of semi-implicit multi-step one-shot methods, that is algorithm (9) with $k \ge 1$. The procedure is quite similar to the case $k = 1$ but with more involved technicalities.

### 4.1 Block iteration matrix and eigenvalue equation

Let $k \ge 1$ be the number of inner iterations for $u$ and $p$. First we express $(\sigma^{n+1}, u^{n+1}, p^{n+1})$ in terms of $(\sigma^n, u^n, p^n)$ in a matrix form as for the case $k = 1$. More precisely, the system (9) can be equivalently written as

$$ \left\{ \begin{array}{l} \sigma^{n+1} = \sigma^n - \tau M^* p^n, \\ u^{n+1} = B^k u^n + T_k M \sigma^n - \tau T_k M M^* p^n + T_k F, \\ p^{n+1} = [(B^*)^k - \tau X_k M M^*] p^n + U_k u^n + X_k M \sigma^n + X_k F - T_k^* H^* g \end{array} \right. \qquad (35) $$

where

$$ T_k := I + B + \dots + B^{k-1} = (I - B)^{-1}(I - B^k), \qquad (36) $$

$$ U_k := (B^*)^{k-1} H^* H + (B^*)^{k-2} H^* H B + \dots + H^* H B^{k-1}, \quad \text{and,} \qquad (37) $$

$$ X_k := \begin{cases} (B^*)^{k-2} H^* H T_1 + (B^*)^{k-3} H^* H T_2 + \dots + H^* H T_{k-1} & \text{if } k \ge 2, \\ 0 & \text{if } k = 1. \end{cases} \qquad (38) $$

Before analyzing recursion (35), we gather in the following lemma some useful properties of $T_k$, $U_k$ and $X_k$.
Using again Lemma A.4 (iii) and (28), we have

$$ \begin{aligned} \frac{G(Q^*, Q)}{\operatorname{Re}(\lambda^2 - \lambda) + \gamma_1 \operatorname{Im}(\lambda^2 - \lambda)} &\le \frac{(\|H\|\|M\| \sin\theta |q_2|^2)}{2\delta_0 |\sin(\theta/2)|} = \frac{2}{\delta_0} \left|\sin\frac{\theta}{2}\right| \cos^2\frac{\theta}{2} \|H\|^2 \|M\|^2 q_2^2 \\ &\le \frac{2}{\delta_0} \|H\|^2 \|M\|^2 q_2^2. \end{aligned} $$

Inserting the two previous inequalities in (34) gives the desired results using definitions (22) and (23) of $q_2$. □

### 3.3 Final result ($k = 1$)

Considering Proposition 3.3 (for $B = 0$), Proposition 3.2 (for real eigenvalues and $B \neq 0$) and taking the bound in Proposition 3.5 (for complex eigenvalues and $B \neq 0$), we obtain a sufficient condition on the descent step $\tau$ to ensure convergence of the semi-implicit one-step one-shot method.

**Theorem 3.9 (Convergence of semi-implicit one-step one-shot).** Under assumption (4), the one-step one-shot method (11) converges for sufficiently small $\tau$. In particular, for $\|B\| < 1$, there exists an explicit piecewise (at most a 4th order) polynomial function $P_1$ and a pure constant $C > 0$ such that $P_1 > 0$ on $[0, 1]$ and it is enough to take

$$ \tau < \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|^4)^4} P_1(\|B\|) + C\alpha \right)^{-1}. $$

We emphasize that the bound on $\tau$ in Theorem 3.9 depends only on $\|B\|, \|M\|, \|H\|$ and the regularization parameter $\alpha$. This bound does not depend on the dimensions of $\sigma, u$ and $g$.

## 4 Convergence of the multi-step one-shot method ($k \ge 1$)

We now tackle the general case of semi-implicit multi-step one-shot methods, that is algorithm (9) with $k \ge 1$. The procedure is quite similar to the case $k = 1$ but with more involved technicalities.

### 4.1 Block iteration matrix and eigenvalue equation

Let $k \ge 1$ be the number of inner iterations for $u$ and $p$. First we express $(\sigma^{n+1}, u^{n+1}, p^{n+1})$ in terms of $(\sigma^n, u^n, p^n)$ in a matrix form as for the case $k = 1$. More precisely, the system (9) can be equivalently written as

$$ \left\{ \begin{aligned} \sigma^{n+1} &= \sigma^n - \tau M^* p^n, \\ u^{n+1} &= B^k u^n + T_k M \sigma^n - \tau T_k M M^* p^n + T_k F, \\ p^{n+1} &= [(B^*)^k - \tau X_k M M^*] p^n + U_k u^n + X_k M \sigma^n + X_k F - T_k^* H^* g \end{aligned} \right. \qquad (35) $$

where

$$ T_k := I + B + \dots + B^{k-1} = (I - B)^{-1}(I - B^k), \qquad (36) $$

$$ U_k := (B^*)^{k-1} H^* H + (B^*)^{k-2} H^* H B + \dots + H^* H B^{k-1}, \quad \text{and,} \qquad (37) $$

$$ X_k := \begin{cases} (B^*)^{k-2} H^* H T_1 + (B^*)^{k-3} H^* H T_2 + \dots + H^* H T_{k-1} & \text{if } k \ge 2, \\ 0 & \text{if } k = 1. \end{cases} \qquad (38) $$

Before analyzing recursion (35), we gather in the following lemma some useful properties of $T_k, U_k$ and $X_k$.
**Lemma 4.1.** We have the following properties.

(i) The matrices $U_k$ and $X_k$ can be respectively rewritten as

$$ 
\begin{aligned}
U_k &= \sum_{i+j=k-1} (B^*)^i H^* H B^j, \quad \forall k \ge 1, \\
X_k &= \sum_{l=0}^{k-2} \sum_{i+j=l} (B^*)^i H^* H B^j = \sum_{l=1}^{k-1} U_l, \quad \forall k \ge 2.
\end{aligned}
$$

(ii) The matrices $U_k$ and $X_k$ are self-adjoint: $U_k^* = U_k$, $X_k^* = X_k$.

(iii) We have the relation

$$ U_k T_k - X_k B^k + X_k = T_k^* H^* H T_k, \quad \forall k \ge 1. \quad (39) $$

*Proof.* Property is easy to check from the definitions (36), (37) and (38). Property (ii) straightforwardly follows from (i).

Now we prove (iii). For $k = 1$, we have $U_1 = H^*H$, $T_1 = I$ and $X_1 = 0$, hence the identity is verified. For $k \ge 2$, we remark that $X_{k+1} = B^*X_k + H^*HT_k$ by definition (38). Thanks to (ii) we know that $X_{k+1}$ is self-adjoint, hence $X_{k+1} = X_{k+1}^* = X_kB + T_k^*H^*H$. On the other hand, from we get that $X_{k+1} = X_k + U_k$. Thus,

$$ X_k + U_k = X_k B + T_k^* H^* H, \quad \text{or equivalently,} \quad U_k = X_k(B - I) + T_k^* H^* H. $$

Finally,

$$ U_k T_k = X_k(B - I)T_k + T_k^* H^* H T_k = X_k(B^k - I) + T_k^* H^* H T_k. $$

□

Now, we consider the errors $(\sigma^n - \sigma_\alpha^{\text{ex}}, u^n - u(\sigma_\alpha^{\text{ex}}), p^n - p(\sigma_\alpha^{\text{ex}}))$ with respect to the regularized solution at the n-th iteration, and, by abuse of notation, we denote them by $(\sigma^n, u^n, p^n)$. We obtain that the errors satisfy

$$ 
\left\{
\begin{aligned}
\sigma^{n+1} &= \frac{1}{1+\tau\alpha}\sigma^n - \frac{\tau}{1+\tau\alpha}M^*p^n, \\
u^{n+1} &= B^k u^n + \frac{1}{1+\tau\alpha}T_k M\sigma^n - \frac{\tau}{1+\tau\alpha}T_k M M^*p^n, \\
p^{n+1} &= \left[ (B^*)^k - \frac{\tau}{1+\tau\alpha}X_k M M^* \right] p^n + U_k u^n + \frac{1}{1+\tau\alpha}X_k M\sigma^n,
\end{aligned}
\right.
\quad (40) $$

or equivalently, by putting in evidence the block iteration matrix

$$ 
\begin{bmatrix} p^{n+1} \\ u^{n+1} \\ \sigma^{n+1} \end{bmatrix} = \begin{bmatrix} (B^*)^k - \frac{\tau}{1+\tau\alpha}X_k M M^* & U_k & \frac{1}{1+\tau\alpha}X_k M \\ -\frac{\tau}{1+\tau\alpha}T_k M M^* & B^k & \frac{1}{1+\tau\alpha}T_k M \\ -\frac{\tau}{1+\tau\alpha}M^* & 0 & \frac{1}{1+\tau\alpha}I \end{bmatrix} \begin{bmatrix} p^n \\ u^n \\ \sigma^n \end{bmatrix}. \quad (41) $$

**Proposition 4.2.** Assume that $\lambda \in \mathbb{C}$, $|\lambda| \ge 1$ is an eigenvalue of the iteration matrix in (41). If $\lambda \in \mathbb{C}$, $\lambda \notin \text{Spec}(B)$ then $\exists y \in \mathbb{C}^{n_\sigma}$, $\|y\| = 1$ such that:

$$ (1 + \tau\alpha)\lambda - 1 + \tau\lambda\langle M^*[\lambda I - (B^*)^k]^{-1}[(\lambda - 1)X_k + T_k^*H^*H T_k](\lambda I - B^k)^{-1}My, y\rangle = 0. \quad (42) $$

In particular, $\lambda = 1$ is not an eigenvalue of the iteration matrix.

The proof is similar to the proof of Proposition 3.1. The slight difference is that in the calculation we use (39) to simplify some terms and exploit the fact that $T_k$ and $(\lambda I - B^k)^{-1}$ commute.
**Lemma 4.1.** We have the following properties.

(i) The matrices $U_k$ and $X_k$ can be respectively rewritten as

$$ 
\begin{aligned}
U_k &= \sum_{i+j=k-1} (B^*)^i H^* H B^j, \quad \forall k \ge 1, \\
X_k &= \sum_{l=0}^{k-2} \sum_{i+j=l} (B^*)^i H^* H B^j = \sum_{l=1}^{k-1} U_l, \quad \forall k \ge 2.
\end{aligned}
$$

(ii) The matrices $U_k$ and $X_k$ are self-adjoint: $U_k^* = U_k$, $X_k^* = X_k$.

(iii) We have the relation

$$ U_k T_k - X_k B^k + X_k = T_k^* H^* H T_k, \quad \forall k \ge 1. \qquad (39) $$

*Proof.* Property is easy to check from the definitions (36), (37) and (38). Property (ii) straightforwardly follows from (i).

Now we prove (iii). For $k = 1$, we have $U_1 = H^*H$, $T_1 = I$ and $X_1 = 0$, hence the identity is verified. For $k \ge 2$, we remark that $X_{k+1} = B^*X_k + H^*HT_k$ by definition (38). Thanks to (ii) we know that $X_{k+1}$ is self-adjoint, hence $X_{k+1} = X_{k+1}^* = X_k B + T_k^* H^* H$. On the other hand, from we get that $X_{k+1} = X_k + U_k$. Thus,

$$ X_k + U_k = X_k B + T_k^* H^* H, \quad \text{or equivalently,} \quad U_k = X_k(B - I) + T_k^* H^* H. $$

Finally,

$$ U_k T_k = X_k(B - I)T_k + T_k^* H^* H T_k = X_k(B^k - I) + T_k^* H^* H T_k. $$

□

Now, we consider the errors $(\sigma^n - \sigma_\alpha^{\text{ex}}, u^n - u(\sigma_\alpha^{\text{ex}}), p^n - p(\sigma_\alpha^{\text{ex}}))$ with respect to the regularized solution at the n-th iteration, and, by abuse of notation, we denote them by $(\sigma^n, u^n, p^n)$. We obtain that the errors satisfy

$$ 
\left\{
\begin{aligned}
\sigma^{n+1} &= \frac{1}{1+\tau\alpha}\sigma^n - \frac{\tau}{1+\tau\alpha}M^*p^n, \\
u^{n+1} &= B^k u^n + \frac{1}{1+\tau\alpha}T_k M\sigma^n - \frac{\tau}{1+\tau\alpha}T_k M M^* p^n, \\
p^{n+1} &= \left[ (B^*)^k - \frac{\tau}{1+\tau\alpha}X_k M M^* \right] p^n + U_k u^n + \frac{1}{1+\tau\alpha}X_k M\sigma^n,
\end{aligned}
\right.
\qquad (40) $$

or equivalently, by putting in evidence the block iteration matrix

$$ 
\begin{bmatrix} p^{n+1} \\ u^{n+1} \\ \sigma^{n+1} \end{bmatrix} = \begin{bmatrix} (B^*)^k - \frac{\tau}{1+\tau\alpha}X_k M M^* & U_k & \frac{1}{1+\tau\alpha}X_k M \\ -\frac{\tau}{1+\tau\alpha}T_k M M^* & B^k & \frac{1}{1+\tau\alpha}T_k M \\ -\frac{\tau}{1+\tau\alpha}M^* & 0 & \frac{1}{1+\tau\alpha}I \end{bmatrix} \begin{bmatrix} p^n \\ u^n \\ \sigma^n \end{bmatrix}. \qquad (41) $$

**Proposition 4.2.** Assume that $\lambda \in \mathbb{C}$, $|\lambda| \ge 1$ is an eigenvalue of the iteration matrix in (41). If $\lambda \in \mathbb{C}$, $\lambda \notin \text{Spec}(B)$ then $\exists y \in \mathbb{C}^{n_\sigma}$, $\|y\| = 1$ such that:

$$ (1 + \tau\alpha)\lambda - 1 + \tau\lambda \langle M^*[\lambda I - (B^*)^k]^{-1}[(\lambda - 1)X_k + T_k^* H^* H T_k](\lambda I - B^k)^{-1}My, y \rangle = 0. \qquad (42) $$

In particular, $\lambda = 1$ is not an eigenvalue of the iteration matrix.

The proof is similar to the proof of Proposition 3.1. The slight difference is that in the calculation we use (39) to simplify some terms and exploit the fact that $T_k$ and $(\lambda I - B^k)^{-1}$ commute.
Proof. Since $\lambda \in \mathbb{C}$ is an eigenvalue of the iteration matrix, there exists a non-zero vector $(\tilde{p}, \tilde{u}, y) \in \mathbb{C}^{n_u+n_u+n_\sigma}$ such that

$$
\left\{
\begin{aligned}
\lambda y &= \frac{1}{1+\tau\alpha}y - \frac{\tau}{1+\tau\alpha}M^*\tilde{p}, \\
\lambda \tilde{u} &= B^k\tilde{u} + \frac{1}{1+\tau\alpha}T_kMy - \frac{\tau}{1+\tau\alpha}T_kMM^*\tilde{p}, \\
\lambda \tilde{p} &= [(B^*)^k - \frac{\tau}{1+\tau\alpha}X_kMM^*]\tilde{p} + U_k\tilde{u} + \frac{1}{1+\tau\alpha}X_kMy.
\end{aligned}
\right.
$$

By inserting the first equation into the second and the third equations, we simplify this system of equations as

$$
\left\{
\begin{aligned}
\lambda y &= \frac{1}{1+\tau\alpha}y - \frac{\tau}{1+\tau\alpha}M^*\tilde{p}, \\
\lambda \tilde{u} &= B^k\tilde{u} + \lambda T_k My, \\
\lambda \tilde{p} &= (B^*)^k\tilde{p} + U_k\tilde{u} + \lambda X_k My.
\end{aligned}
\right.
\qquad (43)
$$

The second equation in (43) gives us directly $\tilde{u}$ in terms of $y$:

$$ \tilde{u} = \lambda(\lambda I - B^k)^{-1}T_k My = \lambda T_k(\lambda I - B^k)^{-1}My. \quad (44) $$

The third equation in (43) gives us $\tilde{p}$ in terms of $\tilde{u}$ and $y$:

$$ \tilde{p} = [\lambda I - (B^*)^k]^{-1} U_k \tilde{u} + \lambda [\lambda I - (B^*)^k]^{-1} X_k My, $$

then by combing with the equation (44), we obtain $\tilde{p}$ in terms of $y$:

$$ \tilde{p} = \lambda(\lambda I - (B^*)^k)^{-1}V(\lambda I - B^k)^{-1}My, \quad (45) $$

where

$$ V := U_k T_k + X_k(\lambda I - B^k) = (\lambda - 1)X_k + T_k^* H^* H T_k $$

thanks to Lemma 4.1. We also see that $y \neq 0$; indeed if $y = 0$ then equations (44) and (45) give $u = 0$ and $p = 0$, that is a contradiction. Next, inserting the expression of $\tilde{p}$ in equation (45) back into the first equation in (43), we get

$$ \lambda y = \frac{1}{1 + \tau \alpha} y - \frac{\tau}{1 + \tau \alpha} \lambda M^* (\lambda I - (B^*)^k)^{-1} V (\lambda I - B^k)^{-1} My, $$

which leads to

$$ [(1 + \alpha\tau)\lambda - 1]y + \tau\lambda M^*[\lambda I - (B^*)^k]^{-1}V(\lambda I - B^k)^{-1}My = 0. \quad (46) $$

Finally, by taking scalar product of (46) with $y$, then dividing by $\|y\|^2$, we obtain (42). Now assume that $\lambda = 1$ is an eigenvalue of the iteration matrix, then (42) yields

$$ \alpha + \|H(I - B)^{-1}My\|^2 = 0, $$

which cannot be true due to the injectivity of $H(I - B)^{-1}M$.

<figure><img src="image_4.png" /></figure>

In the following sections we will show that, for sufficiently small $\tau$, equation (42) admits no solution $|\lambda| \ge 1$, thus algorithm (9) converges. When $\lambda \neq 0$, it is convenient to rewrite (42) as

$$ (1 + \tau\alpha)\lambda^2 - \lambda + \tau \langle M^* [I - (B^*)^k / \lambda]^{-1} [(\lambda - 1)X_k + T_k^* H^* H T_k] (I - B^k / \lambda)^{-1} My, y \rangle = 0. \quad (47) $$

Remark 4.3. The simple scalar case where $n_u, n_\sigma, n_g = 1$ and $\alpha = 0$ is analyzed in [2], for which necessary and sufficient conditions on $\tau$ are derived.
Proof. Since $\lambda \in \mathbb{C}$ is an eigenvalue of the iteration matrix, there exists a non-zero vector $(\tilde{p}, \tilde{u}, y) \in \mathbb{C}^{n_u+n_u+n_\sigma}$ such that

$$
\left\{
\begin{aligned}
\lambda y &= \frac{1}{1+\tau\alpha}y - \frac{\tau}{1+\tau\alpha}M^*\tilde{p}, \\
\lambda \tilde{u} &= B^k\tilde{u} + \frac{1}{1+\tau\alpha}T_kMy - \frac{\tau}{1+\tau\alpha}T_kMM^*\tilde{p}, \\
\lambda \tilde{p} &= [(B^*)^k - \frac{\tau}{1+\tau\alpha}X_kMM^*]\tilde{p} + U_k\tilde{u} + \frac{1}{1+\tau\alpha}X_kMy.
\end{aligned}
\right.
$$

By inserting the first equation into the second and the third equations, we simplify this system of equations as

$$
\left\{
\begin{aligned}
\lambda y &= \frac{1}{1+\tau\alpha}y - \frac{\tau}{1+\tau\alpha}M^*\tilde{p}, \\
\lambda \tilde{u} &= B^k \tilde{u} + \lambda T_k My, \\
\lambda \tilde{p} &= (B^*)^k \tilde{p} + U_k \tilde{u} + \lambda X_k My.
\end{aligned}
\right.
\qquad (43)
$$

The second equation in (43) gives us directly $\tilde{u}$ in terms of $y$:

$$ \tilde{u} = \lambda(\lambda I - B^k)^{-1}T_k My = \lambda T_k(\lambda I - B^k)^{-1}My. \quad (44) $$

The third equation in (43) gives us $\tilde{p}$ in terms of $\tilde{u}$ and $y$:

$$ \tilde{p} = [\lambda I - (B^*)^k]^{-1} U_k \tilde{u} + \lambda [\lambda I - (B^*)^k]^{-1} X_k My, $$

then by combining with the equation (44), we obtain $\tilde{p}$ in terms of $y$:

$$ \tilde{p} = \lambda(\lambda I - (B^*)^k)^{-1}V(\lambda I - B^k)^{-1}My, \quad (45) $$

where

$$ V := U_k T_k + X_k (\lambda I - B^k) = (\lambda - 1) X_k + T_k^* H^* H T_k $$

thanks to Lemma 4.1. We also see that $y \neq 0$; indeed if $y = 0$ then equations (44) and (45) give $u = 0$ and $p = 0$, that is a contradiction. Next, inserting the expression of $\tilde{p}$ in equation (45) back into the first equation in (43), we get

$$ \lambda y = \frac{1}{1 + \tau \alpha} y - \frac{\tau}{1 + \tau \alpha} \lambda M^* (\lambda I - (B^*)^k)^{-1} V (\lambda I - B^k)^{-1} My, $$

which leads to

$$ [(1 + \alpha\tau)\lambda - 1]y + \tau\lambda M^*[\lambda I - (B^*)^k]^{-1}V(\lambda I - B^k)^{-1}My = 0. \quad (46) $$

Finally, by taking scalar product of (46) with $y$, then dividing by $\|y\|^2$, we obtain (42). Now assume that $\lambda = 1$ is an eigenvalue of the iteration matrix, then (42) yields

$$ \alpha + \|H(I - B)^{-1}My\|^2 = 0, $$

which cannot be true due to the injectivity of $H(I - B)^{-1}M$.

<figure><img src="image_4.png" /></figure>

In the following sections we will show that, for sufficiently small $\tau$, equation (42) admits no solution $|\lambda| \ge 1$, thus algorithm (9) converges. When $\lambda \ne 0$, it is convenient to rewrite (42) as

$$ (1 + \tau\alpha)\lambda^2 - \lambda + \tau \langle M^* [I - (B^*)^k / \lambda]^{-1} [(\lambda - 1)X_k + T_k^* H^* H T_k] (I - B^k / \lambda)^{-1} My, y \rangle = 0. \quad (47) $$

Remark 4.3. The simple scalar case where $n_u, n_\sigma, n_g = 1$ and $\alpha = 0$ is analyzed in [2], for which necessary and sufficient conditions on $\tau$ are derived.
Remark 4.4. Note that when $B = 0$ and $k \ge 2$, the semi-implicit $k$-step one-shot (9) is equivalent to the semi-implicit gradient descent method (7), which converges if and only if $(\rho(A^*A) - \alpha)\tau < 2$.

For the analysis we use some auxiliary results proved in Appendix A, and the following bounds for $s(B^k), T_k, X_k$.

**Lemma 4.5.** If $\|B\| < 1$ then for every $k \ge 1$:

$$s(B^k) = s((B^*)^k) \le \frac{1}{1 - \|B\|^k}, \quad \|T_k\| \le \frac{1 - \|B\|^k}{1 - \|B\|}$$

and

$$\|X_k\| \le \frac{\|H\|^2(1 - k\|B\|^{k-1} + (k-1)\|B\|^k)}{(1 - \|B\|)^2}.$$

*Proof.* The bound for $s(B^k)$ is proved using Lemma A.2 and $\|B^k\| \le \|B\|^k$. Next, from (36) we have

$$\|T_k\| \le 1 + \|B\| + \dots + \|B\|^{k-1} = \frac{1 - \|B\|^k}{1 - \|B\|}.$$

From (38), if $k \ge 2$ we have

$$\begin{aligned} \|X_k\| &\le \|H\|^2(\|B\|^{k-2} + \|B\|^{k-3}(1 + \|B\|) + \dots + (1 + \|B\| + \dots + \|B\|^{k-2})) \\ &= \|H\|^2(1 + 2\|B\| + \dots + (k-1)\|B\|^{k-2}) \\ &= \frac{\|H\|^2(1 - k\|B\|^{k-1} + (k-1)\|B\|^k)}{(1 - \|B\|)^2}. \end{aligned}$$

□

## 4.2 Location of eigenvalues in the complex plane

We first establish conditions on the descent step $\tau > 0$ such that the real eigenvalues stay inside the unit disk. Recall that we have already proved that $\lambda = 1$ is not an eigenvalue for any $k$.

**Proposition 4.6 (Real eigenvalues).** Let $0 \ne \rho(B) < 1$ and $\lambda \in \mathbb{R}, \lambda \ne 1, |\lambda| \ge 1$. Then equation (47) cannot hold if $\tau > 0$ and

$$\left( \|M\|^2 \|X_k\| s(B^k)^2 - \frac{1}{2} \alpha \right) \tau < 1.$$

Moreover, if $\|B\| < 1$, the result is also true if $\tau > 0$ and

$$\left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|^2)^2 (1 - \|B\|^k)^2} \left(1 - k\|B\|^{k-1} + (k-1)\|B\|^k\right) - \frac{1}{2}\alpha \right) \tau < 1.$$

*Proof.* When $\lambda \in \mathbb{R}$ equation (47) can be rewritten as

$$(1 + \tau\alpha)\lambda^2 - \lambda + \tau \|H T_k \left(I - \frac{B^k}{\lambda}\right)^{-1} My\|^2 \\ + \tau(\lambda - 1) \langle M^* \left[I - \frac{(B^*)^k}{\lambda}\right]^{-1} X_k \left(I - \frac{B^k}{\lambda}\right)^{-1} My, y \rangle = 0. \quad (48)$$

We show that we can choose $\tau$ so that the left-hand side of the above equation is positive. First, we note that

$$\left| \langle M^* \left[ I - \frac{(B^*)^k}{\lambda} \right]^{-1} X_k \left( I - \frac{B^k}{\lambda} \right)^{-1} My, y \rangle \right| \le \|M\|^2 \|X_k\| s(B^k)^2.$$
Remark 4.4. Note that when $B = 0$ and $k \ge 2$, the semi-implicit $k$-step one-shot (9) is equivalent to the semi-implicit gradient descent method (7), which converges if and only if $(\rho(A^*A) - \alpha)\tau < 2$.

For the analysis we use some auxiliary results proved in Appendix A, and the following bounds for $s(B^k), T_k, X_k$.

**Lemma 4.5.** If $\|B\| < 1$ then for every $k \ge 1$:

$$s(B^k) = s((B^*)^k) \le \frac{1}{1 - \|B\|^k}, \quad \|T_k\| \le \frac{1 - \|B\|^k}{1 - \|B\|}$$

and

$$\|X_k\| \le \frac{\|H\|^2(1 - k\|B\|^{k-1} + (k-1)\|B\|^k)}{(1 - \|B\|)^2}.$$

*Proof.* The bound for $s(B^k)$ is proved using Lemma A.2 and $\|B^k\| \le \|B\|^k$. Next, from (36) we have

$$\|T_k\| \le 1 + \|B\| + \dots + \|B\|^{k-1} = \frac{1 - \|B\|^k}{1 - \|B\|}.$$

From (38), if $k \ge 2$ we have

$$\begin{aligned} \|X_k\| &\le \|H\|^2(\|B\|^{k-2} + \|B\|^{k-3}(1 + \|B\|) + \dots + (1 + \|B\| + \dots + \|B\|^{k-2})) \\ &= \|H\|^2(1 + 2\|B\| + \dots + (k-1)\|B\|^{k-2}) \\ &= \frac{\|H\|^2(1 - k\|B\|^{k-1} + (k-1)\|B\|^k)}{(1 - \|B\|)^2}. \end{aligned}$$

□

## 4.2 Location of eigenvalues in the complex plane

We first establish conditions on the descent step $\tau > 0$ such that the real eigenvalues stay inside the unit disk. Recall that we have already proved that $\lambda = 1$ is not an eigenvalue for any $k$.

**Proposition 4.6 (Real eigenvalues).** Let $0 \ne \rho(B) < 1$ and $\lambda \in \mathbb{R}, \lambda \ne 1, |\lambda| \ge 1$. Then equation (47) cannot hold if $\tau > 0$ and

$$\left( \|M\|^2 \|X_k\| s(B^k)^2 - \frac{1}{2} \alpha \right) \tau < 1.$$

Moreover, if $\|B\| < 1$, the result is also true if $\tau > 0$ and

$$\left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|)^2 (1 - \|B\|^k)^2} \left(1 - k\|B\|^{k-1} + (k-1)\|B\|^k\right) - \frac{1}{2}\alpha \right) \tau < 1.$$

*Proof.* When $\lambda \in \mathbb{R}$ equation (47) can be rewritten as

$$(1 + \tau\alpha)\lambda^2 - \lambda + \tau \|H T_k \left(I - \frac{B^k}{\lambda}\right)^{-1} My\|^2 \\ + \tau(\lambda - 1) \langle M^* \left[I - \frac{(B^*)^k}{\lambda}\right]^{-1} X_k \left(I - \frac{B^k}{\lambda}\right)^{-1} My, y \rangle = 0. \quad (48)$$

We show that we can choose $\tau$ so that the left-hand side of the above equation is positive. First, we note that

$$\left| \langle M^* \left[ I - \frac{(B^*)^k}{\lambda} \right]^{-1} X_k \left( I - \frac{B^k}{\lambda} \right)^{-1} My, y \rangle \right| \le \|M\|^2 \|X_k\| s(B^k)^2.$$
If $\lambda > 1$, we rewrite equation (48) again as

$$ (1 + \tau\alpha)\lambda(\lambda - 1) + \tau\alpha + \tau \|HT_k\left(I - \frac{B^k}{\lambda}\right)^{-1} My\|^2 \\ + \tau(\lambda - 1)\langle M^* \left[I - \frac{(B^*)^k}{\lambda}\right]^{-1} X_k \left(I - \frac{B^k}{\lambda}\right)^{-1} My, y \rangle = 0. $$

Since $\lambda(\lambda - 1) \ge \lambda - 1$, $\|HT_k\left(I - \frac{B^k}{\lambda}\right)^{-1} My\|^2 \ge 0$ and $\tau\alpha \ge 0$, we choose $\tau$ such that

$$ (1 + \tau\alpha)(\lambda - 1) - \tau(\lambda - 1)\|M\|^2 \|X_k\|s(B^k)^2 > 0, $$

or equivalently,

$$ 1 + \tau\alpha - \tau\|M\|^2 \|X_k\|s(B^k) > 0. $$

If $\lambda \le -1$, we consider equation (48). Since $\|HT_k\left(I - \frac{B^k}{\lambda}\right)^{-1} My\|^2 \ge 0$, we choose $\tau$ such that

$$ (1 + \tau\alpha)\lambda^2 - \lambda - \tau(1 - \lambda)\|M\|^2 \|X_k\|s(B^k)^2 > 0, $$

or equivalently,

$$ \frac{(1 + \tau\alpha)\lambda^2 - \lambda}{1 - \lambda} - \tau\|M\|^2 \|X_k\|s(B^k)^2 > 0, $$

Since the function $\lambda \mapsto \frac{(1+\tau\alpha)\lambda^2-\lambda}{\lambda-1}$ is decreasing on $(-\infty, -1]$, it suffices to choose $\tau$ such that

$$ 1 + \frac{\alpha}{2}\tau - \tau\|M\|^2 \|X_k\|s(B^k)^2 > 0, $$

which proves the first statement of the proposition. Finally, the case $\|B\| < 1$ can be deduced using Lemma 4.5. □

For the general case of complex eigenvalues, the study is much more complicated and technical. The following proposition summarizes the results obtained.

**Proposition 4.7 (Complex eigenvalues).** If $0 \ne \rho(B) < 1$, there exists $\tau > 0$ sufficiently small such that equation (47) admits no solution $\lambda \in \mathbb{C}\backslash\mathbb{R}$, $|\lambda| \ge 1$. In particular, if $\|B\| < 1$, given any $\delta_0 > 0$ and $0 < \theta_0 < \frac{\pi}{4}$, one can choose

$$ \tau < \min_{1 \le i \le 3} \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|)^2 (1 - \|B\|^k)^2} \psi_i(k, \|B\|) + C_i \alpha \right)^{-1}, $$

where the functions $\psi_i$ are defined by

$$ \psi_1(k, b) := 4b^{2k} + \sqrt{2}[1 - kb^{k-1} + (k-1)b^k](1 + b^k), $$

$$ \psi_2(k, b) := \left( \frac{1}{2 \sin \frac{\theta_0}{2}} (1 - b^k)^2 + \sqrt{2} (1 - kb^{k-1} + (k-1)b^k) \right) (1 + b^k)^2 $$

$$ \text{and } \psi_3(k, b) := \frac{2c \sin \frac{\theta_0}{2}}{\delta_0} b^{2k} + \frac{\sqrt{c}}{\delta_0} [1 - kb^{k-1} + (k-1)b^k](1 + b^{2k}) \\ + 2 \max \left( \frac{\sqrt{c}}{\delta_0}, \frac{\sqrt{c}}{\cos 2\theta_0} \right) [1 - kb^{k-1} + (k-1)b^k]b^k, $$

and the pure constant $C_i$, $1 \le i \le 3$, are defined by

$$ C_1 := \sqrt{2} - 1, \quad C_2 := \sqrt{2} + \frac{1}{2 \sin \frac{\theta_0}{2}} - 1, \quad \text{and} \quad C_3 := \frac{\sqrt{c}}{\delta_0} - 1 $$
If $\lambda > 1$, we rewrite equation (48) again as

$$ (1 + \tau\alpha)\lambda(\lambda - 1) + \tau\alpha + \tau \|HT_k\left(I - \frac{B^k}{\lambda}\right)^{-1} My\|^2 \\ + \tau(\lambda - 1)\langle M^* \left[I - \frac{(B^*)^k}{\lambda}\right]^{-1} X_k \left(I - \frac{B^k}{\lambda}\right)^{-1} My, y \rangle = 0. $$

Since $\lambda(\lambda - 1) \ge \lambda - 1$, $\|HT_k\left(I - \frac{B^k}{\lambda}\right)^{-1} My\|^2 \ge 0$ and $\tau\alpha \ge 0$, we choose $\tau$ such that

$$ (1 + \tau\alpha)(\lambda - 1) - \tau(\lambda - 1)\|M\|^2\|X_k\|s(B^k)^2 > 0, $$

or equivalently,

$$ 1 + \tau\alpha - \tau\|M\|^2\|X_k\|s(B^k) > 0. $$

If $\lambda \le -1$, we consider equation (48). Since $\|HT_k\left(I - \frac{B^k}{\lambda}\right)^{-1} My\|^2 \ge 0$, we choose $\tau$ such that

$$ (1 + \tau\alpha)\lambda^2 - \lambda - \tau(1 - \lambda)\|M\|^2\|X_k\|s(B^k)^2 > 0, $$

or equivalently,

$$ \frac{(1 + \tau\alpha)\lambda^2 - \lambda}{1 - \lambda} - \tau\|M\|^2\|X_k\|s(B^k)^2 > 0, $$

Since the function $\lambda \mapsto \frac{(1+\tau\alpha)\lambda^2-\lambda}{\lambda-1}$ is decreasing on $(-\infty, -1]$, it suffices to choose $\tau$ such that

$$ 1 + \frac{\alpha}{2}\tau - \tau\|M\|^2\|X_k\|s(B^k)^2 > 0, $$

which proves the first statement of the proposition. Finally, the case $\|B\| < 1$ can be deduced using Lemma 4.5. □

For the general case of complex eigenvalues, the study is much more complicated and technical. The following proposition summarizes the results obtained.

**Proposition 4.7 (Complex eigenvalues).** If $0 \ne \rho(B) < 1$, there exists $\tau > 0$ sufficiently small such that equation (47) admits no solution $\lambda \in \mathbb{C}\backslash\mathbb{R}$, $|\lambda| \ge 1$. In particular, if $\|B\| < 1$, given any $\delta_0 > 0$ and $0 < \theta_0 < \frac{\pi}{4}$, one can choose

$$ \tau < \min_{1 \le i \le 3} \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|)^2 (1 - \|B\|^k)^2} \psi_i(k, \|B\|) + C_i \alpha \right)^{-1}, $$

where the functions $\psi_i$ are defined by

$$ \psi_1(k, b) := 4b^{2k} + \sqrt{2}[1 - kb^{k-1} + (k-1)b^k](1 + b^k), $$

$$ \psi_2(k, b) := \left( \frac{1}{2 \sin \frac{\theta_0}{2}} (1 - b^k)^2 + \sqrt{2} (1 - kb^{k-1} + (k-1)b^k) \right) (1 + b^k)^2 $$

$$ \text{and } \psi_3(k, b) := \frac{2c \sin \frac{\theta_0}{2}}{\delta_0} b^{2k} + \frac{\sqrt{c}}{\delta_0} [1 - kb^{k-1} + (k-1)b^k](1 + b^{2k}) \\ + 2 \max \left( \frac{\sqrt{c}}{\delta_0}, \frac{\sqrt{c}}{\cos 2\theta_0} \right) [1 - kb^{k-1} + (k-1)b^k]b^k, $$

and the pure constant $C_i$, $1 \le i \le 3$, are defined by

$$ C_1 := \sqrt{2} - 1, \quad C_2 := \sqrt{2} + \frac{1}{2 \sin \frac{\theta_0}{2}} - 1, \quad \text{and} \quad C_3 := \frac{\sqrt{c}}{\delta_0} - 1 $$
with

$$c := \frac{1 + 2\delta_0 \sin \frac{3\theta_0}{2} + \delta_0^2}{\cos^2 \frac{3\theta_0}{2}}.$$

**Proof. Step 1.** Rewrite equation (47) by separating real and imaginary parts.

Let $\lambda = R(\cos \theta + i \sin \theta)$ in polar form where $R = |\lambda| \ge 1$ and $\theta \in (-\pi, \pi)$, $\theta \ne 0$. Write $1/\lambda = r(\cos \phi + i \sin \phi)$ in polar form where $r = 1/|\lambda| = 1/R \le 1$ and $\phi = -\theta \in (-\pi, \pi)$. By Lemma A.3 applied to $T = B^k$, we have

$$\left(I - \frac{B^k}{\lambda}\right)^{-1} = P_k(\lambda) + iQ_k(\lambda), \quad \left(I - \frac{(B^*)^k}{\lambda}\right)^{-1} = P_k(\lambda)^* + iQ_k(\lambda)^*$$

where $P_k(\lambda)$ and $Q_k(\lambda)$ are $\mathbb{C}^{n_u \times n_u}$ matrices that satisfy the following bounds in the case $\|B\| < 1$ for all $|\lambda| \ge 1$:

$$\|P_k(\lambda)\| \le p := \frac{1}{1 - \|B\|^k}, \quad (49)$$

$$\|Q_k(\lambda)\| \le q_1 := \frac{\|B\|^k}{1 - \|B\|^k} \quad \text{and} \quad \|Q_k(\lambda)\| \le q_2 |\sin \theta| \quad \text{with} \quad q_2 := \frac{\|B\|^k}{(1 - \|B\|^k)^2}. \quad (50)$$

These bounds still hold in the case $0 \ne \rho(B) < 1$ with

$$p := (1 + \|B^k\|)s(B^k)^2, \quad q_1 := \|B^k\|s(B^k)^2, \quad \text{and} \quad q_2 := \frac{\|B^k\|}{1 - \|B^k\|}. \quad (51)$$

To simplify the notation, we will not explicitly write the dependence of $P_k$ and $Q_k$ on $\lambda$. Now we rewrite (47) as

$$(1 + \tau\alpha)(\lambda^2 - \lambda) + \tau\alpha(\lambda - 1) + \tau\alpha + \tau G_k(P_k^* + iQ_k^*, P_k + iQ_k) + \tau(\lambda - 1)L_k(P_k^* + iQ_k^*, P_k + iQ_k) = 0 \quad (52)$$

where

$$G_k(X, Y) = \langle M^*XT_k^*H^*HT_kYM_y, y \rangle, \forall X, Y \in \mathbb{C}^{n_u \times n_u}$$

and

$$L_k(X, Y) = \langle M^*XX_kYM_y, y \rangle, \forall X, Y \in \mathbb{C}^{n_u \times n_u}.$$

Notice that $G_k$ is a bilinear form and $G_k(X, Y) = G_k(Y^*, X^*)^*$ so that $G_k(X, Y) + G_k(X^*, Y^*)$ is real. Similarly, $L_k$ has the same properties as $G_k$ (note that $X_k^* = X_k$ by Lemma 4.1). With these properties of $G_k$ and $L_k$, we expand (4.2) and take its real and imaginary parts, which yields

$$(1 + \tau\alpha) \operatorname{Re}(\lambda^2 - \lambda) + \tau\alpha \operatorname{Re}(\lambda - 1) + \tau\alpha + \tau G_{1,k} + \tau[\operatorname{Re}(\lambda - 1)L_{1,k} - \operatorname{Im}(\lambda - 1)L_{2,k}] = 0 \quad (53)$$

and

$$(1 + \tau\alpha) \operatorname{Im}(\lambda^2 - \lambda) + \tau\alpha \operatorname{Im}(\lambda - 1) + \tau G_{2,k} + \tau[\operatorname{Im}(\lambda - 1)L_{1,k} + \operatorname{Re}(\lambda - 1)L_{2,k}] = 0 \quad (54)$$

where

$$G_{1,k} := G_k(P_k^*, P_k) - G_k(Q_k^*, Q_k), \quad G_{2,k} := G_k(P_k^*, Q_k) + G_k(Q_k^*, P_k), \\ L_{1,k} := L_k(P_k^*, P_k) - L_k(Q_k^*, Q_k) \quad \text{and} \quad L_{2,k} := L_k(P_k^*, Q_k) + L_k(Q_k^*, P_k).$$

**Step 2.** Use a suitable combination of equations (53) and (54).

Let $\gamma \in \mathbb{R}$. Multiplying equation (54) with $\gamma$ then summing it with equation (53), we obtain:

$$(1 + \tau\alpha)[\operatorname{Re}(\lambda^2 - \lambda) + \gamma \operatorname{Im}(\lambda^2 - \lambda)] + \tau\alpha[\operatorname{Re}(\lambda - 1) + \gamma \operatorname{Im}(\lambda - 1)] + \tau\alpha \\ + \tau G_k(P_k^* + \gamma Q_k^*, P_k + \gamma Q_k) - \tau(1 + \gamma^2)G_k(Q_k^*, Q_k) \\ + \tau[\operatorname{Re}(\lambda - 1) + \gamma \operatorname{Im}(\lambda - 1)]L_{1,k} + \tau[\gamma \operatorname{Re}(\lambda - 1) - \operatorname{Im}(\lambda - 1)]L_{2,k} = 0. \quad (55)$$

Now we consider four cases of $\lambda$ as in the proof for $k = 1$ namely:
with

$$c := \frac{1 + 2\delta_0 \sin \frac{3\theta_0}{2} + \delta_0^2}{\cos^2 \frac{3\theta_0}{2}}.$$

**Proof. Step 1. Rewrite equation (47) by separating real and imaginary parts.**

Let $\lambda = R(\cos \theta + i \sin \theta)$ in polar form where $R = |\lambda| \ge 1$ and $\theta \in (-\pi, \pi)$, $\theta \ne 0$. Write $1/\lambda = r(\cos \phi + i \sin \phi)$ in polar form where $r = 1/|\lambda| = 1/R \le 1$ and $\phi = -\theta \in (-\pi, \pi)$. By Lemma A.3 applied to $T = B^k$, we have

$$\left(I - \frac{B^k}{\lambda}\right)^{-1} = P_k(\lambda) + iQ_k(\lambda), \quad \left(I - \frac{(B^*)^k}{\lambda}\right)^{-1} = P_k(\lambda)^* + iQ_k(\lambda)^*$$

where $P_k(\lambda)$ and $Q_k(\lambda)$ are $\mathbb{C}^{n_u \times n_u}$ matrices that satisfy the following bounds in the case $\|B\| < 1$ for all $|\lambda| \ge 1$:

$$\|P_k(\lambda)\| \le p := \frac{1}{1 - \|B\|^k}, \qquad (49)$$

$$\|Q_k(\lambda)\| \le q_1 := \frac{\|B\|^k}{1 - \|B\|^k} \quad \text{and} \quad \|Q_k(\lambda)\| \le q_2 |\sin \theta| \quad \text{with} \quad q_2 := \frac{\|B\|^k}{(1 - \|B\|^k)^2}. \qquad (50)$$

These bounds still hold in the case $0 \ne \rho(B) < 1$ with

$$p := (1 + \|B^k\|)s(B^k)^2, \quad q_1 := \|B^k\|s(B^k)^2, \quad \text{and} \quad q_2 := \frac{\|B^k\|}{1 - \|B^k\|}. \qquad (51)$$

To simplify the notation, we will not explicitly write the dependence of $P_k$ and $Q_k$ on $\lambda$. Now we rewrite (47) as

$$(1 + \tau\alpha)(\lambda^2 - \lambda) + \tau\alpha(\lambda - 1) + \tau\alpha + \tau G_k(P_k^* + iQ_k^*, P_k + iQ_k) + \tau(\lambda - 1)L_k(P_k^* + iQ_k^*, P_k + iQ_k) = 0 \qquad (52)$$

where

$$G_k(X, Y) = \langle M^*XT_k^*H^*HT_kYM_y, y \rangle, \forall X, Y \in \mathbb{C}^{n_u \times n_u}$$

and

$$L_k(X, Y) = \langle M^*XX_kYM_y, y \rangle, \forall X, Y \in \mathbb{C}^{n_u \times n_u}.$$

Notice that $G_k$ is a bilinear form and $G_k(X, Y) = G_k(Y^*, X^*)^*$ so that $G_k(X, Y) + G_k(X^*, Y^*)$ is real. Similarly, $L_k$ has the same properties as $G_k$ (note that $X_k^* = X_k$ by Lemma 4.1). With these properties of $G_k$ and $L_k$, we expand (4.2) and take its real and imaginary parts, which yields

$$(1 + \tau\alpha) \mathrm{Re}(\lambda^2 - \lambda) + \tau\alpha \mathrm{Re}(\lambda - 1) + \tau\alpha + \tau G_{1,k} + \tau[\mathrm{Re}(\lambda - 1)L_{1,k} - \mathrm{Im}(\lambda - 1)L_{2,k}] = 0 \qquad (53)$$

and

$$(1 + \tau\alpha) \mathrm{Im}(\lambda^2 - \lambda) + \tau\alpha \mathrm{Im}(\lambda - 1) + \tau G_{2,k} + \tau[\mathrm{Im}(\lambda - 1)L_{1,k} + \mathrm{Re}(\lambda - 1)L_{2,k}] = 0 \qquad (54)$$

where

$$G_{1,k} := G_k(P_k^*, P_k) - G_k(Q_k^*, Q_k), \quad G_{2,k} := G_k(P_k^*, Q_k) + G_k(Q_k^*, P_k), \\ L_{1,k} := L_k(P_k^*, P_k) - L_k(Q_k^*, Q_k) \quad \text{and} \quad L_{2,k} := L_k(P_k^*, Q_k) + L_k(Q_k^*, P_k).$$

**Step 2. Use a suitable combination of equations (53) and (54).**

Let $\gamma \in \mathbb{R}$. Multiplying equation (54) with $\gamma$ then summing it with equation (53), we obtain:

$$(1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha[\mathrm{Re}(\lambda - 1) + \gamma \mathrm{Im}(\lambda - 1)] + \tau\alpha \\ + \tau G_k(P_k^* + \gamma Q_k^*, P_k + \gamma Q_k) - \tau(1 + \gamma^2)G_k(Q_k^*, Q_k) \\ + \tau[\mathrm{Re}(\lambda - 1) + \gamma \mathrm{Im}(\lambda - 1)]L_{1,k} + \tau[\gamma \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)]L_{2,k} = 0. \qquad (55)$$

Now we consider four cases of $\lambda$ as in the proof for $k = 1$ namely:
* Case 1. Re($\lambda^2 - \lambda$) $\ge$ 0;

* Case 2. Re($\lambda^2 - \lambda$) < 0 and $\theta \in [\theta_0, \pi - \theta_0] \cup [-\pi + \theta_0, -\theta_0]$ for fixed $0 < \theta_0 < \frac{\pi}{4}$;

* Case 3. Re($\lambda^2 - \lambda$) < 0 and $\theta \in (-\theta_0, \theta_0)$ for fixed $0 < \theta_0 < \frac{\pi}{4}$;

* Case 4. Re($\lambda^2 - \lambda$) < 0 and $\theta \in (\pi - \theta_0, \pi) \cup (-\pi, -\pi + \theta_0)$ for fixed $0 < \theta_0 < \frac{\pi}{4}$.

Cases 1, 2 and 3 are respectively treated in Lemmas 4.8, 4.9 and 4.10 below. Notice that case 4 corresponds to an empty set, according to Lemma A.4 (iv). The statement of the proposition easily follows from the combination of Lemmas 4.8, 4.9 and 4.10. □

In the lemmas below we shall make use of the following obvious properties where $\|y\| = 1$:

$$0 \le G_k(X^*, X) = \|HT_kXM_y\|^2 \le (\|H\|\|M\|\|T_k\|\|X\|)^2, \quad \forall X \in \mathbb{C}^{n_u \times n_u}. \quad (56)$$

$$|L_{1,k}| = |L_k(P_k^*, P_k) - L_k(Q_k^*, Q_k)| \le |L_k(P_k^*, P_k)| + |L_k(Q_k^*, Q_k)| \\ \le \|X_k\|\|M\|^2(\|P_k\|^2 + \|Q_k\|^2) \le \|X_k\|\|M\|^2(p^2 + q_1^2) \quad (57)$$

and

$$|L_{2,k}| = |L_k(P_k^*, Q_k) + L_k(Q_k^*, P_k)| \le |L_k(P_k^*, Q_k)| + |L_k(Q_k^*, P_k)| \\ \le 2\|X_k\|\|M\|^2\|P_k\|\|Q_k\| \le 2\|X_k\|\|M\|^2 pq_1. \quad (58)$$

**Lemma 4.8 (Case 1).** Let $\rho(B) < 1$ and let $|\lambda| \ge 1$ with Re($\lambda^2 - \lambda$) $\ge 0$. Then equation (47) cannot hold if $\tau > 0$ and

$$(4\|H\|^2\|M\|^2\|T_k\|^2\|B^k\|^2s(B^k)^4 + \sqrt{2}\|M\|^2\|X_k\|(1 + 2\|B^k\|)^2s(B^k)^4 + (\sqrt{2} - 1)\alpha)\tau < 1.$$

Moreover, if $\|B\| < 1$, the result is also true if $\tau > 0$ and

$$\left( \frac{\|H\|^2\|M\|^2}{(1 - \|B\|)^2(1 - \|B\|^k)^2} \psi_1(k, \|B\|) + (\sqrt{2} - 1)\alpha \right) \tau < 1.$$

$$\text{where } \psi_1(k, b) := 4b^{2k} + \sqrt{2}(1 - kb^{k-1} + (k-1)b^k)(1 + b^k).$$

*Proof.* Define

$$\gamma_1 = \gamma_1(\lambda) := \begin{cases} 1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) \ge 0, \\ -1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) < 0 \end{cases}$$

as in Lemma A.4. Writing (55) for $\gamma = \gamma_1$ as in Lemma A.4 and using $\gamma_1^2 = 1$ we obtain

$$(1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha[\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)] + \tau\alpha \\ + \tau G_k(P_k^* + \gamma_1 Q_k^*, P_k + \gamma_1 Q_k) - 2\tau G_k(Q_k^*, Q_k) \\ + \tau[\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)]L_{1,k} + \tau[\gamma_1 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)]L_{2,k} = 0. \quad (59)$$

Since $G_k(P_k^* + \gamma_1 Q_k^*, P_k + \gamma_1 Q_k) \ge 0$ by (56) and $\tau\alpha \ge 0$, the left-hand side of (59) is positive if $\tau$ satisfies

$$(1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)] - \tau\alpha[\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)] - 2\tau G_k(Q_k^*, Q_k) \\ - \tau[\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)]|L_{1,k}| - \tau[\gamma_1 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)]|L_{2,k}| > 0,$$
* Case 1. Re($\lambda^2 - \lambda$) $\ge$ 0;

* Case 2. Re($\lambda^2 - \lambda$) < 0 and $\theta \in [\theta_0, \pi - \theta_0] \cup [-\pi + \theta_0, -\theta_0]$ for fixed $0 < \theta_0 < \frac{\pi}{4}$;

* Case 3. Re($\lambda^2 - \lambda$) < 0 and $\theta \in (-\theta_0, \theta_0)$ for fixed $0 < \theta_0 < \frac{\pi}{4}$;

* Case 4. Re($\lambda^2 - \lambda$) < 0 and $\theta \in (\pi - \theta_0, \pi) \cup (-\pi, -\pi + \theta_0)$ for fixed $0 < \theta_0 < \frac{\pi}{4}$.

Cases 1, 2 and 3 are respectively treated in Lemmas 4.8, 4.9 and 4.10 below. Notice that case 4 corresponds to an empty set, according to Lemma A.4 (iv). The statement of the proposition easily follows from the combination of Lemmas 4.8, 4.9 and 4.10. □

In the lemmas below we shall make use of the following obvious properties where $\|y\| = 1$:

$$0 \le G_k(X^*, X) = \|HT_kXM_y\|^2 \le (\|H\|\|M\|\|T_k\|\|X\|)^2, \quad \forall X \in \mathbb{C}^{n_u \times n_u}. \quad (56)$$

$$|L_{1,k}| = |L_k(P_k^*, P_k) - L_k(Q_k^*, Q_k)| \le |L_k(P_k^*, P_k)| + |L_k(Q_k^*, Q_k)| \\ \le \|X_k\|\|M\|^2(\|P_k\|^2 + \|Q_k\|^2) \le \|X_k\|\|M\|^2(p^2 + q_1^2) \quad (57)$$

and

$$|L_{2,k}| = |L_k(P_k^*, Q_k) + L_k(Q_k^*, P_k)| \le |L_k(P_k^*, Q_k)| + |L_k(Q_k^*, P_k)| \\ \le 2\|X_k\|\|M\|^2\|P_k\|\|Q_k\| \le 2\|X_k\|\|M\|^2pq_1. \quad (58)$$

**Lemma 4.8 (Case 1).** Let $\rho(B) < 1$ and let $|\lambda| \ge 1$ with Re($\lambda^2 - \lambda$) $\ge$ 0. Then equation (47) cannot hold if $\tau > 0$ and

$$(4\|H\|^2\|M\|^2\|T_k\|^2\|B^k\|^2s(B^k)^4 + \sqrt{2}\|M\|^2\|X_k\|(1 + 2\|B^k\|)^2s(B^k)^4 + (\sqrt{2} - 1)\alpha)\tau < 1.$$

Moreover, if $\|B\| < 1$, the result is also true if $\tau > 0$ and

$$\left( \frac{\|H\|^2\|M\|^2}{(1 - \|B\|)^2(1 - \|B\|^k)^2} \psi_1(k, \|B\|) + (\sqrt{2} - 1)\alpha \right) \tau < 1.$$

$$\text{where } \psi_1(k, b) := 4b^{2k} + \sqrt{2}(1 - kb^{k-1} + (k-1)b^k)(1 + b^k).$$

*Proof.* Define

$$\gamma_1 = \gamma_1(\lambda) := \begin{cases} 1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) \ge 0, \\ -1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) < 0 \end{cases}$$

as in Lemma A.4. Writing (55) for $\gamma = \gamma_1$ as in Lemma A.4 and using $\gamma_1^2 = 1$ we obtain

$$(1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha[\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)] + \tau\alpha \\ + \tau G_k(P_k^* + \gamma_1 Q_k^*, P_k + \gamma_1 Q_k) - 2\tau G_k(Q_k^*, Q_k) \\ + \tau[\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)]L_{1,k} + \tau[\gamma_1 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)]L_{2,k} = 0. \quad (59)$$

Since $G_k(P_k^* + \gamma_1 Q_k^*, P_k + \gamma_1 Q_k) \ge 0$ by (56) and $\tau\alpha \ge 0$, the left-hand side of (59) is positive if $\tau$ satisfies

$$(1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)] - \tau\alpha[\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)] - 2\tau G_k(Q_k^*, Q_k) \\ - \tau[\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)]|L_{1,k}| - \tau[\gamma_1 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)]|L_{2,k}| > 0,$$
or equivalently,

$$1 + \tau\alpha - \tau\alpha \frac{|\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} - 2\tau \frac{G_k(Q_k^*, Q_k)}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} \\ - \tau \frac{|\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} |L_{1,k}| - \tau \frac{|\gamma_1 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} |L_{2,k}| > 0. \quad (60)$$

Notice that the choice of $\gamma_1$ ensures $\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda) > 0$. In the following we derive upper bounds independent of $\lambda$ for the terms appearing with the negative sign in (60). By Lemma A.4 we have

$$\frac{|\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\sqrt{1 + \gamma_1^2}|\lambda - 1|}{|\lambda(\lambda - 1)|} = \frac{\sqrt{2}}{|\lambda|} \le \sqrt{2}$$

and

$$\frac{|\gamma_1 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\sqrt{1 + \gamma_1^2}|\lambda - 1|}{|\lambda(\lambda - 1)|} = \frac{\sqrt{2}}{|\lambda|} \le \sqrt{2}.$$

Using again Lemma A.4 and (58), we have

$$\frac{G_k(Q_k^*, Q_k)}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\left(\|H\|\|M\|\|T_k\| \sin\theta |q_2|^2\right)^2}{2|\sin(\theta/2)|} = 2 \left|\sin\frac{\theta}{2}\right| \cos^2\frac{\theta}{2} \|H\|^2 \|M\|^2 \|T_k\|^2 q_2^2 \\ \le 2 \|H\|^2 \|M\|^2 \|T_k\|^2 q_2^2.$$

Recall from (57) and (58) that

$$|L_{1,k}| \le \|X_k\| \|M\|^2 (p^2 + q_1^2), \quad |L_{2,k}| \le 2 \|X_k\| \|M\|^2 pq_1.$$

Inserting these previous inequalities in (60) gives the desired result thanks to expressions (49), (50) and (51) of respectively $p, q_1$ and $q_2$, and thanks to Lemma 4.5. $\square$

**Lemma 4.9 (Case 2).** Let $\rho(B) < 1$ and let $|\lambda| \ge 1$, $\mathrm{Re}(\lambda^2 - \lambda) < 0$, $\theta \in [\theta_0, \pi - \theta_0] \cup [-\pi + \theta_0, -\theta_0]$ for given $0 < \theta_0 \le \frac{\pi}{4}$. Then equation (47) cannot hold if $\tau > 0$ and

$$\left( \left( \frac{1}{2 \sin \frac{\theta_0}{2}} \|H\|^2 \|M\|^2 \|T_k\|^2 + \sqrt{2} \|M\|^2 \|X_k\| \right) (1 + 2 \|B^k\|^2)^2 s(B^k)^4 + (\sqrt{2} - 1) \alpha \right) \tau < 1.$$

Moreover, if $\|B\| < 1$, the result is also true if

$$\tau < \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|^2)^2 (1 - \|B\|^k)^2} \psi_2(k, \|B\|) + (\sqrt{2} - 1)\alpha \right)^{-1}$$

$$\text{where } \psi_2(k, b) := \left( \frac{1}{2 \sin \frac{\theta_0}{2}} (1 - b^k)^2 + \sqrt{2} (1 - k b^{k-1} + (k-1)b^k) \right) (1 + b^k)^2.$$

*Proof.* Define

$$\gamma_2 = \gamma_2(\lambda) := \begin{cases} -1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) \ge 0, \\ 1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) < 0 \end{cases}$$

as in Lemma A.4 (ii). Writing (55) for $\gamma = \gamma_2$ as in Lemma A.4 (ii) and using $\gamma_2^2 = 1$, we obtain

$$(1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha[\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)] + \tau\alpha \\ + \tau G_k(P_k^* + \gamma_2 Q_k^*, P_k + \gamma_2 Q_k) - 2\tau G_k(Q_k^*, Q_k) \\ + \tau[\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)]L_{1,k} + \tau[\gamma_2 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)]L_{2,k} = 0. \quad (61)$$
or equivalently,

$$1 + \tau\alpha - \tau\alpha \frac{|\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} - 2\tau \frac{G_k(Q_k^*, Q_k)}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} \\ - \tau \frac{|\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} |L_{1,k}| - \tau \frac{|\gamma_1 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} |L_{2,k}| > 0. \quad (60)$$

Notice that the choice of $\gamma_1$ ensures $\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda) > 0$. In the following we derive upper bounds independent of $\lambda$ for the terms appearing with the negative sign in (60). By Lemma A.4 we have

$$\frac{|\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\sqrt{1 + \gamma_1^2}|\lambda - 1|}{|\lambda(\lambda - 1)|} = \frac{\sqrt{2}}{|\lambda|} \le \sqrt{2}$$

and

$$\frac{|\gamma_1 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\sqrt{1 + \gamma_1^2}|\lambda - 1|}{|\lambda(\lambda - 1)|} = \frac{\sqrt{2}}{|\lambda|} \le \sqrt{2}.$$

Using again Lemma A.4 and (58), we have

$$\frac{G_k(Q_k^*, Q_k)}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_1 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\left(\|H\|\|M\|\|T_k\| \sin\theta |q_2|^2\right)^2}{2|\sin(\theta/2)|} = 2 \left|\sin\frac{\theta}{2}\right| \cos^2\frac{\theta}{2} \|H\|^2 \|M\|^2 \|T_k\|^2 q_2^2 \\ \le 2\|H\|^2 \|M\|^2 \|T_k\|^2 q_2^2.$$

Recall from (57) and (58) that

$$|L_{1,k}| \le \|X_k\| \|M\|^2 (p^2 + q_1^2), \quad |L_{2,k}| \le 2 \|X_k\| \|M\|^2 pq_1.$$

Inserting these previous inequalities in (60) gives the desired result thanks to expressions (49), (50) and (51) of respectively $p, q_1$ and $q_2$, and thanks to Lemma 4.5. $\square$

**Lemma 4.9 (Case 2).** Let $\rho(B) < 1$ and let $|\lambda| \ge 1$, $\mathrm{Re}(\lambda^2 - \lambda) < 0$, $\theta \in [\theta_0, \pi - \theta_0] \cup [-\pi + \theta_0, -\theta_0]$ for given $0 < \theta_0 \le \frac{\pi}{4}$. Then equation (47) cannot hold if $\tau > 0$ and

$$\left( \left( \frac{1}{2 \sin \frac{\theta_0}{2}} \|H\|^2 \|M\|^2 \|T_k\|^2 + \sqrt{2} \|M\|^2 \|X_k\| \right) (1 + 2\|B^k\|^2)^2 s(B^k)^4 + (\sqrt{2} - 1)\alpha \right) \tau < 1.$$

Moreover, if $\|B\| < 1$, the result is also true if

$$\tau < \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|^2)^2 (1 - \|B\|^k)^2} \psi_2(k, \|B\|) + (\sqrt{2} - 1)\alpha \right)^{-1}$$

$$\text{where } \psi_2(k, b) := \left( \frac{1}{2 \sin \frac{\theta_0}{2}} (1 - b^k)^2 + \sqrt{2} (1 - k b^{k-1} + (k-1)b^k) \right) (1 + b^k)^2.$$

*Proof.* Define

$$\gamma_2 = \gamma_2(\lambda) := \begin{cases} -1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) \ge 0, \\ 1 & \text{if } \mathrm{Im}(\lambda^2 - \lambda) < 0 \end{cases}$$

as in Lemma A.4 (ii). Writing (55) for $\gamma = \gamma_2$ as in Lemma A.4 (ii) and using $\gamma_2^2 = 1$, we obtain

$$(1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha[\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)] + \tau\alpha \\ + \tau G_k(P_k^* + \gamma_2 Q_k^*, P_k + \gamma_2 Q_k) - 2\tau G_k(Q_k^*, Q_k) \\ + \tau[\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)]L_{1,k} + \tau[\gamma_2 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)]L_{2,k} = 0. \quad (61)$$
Since $G_k(Q_k^*, Q_k) \ge 0$ by (56), the left-hand side of (61) is negative if $\tau$ satisfies

$$ (1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha|\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)| + \tau\alpha \\ + \tau G_k(P_k^* + \gamma_2 Q_k^*, P_k + \gamma_2 Q_k) \\ + \tau|\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)||L_{1,k}| + \tau|\gamma_1 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)||L_{2,k}| < 0. $$

or equivalently,

$$ -1 - \tau\alpha + \tau\alpha \frac{|\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)|}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} + \frac{\tau\alpha}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \\ + \tau \frac{G_k(P_k^* + \gamma_2 Q_k^*, P_k + \gamma_2 Q_k)}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} + \tau \frac{|\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)|}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} |L_{1,k}| \\ + \tau \frac{|\gamma_2 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)|}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} |L_{2,k}| < 0. \quad (62) $$

Notice that the choice of $\gamma_2$ ensures $\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda) < 0$. In the following we derive upper bounds independent of $\lambda$ for the terms appearing with the positive sign in (62). By Lemma A.4 (ii) we have

$$ \frac{|\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)|}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \le \frac{\sqrt{1 + \gamma_2^2}|\lambda - 1|}{|\lambda(\lambda - 1)|} = \frac{\sqrt{2}}{|\lambda|} \le \sqrt{2}, \\ \frac{|\gamma_2 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)|}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \le \frac{\sqrt{1 + \gamma_2^2}|\lambda - 1|}{|\lambda(\lambda - 1)|} = \frac{\sqrt{2}}{|\lambda|} \le \sqrt{2}, $$

and

$$ \frac{1}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \le \frac{1}{2 \sin \frac{\theta_0}{2}}. $$

Using again Lemma A.4 (ii) and (56), we have

$$ \frac{G_k(P_k^* + \gamma_2 Q_k^*, P_k + \gamma_2 Q_k)}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \le \frac{\|H\|^2 \|M\|^2 \|T_k\|^2 (p + q_1)^2}{2 \sin \frac{\theta_0}{2}}. $$

Recall from (57) and (58) that

$$ |L_{1,k}| \le \|X_k\| \|M\|^2 (p^2 + q_1^2), \quad |L_{2,k}| \le 2 \|X_k\| \|M\|^2 pq_1. $$

Inserting these previous inequalities in (62) gives the desired result thanks to expressions (49), (50) and (51) of respectively $p$ and $q_1$, and thanks to Lemma 4.5. $\square$

**Lemma 4.10** (Case 3). Let $\rho(B) < 1$ and let $|\lambda| \ge 1$, $\mathrm{Re}(\lambda^2 - \lambda) < 0$, $\theta \in (-\theta_0, \theta_0)$ for given $0 < \theta_0 \le \frac{\pi}{4}$. For any $\delta_0 > 0$, equation (47) cannot hold if $\tau > 0$ and

$$ \left( \left[ \frac{2c \sin \frac{\theta_0}{2}}{\delta_0} \|H\|^2 \|M\|^2 \|T_k\|^2 \|B^k\|^2 + \frac{\sqrt{c}}{\delta_0} \|M\|^2 \|X_k\| (1 + 2\|B^k\| + 2\|B^k\|^2) \right. \right. \\ \left. \left. + 2 \max \left( \frac{\sqrt{c}}{\delta_0}, \frac{\sqrt{c}}{\cos 2\theta_0} \right) \|M\|^2 \|X_k\| (\|B^k\| + \|B^k\|^2) \right] s(B^k)^2 + (\sqrt{2} - 1)\alpha \right) \tau < 1 $$

where $c = c(\theta_0, \delta_0) := (1 + 2\delta_0 \sin \frac{3\theta_0}{2} + \delta_0^2) / \cos^2 \frac{3\theta_0}{2}$. Moreover, if $\|B\| < 1$, the result is also true if $\tau > 0$ and

$$ \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|^2)(1 - \|B\|^k)^2} \psi_3(k, \|B\|) + \left( \frac{\sqrt{c}}{\delta_0} - 1 \right) \alpha \right) \tau < 1 $$

where

$$ \psi_3(k, b) := \frac{2c \sin \frac{\theta_0}{2}}{\delta_0} b^{2k} + \left( \frac{\sqrt{c}}{\delta_0} (1 + b^{2k}) + 2 \max \left( \frac{\sqrt{c}}{\delta_0}, \frac{\sqrt{c}}{\cos 2\theta_0} \right) b^k \right) (1 - k b^{k-1} + (k-1)b^k)b^k. $$
Since $G_k(Q_k^*, Q_k) \ge 0$ by (56), the left-hand side of (61) is negative if $\tau$ satisfies

$$ (1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha|\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)| + \tau\alpha \\ + \tau G_k(P_k^* + \gamma_2 Q_k^*, P_k + \gamma_2 Q_k) \\ + \tau|\mathrm{Re}(\lambda - 1) + \gamma_1 \mathrm{Im}(\lambda - 1)||L_{1,k}| + \tau|\gamma_1 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)||L_{2,k}| < 0. $$

or equivalently,

$$ -1 - \tau\alpha + \tau\alpha \frac{|\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)|}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} + \frac{\tau\alpha}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \\ + \tau \frac{G_k(P_k^* + \gamma_2 Q_k^*, P_k + \gamma_2 Q_k)}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} + \tau \frac{|\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)|}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} |L_{1,k}| \\ + \tau \frac{|\gamma_2 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)|}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} |L_{2,k}| < 0. \quad (62) $$

Notice that the choice of $\gamma_2$ ensures $\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda) < 0$. In the following we derive upper bounds independent of $\lambda$ for the terms appearing with the positive sign in (62). By Lemma A.4 (ii) we have

$$ \frac{|\mathrm{Re}(\lambda - 1) + \gamma_2 \mathrm{Im}(\lambda - 1)|}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \le \frac{\sqrt{1 + \gamma_2^2}|\lambda - 1|}{|\lambda(\lambda - 1)|} = \frac{\sqrt{2}}{|\lambda|} \le \sqrt{2}, $$

$$ \frac{|\gamma_2 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)|}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \le \frac{\sqrt{1 + \gamma_2^2}|\lambda - 1|}{|\lambda(\lambda - 1)|} = \frac{\sqrt{2}}{|\lambda|} \le \sqrt{2}, $$

and

$$ \frac{1}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \le \frac{1}{2 \sin \frac{\theta_0}{2}}. $$

Using again Lemma A.4 (ii) and (56), we have

$$ \frac{G_k(P_k^* + \gamma_2 Q_k^*, P_k + \gamma_2 Q_k)}{|\mathrm{Re}(\lambda^2 - \lambda) + \gamma_2 \mathrm{Im}(\lambda^2 - \lambda)|} \le \frac{\|H\|^2 \|M\|^2 \|T_k\|^2 (p + q_1)^2}{2 \sin \frac{\theta_0}{2}}. $$

Recall from (57) and (58) that

$$ |L_{1,k}| \le \|X_k\| \|M\|^2 (p^2 + q_1^2), \quad |L_{2,k}| \le 2 \|X_k\| \|M\|^2 pq_1. $$

Inserting these previous inequalities in (62) gives the desired result thanks to expressions (49), (50) and (51) of respectively $p$ and $q_1$, and thanks to Lemma 4.5. $\square$

**Lemma 4.10 (Case 3).** Let $\rho(B) < 1$ and let $|\lambda| \ge 1$, $\mathrm{Re}(\lambda^2 - \lambda) < 0$, $\theta \in (-\theta_0, \theta_0)$ for given $0 < \theta_0 \le \frac{\pi}{4}$. For any $\delta_0 > 0$, equation (47) cannot hold if $\tau > 0$ and

$$ \left( \left[ \frac{2c \sin \frac{\theta_0}{2}}{\delta_0} \|H\|^2 \|M\|^2 \|T_k\|^2 \|B^k\|^2 + \frac{\sqrt{c}}{\delta_0} \|M\|^2 \|X_k\| (1 + 2\|B^k\| + 2\|B^k\|^2) \right. \right. \\ \left. \left. + 2 \max \left( \frac{\sqrt{c}}{\delta_0}, \frac{\sqrt{c}}{\cos 2\theta_0} \right) \|M\|^2 \|X_k\| (\|B^k\| + \|B^k\|^2) \right] s(B^k)^2 + (\sqrt{2} - 1)\alpha \right) \tau < 1 $$

where $c = c(\theta_0, \delta_0) := (1 + 2\delta_0 \sin \frac{3\theta_0}{2} + \delta_0^2) / \cos^2 \frac{3\theta_0}{2}$. Moreover, if $\|B\| < 1$, the result is also true if $\tau > 0$ and

$$ \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|^2)(1 - \|B\|^k)^2} \psi_3(k, \|B\|) + \left( \frac{\sqrt{c}}{\delta_0} - 1 \right) \alpha \right) \tau < 1 $$

where

$$ \psi_3(k, b) := \frac{2c \sin \frac{\theta_0}{2}}{\delta_0} b^{2k} + \left( \frac{\sqrt{c}}{\delta_0} (1 + b^{2k}) + 2 \max \left( \frac{\sqrt{c}}{\delta_0}, \frac{\sqrt{c}}{\cos 2\theta_0} \right) b^k \right) (1 - k b^{k-1} + (k-1)b^k)b^k. $$
Proof. Define

$$ \gamma_3 = \gamma_3(\mathrm{sign}(\theta)) := \begin{cases} \left(\delta_0 + \sin \frac{3\theta_0}{2}\right) / \cos \frac{3\theta_0}{2} & \text{if } \theta > 0, \\ -\left(\delta_0 + \sin \frac{3\theta_0}{2}\right) / \cos \frac{3\theta_0}{2} & \text{if } \theta < 0 \end{cases} $$

as in Lemma A.4 (iii). Writing (55) for $\gamma = \gamma_3$ we obtain

$$ (1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha[\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)] + \tau\alpha \\ + \tau G_k(P_k^* + \gamma_3 Q_k^*, P_k + \gamma_3 Q_k) - (1 + \gamma_3^2)\tau G_k(Q_k^*, Q_k) \\ + \tau ([\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)]L_{1,k} + [\gamma_3 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)]L_{2,k}) = 0. \quad (63) $$

Since $G(P^* + \gamma_3 Q^*, P + \gamma_3 Q) \ge 0$ and $\tau\alpha \ge 0$, the left-hand side of (63) is positive if $\tau$ satisfies

$$ (1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)] - \tau\alpha|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)| - (1 + \gamma_3^2)\tau G_k(Q_k^*, Q_k) \\ - \tau|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)||L_{1,k}| + |\gamma_3 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)||L_{2,k}| > 0, $$

or equivalently,

$$ 1 + \tau\alpha - \tau\alpha \frac{|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} - \tau(1 + \gamma_3^2) \frac{G_k(Q_k^*, Q_k)}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} \\ - \tau \frac{|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} |L_{1,k}| - \tau \frac{|\gamma_3 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} |L_{2,k}| > 0. \quad (64) $$

Notice that the choice of $\gamma_3$ ensures $\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda) > 0$, also

$$ 1 + \gamma_3^2 = 1 + \frac{\left(\delta_0 + \sin \frac{3\theta_0}{2}\right)^2}{\cos^2 \frac{3\theta_0}{2}} = \frac{1 + 2\delta_0 \sin \frac{3\theta_0}{2} + \delta_0^2}{\cos^2 \frac{3\theta_0}{2}} =: c $$

is a constant greater than $\delta_0^2$. In the following we derive upper bounds independent of $\lambda$ for the terms appearing with the negative sign in (64). By Lemma A.4 (iii) we have

$$ \frac{|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\sqrt{1 + \gamma_3^2}}{\delta_0} = \frac{\sqrt{c}}{\delta_0} $$

and

$$ \frac{|\gamma_3 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} \le \max \left( \frac{\sqrt{1 + \gamma_3^2}}{\delta_0}, \frac{\sqrt{1 + \gamma_3^2}}{\cos 2\theta_0} \right) = \max \left( \frac{\sqrt{c}}{\delta_0}, \frac{\sqrt{c}}{\cos 2\theta_0} \right). $$

Using again Lemma A.4 (iii) and (56), we have

$$ \frac{G_k(Q_k^*, Q_k)}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\left(\|H\|\|M\|\|T_k\| \sin\theta |q_2|^2\right)^2}{2\delta_0 |\sin(\theta/2)|} = \frac{2}{\delta_0} \left|\sin\frac{\theta}{2}\right| \cos^2\frac{\theta}{2} \|H\|^2 \|M\|^2 \|T_k\|^2 q_2^2 \\ \le \frac{2}{\delta_0} \|H\|^2 \|M\|^2 \|T_k\|^2 q_2^2. $$

Recall from (57) and (58) that

$$ |L_{1,k}| \le \|X_k\| \|M\|^2(p^2 + q_1^2), \quad |L_{2,k}| \le 2 \|X_k\| \|M\|^2 pq_1. $$

Inserting these previous inequalities in (64) gives the desired result thanks to expressions (49), (50) and (51) of respectively $p, q_1$ and $q_2$, and thanks to Lemma 4.5. $\square$
Proof. Define

$$ \gamma_3 = \gamma_3(\mathrm{sign}(\theta)) := \begin{cases} \left(\delta_0 + \sin \frac{3\theta_0}{2}\right) / \cos \frac{3\theta_0}{2} & \text{if } \theta > 0, \\ -\left(\delta_0 + \sin \frac{3\theta_0}{2}\right) / \cos \frac{3\theta_0}{2} & \text{if } \theta < 0 \end{cases} $$

as in Lemma A.4 (iii). Writing (55) for $\gamma = \gamma_3$ we obtain

$$ (1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)] + \tau\alpha[\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)] + \tau\alpha \\ + \tau G_k(P_k^* + \gamma_3 Q_k^*, P_k + \gamma_3 Q_k) - (1 + \gamma_3^2)\tau G_k(Q_k^*, Q_k) \\ + \tau ([\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)]L_{1,k} + [\gamma_3 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)]L_{2,k}) = 0. \quad (63) $$

Since $G(P^* + \gamma_3 Q^*, P + \gamma_3 Q) \ge 0$ and $\tau\alpha \ge 0$, the left-hand side of (63) is positive if $\tau$ satisfies

$$ (1 + \tau\alpha)[\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)] - \tau\alpha|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)| - (1 + \gamma_3^2)\tau G_k(Q_k^*, Q_k) \\ - \tau|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)||L_{1,k}| + |\gamma_3 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)||L_{2,k}| > 0, $$

or equivalently,

$$ 1 + \tau\alpha - \tau\alpha \frac{|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} - \tau(1 + \gamma_3^2) \frac{G_k(Q_k^*, Q_k)}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} \\ - \tau \frac{|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} |L_{1,k}| - \tau \frac{|\gamma_3 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} |L_{2,k}| > 0. \quad (64) $$

Notice that the choice of $\gamma_3$ ensures $\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda) > 0$, also

$$ 1 + \gamma_3^2 = 1 + \frac{\left(\delta_0 + \sin \frac{3\theta_0}{2}\right)^2}{\cos^2 \frac{3\theta_0}{2}} = \frac{1 + 2\delta_0 \sin \frac{3\theta_0}{2} + \delta_0^2}{\cos^2 \frac{3\theta_0}{2}} =: c $$

is a constant greater than $\delta_0^2$. In the following we derive upper bounds independent of $\lambda$ for the terms appearing with the negative sign in (64). By Lemma A.4 (iii) we have

$$ \frac{|\mathrm{Re}(\lambda - 1) + \gamma_3 \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\sqrt{1 + \gamma_3^2}}{\delta_0} = \frac{\sqrt{c}}{\delta_0} $$

and

$$ \frac{|\gamma_3 \mathrm{Re}(\lambda - 1) - \mathrm{Im}(\lambda - 1)|}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} \le \max \left( \frac{\sqrt{1 + \gamma_3^2}}{\delta_0}, \frac{\sqrt{1 + \gamma_3^2}}{\cos 2\theta_0} \right) = \max \left( \frac{\sqrt{c}}{\delta_0}, \frac{\sqrt{c}}{\cos 2\theta_0} \right). $$

Using again Lemma A.4 (iii) and (56), we have

$$ \frac{G_k(Q_k^*, Q_k)}{\mathrm{Re}(\lambda^2 - \lambda) + \gamma_3 \mathrm{Im}(\lambda^2 - \lambda)} \le \frac{\left(\|H\|\|M\|\|T_k\| \sin\theta |q_2|^2\right)^2}{2\delta_0 |\sin(\theta/2)|} = \frac{2}{\delta_0} \left|\sin\frac{\theta}{2}\right| \cos^2\frac{\theta}{2} \|H\|^2 \|M\|^2 \|T_k\|^2 q_2^2 \\ \le \frac{2}{\delta_0} \|H\|^2 \|M\|^2 \|T_k\|^2 q_2^2. $$

Recall from (57) and (58) that

$$ |L_{1,k}| \le \|X_k\| \|M\|^2(p^2 + q_1^2), \quad |L_{2,k}| \le 2 \|X_k\| \|M\|^2 pq_1. $$

Inserting these previous inequalities in (64) gives the desired result thanks to expressions (49), (50) and (51) of respectively $p, q_1$ and $q_2$, and thanks to Lemma 4.5. $\square$
## 4.3 Final result ($k \ge 1$)

Considering Proposition 4.6 (for real eigenvalues) and taking the bound in Proposition 4.7 (for complex eigenvalues), we finally obtain a sufficient condition on the descent step $\tau$ to ensure convergence of the multi-step one-shot method.

**Theorem 4.11** (Convergence of semi-implicit $k$-step one-shot, $k \ge 1$). Under assumption (4), the $k$-step one-shot method (9), $k \ge 1$, converges for sufficiently small $\tau$. In particular, for $\|B\| < 1$, there exists an explicit piecewise (at most a (4k)th order) polynomial function $P_k$ and a pure constant $C > 0$ independent of $k$ such that $P_k > 0$ on $[0, 1)$ and it is enough to take $\tau > 0$ and

$$ \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|)^2 (1 - \|B\|^k)^2} P_k(\|B\|) + C\alpha \right) \tau < 1. $$

We emphasize that the bound of $\tau$ in Theorem 4.11 depends only on $\|B\|$, $\|M\|$, $\|H\|$, the number of inner iterations $k$ and the regularization parameter $\alpha$. This bound in fact does not depend on the dimensions of $\sigma$, $u$ and $g$. Also, Theorem 4.11 includes the case $B = 0$, but the bound for $\tau$ in this case is quite far from optimal. Note that the optimal bound for $\tau$ in the case $B = 0$ can be found in Proposition 3.3 (for $k = 1$) and Remark 4.4 (for $k \ge 2$).

# 5 Numerical experiments on a toy problem

In order to compare the performance of classical gradient descent algorithm (7) with the $k$-step one-shot algorithms (9), we propose the following toy model related to the inverse conductivity problem in a cavity. Given $\Omega \subset \mathbb{R}^2$ an open bounded regular domain, we consider a system of $M$ direct problems, each of which is the Helmholtz equation for the linearized scattered field $u \in H^1(\Omega)$ given by

$$ \left\{ \begin{array}{ll} \text{div}(\sigma_0 \nabla u) + \omega^2 u = \text{div}(\sigma \nabla u_0), & \text{in } \Omega, \\ u = 0, & \text{on } \partial\Omega, \end{array} \right. \qquad (65) $$

where the incident field $u_0 \in H^1(\Omega)$ satisfies

$$ \left\{ \begin{array}{ll} \text{div}(\sigma_0 \nabla u_0) + \omega^2 u_0 = 0, & \text{in } \Omega, \\ u_0 = f, & \text{on } \partial\Omega, \end{array} \right. \qquad (66) $$

with the boundary data $f = f_i \in H^{1/2}(\partial\Omega)$, $1 \le i \le M$. Here $\sigma$ and $\sigma_0 \in L^\infty(\Omega)$ are supposed to be positive functions such that the support $\bar{\Omega}_0$ of $\sigma$ is strictly included inside $\Omega$. Equation (65) is obtained by linearizing the scattered field problem using the Born approximation, and $\sigma$ is the conductivity contrast with respect to $\sigma_0$. The variational formulation of (65) is: find $u \in H_0^1(\Omega)$ such that

$$ \int_{\Omega} \sigma_0 \nabla u \cdot \nabla v \, dx - \int_{\Omega} \omega^2 uv \, dx = \int_{\Omega_0} \sigma \nabla u_0 \cdot \nabla v \, dx, \quad \forall v \in H_0^1(\Omega). \qquad (67) $$

We consider the inverse problem of retrieving $\sigma$ from measurements $g := \sigma_0 \frac{\partial u}{\partial \nu}|_{\partial\Omega}$ (we indeed collect $M$ different measurements).

By discretizing $u$ using $\mathbb{P}^1$-Lagrange finite elements on a mesh $\mathcal{T}_h(\Omega)$ of $\Omega$, and $\sigma$ by $\mathbb{P}^0$-Lagrange finite elements on a coarser mesh $\tilde{\mathcal{T}}_{h'}(\Omega_0)$ of $\Omega_0$, the discretization of (67) leads to a linear system of the form

$$ A_1 \vec{u} = A_2 \vec{\sigma}, \qquad (68) $$

where $\vec{u} \in \mathbb{R}^{n_u}$, $\vec{\sigma} \in \mathbb{R}^{n_\sigma}$ with $n_u$ denoting the number of inner nodes of $\mathcal{T}_h(\Omega)$ and $n_\sigma$ denoting the number of triangles in $\tilde{\mathcal{T}}_{h'}(\Omega_0)$. In order to rewrite the linear system in the form (1) with a controllable norm for the matrix $B$, we choose to parameterize $\sigma_0$ as $\sigma_0 = \tilde{\sigma}_0 + \delta\sigma_r$ where $\delta > 0$
## 4.3 Final result ($k \ge 1$)

Considering Proposition 4.6 (for real eigenvalues) and taking the bound in Proposition 4.7 (for complex eigenvalues), we finally obtain a sufficient condition on the descent step $\tau$ to ensure convergence of the multi-step one-shot method.

**Theorem 4.11** (Convergence of semi-implicit $k$-step one-shot, $k \ge 1$). Under assumption (4), the $k$-step one-shot method (9), $k \ge 1$, converges for sufficiently small $\tau$. In particular, for $\|B\| < 1$, there exists an explicit piecewise (at most a (4k)th order) polynomial function $P_k$ and a pure constant $C > 0$ independent of $k$ such that $P_k > 0$ on $[0, 1)$ and it is enough to take $\tau > 0$ and

$$ \left( \frac{\|H\|^2 \|M\|^2}{(1 - \|B\|)^2 (1 - \|B\|^k)^2} P_k(\|B\|) + C\alpha \right) \tau < 1. $$

We emphasize that the bound of $\tau$ in Theorem 4.11 depends only on $\|B\|$, $\|M\|$, $\|H\|$, the number of inner iterations $k$ and the regularization parameter $\alpha$. This bound in fact does not depend on the dimensions of $\sigma$, $u$ and $g$. Also, Theorem 4.11 includes the case $B = 0$, but the bound for $\tau$ in this case is quite far from optimal. Note that the optimal bound for $\tau$ in the case $B = 0$ can be found in Proposition 3.3 (for $k = 1$) and Remark 4.4 (for $k \ge 2$).

# 5 Numerical experiments on a toy problem

In order to compare the performance of classical gradient descent algorithm (7) with the $k$-step one-shot algorithms (9), we propose the following toy model related to the inverse conductivity problem in a cavity. Given $\Omega \subset \mathbb{R}^2$ an open bounded regular domain, we consider a system of $M$ direct problems, each of which is the Helmholtz equation for the linearized scattered field $u \in H^1(\Omega)$ given by

$$ \left\{ \begin{array}{ll} \text{div}(\sigma_0 \nabla u) + \omega^2 u = \text{div}(\sigma \nabla u_0), & \text{in } \Omega, \\ u = 0, & \text{on } \partial\Omega, \end{array} \right. \qquad (65) $$

where the incident field $u_0 \in H^1(\Omega)$ satisfies

$$ \left\{ \begin{array}{ll} \text{div}(\sigma_0 \nabla u_0) + \omega^2 u_0 = 0, & \text{in } \Omega, \\ u_0 = f, & \text{on } \partial\Omega, \end{array} \right. \qquad (66) $$

with the boundary data $f = f_i \in H^{1/2}(\partial\Omega)$, $1 \le i \le M$. Here $\sigma$ and $\sigma_0 \in L^\infty(\Omega)$ are supposed to be positive functions such that the support $\bar{\Omega}_0$ of $\sigma$ is strictly included inside $\Omega$. Equation (65) is obtained by linearizing the scattered field problem using the Born approximation, and $\sigma$ is the conductivity contrast with respect to $\sigma_0$. The variational formulation of (65) is: find $u \in H_0^1(\Omega)$ such that

$$ \int_{\Omega} \sigma_0 \nabla u \cdot \nabla v \, dx - \int_{\Omega} \omega^2 uv \, dx = \int_{\Omega_0} \sigma \nabla u_0 \cdot \nabla v \, dx, \quad \forall v \in H_0^1(\Omega). \qquad (67) $$

We consider the inverse problem of retrieving $\sigma$ from measurements $g := \sigma_0 \frac{\partial u}{\partial \nu}|_{\partial\Omega}$ (we indeed collect $M$ different measurements).

By discretizing $u$ using $\mathbb{P}^1$-Lagrange finite elements on a mesh $\mathcal{T}_h(\Omega)$ of $\Omega$, and $\sigma$ by $\mathbb{P}^0$-Lagrange finite elements on a coarser mesh $\tilde{\mathcal{T}}_{h'}(\Omega_0)$ of $\Omega_0$, the discretization of (67) leads to a linear system of the form

$$ A_1 \vec{u} = A_2 \vec{\sigma}, \qquad (68) $$

where $\vec{u} \in \mathbb{R}^{n_u}$, $\vec{\sigma} \in \mathbb{R}^{n_\sigma}$ with $n_u$ denoting the number of inner nodes of $\mathcal{T}_h(\Omega)$ and $n_\sigma$ denoting the number of triangles in $\tilde{\mathcal{T}}_{h'}(\Omega_0)$. In order to rewrite the linear system in the form (1) with a controllable norm for the matrix $B$, we choose to parameterize $\sigma_0$ as $\sigma_0 = \tilde{\sigma}_0 + \delta\sigma_r$ where $\delta > 0$
<figure><img src="image_5.png" /><figcaption>Figure 2: The configuration for the experiment in Section 5.1. The small circles outside the mesh for $u$ indicate the source locations $y_i$, $1 \le i \le 6$.</figcaption></figure>

is a small parameter and $\sigma_r \le 1$ is a random function. With this choice of $\sigma_0$, we can write the matrix $A_1$ as $A_1 = A_{11} + \delta A_{12}$ where $A_{11}$ corresponds to the discretization of the bilinear form $\int_{\Omega}(\tilde{\sigma}_0 \nabla u \cdot \nabla v - \omega^2 uv) dx$ on $H_0^1(\Omega) \times H_0^1(\Omega)$ and $A_{12}$ corresponds to the discretization of the bilinear form $\int_{\Omega} \sigma_r \nabla u \cdot \nabla v dx$ on $H_0^1(\Omega) \times H_0^1(\Omega)$. For $\omega^2$ not a Dirichlet eigenvalue of $-\text{div}(\tilde{\sigma}_0 \nabla u)$ in $\Omega$, the matrix $A_{11}$ is invertible and we can equivalently write (68) as

$$ \vec{u} = A_{11}^{-1}(-\delta A_{12}\vec{u} + A_2\vec{\sigma}) $$

which is in the form (1) with $B = -\delta A_{11}^{-1}A_{12}$, $M = A_{11}^{-1}A_2$ and $F = 0$. Indeed $\|B\| < 1$ for sufficiently small $\delta$. We employed the finite element library FreeFEM [14] to generate the matrices $A_{11}$, $A_{12}$, $A_2$ and the measurement operator $H \in \mathbb{R}^{n_g \times n_u}$, which is the discretization of the normal trace operator $\sigma_0 \frac{\partial u}{\partial \nu}|_{\partial \Omega}$, where $n_g$ denotes the number of nodes on $\partial \Omega$.

For the numerical tests below, we set $\omega = 2\pi$, $\tilde{\sigma}_0 = 1$, $\delta = 0.01$ and the mesh size $h = \lambda/20 = 0.05$ where $\lambda = \sqrt{\tilde{\sigma}_0}2\pi/\omega = 1$. The domain $\Omega$ is the disk of radius $R = 2\lambda$ and $\Omega_0$ is formed by three squares as in Figures 2 and 5. To generate measurements, we use $M = 6$ resulting in six different incident fields $u_0$ corresponding to imposing boundary data $f = f_i$, $1 \le i \le 6$, where $f_i(x) = Y_0(\omega|x - y_i|)$ with $y_i$ located (outside $\Omega$) on the circle of radius $R + 0.25\lambda$ (see Figure 2). The function $Y_0$ is the Bessel function of the second kind of order zero. The cost functional is then the sum of the cost functionals associated with each of the incident fields. We take as exact solution $\sigma^{\text{ex}} = 10$ in each square, and as initial guess $\sigma^0 = 12$ in each square.

## 5.1 The case of noise-free data

We first consider the case of noise-free data, without regularization ($\alpha = 0$). Here, the domain $\Omega_0$ is formed by three squares of size $\lambda/4$ distributed as shown in Figure 2. We set the mesh size $h' = \lambda/4$ so that each square is divided into two triangles (see Figure 2b), which gives $n_\sigma = 6$. The mesh used for $\Omega$ (see Figure 2a) leads to $n_u = 5849$. The boundary mesh used for generating the data $g$ (see Figure 2c) is two times coarser than the mesh for $u$ and this gives $n_g = 125$. We compare the performances of $k$-step one-shot methods (9) (which coincide with (8) in the present case $\alpha = 0$). Recall that $k$ is the number of inner iterations on the direct and adjoint problems. We consider two series of experiments.

In the first one, we study the dependence on the descent step $\tau$. In Figure 3a–3b and 3c–3d we respectively fix $k = 1$ and $k = 2$ and compare $k$-step one-shot methods with the gradient descent method. We plot in semi-log scale the value of the cost functional in terms of the (outer) iteration number $n$ in (6) and (9). We can verify that for sufficiently small $\tau$, the $k$-step one-shot
<figure><img src="image_5.png" /><figcaption>Figure 2: The configuration for the experiment in Section 5.1. The small circles outside the mesh for $u$ indicate the source locations $y_i$, $1 \le i \le 6$.</figcaption></figure>

is a small parameter and $\sigma_r \le 1$ is a random function. With this choice of $\sigma_0$, we can write the matrix $A_1$ as $A_1 = A_{11} + \delta A_{12}$ where $A_{11}$ corresponds to the discretization of the bilinear form $\int_{\Omega}(\tilde{\sigma}_0 \nabla u \cdot \nabla v - \omega^2 uv) \mathrm{d}x$ on $H_0^1(\Omega) \times H_0^1(\Omega)$ and $A_{12}$ corresponds to the discretization of the bilinear form $\int_{\Omega} \sigma_r \nabla u \cdot \nabla v \mathrm{d}x$ on $H_0^1(\Omega) \times H_0^1(\Omega)$. For $\omega^2$ not a Dirichlet eigenvalue of $-\mathrm{div}(\tilde{\sigma}_0 \nabla u)$ in $\Omega$, the matrix $A_{11}$ is invertible and we can equivalently write (68) as

$$ \vec{u} = A_{11}^{-1}(-\delta A_{12}\vec{u} + A_2\vec{\sigma}) $$

which is in the form (1) with $B = -\delta A_{11}^{-1}A_{12}$, $M = A_{11}^{-1}A_2$ and $F = 0$. Indeed $\|B\| < 1$ for sufficiently small $\delta$. We employed the finite element library FreeFEM [14] to generate the matrices $A_{11}$, $A_{12}$, $A_2$ and the measurement operator $H \in \mathbb{R}^{n_g \times n_u}$, which is the discretization of the normal trace operator $\sigma_0 \frac{\partial u}{\partial \nu}|_{\partial \Omega}$, where $n_g$ denotes the number of nodes on $\partial \Omega$.

For the numerical tests below, we set $\omega = 2\pi$, $\tilde{\sigma}_0 = 1$, $\delta = 0.01$ and the mesh size $h = \lambda/20 = 0.05$ where $\lambda = \sqrt{\tilde{\sigma}_0}2\pi/\omega = 1$. The domain $\Omega$ is the disk of radius $R = 2\lambda$ and $\Omega_0$ is formed by three squares as in Figures 2 and 5. To generate measurements, we use $M = 6$ resulting in six different incident fields $u_0$ corresponding to imposing boundary data $f = f_i$, $1 \le i \le 6$, where $f_i(x) = Y_0(\omega|x - y_i|)$ with $y_i$ located (outside $\Omega$) on the circle of radius $R + 0.25\lambda$ (see Figure 2). The function $Y_0$ is the Bessel function of the second kind of order zero. The cost functional is then the sum of the cost functionals associated with each of the incident fields. We take as exact solution $\sigma^\mathrm{ex} = 10$ in each square, and as initial guess $\sigma^0 = 12$ in each square.

## 5.1 The case of noise-free data

We first consider the case of noise-free data, without regularization ($\alpha = 0$). Here, the domain $\Omega_0$ is formed by three squares of size $\lambda/4$ distributed as shown in Figure 2. We set the mesh size $h' = \lambda/4$ so that each square is divided into two triangles (see Figure 2b), which gives $n_\sigma = 6$. The mesh used for $\Omega$ (see Figure 2a) leads to $n_u = 5849$. The boundary mesh used for generating the data $g$ (see Figure 2c) is two times coarser than the mesh for $u$ and this gives $n_g = 125$. We compare the performances of $k$-step one-shot methods (9) (which coincide with (8) in the present case $\alpha = 0$). Recall that $k$ is the number of inner iterations on the direct and adjoint problems. We consider two series of experiments.

In the first one, we study the dependence on the descent step $\tau$. In Figure 3a–3b and 3c–3d we respectively fix $k = 1$ and $k = 2$ and compare $k$-step one-shot methods with the gradient descent method. We plot in semi-log scale the value of the cost functional in terms of the (outer) iteration number $n$ in (6) and (9). We can verify that for sufficiently small $\tau$, the $k$-step one-shot
<figure><img src="image_6.png" /><figcaption>Figure 3: Convergence curves of gradient descent and $k$-step one-shot: dependence on the descent step $\tau$.</figcaption></figure>

methods converge. This is not always the case for larger value of $\tau$. In particular, for $\tau = 2$, while gradient descent and 2-step one-shot converge, 1-step one-shot diverges. Oscillations may appear on the convergence curve for certain values of $\tau$, but they gradually disappear when $\tau$ gets smaller. For sufficiently small $\tau$, the convergence curves of one-shot methods are comparable to the one of gradient descent.

In the second series of experiments, we study the dependence on the number of inner iterations $k$, for fixed $\tau$. First (Figures 4a and 4c), we investigate for which $k$ the convergence curve of $k$-step one-shot is comparable with the one of gradient descent, and, as in the previous figures, on the horizontal axis we indicate the (outer) iteration number $n$. For $\tau = 2$ (see Figure 4a), we observe that for $k = 3, 4$ the convergence curves of $k$-step one-shot are close to the one of gradient descent. Note that with 3 inner iterations the $L^2$ error between $u^n$ and the exact solution to the forward problem ranges between $59 \cdot 10^{-9}$ and $0.48334$ for different $n$ in (9); in fact, this error is rather significant at the beginning of the iteration, then it reduces as we get closer to the convergence for the parameter $\sigma$. Therefore, incomplete inner iterations on the forward problem are enough to have good precision on the solution of the inverse problem. In the particular case $\tau = 2.5$ (see Figure 4c), we observe an interesting phenomenon: when $k = 3, 4, 10$, with $k$-step one-shot the cost functional decreases even faster than with gradient descent. For larger values of $k$, for example $k = 14$, the convergence curve of one-shot method gets closer to the one of gradient descent as one may expect. Next, in Figures 4b and 4d, we display the results of the same experiment as in Figures 4a and 4c, but this time on the horizontal axis we indicate the accumulated inner iteration number, which is equal to $kn$ where $n$ is the number of outer iterations. This would allow us to compare the overall
<figure><img src="image_6.png" /><figcaption>Figure 3: Convergence curves of gradient descent and $k$-step one-shot: dependence on the descent step $\tau$.</figcaption></figure>

methods converge. This is not always the case for larger value of $\tau$. In particular, for $\tau = 2$, while gradient descent and 2-step one-shot converge, 1-step one-shot diverges. Oscillations may appear on the convergence curve for certain values of $\tau$, but they gradually disappear when $\tau$ gets smaller. For sufficiently small $\tau$, the convergence curves of one-shot methods are comparable to the one of gradient descent.

In the second series of experiments, we study the dependence on the number of inner iterations $k$, for fixed $\tau$. First (Figures 4a and 4c), we investigate for which $k$ the convergence curve of $k$-step one-shot is comparable with the one of gradient descent, and, as in the previous figures, on the horizontal axis we indicate the (outer) iteration number $n$. For $\tau = 2$ (see Figure 4a), we observe that for $k = 3, 4$ the convergence curves of $k$-step one-shot are close to the one of gradient descent. Note that with 3 inner iterations the $L^2$ error between $u^n$ and the exact solution to the forward problem ranges between $59 \cdot 10^{-9}$ and $0.48334$ for different $n$ in (9); in fact, this error is rather significant at the beginning of the iteration, then it reduces as we get closer to the convergence for the parameter $\sigma$. Therefore, incomplete inner iterations on the forward problem are enough to have good precision on the solution of the inverse problem. In the particular case $\tau = 2.5$ (see Figure 4c), we observe an interesting phenomenon: when $k = 3, 4, 10$, with $k$-step one-shot the cost functional decreases even faster than with gradient descent. For larger values of $k$, for example $k = 14$, the convergence curve of one-shot method gets closer to the one of gradient descent as one may expect. Next, in Figures 4b and 4d, we display the results of the same experiment as in Figures 4a and 4c, but this time on the horizontal axis we indicate the accumulated inner iteration number, which is equal to $kn$ where $n$ is the number of outer iterations. This would allow us to compare the overall
<figure><img src="image_7.png" /><figcaption>(a) Gradient descent and $k$-step one-shot with $\tau = 2$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) $k$-step one-shot with $\tau = 2$ in terms of the accumulated inner iteration number.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Gradient descent and $k$-step one-shot with $\tau = 2.5$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) $k$-step one-shot with $\tau = 2.5$ in terms of the accumulated inner iteration number.</figcaption></figure>

Figure 4: Convergence curves of gradient descent and $k$-step one-shot: dependence on the number of inner iterations $k$.

speed of convergence among $k$-step one-shot methods. For $\tau = 2$ (respectively for $\tau = 2.5$), $k = 3$ and $k = 4$ (respectively $k = 3$) appear to provide the fastest rate of convergence. This confirms the potential interest that this type of method would have in solving large-scale inverse problem since few inner iterations are capable of providing fast convergence.

## 5.2 The case of noisy data

We consider now the case where the measurements $g$ are corrupted with noise. More specifically, we replace the vector $g = Hu(\sigma^{\text{ex}})$ in the cost functional by the vector $g^{\varepsilon} \in \mathbb{R}^{n_g}$ with $g_i^{\varepsilon} := g_i + \varepsilon_i g_i$, where the $\varepsilon_i$ are random numbers uniformly distributed between $-\varepsilon$ and $\varepsilon$ for a noise level $\varepsilon$. In order to hit the ill-posedness of the inverse problem, we artificially increase the size of the discretization space for the parameter $\sigma$ to $n_{\sigma} = 698$ by enlarging the size of the squares (now the size of their edges equals $2\lambda/3$ and the distance of their center from the boundary equals $\lambda$, see Figure 5b). We also show in Figure 5a the mesh for $u$, which is about 20 times finer than the boundary mesh used for generating the 6 data $g^{\varepsilon}$ (see Figure 5c). In this new configuration, $n_u = 5582$ and $n_g = 12$. Notice also that the meshes of the squares in Figure 5a and 5b do not coincide.

We perform several numerical tests with different noise levels: $\varepsilon = 1\%, 3\%$ and $5\%$. We choose the regularization parameter $\alpha$ depending on the noise level $\varepsilon$ in order to optimize the accuracy of the reconstruction. This choice, which is made by trial and error, does not affect much the convergence of the algorithms (see Figures 6 and 7) and mainly reduces the size of the oscillations
<figure><img src="image_7.png" /><figcaption>Figure 4: Convergence curves of gradient descent and *k*-step one-shot: dependence on the number of inner iterations *k*.</figcaption></figure>

speed of convergence among *k*-step one-shot methods. For $\tau = 2$ (respectively for $\tau = 2.5$), $k = 3$ and $k = 4$ (respectively $k = 3$) appear to provide the fastest rate of convergence. This confirms the potential interest that this type of method would have in solving large-scale inverse problem since few inner iterations are capable of providing fast convergence.

## 5.2 The case of noisy data

We consider now the case where the measurements *g* are corrupted with noise. More specifically, we replace the vector *g* = *Hu*(*σ*<sup>ex</sup>) in the cost functional by the vector *g*<sup>ε</sup> ∈ R<sup>n<sub>g</sub></sup> with *g*<sub>i</sub><sup>ε</sup> := *g*<sub>i</sub> + *ε*<sub>i</sub>*g*<sub>i</sub>, where the *ε*<sub>i</sub> are random numbers uniformly distributed between −ε and ε for a noise level ε. In order to hit the ill-posedness of the inverse problem, we artificially increase the size of the discretization space for the parameter σ to *n*<sub>σ</sub> = 698 by enlarging the size of the squares (now the size of their edges equals 2λ/3 and the distance of their center from the boundary equals λ, see Figure 5b). We also show in Figure 5a the mesh for *u*, which is about 20 times finer than the boundary mesh used for generating the 6 data *g*<sup>ε</sup> (see Figure 5c). In this new configuration, *n*<sub>*u*</sub> = 5582 and *n*<sub>*g*</sub> = 12. Notice also that the meshes of the squares in Figure 5a and 5b do not coincide.

We perform several numerical tests with different noise levels: ε = 1%, 3% and 5%. We choose the regularization parameter α depending on the noise level ε in order to optimize the accuracy of the reconstruction. This choice, which is made by trial and error, does not affect much the convergence of the algorithms (see Figures 6 and 7) and mainly reduces the size of the oscillations
<figure><img src="image_11.png" /><figcaption>Figure 5: The configuration for the experiment in Section 5.2. The small circles outside the mesh for $u$ indicate the source locations $y_i$, $1 \le i \le 6$.</figcaption></figure>

for the reconstructed $\sigma$. The convergence curves displayed in Figure 6 for different noise levels show that the semi-implicit $k$-step one-shot methods ($k = 3, 4$) require a similar number of outer iterations as the semi-implicit gradient descent algorithm to achieve the same precision. We also see from Figure 6 that the convergence curves with $\alpha \neq 0$ look less steep than those with $\alpha = 0$. Since in the case of noisy data the convergence for the cost functional does not imply in general the accuracy of the reconstructed parameter $\sigma$, we also plot convergence curves for the relative error on $\sigma$ in Figure 7 to check the quality of the reconstruction. Also in these plots we see that the semi-implicit $k$-step one-shot methods ($k = 3, 4$) require a similar number of outer iterations as the semi-implicit gradient descent algorithm to achieve the same accuracy. Indeed, this proves the potential of these methods since only a few inner iterations are used. Moreover, the chosen values of $\alpha \neq 0$ adapted to the noise level give fairly better relative error curves. Finally, to better show the effect of the regularization, Figures 8a–8f display the final reconstructions for $\sigma$ by semi-implicit 3-step one-shot, whose relative error curves were presented in Figures 7a–7f. On the left of the color scale we also indicate the actual maximum and minimum values attained for each reconstruction; note that this range gets wider for higher noise level. These plots confirm the benefit of regularization for treating noisy data.

## 5.3 Robustness with respect to the size of the discretized problem

In the following experiment, we would like to confirm that the convergence rate of the $k$-step one-shot methods is asymptotically independent of the size of the discretized problem similarly to what is known for the gradient descent algorithms. In order to modify the number of discretization points without modifying the structure of the background media, we replace the random contribution $\sigma_r$ with the constant function $\sigma_r = 1$. This means that we choose $\sigma_0 = \tilde{\sigma}_0 + \delta$. (we fix $\delta = 0.01$). We keep the same configuration as in the experiment presented in Figure 5 in terms of number of sources, number of measurements, and values of the exact conductivity contrast and of the initial guess (i.e. $\sigma^{\text{ex}} = 10$ and $\sigma^0 = 12$ in each square), but we change the configuration of the support of the unknown conductivity contrast to investigate another test case. The exact data $g_i$, $1 \le i \le 6$ are generated using a fine mesh corresponding to a mesh size $h = \lambda/60$ and a number of degrees of freedom 157844. The data are then corrupted with $\varepsilon = 3\%$ relative random noise. For the inversion, we use two different meshes respectively corresponding to mesh sizes $h = \lambda/40$ ($n_u = 70093$) and $h = \lambda/55$ ($n_u = 128490$). This experiment is illustrated by Figure 9, where the mesh with $h = \lambda/40$ is displayed. We also show the mesh for the unknown $\sigma$, that results in $n_\sigma = 10$ unknowns. For the inversion, we choose to show the result for the semi-implicit 2-step one-shot method with fixed
<figure><img src="image_8.png" /><figcaption>Figure 5: The configuration for the experiment in Section 5.2. The small circles outside the mesh for $u$ indicate the source locations $y_i$, $1 \le i \le 6$.</figcaption></figure>

for the reconstructed $\sigma$. The convergence curves displayed in Figure 6 for different noise levels show that the semi-implicit $k$-step one-shot methods ($k = 3, 4$) require a similar number of outer iterations as the semi-implicit gradient descent algorithm to achieve the same precision. We also see from Figure 6 that the convergence curves with $\alpha \neq 0$ look less steep than those with $\alpha = 0$. Since in the case of noisy data the convergence for the cost functional does not imply in general the accuracy of the reconstructed parameter $\sigma$, we also plot convergence curves for the relative error on $\sigma$ in Figure 7 to check the quality of the reconstruction. Also in these plots we see that the semi-implicit $k$-step one-shot methods ($k = 3, 4$) require a similar number of outer iterations as the semi-implicit gradient descent algorithm to achieve the same accuracy. Indeed, this proves the potential of these methods since only a few inner iterations are used. Moreover, the chosen values of $\alpha \neq 0$ adapted to the noise level give fairly better relative error curves. Finally, to better show the effect of the regularization, Figures 8a–8f display the final reconstructions for $\sigma$ by semi-implicit 3-step one-shot, whose relative error curves were presented in Figures 7a–7f. On the left of the color scale we also indicate the actual maximum and minimum values attained for each reconstruction; note that this range gets wider for higher noise level. These plots confirm the benefit of regularization for treating noisy data.

## 5.3 Robustness with respect to the size of the discretized problem

In the following experiment, we would like to confirm that the convergence rate of the $k$-step one-shot methods is asymptotically independent of the size of the discretized problem similarly to what is known for the gradient descent algorithms. In order to modify the number of discretization points without modifying the structure of the background media, we replace the random contribution $\sigma_r$ with the constant function $\sigma_r = 1$. This means that we choose $\sigma_0 = \tilde{\sigma}_0 + \delta$. (we fix $\delta = 0.01$). We keep the same configuration as in the experiment presented in Figure 5 in terms of number of sources, number of measurements, and values of the exact conductivity contrast and of the initial guess (i.e. $\sigma^{\text{ex}} = 10$ and $\sigma^0 = 12$ in each square), but we change the configuration of the support of the unknown conductivity contrast to investigate another test case. The exact data $g_i$, $1 \le i \le 6$ are generated using a fine mesh corresponding to a mesh size $h = \lambda/60$ and a number of degrees of freedom 157844. The data are then corrupted with $\varepsilon = 3\%$ relative random noise. For the inversion, we use two different meshes respectively corresponding to mesh sizes $h = \lambda/40$ ($n_u = 70093$) and $h = \lambda/55$ ($n_u = 128490$). This experiment is illustrated by Figure 9, where the mesh with $h = \lambda/40$ is displayed. We also show the mesh for the unknown $\sigma$, that results in $n_\sigma = 10$ unknowns. For the inversion, we choose to show the result for the semi-implicit 2-step one-shot method with fixed
<figure><img src="image_12.png" /><figcaption>Figure 6: Convergence curves of semi-implicit gradient descent and k-step one-shot with different noise levels $\epsilon$ and regularization parameters $\alpha$. The descent step is $\tau = 4.7$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 7: Convergence curves for the relative error on the parameter $\sigma$ of semi-implicit gradient descent and k-step one-shot with different noise levels $\epsilon$ and regularization parameters $\alpha$. The descent step is $\tau = 4.7$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 6: Convergence curves of semi-implicit gradient descent and k-step one-shot with different noise levels $\epsilon$ and regularization parameters $\alpha$. The descent step is $\tau = 4.7$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: Convergence curves for the relative error on the parameter $\sigma$ of semi-implicit gradient descent and k-step one-shot with different noise levels $\epsilon$ and regularization parameters $\alpha$. The descent step is $\tau = 4.7$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 8: Reconstructed $\sigma$ by semi-implicit 3-step one-shot with different noise levels $\varepsilon$. The exact $\sigma$ is a constant function in each square with the value $\sigma^{\text{ex}} = 10$.</figcaption></figure>

descent step $\tau = 2$ and regularization parameter $\alpha = 2 \cdot 10^{-4}$. As shown in Figure 10, the number of outer iterations is not much affected by the size of the discrete system. We observe the same behavior for other semi-implicit $k$-step one-shot methods with $k \ge 3$.

## 5.4 Dependence of the number of outer iterations on the norm of B

The purpose of this experiment is to illustrate how the convergence speed would depend on the norm of $B$. Although we did not analyze the convergence rate, one indeed expects that the smaller the norm of $B$, the faster the convergence. This is what is observed in Figure 11, where we repeat the same experiment as in Figure 10 for $h = \lambda/40$, but we vary the values of $\delta$ ($\delta = 0.01$ and $\delta = 0.02$). For a given accuracy, the number of outer iterations for $\delta = 0.01$ is less than the one for $\delta = 0.02$.

## 6 Conclusion

We have proved sufficient conditions on the descent step for the convergence of semi-implicit multi-step one-shot methods, with a regularization parameter $\alpha \ge 0$. This complements the results obtained in our research report [2] for explicit schemes with no regularization. Although these bounds on the descent step are not optimal, to our knowledge no other bounds, explicit in the number of inner iterations, are available in the literature for multi-step one-shot methods. Furthermore, we have shown in the numerical experiments that very few inner iterations on the forward and adjoint problems may be enough to guarantee similar results as for the classical gradient descent algorithm.

These preliminary numerical results are just proof of concept since the size of the direct problem is not very large. In our future work, we shall carry out more realistic numerical investigation where the iterative solvers are based on domain decomposition methods (see e.g. [4]), which are well adapted to large-scale problems. In addition, the inner fixed point iterations could be replaced by more efficient Krylov subspace methods, such as conjugate gradient or GMRES, and one could use L-BFGS instead of gradient descent as optimization algorithm. Another interesting issue is how to
<figure><img src="image_11.png" /><figcaption>Figure 8: Reconstructed $\sigma$ by semi-implicit 3-step one-shot with different noise levels $\varepsilon$. The exact $\sigma$ is a constant function in each square with the value $\sigma^{\text{ex}} = 10$.</figcaption></figure>

descent step $\tau = 2$ and regularization parameter $\alpha = 2 \cdot 10^{-4}$. As shown in Figure 10, the number of outer iterations is not much affected by the size of the discrete system. We observe the same behavior for other semi-implicit $k$-step one-shot methods with $k \ge 3$.

## 5.4 Dependence of the number of outer iterations on the norm of B

The purpose of this experiment is to illustrate how the convergence speed would depend on the norm of $B$. Although we did not analyze the convergence rate, one indeed expects that the smaller the norm of $B$, the faster the convergence. This is what is observed in Figure 11, where we repeat the same experiment as in Figure 10 for $h = \lambda/40$, but we vary the values of $\delta$ ($\delta = 0.01$ and $\delta = 0.02$). For a given accuracy, the number of outer iterations for $\delta = 0.01$ is less than the one for $\delta = 0.02$.

## 6 Conclusion

We have proved sufficient conditions on the descent step for the convergence of semi-implicit multi-step one-shot methods, with a regularization parameter $\alpha \ge 0$. This complements the results obtained in our research report [2] for explicit schemes with no regularization. Although these bounds on the descent step are not optimal, to our knowledge no other bounds, explicit in the number of inner iterations, are available in the literature for multi-step one-shot methods. Furthermore, we have shown in the numerical experiments that very few inner iterations on the forward and adjoint problems may be enough to guarantee similar results as for the classical gradient descent algorithm.

These preliminary numerical results are just proof of concept since the size of the direct problem is not very large. In our future work, we shall carry out more realistic numerical investigation where the iterative solvers are based on domain decomposition methods (see e.g. [4]), which are well adapted to large-scale problems. In addition, the inner fixed point iterations could be replaced by more efficient Krylov subspace methods, such as conjugate gradient or GMRES, and one could use L-BFGS instead of gradient descent as optimization algorithm. Another interesting issue is how to
<figure><img src="image_15.png" /><figcaption>Figure 9: The configuration for the experiment of Sections 5.3 and 5.4. The values of the exact conductivity contrast is $\sigma^{\text{ex}} = 10$ in each square. The small circles outside the mesh for $u$ indicate the source locations $y_i$, $1 \le i \le 6$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 10: Result for the experiment in Section 5.3 displaying the convergence rate for semi-implicit 2-step one-shot and gradient descent with different mesh sizes for the state equations. The descent step is $\tau = 2$, the noise level is $\varepsilon = 3\%$ and the regularization parameter is $\alpha = 2 \cdot 10^{-4}$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 11: Result for the experiment in Section 5.4 displaying the dependence of the convergence rate for semi-implicit 2-step one-shot with respect to $\delta$. The descent step is $\tau = 2$, the noise level is $\varepsilon = 3\%$ and the regularization parameter is $\alpha = 2 \cdot 10^{-4}$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 9: The configuration for the experiment of Sections 5.3 and 5.4. The values of the exact conductivity contrast is $\sigma^{\text{ex}} = 10$ in each square. The small circles outside the mesh for $u$ indicate the source locations $y_i$, $1 \le i \le 6$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 10: Result for the experiment in Section 5.3 displaying the convergence rate for semi-implicit 2-step one-shot and gradient descent with different mesh sizes for the state equations. The descent step is $\tau = 2$, the noise level is $\varepsilon = 3\%$ and the regularization parameter is $\alpha = 2 \cdot 10^{-4}$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 11: Result for the experiment in Section 5.4 displaying the dependence of the convergence rate for semi-implicit 2-step one-shot with respect to $\delta$. The descent step is $\tau = 2$, the noise level is $\varepsilon = 3\%$ and the regularization parameter is $\alpha = 2 \cdot 10^{-4}$.</figcaption></figure>
adapt the number of inner iterations in the course of the outer iterations. Moreover, based on this
linear inverse problem study, we plan to tackle the challenging case of non-linear inverse problems,
as studied in the seminal paper [12].

References

[1] L. Audibert, H. Girardon, H. Haddar, and P. Jolivet. Inversion of eddy-current signals using a level-set method and block Krylov solvers. *SIAM Journal on Scientific Computing*, 45(3):B366-B389, 2023.

[2] M. Bonazzoli, H. Haddar, and T.A. Vu. Convergence analysis of multi-step one-shot methods for linear inverse problems. Research Report RR-9477, Inria Saclay; ENSTA ParisTech, 2022.

[3] M. Burger and W. Mühlhuber. Iterative regularization of parameter identification problems by sequential quadratic programming methods. *Inverse Problems*, 18:943–969, 2002.

[4] V. Dolean, P. Jolivet, and F. Nataf. *An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation*. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2015.

[5] N. Gauger, A. Griewank, A. Hamdi, C. Kratzenstein, E. Özkaya, and T. Slawig. Automated extension of fixed point PDE solvers for optimal design with bounded retardation. In *Constrained Optimization and Optimal Control for Partial Differential Equations, International Series of Numerical Mathematics*, pages 99–122. Springer Basel, 2012.

[6] A. Greenbaum. *Iterative Methods for Solving Linear Systems*. Number 17 in Frontiers in Applied Mathematics. Soc. for Industrial and Applied Math, Philadelphia, 1997.

[7] A. Griewank. Projected Hessians for Preconditioning in One-Step One-Shot Design Optimization. In *Large-Scale Nonlinear Optimization*, volume 83, pages 151–171. Springer US, Boston, MA, 2006. Series Title: Nonconvex Optimization and Its Applications.

[8] S. Günther, N. R. Gauger, and Q. Wang. Simultaneous single-step one-shot optimization with unsteady PDEs. *Journal of Computational and Applied Mathematics*, 294:12–22, 2016.

[9] E. Haber and U. M. Ascher. Preconditioned all-at-once methods for large, sparse parameter estimation problems. *Inverse Problems*, 17(6):1847–1864, 2001.

[10] A. Hamdi and A. Griewank. Reduced quasi-Newton method for simultaneous design and optimization. *Computational Optimization and Applications*, 49(3):521–548, 2009.

[11] A. Hamdi and A. Griewank. Properties of an augmented Lagrangian for design optimization. *Optimization Methods and Software*, 25(4):645–664, 2010.

[12] Martin Hanke, Andreas Neubauer, and Otmar Scherzer. A convergence analysis of the Landweber iteration for nonlinear ill-posed problems. *Numerische Mathematik*, 72(1):21–37, 1995.

[13] S.B. Hazra, V. Schulz, J. Brezillon, and N.R. Gauger. Aerodynamic shape optimization using simultaneous pseudo-timestepping. *Journal of Computational Physics*, 204(1):46–64, 2005.

[14] F. Hecht. New development in FreeFem++. *J. Numer. Math.*, 20(3-4):251–265, 2012.

[15] B. Kaltenbacher. Regularization based on all-at-once formulations for inverse problems. *SIAM Journal on Numerical Analysis*, 54(4):2594–2618, 2016.
adapt the number of inner iterations in the course of the outer iterations. Moreover, based on this
linear inverse problem study, we plan to tackle the challenging case of non-linear inverse problems,
as studied in the seminal paper [12].

References

[1] L. Audibert, H. Girardon, H. Haddar, and P. Jolivet. Inversion of eddy-current signals using a level-set method and block Krylov solvers. *SIAM Journal on Scientific Computing*, 45(3):B366-B389, 2023.

[2] M. Bonazzoli, H. Haddar, and T.A. Vu. Convergence analysis of multi-step one-shot methods for linear inverse problems. Research Report RR-9477, Inria Saclay; ENSTA ParisTech, 2022.

[3] M. Burger and W. Mühlhuber. Iterative regularization of parameter identification problems by sequential quadratic programming methods. *Inverse Problems*, 18:943–969, 2002.

[4] V. Dolean, P. Jolivet, and F. Nataf. *An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation*. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2015.

[5] N. Gauger, A. Griewank, A. Hamdi, C. Kratzenstein, E. Özkaya, and T. Slawig. Automated extension of fixed point PDE solvers for optimal design with bounded retardation. In *Constrained Optimization and Optimal Control for Partial Differential Equations, International Series of Numerical Mathematics*, pages 99–122. Springer Basel, 2012.

[6] A. Greenbaum. *Iterative Methods for Solving Linear Systems*. Number 17 in Frontiers in Applied Mathematics. Soc. for Industrial and Applied Math, Philadelphia, 1997.

[7] A. Griewank. Projected Hessians for Preconditioning in One-Step One-Shot Design Optimization. In *Large-Scale Nonlinear Optimization*, volume 83, pages 151–171. Springer US, Boston, MA, 2006. Series Title: Nonconvex Optimization and Its Applications.

[8] S. Günther, N. R. Gauger, and Q. Wang. Simultaneous single-step one-shot optimization with unsteady PDEs. *Journal of Computational and Applied Mathematics*, 294:12–22, 2016.

[9] E. Haber and U. M. Ascher. Preconditioned all-at-once methods for large, sparse parameter estimation problems. *Inverse Problems*, 17(6):1847–1864, 2001.

[10] A. Hamdi and A. Griewank. Reduced quasi-Newton method for simultaneous design and optimization. *Computational Optimization and Applications*, 49(3):521–548, 2009.

[11] A. Hamdi and A. Griewank. Properties of an augmented Lagrangian for design optimization. *Optimization Methods and Software*, 25(4):645–664, 2010.

[12] Martin Hanke, Andreas Neubauer, and Otmar Scherzer. A convergence analysis of the Landweber iteration for nonlinear ill-posed problems. *Numerische Mathematik*, 72(1):21–37, 1995.

[13] S.B. Hazra, V. Schulz, J. Brezillon, and N.R. Gauger. Aerodynamic shape optimization using simultaneous pseudo-timestepping. *Journal of Computational Physics*, 204(1):46–64, 2005.

[14] F. Hecht. New development in FreeFem++. *J. Numer. Math.*, 20(3-4):251–265, 2012.

[15] B. Kaltenbacher. Regularization based on all-at-once formulations for inverse problems. *SIAM Journal on Numerical Analysis*, 54(4):2594–2618, 2016.
[16] B Kaltenbacher. All-at-once versus reduced iterative methods for time dependent inverse problems. *Inverse Problems*, 33(6):064002, 2017.

[17] B. Kaltenbacher, A. Kirchner, and B. Vexler. Goal oriented adaptivity in the IRGNM for parameter identification in PDEs II: all-at-once formulations. *Inverse Problems*, 30:045002, 2014.

[18] M. Marden. *Geometry of Polynomials*. Number 3 in Mathematical Surveys and Monographs. American Math. Soc, Providence, RI, 2nd edition, 1966.

[19] Tram Thi Ngoc Nguyen. Landweber–Kaczmarz for parameter identification in time-dependent inverse problems: All-at-once versus reduced version. *Inverse Problems*, 35(3):035009, 2019.

[20] Tram Thi Ngoc Nguyen. Bi-level iterative regularization for inverse problems in nonlinear PDEs. *Inverse Problems*, 40(4):045020, 2024.

[21] E. Özkaya and N. R. Gauger. Single-step One-shot Aerodynamic Shape Optimization. In *Optimal Control of Coupled Systems of Partial Differential Equations*, volume 158, pages 191–204. Birkhäuser Basel, Basel, 2009. Series Title: International Series of Numerical Mathematics.

[22] V. Schulz and I. Gherman. One-Shot Methods for Aerodynamic Shape Optimization. In *MEGADESIGN and MegaOpt - German Initiatives for Aerodynamic Simulation and Optimization in Aircraft Design*, volume 107, pages 207–220. Springer Berlin Heidelberg, Berlin, Heidelberg, 2009. Series Title: Notes on Numerical Fluid Mechanics and Multidisciplinary Design.

[23] A. Shenoy, M. Heinkenschloss, and E. M. Cliff. Airfoil design by an all-at-once method. *International Journal of Computational Fluid Dynamics*, 11(1-2):3–25, 1997.

[24] S. Ta'asan. "One Shot" Methods for Optimal Control of Distributed Parameter Systems I: Finite Dimensional Control. Technical Report 91-2, ICASE, Hampton, 1991.

[25] S. Ta'asan, G. Kuruvila, and M. Salas. Aerodynamic design and optimization in one shot. In *30th Aerospace Sciences Meeting and Exhibit*, Reno, NV, U.S.A., 1992. American Institute of Aeronautics and Astronautics.

[26] A. Tarantola and B. Valette. Generalized nonlinear inverse problems solved using the least squares criterion. *Reviews of Geophysics*, 20(2):219–232, 1982.

[27] P.-H. Tournier, I. Aliferis, M. Bonazzoli, M. de Buhan, M. Darbas, V. Dolean, F. Hecht, P. Jolivet, I. El Kanfoud, C. Migliaccio, F. Nataf, C. Pichot, and S. Semenov. Microwave tomographic imaging of cerebrovascular accidents by using high-performance computing. *Parallel Computing*, 85:88–97, 2019.

[28] T. van Leeuwen and F. J. Herrmann. Mitigating local minima in full-waveform inversion by expanding the search space. *Geophysical Journal International*, 195(1):661–667, 2013.

[29] T. van Leeuwen and F. J. Herrmann. A penalty method for PDE-constrained optimization in inverse problems. *Inverse Problems*, 32(1):015007, 2015.

# A Some useful lemmas

We state auxiliary results about matrices like those appearing in the eigenvalue equations (20) and (47).
[16] B Kaltenbacher. All-at-once versus reduced iterative methods for time dependent inverse problems. *Inverse Problems*, 33(6):064002, 2017.

[17] B. Kaltenbacher, A. Kirchner, and B. Vexler. Goal oriented adaptivity in the IRGNM for parameter identification in PDEs II: all-at-once formulations. *Inverse Problems*, 30:045002, 2014.

[18] M. Marden. *Geometry of Polynomials*. Number 3 in Mathematical Surveys and Monographs. American Math. Soc, Providence, RI, 2nd edition, 1966.

[19] Tram Thi Ngoc Nguyen. Landweber–Kaczmarz for parameter identification in time-dependent inverse problems: All-at-once versus reduced version. *Inverse Problems*, 35(3):035009, 2019.

[20] Tram Thi Ngoc Nguyen. Bi-level iterative regularization for inverse problems in nonlinear PDEs. *Inverse Problems*, 40(4):045020, 2024.

[21] E. Özkaya and N. R. Gauger. Single-step One-shot Aerodynamic Shape Optimization. In *Optimal Control of Coupled Systems of Partial Differential Equations*, volume 158, pages 191–204. Birkhäuser Basel, Basel, 2009. Series Title: International Series of Numerical Mathematics.

[22] V. Schulz and I. Gherman. One-Shot Methods for Aerodynamic Shape Optimization. In *MEGADESIGN and MegaOpt - German Initiatives for Aerodynamic Simulation and Optimization in Aircraft Design*, volume 107, pages 207–220. Springer Berlin Heidelberg, Berlin, Heidelberg, 2009. Series Title: Notes on Numerical Fluid Mechanics and Multidisciplinary Design.

[23] A. Shenoy, M. Heinkenschloss, and E. M. Cliff. Airfoil design by an all-at-once method. *International Journal of Computational Fluid Dynamics*, 11(1-2):3–25, 1997.

[24] S. Ta'asan. "One Shot" Methods for Optimal Control of Distributed Parameter Systems I: Finite Dimensional Control. Technical Report 91-2, ICASE, Hampton, 1991.

[25] S. Ta'asan, G. Kuruvila, and M. Salas. Aerodynamic design and optimization in one shot. In *30th Aerospace Sciences Meeting and Exhibit*, Reno, NV, U.S.A., 1992. American Institute of Aeronautics and Astronautics.

[26] A. Tarantola and B. Valette. Generalized nonlinear inverse problems solved using the least squares criterion. *Reviews of Geophysics*, 20(2):219–232, 1982.

[27] P.-H. Tournier, I. Aliferis, M. Bonazzoli, M. de Buhan, M. Darbas, V. Dolean, F. Hecht, P. Jolivet, I. El Kanfoud, C. Migliaccio, F. Nataf, C. Pichot, and S. Semenov. Microwave tomographic imaging of cerebrovascular accidents by using high-performance computing. *Parallel Computing*, 85:88–97, 2019.

[28] T. van Leeuwen and F. J. Herrmann. Mitigating local minima in full-waveform inversion by expanding the search space. *Geophysical Journal International*, 195(1):661–667, 2013.

[29] T. van Leeuwen and F. J. Herrmann. A penalty method for PDE-constrained optimization in inverse problems. *Inverse Problems*, 32(1):015007, 2015.

## A Some useful lemmas

We state auxiliary results about matrices like those appearing in the eigenvalue equations (20) and (47).
**Lemma A.1.** Let $(\mathbb{C}^{n \times n}, \|\cdot\|)$ be a normed space and $T \in \mathbb{C}^{n \times n}$. If $\rho(T) < 1$, then

$$ \sum_{k=0}^{\infty} T^k \text{ converges and } \sum_{k=0}^{\infty} T^k = (I - T)^{-1}. $$

Moreover, if $\|T\| < 1$, $\|(I - T)^{-1}\| \le \frac{1}{1 - \|T\|}$.

**Lemma A.2.** Let $T \in \mathbb{C}^{n \times n}$ such that $\rho(T) < 1$. Set

$$ s(T) := \sup_{z \in \mathbb{C}, |z| \ge 1} \|(I - T/z)^{-1}\| \quad (69) $$

then $0 < \|(I - T)^{-1}\| \le s(T) = s(T^*) < +\infty$. Moreover, if $\|T\| < 1$, $0 < s(T) \le \frac{1}{1 - \|T\|}$.

*Proof.* The existence of $s(T)$ (and also $s(T^*)$) is deduced from the fact that the functional $z \mapsto \|(I - T/z)^{-1}\|$, with $z \in \mathbb{C}, |z| \ge 1$, is well-defined and continuous. For every $z \in \mathbb{C}, |z| \ge 1$ we have

$$ \|(I - T/z)^{-1}\| = \|((I - T/z)^{-1})^*\| = \|(I - T^*/z^*)^{-1}\| \le s(T^*) $$

and

$$ \|(I - T^*/z)^{-1}\| = \|((I - T^*/z)^{-1})^*\| = \|(I - T/z^*)^{-1}\| \le s(T), $$

thus $s(T) = s(T^*)$. The second part of conclusion is obtained by Lemma A.1. □

The following lemma says that, for $T \in \mathbb{C}^{n \times n}$ and $\lambda \in \mathbb{C}, |\lambda| \ge 1$, we can decompose $(I - \frac{T}{\lambda})^{-1} = P(\lambda) + iQ(\lambda)$ and $(I - \frac{T^*}{\lambda})^{-1} = P(\lambda)^* + iQ(\lambda)^*$, and gives bounds for $P(\lambda)$ and $Q(\lambda)$.

**Lemma A.3.** Let $T \in \mathbb{C}^{n \times n}$ such that $\rho(T) < 1$ and $\lambda \in \mathbb{C}, |\lambda| \ge 1$. Write $\frac{1}{\lambda} = r(\cos \phi + i \sin \phi)$ in polar form, where $0 < r \le 1$ and $\phi \in [-\pi, \pi]$. Then

$$ \left(I - \frac{T}{\lambda}\right)^{-1} = P(\lambda) + iQ(\lambda) \quad \text{and} \quad \left(I - \frac{T^*}{\lambda}\right)^{-1} = P(\lambda)^* + iQ(\lambda)^* $$

where

$$ P(\lambda) = (I - r \cos \phi T)(I - 2r \cos \phi T + r^2 T^2)^{-1} $$

and

$$ Q(\lambda) = r \sin \phi T (I - 2r \cos \phi T + r^2 T^2)^{-1} $$

are $\mathbb{C}^{n \times n}$-valued functions. We also have the following properties:

$$ (i) \quad \|P(\lambda)\| \le (1 + \|T\|) s(T)^2 \quad \text{and} \quad \|Q(\lambda)\| \le |\sin \phi| \|T\| s(T)^2 \le \|T\| s(T)^2. $$

(ii) Moreover if $\|T\| < 1$ then

$$ \|P(\lambda)\| \le \frac{1}{1 - \|T\|} \quad \text{and} \quad \|Q(\lambda)\| \le \frac{\|T\|}{1 - \|T\|}. $$

*Proof.* The first part of the lemma is verified by direct computation, using

$$ \begin{align*} (I - T/\lambda)^{-1} &= (I - T/\lambda^*) [(I - T/\lambda)(I - T/\lambda^*)]^{-1}, \\ (I - T^*/\lambda)^{-1} &= [(I - T^*/\lambda^*)(I - T^*/\lambda)]^{-1} (I - T^*/\lambda^*) \\ \text{and } (I - T/\lambda)(I - T/\lambda^*) &= I - 2r \cos \phi T + r^2 T^2. \end{align*} $$