Bi-level iterative regularization for inverse problems in nonlinear PDEs

We investigate the inverse problem of recovering the spatially dependent parameter θ in the nonlinear evolution system

$$\begin{align} \begin{aligned} &\dot{u}\left(t\right)+f\left(t,\theta,u\left(t\right)\right) = 0\qquad t\in\left(0,T\right)\\ &u\left(0\right) = A\theta. \end{aligned} \end{align} \tag{ 1 }$$

with nonlinear model f and linear initial condition A; here, $\dot{u}$ denotes time derivative of u. The parameter identification is carried out with additional linear measurement of the state u

$$\begin{align} y = Lu, \end{align} \tag{ 2 }$$

which is usually contaminated by noise of some level δ, resulting in y^δ with $\|y^\delta-y\|\unicode{x2A7D}\delta$.

The parameter identification (1) and (2) can be formulated in a reduced setting with the parameter-to-observation map G as the forward operator

$$\begin{align} G\left(\theta\right) = y\qquad{\text{with}}\qquad G = L\circ S:X\to\mathcal{Y}, \quad \theta\mapsto y. \end{align} \tag{ 3 }$$

Essentially, S is a well-defined parameter-to-solution map whose output is the unique solution of the nonlinear PDE (1)

$$\begin{align} S:X\to \mathcal{V},\quad \theta\mapsto u\quad\text{solving PDE }\left(1\right), \end{align} \tag{ 4 }$$

while L is the linear observation operator

$$\begin{align} L:\left(\mathcal{V}\subseteq\right)\,\mathcal{U}\to\mathcal{Y}, \quad u\mapsto y. \end{align} \tag{ 5 }$$

This abstract setting is the framework considered throughout the paper.

1.1.1. Regularization with state approximation error.

1.1.2. Introducing the lower-level.

Adapting our notation to the framework of iterative methods, we call the iteration for estimating u the lower-level iteration, and the iteration for reconstructing θ the upper-level iteration. To this end, we first recast the lower level problem as an inverse problem in which the unknown is the state u. The forward map, at a given parameter θ, is the PDE model (1)

$$\begin{align} F\left(\theta\right):\mathcal{V}\to \mathcal{U}^*\times H \qquad F\left(\theta\right)\left(u\right): = \left(\dot{u}+f\left(t,\theta,u\left(t\right)\right);A\theta-u|_{t = 0}\right) \end{align} \tag{ 6 }$$

with the nonlinear model $f:X\times\mathcal{V}\to \mathcal{U}^*$ and the initial condition $A:X\to H$. We aim at minimizing the PDE residual towards zero, that is solving

$$\begin{align} F\left(\theta\right)\left(u\right) = \varphi: = \left(0,0\right) \end{align} \tag{ 7 }$$

with the image $\varphi\in\mathcal{U}^*\times H$ simply being zero. Seeking a good approximation $\widetilde{S}(\theta)$ for the exact state $S(\theta)$ in (4) is equivalent to finding an $\tilde{u}\in\mathcal{V}$ such that the PDE residual (7) is sufficiently small. Here, the choice of function space plays an important role and will be discussed in detail later.

The state approximation error $\epsilon: = \|\tilde{u}-u\|$ in the lower-level iteration must be handled with care: Allowing it to be large saves computational time, yet at the same time, it must be kept small enough such that when it enters the upper-level iteration, one can still obtain an acceptable reconstruction result for θ.

1.1.3. Bi-level as a mixed approach.

The suggested approach could be seen as a combined version of the reduced formulation (3)–(5) and an all-at-once approach. All-at-once formulations is a recent approach introduced in optimization with PDE constraint. This concept was subsequently studied in the context of inverse problems [6, 7, 16, 27, 28, 48], and more recently specifically for time dependent inverse problems [32, 41]. In this setting, one constructs, instead of the reduced map G (3), the higher-dimensional all-at-once map $\mathbb{G}$ and

$$\begin{align} \begin{aligned} &\text{solve:}\quad \mathbb{G}\left(\theta,u\right) = \left(0,0,y\right)\\ &\text{with}\quad\mathbb{G}:X\times\mathcal{V}\to \mathcal{U}^*\times H\times\mathcal{Y}, \quad \mathbb{G}\left(\theta,u\right): = \left(\dot{u}+f\left(t,\theta,u\left(t\right)\right);A\theta-u|_{t = 0};Lu\right). \end{aligned} \end{align} \tag{ 8 }$$

By treating u as an independent unknown, one jointly recovers the parameter and the state $(\theta,u)$, bypassing the step of solving the nonlinear PDE (1) exactly. This feature makes the all-at-once approach powerful in practice. A comparison between these formulations can be found in [32, 41].

We remark that the classical reduced approach usually brings in certain structural conditions on the model f to induce unique existence for the PDE (1). The new all-at-once approach (8), as explained, bypasses the need to study unique existence as well as solving the PDE exactly. This raises the question: Can one utilize accessible properties of f, and still bypass the exact solution for any nonlinear PDE? This is the viewpoint that we take in this study. In particular, we employ coercivity of the PDE model to ensure convergence of the lower-level iteration. Accordingly, as with the all-at-once approach, no nonlinear PDE needs to be solved exactly. Essentially, bi-level as a mixed approach combining the advantages of the reduced setting and the all-at-once setting is the new perspective of this work.

1.1.4. Connection and contrast to other directions.

1.1.5. Contribution.

1.1.6. Structure.

1.2.1. Setting

1.2.2. Notation

We now establish the bi-level Landweber algorithm 1. The main idea is to first perform upper-level iteration for θ by the standard Landweber iteration with the exact forward map $G = L\circ S$. We then replace the parameter-to-state map S by an approximation $\widetilde{S}_{K(j)}$ obtained through a lower-level Landweber iteration. $\widetilde{S}_{K(j)}$ enters the forward and backward (adjoint) evaluation in the upper-level iteration, yielding a gradient update step for θ. After Algorithm 1 is constructed, we perform some preliminary error analysis in section 2.2.

First, recalling the set up of the forward map $G(\theta) = L\circ S(\theta)$ in (3)–(5), the standard Landweber iteration for θ in this case reads as

$$\begin{align*} \theta^{\delta,j+1} = \theta^{\delta,j} -S^{{\prime}}\left(\theta^{\delta,j}\right)^*L^*\left(LS\left(\theta^{\delta,j}\right)-y^\delta\right) \qquad j\unicode{x2A7D} j^*\left(\delta\right). \end{align*}$$

Replacing S by its approximation $\widetilde{S}_{K(j)}$, the Landweber iteration with $\widetilde{S}_{K(j)}$ becomes

$$\begin{align} \widetilde{\theta}^{\delta,j+1} = \widetilde{\theta}^{\delta,j} -S^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)^*L^*\left(L\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)-y^\delta\right) \qquad j\unicode{x2A7D}\widetilde{j}^*\left(\delta\right), \end{align} \tag{ 9 }$$

where at each j, that is, each $\widetilde{\theta}^{\delta,j}$, the state approximation $\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j}): = u^{\,j}_{K(j)}$ is the outcome of the lower-level iteration

$$\begin{align} u^{\,j}_{k+1} = u^{\,j}_k-F_u^{{\prime}}\left(\widetilde{\theta}^{\delta,j},u^{\,j}_k\right)^*\left(F\left(\widetilde{\theta}^{\delta,j},u^{\,j}_k\right)-\varphi\right) \qquad k< K\left(j\right) \end{align} \tag{ 10 }$$

with the PDE map $F(\widetilde{\theta}^{\delta,j})$ described in (6) and (7). The stopping rule K(j) for the lower-level iteration as well as $\widetilde{j}^*$ for the upper-level iteration will be the subject of the later sections.

In order to see how the state approximation enters the adjoint $S^{{\prime}}(\widetilde{\theta}^{\delta,j})^*$, we first of all derive the adjoint process.

Proposition 1 (Adjoint). Let $D_U:U\to U^*$ and $I_X:X^*\to X$ be isomorphisms.

Then, the Hilbert space adjoint of Gʹ with G defined in (3)–(5) is given by

$$\begin{align} &G^{{\prime}}\left(\theta\right)^*:\mathcal{Y}\to X && G^{{\prime}}\left(\theta\right)^* = S^{{\prime}}\left(\theta\right)^*L^* \end{align} \tag{ 11 }$$

$$\begin{align} &L^*:\mathcal{Y} \to \mathcal{U} \nonumber\\ &S^{{\prime}}\left(\theta\right)^*:\mathcal{U}\to X && S^{{\prime}}\left(\theta\right)^*v = -\int_0^T I_X\;f^{{\prime}}_\theta\left(\theta,u\right)\left(t\right)^\star z\left(t\right)\,\textrm{d}t +I_XA^\star z\left(0\right), \end{align} \tag{ 12 }$$

where z is the solution to the final value problem at source D_U v

$$\begin{align} \begin{cases} &-\dot{z}\left(t\right)+f^{{\prime}}_u\left(\theta,u\right)\left(t\right)^\star z\left(t\right) = D_Uv\left(t\right) \qquad t\in\left(0,T\right)\\ &z\left(T\right) = 0. \end{cases} \end{align} \tag{ 13 }$$

Here, the Banach space adjoints are

$$\begin{align*} &f^{{\prime}}_u\left(t,\theta,u\left(t\right)\right):U\to U^*,\qquad &&f^{{\prime}}_\theta\left(t,\theta,u\left(t\right)\right):X\to U^*, \qquad && A:X\to H = L^2\left(\Omega\right), \\ &f^{{\prime}}_u\left(t,\theta,u\left(t\right)\right)^\star:U^{**} = U\to U^*,&& f^{{\prime}}_\theta\left(t,\theta,u\left(t\right)\right)^\star:U\to X^*, && A^\star:H^* = H\to X^*. \end{align*}$$

Proof. As the parameter-to-state map $S:X\to\mathcal{V}$ in (4) is assumed to differentiable, $S^{{\prime}}(\theta)\xi = :p\in\mathcal{V}$ solves the sensitivity (or linearized) equation

$$\begin{align} \begin{aligned} &\dot{p}\left(t\right)+f^{{\prime}}_u\left(\theta,u\right)\left(t\right)p\left(t\right) = -f^{{\prime}}_\theta\left(\theta,u\right)\left(t\right)\xi \qquad t\in\left(0,T\right)\\ &p\left(0\right) = A\xi, \end{aligned} \end{align} \tag{ 14 }$$

where $u = S(\theta)$ is the solution to the original PDE (1). Note that $p\in\mathcal{V}$, with $\mathcal{V}$ being smooth as defined in section 1.2.1, along with the Nemytskii assumption on f ensure $\dot{p}+f^{{\prime}}_u(\theta,u)p\in {\mathcal{U}^*}$ and $f^{{\prime}}_\theta(\theta,u)\xi\in {\mathcal{U}^*}$, confirming well-definedness of the sensitivity equation.

In the following computation, we will respectively use (13), integration by parts and (14) with the helps of isomorphisms to translate inner products to dual pairings. For every $\xi\in X, v\in\mathcal{U}$, one has

$$\begin{align*} \left(S^{{\prime}}\left(\theta\right)\xi,v\right)_\mathcal{U} & = \int_0^T\left( p, v\right)_U\textrm{d}t = \int_0^T\langle\, p, D_Uv\rangle_{U,U^*}\textrm{d}t\\ & = \int_0^T\int_\Omega \langle\, p,-\dot{z}+f^{{\prime}}_u\left(\theta,u\right)^\star z \rangle_{U,U^*}\,\textrm{d}t\,\textrm{d}x\\ & = \int_0^T\int_\Omega \left(\dot{p}+f^{{\prime}}_u\left(\theta,u\right)p \right)z\,\textrm{d}t\,\textrm{d}x + \int_\Omega A\xi\, z\left(0\right)\,\textrm{d}x\\ & = \int_0^T \langle \dot{p}+f^{{\prime}}_u\left(\theta,u\right)p,z\rangle_{U^*,U} \textrm{d}t + \left(A\xi, z\left(0\right)\right)_H\\ & = \int_0^T \langle-f^{{\prime}}_\theta\left(\theta,u\right)\xi, z\rangle_{U^*,U}\,\textrm{d}t + \langle \xi, A^\star z\left(0\right)\rangle_{X,X^*}\\ & = \int_0^T \langle\xi, -f^{{\prime}}_\theta\left(\theta,u\right)^\star z\rangle_{X,X^*}\,\textrm{d}t + \langle \xi, A^\star z\left(0\right)\rangle_{X,X^*}\\ & = \left(\xi, I_X\int_0^T -f^{{\prime}}_\theta\left(\theta,u\right)^\star z + I_XA^\star z\left(0\right)\right)_X = :\left(\xi, S^{{\prime}}\left(\theta\right)^*v\right)_X. \end{align*}$$

This proves that the adjoints of $S^{{\prime}}(\theta)$ are as shown in (12). Applying the chain rule, we obtain the adjoint of the composition map

$$\begin{align*}G^{{\prime}}\left(\theta\right)^* = \left(L\circ S^{{\prime}}\left(\theta\right)\right)^* = S^{{\prime}}\left(\theta\right)^*L^*\end{align*}$$

as claimed in (11). □

With the adjoints detailed, we can now state the Landweber iterations for the exact problem (single-level) and the approximated problem (bi-level).

Remark 1. In the bi-level algorithm, no exact S appears; equivalently, no nonlinear PDE needs to be exactly solved. The state approximation $\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j})$ enters the residual $(L\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j})-y^\delta)$ as well as the adjoint (20). We emphasize that the term $S^{{\prime}}(\widetilde{\theta}^{\delta,j})^*$ appearing here is not, in fact, the adjoint of the derivative of the exact S evaluated at $\widetilde{\theta}^{\delta,j}$, but an approximation thereof.

At this stage, the Landweber algorithms for single-level and bi-level have been established. The two algorithms differ in the state and the adjoint as mentioned in remark 1. We thus examine the system output error and the adjoint error with respect to a given approximation error $\epsilon_{K(j)}$; this is the focus of this section.

The state approximation error $\epsilon_{K(j)}$ will be quantified in the succeeding section. The subscript K(j) alludes to the fact that we will eventually derive a stopping rule K of the lower-level with respect to the upper-level iterate j. For now, we assume that K(j) is given.

Lemma 1 (System output error). Denote the state approximation error $\epsilon_{K(j)}$ at the input $\widetilde{\theta}^{\delta,j}\in B_R(\theta^\dagger)$ by

$$\begin{align} \epsilon_{K\left(j\right)}: = \left\|S\left(\widetilde{\theta}^{\delta,j}\right)-\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)\right\|_\mathcal{U}, \end{align} \tag{ 21 }$$

and assume that there exists a constant $M_S\gt0$ such that

$$\begin{align} \|S^{{\prime}}\left(\theta\right)\|_{X\to\mathcal{U}}\unicode{x2A7D} M_S \qquad \forall \theta\in B_R\left(\theta^\dagger\right). \end{align} \tag{ 22 }$$

Given any other input $\theta^{\,j}\in B_R(\theta^\dagger)$, the system output error is

$$\begin{align} \|LS\left(\theta^{\,j}\right)-L\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)\|\unicode{x2A7D} M_S\|L\|\|\theta^{\delta,j}-\widetilde{\theta}^{\delta,j}\| + \|L\|\epsilon_{K\left(j\right)}. \end{align} \tag{ 23 }$$

Proof. Differentiability of $S:X \to \mathcal{V}$, thus of $S:X \to \mathcal{U}$ due to $\mathcal{V}\hookrightarrow\mathcal{U}$, and application of the mean value theorem imply

$$\begin{align*} S\left(y\right)-S\left(x\right) = \int_0^1S^{{\prime}}\left(x+\lambda\left(y-x\right)\right)\,\textrm{d}\lambda\left(y-x\right) \end{align*}$$

meaning

$$\begin{align*} \left\|LS\left(\theta^{\,j}\right)-L\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)\right\|_\mathcal{Y}&\unicode{x2A7D} \left\|LS\left(\theta^{\,j}\right)-LS\left(\widetilde{\theta}^{\delta,j}\right)\right\|_\mathcal{Y}+\left\|LS\left(\widetilde{\theta}^{\delta,j}\right)-L\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)\right\|_\mathcal{Y}\\ &\unicode{x2A7D}\left\|L\int_0^1 S^{{\prime}}\left(\widetilde{\theta}^{\delta,j}-\lambda\left(\theta^{\,j}-\widetilde{\theta}^{\delta,j}\right)\right)\,d\lambda\left(\theta^{\,j}-\widetilde{\theta}^{\delta,j}\right) \right\|_\mathcal{Y} + \|L\|_{\mathcal{U}\to\mathcal{Y}}\epsilon_{K\left(j\right)}\\ &\unicode{x2A7D}\|L\|_{\mathcal{U}\to\mathcal{Y}}\sup_{\theta\in B_R\left(\theta^\dagger\right)}\left\|S^{{\prime}}\left(\theta\right)\right\|_{X\to\mathcal{U}}\left\|\theta^{\,j}-\widetilde{\theta}^{\delta,j} \right\|_X + \|L\|_{\mathcal{U}\to\mathcal{Y}}\epsilon_{K\left(j\right)}\\ &\unicode{x2A7D} M_S\|L\|\left\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\right\| + \|L\|\epsilon_{K\left(j\right)}. \end{align*}$$

 □

Remark 2.  M_S in (22) is a locally uniform bound for solutions of the linearized PDE (14). A sufficient condition for (22) to hold is Lipschitz continuity of the parameter-to-state map S. When θ plays the role of a source term or initial data, results regarding Lipschitz continuity of S for general classes of nonlinear PDE model f exist in the literature [44, theorem 8.32].

We now turn our attention to the discrepancy between the exact adjoint and the approximate adjoint with respect to the state approximation error $\epsilon_{K(j)}$.

Lemma 2 (Adjoint state error). Let the assumptions in lemma 1 hold. Suppose that the solution z of the adjoint PDE

$$\begin{align*} \begin{cases} &-\dot{z}+f^{{\prime}}_u\left(\theta,u\right)^\star z = h \\ &z\left(T\right) = 0, \end{cases} \end{align*}$$

depends continuously on the source term

$$\begin{align} \|z\|_\mathcal{U}+\|z\left(0\right)\|_H\unicode{x2A7D} C_{f_u} \|h\|_{\mathcal{U}^*} \qquad \forall \left(\theta,u\right)\in B_R\left(\theta^\dagger\right)\times B_r\left(u^\dagger\right). \end{align} \tag{ 24 }$$

Furthermore, assume that $f^{{\prime}}_\theta,f^{{\prime}}_u$ are Lipschitz continuous, with

$$\begin{align} \begin{aligned} &\|\,f^{{\prime}}_\theta\left(\theta_1,u_1\right)-f^{{\prime}}_\theta\left(\theta_2,u_2\right)\|_{X\to{\mathcal{U}^*}}+\|\,f^{{\prime}}_u\left(\theta_1,u_1\right)-f^{{\prime}}_u\left(\theta_2,u_2\right)\|_{\mathcal{U}\to{\mathcal{U}^*}}\\ & \quad \unicode{x2A7D} L_{\nabla f} \left(\|\theta_1-\theta_2\|_X+\|u_1-u_2\|_\mathcal{U}\right), \end{aligned} \end{align} \tag{ 25 }$$

and denote $L_{R,r}: = L_{\nabla f}\left(R+rC_{\mathcal{V}\to\mathcal{U}}\right)+\|\,f^{{\prime}}_\theta(\theta^\dagger,u^\dagger)\|_{X\to{\mathcal{U}^*}}+\|\,f^{{\prime}}_u(\theta^\dagger,u^\dagger)\|_{\mathcal{U}\to{\mathcal{U}^*}}$.

Then, given any other input $\theta^{\,j}\in B_R(\theta^\dagger)$, we obtain the adjoint error

$$\begin{align} \begin{aligned} \|S^{{\prime}}&\left(\theta^{\,j}\right)^*-S^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)^*\|_{\mathcal{U}\to X}\\ &\unicode{x2A7D} {L_{\nabla f}} C_{f_u}\|I_X\|\|D_U\|\left[L_{R,r} C_{f_u}+1+\|A\|C_{f_u}\right]\left(\left(1+M_S\right)\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|+ \epsilon_{K\left(j\right)}\right)\\ & = :C_{\nabla f,A}\left(\left(1+M_S\right)\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|+ \epsilon_{K\left(j\right)}\right). \end{aligned} \end{align} \tag{ 26 }$$

Proof. The Lipschitz condition (25) implies boundedness of the derivatives, due to

$$\begin{align} &\|\,f^{{\prime}}_\theta\left(\theta,u\right)\|_{X\to{\mathcal{U}^*}}+\|\,f^{{\prime}}_u\left(\theta,u\right)\|_{\mathcal{U}\to{\mathcal{U}^*}}\nonumber\\ & \quad \unicode{x2A7D}\|\,f^{{\prime}}_\theta\left(\theta,u\right)-f^{{\prime}}_\theta\left(\theta^\dagger,u^\dagger\right)\|+\|\,f^{{\prime}}_u\left(\theta,u\right)-f^{{\prime}}_u\left(\theta^\dagger,u^\dagger\right)\| +\|\,f^{{\prime}}_\theta\left(\theta^\dagger,u^\dagger\right)\|+\|\,f^{{\prime}}_u\left(\theta^\dagger,u^\dagger\right)\| \nonumber\\ & \quad \unicode{x2A7D} L_{\nabla f}\left(\|\theta-\theta^\dagger\|_X+C_{\mathcal{V}\to\mathcal{U}}\|u-u^\dagger\|_\mathcal{V}\right)+\|\,f^{{\prime}}_\theta\left(\theta^\dagger,u^\dagger\right)\|+\|\,f^{{\prime}}_u\left(\theta^\dagger,u^\dagger\right)\| \nonumber\\ & \quad \unicode{x2A7D} L_{\nabla f}\left(R+rC_{\mathcal{V}\to\mathcal{U}}\right)+\|\,f^{{\prime}}_\theta\left(\theta^\dagger,u^\dagger\right)\|+\|\,f^{{\prime}}_u\left(\theta^\dagger,u^\dagger\right)\| \nonumber\\ & \quad = L_{R,r}, \end{align} \tag{ 27 }$$

recalling that R and r are the radii of the balls $B_R(\theta^\dagger)\subset X$ and $B_r(u^\dagger)\subset \mathcal{V}$. We will use this bound in evaluating the term Q₂ later in the proof.

Consulting Algorithm 1, subtracting the approximate adjoint state (20) from the exact adjoint state (17), then inserting the mixed term $f^{{\prime}}_\theta(\widetilde{\theta}^{\delta,j},\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j}))^*z$, we see that for any $v\in \mathcal{U}$,

$$\begin{align*} &\left(S^{{\prime}}\left(\theta^{\,j}\right)-S^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)\right)^*v\\ & \quad = -I_X\int_0^T \left(f^{{\prime}}_\theta\left(\theta^{\,j},S\left(\theta^{\,j}\right)\right)\left(t\right)^\star z\left(t\right)- f^{{\prime}}_\theta\left(\widetilde{\theta}^{\delta,j},\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right) \right)\left(t\right)^\star \widetilde{z}\left(t\right)\right)\textrm{d}t+I_XA^*\left(z\left(0\right)-\widetilde{z}\left(0\right)\right)\\ & \quad = -I_X\int_0^T \left(f^{{\prime}}_\theta\left(\theta^{\,j},S\left(\theta^{\,j}\right)\right)\left(t\right)^\star -f^{{\prime}}_\theta\left(\widetilde{\theta}^{\delta,j},\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right) \right)\left(t\right)^\star\right)z(t)\,\textrm{d}t\\ & \qquad- I_X\int_0^T f^{{\prime}}_\theta(\widetilde{\theta}^{\delta,j},\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j}))(t)^\star(z(t)-\widetilde{z}(t))\,\textrm{d}t\\ & \qquad+I_XA^\star(z(0)-\widetilde{z}(0)) = :Q_1+Q_2+Q_3. \end{align*}$$

We recall that z and $\widetilde{z}$ are, respectively, the solutions of

$$\begin{align*} \begin{cases} &-\dot{z}+f^{{\prime}}_u\left(\theta^{\,j},S\left(\theta^{\,j}\right)\right)^\star z = D_Uv \\ &z\left(T\right) = 0, \end{cases} \text{and } \begin{cases} &-\dot{\widetilde{z}}+f^{{\prime}}_u\left(\widetilde{\theta}^{\delta,j},\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right) \right)^\star \widetilde{z} = D_Uv\\ &\widetilde{z}\left(T\right) = 0. \end{cases} \end{align*}$$

Inspecting Q₂, we see that the difference $z-\widetilde{z}$ solves

$$\begin{align} \begin{cases} &-(\dot{z}-\dot{\widetilde{z}})+f^{{\prime}}_u(\widetilde{\theta}^{\delta,j},\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j}) )^\star (z-\widetilde{z}) = \left(-f^{{\prime}}_u(\theta^{\,j},S(\theta^{\,j})+f^{{\prime}}_u(\widetilde{\theta}^{\delta,j},\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j})\right)^\star z\\ &\widetilde{z}(T) = 0. \end{cases} \end{align} \tag{ 28 }$$

Now, using (24), we can bound z by the source term D_U v. In a similar fashion, we bound $z-\widetilde{z}$ by the source term $(-f^{{\prime}}_u(\theta^{\,j},S(\theta^{\,j})+f^{{\prime}}_u(\widetilde{\theta}^{\delta,j},\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j}))^\star z$, then estimate this source further via the Lipschitz condition (25). More precisely, we have

$$\begin{align*} &\|z\|_\mathcal{U}\unicode{x2A7D} C_{f_u}\|D_Uv\|_{\mathcal{U}^*}\unicode{x2A7D} C_{f_u}\|D_U\|_{U\to U^*}\|v\|_\mathcal{U},\\ &\|z-\widetilde{z}\|_\mathcal{U}+\|z(0)-\widetilde{z}(0)\|_H\\ &\qquad\unicode{x2A7D} C_{f_u} \|\,f^{{\prime}}_u(\theta^{\,j},S(\theta^{\,j})-f^{{\prime}}_u(\widetilde{\theta}^{\delta,j},\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j})\|_{\mathcal{U}\to{\mathcal{U}^*}}\|z\|_{\mathcal{U}}\\ &\qquad\unicode{x2A7D} L_{\nabla f} (C_{f_u})^2\|D_U\|_{U\to U^*}\left(\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|_X+\|S(\theta^{\,j})-\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j})\|_\mathcal{U}\right)\|v\|_\mathcal{U}. \end{align*}$$

Using these and the boundedness of the derivative from (27), we can estimate Q₁, Q₂ and Q₃ by the parameter difference $\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|$ and their difference propagating through the state approximation $\|S(\theta^{\,j})-\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j})\|$, i.e.

$$\begin{align*} \|Q_2\|_X&\unicode{x2A7D} \|I_X\|_{X^*\to X}\left\|\,f^{{\prime}}_\theta\left(\widetilde{\theta}^{\delta,j},\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)\right)^\star\right\|_{\mathcal{U}\to X^*}\|z-\widetilde{z}\|_\mathcal{U}\\ &\unicode{x2A7D} \|I_X\|L_{R,r}L_{\nabla f} \left(C_{f_u}\right)^2\|D_U\| \left(\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|_X+\|S\left(\theta^{\,j}\right)-\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)\|_\mathcal{U}\right) \|v\|_\mathcal{U},\\ \|Q_1\|_X&\unicode{x2A7D} \|I_X\|L_{\nabla f} \left(\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|_X+\|S\left(\theta^{\,j}\right)-\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)\|_\mathcal{U}\right) C_{f_u}\|D_U\|_{U\to U^*}\|v\|_\mathcal{U}\\ & = \|I_X\|{L_{\nabla f}} C_{f_u}\|D_U\| \left(\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|_X+\|S\left(\theta^{\,j}\right)-\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)\|_\mathcal{U}\right)\|v\|_\mathcal{U},\\ \|Q_3\|_X&\unicode{x2A7D} \|I_X\|_{X^*\to X}\|A^\star\|_{H\to X^*}\|z\left(0\right)-\widetilde{z}\left(0\right)\|_H\\ &\unicode{x2A7D} \|I_X\|\|A\| L_{\nabla f} \left(C_{f_u}\right)^2\|D_U\|\left(\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|_X+\|S\left(\theta^{\,j}\right)-\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)\|_\mathcal{U}\right)\|v\|_\mathcal{U}. \end{align*}$$

In Q₂, we have applied the estimate $\|\,f^{{\prime}}_\theta(\widetilde{\theta}^{\delta,j},\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j}))^\star\|_{\mathcal{U}\to X^*}\unicode{x2A7D} L_{R,r}$ derived in (27).

In the last step, combining Q₁, Q₂ and Q₃, then applying lemma 1 for $L = \text{Id}$ to further estimate $\|S(\theta^{\,j})-\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j})\|_\mathcal{U}$, we arrive at

$$\begin{align*} &\left\|\left(S^{{\prime}}\left(\theta^{\,j}\right)-S^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)\right)^*\right\|_{\mathcal{U}\to X}\\ & \quad \unicode{x2A7D}\|I_X\| {L_{\nabla f}} C_{f_u}\|D_U\|\left[L_{R,r} C_{f_u}+1+\|A\|C_{f_u}\right]\left(\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|_X+\|S\left(\theta^{\,j}\right)-\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)\|_\mathcal{U}\right)\\ & \quad \unicode{x2A7D} \|I_X\| {L_{\nabla f}} C_{f_u}\|D_U\|\left[L_{R,r} C_{f_u}+1+\|A\|C_{f_u}\right]\left(\left(1+M_S\right)\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|+ \epsilon_{K\left(j\right)}\right), \end{align*}$$

which is (26); this completes the proof. □

At this point, we have analyzed the system output error and adjoint error in the bi-level algorithm 1 via the state approximation error $\epsilon_{K(j)}$. We now study $\epsilon_{K(j)}$ itself.

In this section, we analyze an a posteriori choice for the stopping index $\widetilde{j}^*$ of the upper-level iteration. The analysis also yields a stopping rule K(j) for the lower-level problem.

Proposition 3 (Boundedness—Fejér monotnicity). Let assumptions 1–3 hold, and let $\widetilde{\theta}^{\delta,j}\in B_R(\theta^\dagger)$ be given. Then $\widetilde{\theta}^{\delta,j+1}$ in (49) is a better approximation, i.e.

$$\begin{align*}\|\widetilde{\theta}^{\delta,j+1}-\theta^\dagger\|<\|\widetilde{\theta}^{\delta,j}-\theta^\dagger\|\end{align*}$$

(in particular, the whole sequence remains in the ball $B_R(\theta^\dagger)$) if, for some $\gamma(j)\unicode{x2A7E} 0$, the two following conditions are satisfied:

• 1.  
  The lower-level (19) is iterated until step K(j), such that the state approximation error satisfies

  $$\begin{align} \mathcal{O}\left(\frac{1}{K\left(j\right)^{1/\left(2\alpha\right)}}\right) = \epsilon_{K\left(j\right)}\unicode{x2A7D}\gamma\left(j\right)\|LS\left(\widetilde{\theta}^{\delta,j}\right)-y^\delta\| \end{align} \tag{ 50 }$$

  with some nonnegative and sufficiently small $\gamma(j)$.
• 2.  
  In the upper-level, the residual satisfies the relation

  $$\begin{align} & \|LS\left(\widetilde{\theta}^{\delta,j}\right)-y^\delta\|>\delta \,\Gamma\left(j\right) > 0,\\ \Gamma\left(j\right) & : = 2\dfrac{1+C_{tc}+\|L\|K_R{\gamma\left(j\right)}}{2-\left(1+\epsilon\right)M_R^2-2C_{tc}- \|L\|\left( 2K_R{\gamma\left(j\right)}+M_R^2\|L\|{\gamma\left(j\right)^2} \right)\left(1+\frac{1}{\epsilon}\right)} \nonumber \end{align} \tag{ 51 }$$

  for some ε > 0 with C_tc as in (48) and M_R , µ_R , K_R as in (47).

Proof. Starting from (49), by inserting $LS(\widetilde{\theta}^{\delta,j+1})$ and applying Young's inequality $(a+b)^2\unicode{x2A7D} (1+\epsilon) a^2+(1+1/\epsilon)b^2$, ε > 0, one estimates the error discrepancy between two consecutive iterations as

$$\begin{align} e_{j}^\delta: & = \|\widetilde{\theta}^{\delta,j+1}-\theta^\dagger\|^2-\|\widetilde{\theta}^{\delta,j}-\theta^\dagger\|^2\nonumber\\ & = 2\left( \widetilde{\theta}^{\delta,j+1}-\widetilde{\theta}^{\delta,j},\widetilde{\theta}^{\delta,j}-\theta^\dagger \right)+\|\widetilde{\theta}^{\delta,j+1}-\widetilde{\theta}^{\delta,j}\|^2 \nonumber\\ & = -2\left( L\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)-y^\delta,LS^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)\left(\widetilde{\theta}^{\delta,j}-\theta^\dagger\right) \right) \nonumber\\ & \qquad +\|S^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)^*L^*\left(L\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)-y^\delta\right)\|^2 \nonumber\\ &\unicode{x2A7D}-2\left( LS\left(\widetilde{\theta}^{\delta,j}\right)-y^\delta,LS^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)\left(\widetilde{\theta}^{\delta,j}-\theta^\dagger\right) \right)\nonumber\\ & \quad -2\left( L\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)-LS\left(\widetilde{\theta}^{\delta,j}\right),LS^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)\left(\widetilde{\theta}^{\delta,j}-\theta^\dagger\right) \right) \nonumber\\ & \quad +\left(1+\epsilon\right)\|S^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)^*L^*\left(LS\left(\widetilde{\theta}^{\delta,j}\right)-y^\delta\right)\|^2\nonumber\\ & \quad +\left(1+\frac{1}{\epsilon}\right)\|S^{{\prime}}(\widetilde{\theta}^{\delta,j})^*L^*(L\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j})-LS(\widetilde{\theta}^{\delta,j}))\|^2 \nonumber\\ & = : E_1+E_2+E_3+E_4, \end{align} \tag{ 52 }$$

recalling that $S(\theta)$ is the exact evaluation, while $\widetilde{S}_{K(j)}(\theta)$ is the approximation via the lower-level iteration (19).

For E₁, we insert $y^\delta-LS(\widetilde{\theta}^{\delta,j})$. For E₂, we recall the state approximation error $\epsilon_{K(j)}: = \|S(\widetilde{\theta}^{\delta,j})-\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j})\|$, and use the second part of the tangential cone condition (47). For E₃ and E₄, we employ boundedness of the derivative from (46), as well as the definition of $\epsilon_{K(j)}$ for E₄. Explicitly, one has

$$\begin{align*} E_1& = -2\|LS(\widetilde{\theta}^{\delta,j})-y^\delta\|^2-2\left(LS(\widetilde{\theta}^{\delta,j})-y^\delta, y^\delta-LS(\widetilde{\theta}^{\delta,j})+LS^{{\prime}}(\widetilde{\theta}^{\delta,j})(\widetilde{\theta}^{\delta,j}-\theta^\dagger\right),\\ E_2&\unicode{x2A7D} 2\|L\|\epsilon_{K(j)} K_R\|LS(\widetilde{\theta}^{\delta,j})-LS(\theta^\dagger)\| = 2\|L\|\epsilon_{K(j)} K_R\|LS(\widetilde{\theta}^{\delta,j})-y\|,\\ E_3&\unicode{x2A7D} (1+\epsilon)M_R^2\|LS(\widetilde{\theta}^{\delta,j})-y^\delta\|^2,\\ E_4&\unicode{x2A7D} \left(1+\frac{1}{\epsilon}\right)M_R^2\|L\|^2 \epsilon_{K(j)}^2. \end{align*}$$

Combining these estimates then inserting y and applying tangential cone condition in the equivalent form (48), we get

$$\begin{align*} e_{j}^\delta & \unicode{x2A7D} \big(E_1+E_3\big)+\big(E_2+E_4\big)\\ &\unicode{x2A7D} -\big(2-(1+\epsilon)M_R^2\big)\|LS(\widetilde{\theta}^{\delta,j})-y^\delta\|^2 -2\left( LS(\widetilde{\theta}^{\delta,j})-y^\delta,y^\delta-LS(\widetilde{\theta}^{\delta,j})+LS^{{\prime}}(\widetilde{\theta}^{\delta,j})(\widetilde{\theta}^{\delta,j}-\theta^\dagger) \right)\\ &\hspace{1cm}+\|L\|\epsilon_{K(j)}\left[ 2K_R\|LS(\widetilde{\theta}^{\delta,j})-y\|+\left(1+\frac{1}{\epsilon}\right)M_R^2\|L\| \epsilon_{K(j)}\right]\\ &\unicode{x2A7D} -\big(2-(1+\epsilon)M_R^2\big)\|LS(\widetilde{\theta}^{\delta,j})-y^\delta\|^2+ 2\|LS(\widetilde{\theta}^{\delta,j})-y^\delta\|\left( \delta+C_{tc}(\delta+\|LS(\widetilde{\theta}^{\delta,j})-y^\delta\|) \right)\\ & \hspace{1cm} +\|L\|\epsilon_{K(j)}\left[ 2K_R(\delta+\|LS(\widetilde{\theta}^{\delta,j})-y^\delta\|) +\left(1+\frac{1}{\epsilon}\right)M_R^2\|L\| \epsilon_{K(j)}) \right]\\ & = \|LS(\widetilde{\theta}^{\delta,j})-y^\delta\| \left[ -\Big(2-(1+\epsilon)M_R^2-2C_{tc}\Big)\|LS(\widetilde{\theta}^{\delta,j})-y^\delta\| + 2\delta(1+C_{tc})\right]\\ & \hspace{1cm} +\|L\|\epsilon_{K(j)}\left[ 2K_R\delta+2K_R\|LS(\widetilde{\theta}^{\delta,j})-y^\delta\|) +\left(1+\frac{1}{\epsilon}\right)M_R^2\|L\| \epsilon_{K(j)}) \right]. \end{align*}$$

This estimate reveals that if the state approximation error $\epsilon_{K(j)}$ is sufficiently small in comparison to $\|LS(\widetilde{\theta}^{\delta,j})-y^\delta\|$, we obtain negativity of $e_{j}^\delta$. This motivates us to control the lower-level error $\epsilon_{K(j)}$ by the upper-level residual $\gamma(j)\|LS(\widetilde{\theta}^{\delta,j})-y^\delta\|$ as in (50). Indeed, (50) leads us to

$$\begin{align*} e_{j}^\delta\unicode{x2A7D} &\ \|LS\left(\widetilde{\theta}^{\delta,j}\right)-y^\delta\|\,\\ & \cdot \left[ -\left(2-\left(1+\epsilon\right)M_R^2-2C_{tc}- \|L\|\left( 2K_R{\gamma\left(j\right)}+M_R^2\|L\|{\gamma\left(j\right)^2} \right)\left(1+\frac{1}{\epsilon}\right) \right)\left\|LS\left(\widetilde{\theta}^{\delta,j}\right)-y^\delta\right\|\right.\\ & \left. +\ 2\delta\left(1+C_{tc}+\|L\|K_R\gamma\left(j\right)\right) \right]\quad <0, \end{align*}$$

where negativity is achieved if we stop the upper-level iteration according to the rule (51). Fejér monotonicity follows. □

In view of (51), we suggest the posterior upper-level stopping rule

$$\begin{align} \boxed{\left\|LS\left(\widetilde{\theta}^{\delta,\widetilde{j}^*}\right)-y^\delta\right\|\unicode{x2A7D}\tau \delta< \left\|LS\left(\widetilde{\theta}^{\delta,j}\right)-y^\delta\right\|, \quad \tau > \Gamma\left(j\right), \quad 0\unicode{x2A7D} j<\widetilde{j}^* } \end{align} \tag{ 53 }$$

with $\Gamma(j)$ as in (51). By Fejér monotonicity, $\widetilde{\theta}^{\delta,j+1}$ is a better approximation of $\theta^\dagger$ than $\widetilde{\theta}^{\delta,j}$, as long as the discrepancy (53) is respected.

Remark 8 (Stopping rule—consistency check). We wish to emphasize that when S is solved exactly or an analytic solution of the PDE (1) can be derived, then in (51) one can set $K(j) = \infty$, $\gamma(j) = \epsilon = 0$ and bypass the lower-level problem. In this event, (51) reduces to

$$\begin{align*} \left\|LS\left(\widetilde{\theta}^{\delta,j}\right)-y^\delta\right\|>2\dfrac{1+C_{tc}}{1-2C_{tc}}\,\delta \qquad \text{for }M_R = 1, C_{tc}<\frac{1}{2}, \end{align*}$$

coinciding with the known discrepancy principle in [29, section 2.2, (2.12)].

At this point, it is natural to ask how one may obtain the exact evaluation $S(\widetilde{\theta}^{\delta,j})$ required to compare the residual $\|LS(\widetilde{\theta}^{\delta,j})-y^\delta\|$ with τδ as the discrepancy principle (53) suggests, and how terminate the lower-level iteration w.r.t the upper-level residual (50). A possible scenario is that one may have access to the system output $LS:\widetilde{\theta}^{\delta,j}\mapsto LS(\widetilde{\theta}^{\delta,j})$, despite $S(\widetilde{\theta}^{\delta,j})$ not being directly accessible. In this event, the proposed posterior stopping rule (53) can be verified. If this is not the case, an a priori stopping rule $\widetilde{j}^*$ via the noise level δ is desirable.

Proposition 4 (Uniform boundedness). Let assumptions 1–3 hold. All the iterates $\widetilde{\theta}^{\delta,j}$ with $0\unicode{x2A7D} j\unicode{x2A7D}\widetilde{j}^*(\delta)$ remain in the ball $B_R(\theta^\dagger)$ if, for some $\gamma(j)\unicode{x2A7E} 0$, the two following conditions are satisfied:

• 1.  
  The lower-level problem (19) is iterated until step K(j), such that the state approximation error satisfies

  $$\begin{align} \epsilon_{K\left(j\right)} = \mathcal{O}\left(\frac{1}{K\left(j\right)^{1/\left(2\alpha\right)}}\right)\unicode{x2A7D} \delta\gamma\left(j\right) \end{align} \tag{ 54 }$$
• 2.  
  In the upper-level, the initialization and the prior stopping index satisfy

  $$\begin{align} \begin{aligned} &\|\theta^0-\theta^\dagger\|\unicode{x2A7D} R_0 \lt R,\\ &\widetilde{j}^*\left(\delta\right)\delta^2\left(\frac{K_R^2}{\mu_R}+5M_R^2+ 2\frac{K_R^2}{\mu_R}\|L\|^2\gamma\left(j\right)^2 +\frac{7}{2}M_R^2\|L\|^2\gamma\left(j\right)^2\right)\unicode{x2A7D} R^2-R_0^2. \end{aligned} \end{align} \tag{ 55 }$$

Proof. Similarly to proposition 3, we will evaluate the error $e_{j}^\delta$ between two consecutive iterations. We begin with the estimate (52), where in the last line, we further insert y into E₁ and E₃ (that is, the terms containing y^δ ). So

$$\begin{align*} &e_{j}^\delta = \|\widetilde{\theta}^{\delta,j+1}-\theta^\dagger\|^2-\|\widetilde{\theta}^{\delta,j}-\theta^\dagger\|^2\unicode{x2A7D} A+B. \end{align*}$$

where A has the same form as as $(E_1+E_2+E_3+E_4)$ in (52), with y replacing y^δ in E₁, E₃, and B is the residual generated by inserting y.

In A, we can directly apply the first part of the tangential cone condition (47) to E₁; the rest is estimated in a similar manner as in (52). That is,

$$\begin{align*} A& = -2\left( LS\left(\widetilde{\theta}^{\delta,j}\right)-y,LS^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)\left(\widetilde{\theta}^{\delta,j}-\theta^\dagger\right) \right)\\ & \quad -2\left( L\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)-LS\left(\widetilde{\theta}^{\delta,j}\right),LS^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)\left(\widetilde{\theta}^{\delta,j}-\theta^\dagger\right) \right) \nonumber\\ & \quad +\left(1+\epsilon\right)\left\|S^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)^*L^*\left(LS\left(\widetilde{\theta}^{\delta,j}\right)-y\right)\right\|^2\\ & \quad +\left(1+\frac{1}{\epsilon}\right)\left\|S^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)^*L^*\left(L\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)-LS\left(\widetilde{\theta}^{\delta,j}\right)\right)\right\|^2\\ &\unicode{x2A7D}-2\left(M_R^2+\mu_R\right)\left\|LS\left(\widetilde{\theta}^{\delta,j}\right)-y\right\|^2+2\|L\|\epsilon_{K\left(j\right)} K_R\left\|LS\left(\widetilde{\theta}^{\delta,j}\right)-y\right\|\\ & \quad +(1+\epsilon)(1+\epsilon^{{\prime}})M_R^2\left\|LS(\widetilde{\theta}^{\delta,j})-y\right\|^2+\left(1+\frac{1}{\epsilon}\right) M_R^2 \|L\|^2\epsilon_{K(j)}^2\\ & = : A_1+A_2+A_3+A_4. \end{align*}$$

Noticing that in comparison to E₃, A₃ has the extra constant $(1+\epsilon^{{\prime}})$ generated from splitting the square term $(a+b)^2\unicode{x2A7D} (1+\epsilon^{{\prime}})a^2+(1+\frac{1}{\epsilon^{{\prime}}})b^2$, $\epsilon^{{\prime}}\gt0$ when inserting y.

In B, we use the second part of the tangential cone condition (47), boundedness of the derivative as in (46), then finally Young's inequality $2ab\unicode{x2A7D} \mu_R a^2+b^2/\mu_R$ for $a: = \|LS(\widetilde{\theta}^{\delta,j})-y\|$, $b: = \delta K_R$. This leads to

$$\begin{align*} B& = -2\left( y-y^\delta,LS^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)\left(\widetilde{\theta}^{\delta,j}-\theta^\dagger\right) \right)+ \left(1+\epsilon\right)\left(1+\frac{1}{\epsilon^{{\prime}}}\right)\left\|S^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)^*L^*\left(y-y^\delta\right)\right\|^2\\ &\unicode{x2A7D} 2\delta K_R\left\|LS\left(\widetilde{\theta}^{\delta,j}\right)-y\right\|+\left(1+\epsilon\right)\left(1+\frac{1}{\epsilon^{{\prime}}}\right)\delta^2 M_R^2\\ &\unicode{x2A7D} \mu_R\left\|LS\left(\widetilde{\theta}^{\delta,j}\right)-y\right\|^2+\delta^2\left(\frac{K_R^2}{\mu_R}+\left(1+\epsilon\right)\left(1+\frac{1}{\epsilon^{{\prime}}}\right)M_R^2\right)\\ & = :B_1+B_2 \end{align*}$$

We now set $\epsilon = \epsilon^{{\prime}} = \sqrt{2}-1$ in both A and B. In A₃, we thus have the constant $(1+\epsilon)(1+\epsilon^{{\prime}}) = 2$; in A₄, we have $(1+\frac{1}{\epsilon}) = \sqrt{2}/(\sqrt{2}-1)\lt7/2$, while in B₂, we have $(1+\epsilon)(1+\frac{1}{\epsilon^{{\prime}}})\lt5$.

Next, we bound $\epsilon_{K(j)}$ in A₂, A₄ by $\delta \gamma(j)$ as the rule (54) suggests. In A₂, we further apply Young's inequality $2ab\unicode{x2A7D} \mu_Ra^2/2 +2b^2/\mu_R$ for $a: = \|LS(\widetilde{\theta}^{\delta,j})-y\|$, $b: = \|L\|\epsilon_{K(j)} K_R$.

Now summing A and B, we see that $(A_3+B_1+\mu_Ra^2/2)$ is absorbed by part of A₁, that is

$$\begin{align} e_{j}^\delta&\unicode{x2A7D} A +B \unicode{x2A7D} \left(A_1+A_3+B_1+\mu_Ra^2/2\right) + \left(B_2+2b^2/\mu_R+A_4\right)\nonumber\\ &\unicode{x2A7D} -\frac{\mu_R}{2}\left\|LS\left(\widetilde{\theta}^{\delta,j}\right)-y\right\|^2+\delta^2\left(\frac{K_R^2}{\mu_R}+5M_R^2+ 2\frac{K_R^2}{\mu_R}\|L\|^2\gamma\left(j\right)^2 +\frac{7}{2}M_R^2\|L\|^2\gamma\left(j\right)^2\right). \end{align} \tag{ 56 }$$

Next, summing up $e_{j}^\delta$ for $j = 0\ldots\bar{j}$ for any $\bar{j}\unicode{x2A7D}\widetilde{j}^*(\delta)$ with the upper-stopping index $\widetilde{j}^*$ proposed in (55), we arrive at

$$\begin{align} \frac{\mu_R}{2}&\sum_{j = 0}^{\bar{j}}\left\|LS\left(\widetilde{\theta}^{\delta,j}\right)-y\right\|^2+\|\widetilde{\theta}^{\delta,\bar{j}}-\theta^\dagger\|^2 \nonumber\\ &\unicode{x2A7D} \|\theta^0-\theta^\dagger\|^2+\widetilde{j}^*\left(\delta\right)\delta^2\left(\frac{K_R^2}{\mu_R}+5M_R^2+ 2\frac{K_R^2}{\mu_R}\|L\|^2\gamma\left(j\right)^2 +\frac{7}{2}M_R^2\|L\|^2\gamma\left(j\right)^2\right)\nonumber\\ &\unicode{x2A7D} R^2. \end{align} \tag{ 57 }$$

This shows that all the iterates until $\widetilde{j}^*(\delta)$ remain in the ball

$$\begin{align*}\widetilde{\theta}^{\delta,j}\in B_R\left(\theta^\dagger\right) \quad \forall j\unicode{x2A7D}\widetilde{j}^*\left(\delta\right),\end{align*}$$

and the proof is complete. □

In comparison to proposition 3, a prior stopping rule cannot yield Fejér monotonicity of the iterate with noisy data. As monotonicity is essential in ensuring the regularizing property of the Landweber iteration, and also due to the bi-level factor, we need to take particular care regarding stability of the upper-level.

The following Proposition is crucial to derive a more concrete stopping rule, which can ensure controllability of the approximation errors propagating to the upper-level. At this point, the preliminary error analysis in section 2.2 comes into play. We remind ourselves that the circumflex $\widetilde{(\cdot)}$ denotes the approximation/bi-level scheme; notation without the circumflex refers to the exact/single-level scheme.

Proposition 5 (Noise propagation-approximation error). Let the assumption in proposition 4 hold. Assume moreover that the lower-level precision (54) satisfies

$$\begin{align} \epsilon_{K\left(j\right)}\unicode{x2A7D}\delta\gamma\left(j\right)\unicode{x2A7D} \frac{\delta}{q^{\,j}} \qquad \text{for some }q\unicode{x2A7E}1. \end{align} \tag{ 58 }$$

Then, we achieve the estimate

$$\begin{align} &\|\widetilde{\theta}^{\delta,j}-\theta^{\,j}\| \unicode{x2A7D} \delta\,\hat\Gamma\left(j\right),\\ & \hat{\Gamma}\left(j\right): = M_R\frac{\left[1+M_R^2+\left(1+M_S\right)\mathbb{D} \right]^j-1}{M_R^2+\left(1+M_S\right)\mathbb{D}} +\frac{\left(M_R\|L\|+ \mathbb{D}\right)}{q} \frac{\left[1+M_R^2+\left(1+M_S\right)\mathbb{D} \right]^j-\left(1/q\right)^j}{\left[1+M_R^2+\left(1+M_S\right)\mathbb{D} \right]-\left(1/q\right) } \nonumber \end{align} \tag{ 59 }$$

for all $j\unicode{x2A7D}\widetilde{j}^*(\delta)$, where

$$\begin{align*} \delta\unicode{x2A7D}\bar{\delta}, \qquad\gamma\left(j\right)\unicode{x2A7D} \bar{\gamma},\qquad \mathbb{D}: = C_{\nabla f,A}\|L\|\left(\bar{\delta}+M_R R+\|L\|\bar{\delta}\bar{\gamma}\right),\qquad C_{\nabla f,A} \text{ as in } \left(26\right). \end{align*}$$

Proof. We begin with comparing the $(j+1)$-iteration of the single-level scheme with exact data, and of the bi-level scheme with noisy data

$$\begin{align*} &\theta^{j+1} = \theta^{\,j} -S^{{\prime}}(\theta^{\,j})^*L^*\left(LS(\theta^{j}-y\right),\\ &\widetilde{\theta}^{\delta,j+1} = \widetilde{\theta}^{\delta,j} -S^{{\prime}}(\widetilde{\theta}^{\delta,j})^*L^*\left(L\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j})-y^\delta\right). \end{align*}$$

Inserting $S^{{\prime}}(\theta^{\,j})^*L^*(L\widetilde{S}_{K(j)}(\widetilde{\theta}^{\delta,j})-y^\delta)$ (first estimate below) and using the analysis of the state approximation error derived in lemma 1 (second and third estimates) and of the adjoint approximation error derived lemma 2 (second estimate), we deduce

$$\begin{align*} \|\theta^{j+1}-\widetilde{\theta}^{\delta,j+1}\| & \unicode{x2A7D}\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|\\ & \quad+\left\|S^{{\prime}}\left(\theta^{\,j}\right)^*L^*\left(LS\left(\theta^{\,j}\right)-L\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)\right)\right\|+\|S^{{\prime}}\left(\theta^{\,j}\right)^*L^*\|\delta\\ &\quad+\left\|\left(S^{{\prime}}\left(\theta^{\,j}\right)^*-S^{{\prime}}\left(\widetilde{\theta}^{\delta,j}\right)^*\right)L^*\left(L\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)-y^\delta\right)\right\|\\ &\unicode{x2A7D}\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\| + M_R\left(M_R\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|+\|L\|\epsilon_{K\left(j\right)}\right)+M_R\delta\\ & \quad+ C_{\nabla f,A}\left(\left(1+M_S\right)\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|+ \epsilon_{K\left(j\right)}\right)\|L\|\left(\delta+\|L\widetilde{S}_{K\left(j\right)}\left(\widetilde{\theta}^{\delta,j}\right)-LS\left(\theta^\dagger\right)\|\right)\\ & \unicode{x2A7D}\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\| + M_R\left(M_R\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|+\|L\|\epsilon_{K\left(j\right)}\right)+M_R\delta\\ & \quad+ C_{\nabla f,A}\left(\left(1+M_S\right)\|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|+ \epsilon_{K(j)}\right)\|L\|\left(\delta+M_R\|\widetilde{\theta}^{\delta,j}-\theta^\dagger\|+\|L\|\epsilon_{K(j)}\right). \end{align*}$$

Simplify the notation by denoting

$$\begin{align*} &T^{\,j}: = \|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\|,\\ &\|L\|\left(\delta+M_R\|\widetilde{\theta}^{\delta,j}-\theta^\dagger\|+\|L\|\epsilon_{K\left(j\right)}\right)\unicode{x2A7D}\|L\|\left(\delta+M_R R+\|L\|\epsilon_{K\left(j\right)}\right) = :B_{\delta,\epsilon_{K\left(j\right)}}, \end{align*}$$

where in $B_{\delta,\epsilon_{K(j)}}$, we have invoked uniform boundedness of the iterates in the ball $B_R(\theta^\dagger)$ proved in proposition 4. We now group the estimate in terms of $T^{\,j}, \delta$ and $\epsilon_{K(j)}$, and obtain

$$\begin{align} &T^{j+1}\\ &\quad\unicode{x2A7D} T^{\,j} \left[1+M_R^2+C_{\nabla f,A}\left(1+M_S\right)B_{\delta,\epsilon_{K\left(j\right)}} \right] + \delta M_R +\epsilon_{K\left(j\right)}\left[ M_R\|L\|+ C_{\nabla f,A}B_{\delta,\epsilon_{K\left(j\right)}} \right].\nonumber \end{align} \tag{ 60 }$$

The estimate (60) exposes an important fact: if $\epsilon_{K(j)}$ decays sufficiently fast with respect to δ then we can bound the term $T^{j+1}$ from above. Indeed, (54) suggests running the lower-level iterate until $\epsilon_{K(j)}\unicode{x2A7D}\delta\gamma(j)$. Employing the upper bounds $\bar{\delta}$ and $\bar{\gamma}$ from the statement of the proposition, we have

$$\begin{align*}C_{\nabla f,A}B_{\delta,\epsilon_{K\left(j\right)}}\unicode{x2A7D} C_{\nabla f,A}\|L\|\left(\bar{\delta}+M_R R+\|L\|\bar{\delta}\bar{\gamma}\right) = :\mathbb{D}.\end{align*}$$

This results in the recursive relation

$$\begin{align} T^{j+1}\unicode{x2A7D} T^{\,j} \left[1+M_R^2+\left(1+M_S\right)\mathbb{D} \right] + \delta M_R + \delta\gamma\left(j\right) \left[ M_R\|L\|+ \mathbb{D}\right]. \end{align} \tag{ 61 }$$

By induction and choosing the same initialization $\widetilde{\theta}^{\delta,0} = \theta^{0}$, i.e. T⁰ = 0, we arrive at the important estimate

$$\begin{align} &T^{\,j}\unicode{x2A7D} \delta\left( M_R\frac{\left[1+M_R^2+\left(1+M_S\right)\mathbb{D} \right]^j-1}{M_R^2+\left(1+M_S\right)\mathbb{D}}\right.\\ & \left. \qquad\qquad\qquad+\left( M_R\|L\|+ \mathbb{D}\right) \sum_{i = 1}^j \left[1+M_R^2+\left(1+M_S\right)\mathbb{D} \right]^i\gamma\left(j-i\right)\right) \qquad\text{for all }j\unicode{x2A7D}\widetilde{j}^*\left(\delta\right).\nonumber \end{align} \tag{ 62 }$$

By (58), $\gamma(j)\unicode{x2A7D}\frac{1}{q^{\,j}}$; we can thus collapse the sum in (62), and arrive at the key estimate

$$\begin{align*} T^{\,j} = \|\theta^{\,j}-\widetilde{\theta}^{\delta,j}\| \unicode{x2A7D} \delta\hat{\Gamma}\left(j\right) \end{align*}$$

for all $j\unicode{x2A7D}\widetilde{j}^*(\delta)$, as required. □

This immediately suggests a stopping rule:

Assumption 4 (prior stopping rule). The upper-stopping index $\widetilde{j}^*(\delta)$ satisfies

$$\begin{align} \begin{aligned} & \delta \,\hat{\Gamma}\left(\widetilde{j}^*\left(\delta\right)\right) \to 0 \,\,\text{as}\,\, \delta\to 0, \end{aligned} \end{align} \tag{ 63 }$$

with $\hat{\Gamma}$ as in (59).

With all the auxiliary results prepared, we now state the main result.

Theorem 2 (Convergence of bi-level iteration). Assuming that the following both hold:

• (a)  
  The upper-level stopping index $\widetilde{j}^*(\delta)$ is chosen according to (55) and (63).
• (b)  
  The lower-level stopping index K(j) is chosen according to (58) at each $j\unicode{x2A7D}\widetilde{j}^*(\delta)$.

Then, the bi-level iterates $\widetilde{\theta}^{\delta,\widetilde{j}^*(\delta)}$ in algorithm 1 converge to a solution of the inverse problem (1) and (2) as δ → 0 .

Proof. We begin by claiming that for exact data y, the whole sequence θ^j run with the exact/single-level scheme converges in the weak sense to a solution of the inverse problem (1) and (2). Indeed, looking at the proof of proposition 4, in particular the expression (56), by setting δ = 0, $\gamma(j) = 0$ we have

$$\begin{align*}e_j^{\delta = 0} = \|\theta^{j+1}-\theta^\dagger\|^2-\|\theta^{\,j}-\theta^\dagger\|^2<0 \qquad j\in\mathbb{N}\end{align*}$$

showing Fejér monotonicity. (57) with $\bar{j} = \infty$, $G = L\circ S$ becomes

$$\begin{align*} \frac{\mu_R}{2}&\sum_{j = 0}^\infty\|G\left(\theta^{\,j}\right)-y\|^2\unicode{x2A7D} R^2. \end{align*}$$

These two properties are exactly the same as (35) and (36) in theorem 1, the key steps to prove convergence of the iteration, with θ in place of u, and G in place of ${_{\theta}{\mathbf{F}}}$. Therefore, following the line of the proof of theorem 1, part (i), noting that the tangential cone condition (even the strong version) is also assumed here, we can assert that the whole iterate sequence of the single-level algorithm with noise-free data converges to a solution to the inverse problem problem, i.e.

$$\begin{align*} \theta^{\,j}\overset{j\to\infty}{\rightharpoonup} \theta^\dagger \qquad\text{with}\qquad G\left(\theta^\dagger\right) = y. \end{align*}$$

Moreover, one can claim strong convergence. This is done in [29, theorem 2.4] (which we do not repeat here) by showing that $e_j^{\delta = 0}$ in fact forms a Cauchy sequence; thus,

$$\begin{align} \theta^{\,j}\overset{j\to\infty}{\to} \theta^\dagger \qquad\text{with}\qquad G\left(\theta^\dagger\right) = y. \end{align} \tag{ 64 }$$

This, together with proposition 5 and the stopping rule (63), allows us to conclude convergence of the sequence $\widetilde{\theta}^{\delta_n,\widetilde{j}^*(\delta_n)}$ via the decomposition

$$\begin{align*} \widetilde{\theta}^{\delta_n,\widetilde{j}^*\left(\delta_n\right)}-\theta^\dagger = \left(\widetilde{\theta}^{\delta_n,\widetilde{j}^*\left(\delta_n\right)} - \theta^{\widetilde{j}^*\left(\delta_n\right)}\right) + \left(\theta^{\widetilde{j}^*\left(\delta_n\right)}-\theta^\dagger\right) \, \overset{n\to\infty}{\to} \, 0, \end{align*}$$

where convergence of the first term follows from proposition 5 with the rule (63), and convergence of the second term is (64). The proof ends here. □

As presented in remark 4 of section 3, coercivity (34) is the essential property that we extract from well-posedness of S to achieve the convergence rate for the lower-level iterations. In what follows, we present two classes of PDEs where coercivity can be established: Lipschitz nonlineartiy and monotone-type nonlinearity (66).

Let us consider the PDE (29), where the unknowns $\theta: = (a,c,\phi,u_0)\in X$ consist of the diffusion coefficient a, potential c, source term φ and initial state u₀. These unknown parameters are embedded in the PDE model via

$$\begin{align} {_{\theta}{\mathbf{F}}}\left(u\right):\mathcal{V}\mapsto\mathcal{U}^*\times H\qquad {_{\theta}{\mathbf{F}}}\left(u\right) = \begin{pmatrix} \dot{u}-\nabla\cdot\left(a\nabla u\right)+ c u + \Phi\left(u\right)- \phi\\ u_0-u|_{t = 0} \end{pmatrix} \end{align} \tag{ 65 }$$

in the Hilbert space framework (cf section 1.2)

$$\begin{align} \begin{aligned} &U = H^1_0(\Omega), \quad\qquad H = L^2(\Omega), \qquad U^* = H^{-1}(\Omega),\nonumber\\ &\mathcal{V} = L^2(0,T;U)\cap H^1(0,T;U^*), \qquad\mathcal{U}^* = L^2(0,T;H^{-1}(\Omega),\\ &X = (X_a,X_c,X_\phi,X_0),\\ &X_a = \left\{a\in H^2(\Omega):a\unicode{x2A7E}\underline{a}>0\right\},\qquad X_c = X_0 = L^2(\Omega),\qquad X_\phi = H^{-1}(\Omega) \end{aligned} \end{align}$$

with bounded Lipschitz domain $\Omega\subset\mathbb{R}^3$.

In (65), the nonlinearity $\Phi:\mathcal{V}\to\mathcal{U}^*$ is a Nemytskii operator induced by $\Phi:U\to U^*$ with (by abuse of notation) $[\Phi(u)](t): = \Phi(u(t))$. We consider

$$\begin{align} \Phi \begin{cases} \text{either (a) Lipschitz continuous:} &\|\Phi\left(u\right)-\Phi\left(v\right)\|_{U^*}\unicode{x2A7D} L_\Phi\|u-v\|_H,\\ \text{or (b) of monotone-type:} &\langle \Phi\left(u\right)-\Phi\left(v\right),u-v\rangle_{U^*,U} +M_\Phi\|u-v\|_H^2\unicode{x2A7E} 0 \end{cases} \end{align} \tag{ 66 }$$

for all u, $v\in U$ with $L_\Phi$, $M_\Phi\unicode{x2A7E}0$. Clearly, case (b) covers the class of standard monotonic functions.

Proposition 6 (Coercivity verification). The classes of PDEs with Lipschitz or monotone-type nonlinearities (65) and (66) satisfy the coercivity property listed in Assumption 2 with α = 1 and in $\mathcal{U}$-norm. That is, the relation

$$\begin{align} C_{coe}\|{_{\theta}{\mathbf{F}}}\left(u\right)-{_{\theta}{\mathbf{F}}}\left(v\right)\|_{\mathcal{U}^*\times H}\unicode{x2A7E} \|u-v\|_{\mathcal{U}} \qquad \forall u,\,v\in\mathcal{V} \end{align} \tag{ 67 }$$

holds uniformly for all $\theta: = (a,c,\phi,u_0)\in B_R(\theta^\dagger)$ with a positive constant C_coe .

The property (67) is also referred to as Lipschitz stability.

Proof. We begin by testing ${_{\theta}{\mathbf{F}}}(u)-{_{\theta}{\mathbf{F}}}(v)$ with $(u-v; 0)\in\mathcal{V}\times H$. Setting $\|u\|_{H_0^1}: = \|\nabla u\|_{L^2}$ due to the Gagliardo-Nirenberg-Sobolev inequality and recalling $a\unicode{x2A7E} \underline{a}\gt0$, we write

$$\begin{align} &\langle {_{\theta}{\mathbf{F}}}\left(u\right)-{_{\theta}{\mathbf{F}}}\left(v\right),\left(u-v; 0\right) \rangle_{U^*\times H,U\times H} \nonumber\\ & \quad = \langle \dot{u}-\dot{v}-\nabla\cdot\left(a \nabla\left(u-v\right)\right)+ c \left(u-v\right) + \Phi\left(u\right)-\Phi\left(v\right),u-v\rangle_{U^*,U}\\ & \quad \unicode{x2A7E}\frac{1}{2}\frac{d}{\textrm{d}t}\|u-v\|^2_{L^2} + \underline{a}\|u-v\|_{H^1}^2 + \langle c\left(u-v\right),u-v\rangle_{U^*,U}+\langle \Phi\left(u\right)-\Phi\left(v\right),u-v\rangle_{U^*,U}\nonumber \end{align} \tag{ 68 }$$

Estimating the c-term by using Hölder's inequaltiy with $p = q = 2$, then by Young's inequality $xy\unicode{x2A7D} x^p/p+y^q/q, p = 4/3,q = 4$, we get

$$\begin{align*} &|\langle c\left(u-v\right),u-v\rangle_{U^*,\,U}|\\ &\qquad\unicode{x2A7D} \|c\|_{L^2}\|\left(u-v\right)^{3/2}\left(u-v\right)^{1/2}\|_{L^2}\unicode{x2A7D} \left(C_{H^1\to L^6}^\Omega\right)^{3/2}\|c\|_{L^2}\|u-v\|^{3/2}_{H^1}\|u-v\|^{1/2}_{L^2} \\ &\qquad\unicode{x2A7D} \frac{3}{4}\left(\underline{a}^{3/4}\|u-v\|^{3/2}_{H^1}\right)^{4/3}+\frac{1}{4} \left(\frac{\left(C_{H^1\to L^6}^\Omega\right)^{3/2}}{\underline{a}^{3/4} }\|c\|_{L^2}\|u-v\|^{1/2}\right)^4\\ &\qquad = \frac{3\underline{a}}{4}\|u-v\|^2_{H^1}+\frac{\left(C_{H^1\to L^6}^\Omega\right)^6}{4\underline{a}^3}\|c\|^4_{L^2}\|u-v\|^2_{L^2}, \end{align*}$$

where $C_{H^1\to L^6}^\Omega$ indicates the application of the continuous embedding $H^1(\Omega)\hookrightarrow L^6(\Omega), \Omega\subset\mathbb{R}^3$. For the nonlinear term Φ in (66), one observes

$$\begin{align*} \langle\Phi\left(u\right)-\Phi\left(v\right),u-v\rangle \begin{cases} \unicode{x2A7D} \frac{\underline{a}}{8}\|u-v\|^2_{H^1}+\frac{8\left(L_\Phi\right)^2}{\underline{a}}\|u-v\|^2_{L^2} \qquad &\text{case (a)}\\ \unicode{x2A7E} -M_\Phi\|u-v\|_{L^2}^2 &\text{case (b)} \end{cases} \end{align*}$$

where for the case (a), we have used Young's inequality $xy\unicode{x2A7D} x^2/\epsilon+\epsilon y^2, \epsilon = \underline{a}/8$. Combining the above estimates of the nonlinearity Φ, the c-term to (68), we have that, for all $t\in[0,T]$

$$\begin{align} &\frac{1}{2}\frac{\textrm{d}}{\textrm{d}t}\|u\left(t\right)-v\left(t\right)\|^2_{L^2} + \frac{\underline{a}}{8}\|u-v\|_{H^1}^2 \unicode{x2A7D} \langle {_{\theta}{\mathbf{F}}}\left(u\right)-{_{\theta}{\mathbf{F}}}\left(v\right),\left(u-v; 0\right)\rangle_{U^*\times H,U\times H} +C_{\Phi,R}\|u-v\|^2_{L^2} \nonumber\\ &\frac{\textrm{d}}{\textrm{d}t}\left(\frac{1}{2}\|u-v\|^2_{L^2}\right) + \left(\frac{\underline{a}}{8}-\epsilon\right)\|u\left(s\right)-v\left(s\right)\|_{H^1}^2\unicode{x2A7D} \frac{1}{\epsilon}\|{_{\theta}{\mathbf{F}}}\left(u\right)-{_{\theta}{\mathbf{F}}}\left(v\right)\|_{H^{-1}}^2+C_{\Phi,R}\|u-v\|^2_{L^2}, \end{align} \tag{ 69 }$$

where in the second line, we have applied one more time Young's inequality $ab\unicode{x2A7D} \epsilon a^2+b^2/\epsilon, \epsilon\lt\underline{a}/8$ and grouped the L²-terms with the nonnegative constant

$$\begin{align*}C_{\Phi,R}: = \left(C_{H^1\to L^6}^\Omega\right)^6\|c\|^4_{L^2}/\left(4\underline{a}^3\right)+8\left(L_\Phi\right)^2/\underline{a}+M_\Phi.\end{align*}$$

In the second step, taking the integral of (69) over the time interval $(0,t)$, $t\unicode{x2A7D} T$ leads to

$$\begin{align} &\left(\frac{\underline{a}}{8}-\epsilon\right)\|u-v\|_{L^2\left(0,t;H^1\left(\Omega\right)\right)}^2+\frac{1}{2}\|u\left(t\right)-v\left(t\right)\|^2_{L^2} \nonumber\\ &\qquad\unicode{x2A7D} \frac{1}{\epsilon}\|{_{\theta}{\mathbf{F}}}\left(u\right)-{_{\theta}{\mathbf{F}}}\left(v\right)\|_{L^2\left(0,t;H^{-1}\left(\Omega\right)\right)}^2+\frac{1}{2}\|\left(u-v\right)|_{t = 0}\|^2_{L^2}+C_{\Phi,R}\int_0^t\|u\left(s\right)-v\left(s\right)\|^2_{L^2}\,\textrm{d}s \nonumber\\ &\qquad\unicode{x2A7D} \left(\frac{1}{\epsilon}+\frac{1}{2}\right)\|{_{\theta}{\mathbf{F}}}\left(u\right)-{_{\theta}{\mathbf{F}}}\left(v\right)\|_{\mathcal{U}^*\times H}^2+C_{\Phi,R}\int_0^t\|u\left(s\right)-v\left(s\right)\|^2_{L^2}\,\textrm{d}s. \end{align} \tag{ 70 }$$

Now, we apply Grönwall's inequality in the form

$$\begin{align*}x\left(t\right)\unicode{x2A7D} \Pi+C_{\Phi,R}\int_0^t x\left(s\right)\,\textrm{d}s\quad\Rightarrow\quad x\left(t\right)\unicode{x2A7D} \Pi + C_{\Phi,R}\int_0^t\Pi \text{exp}^{C_{\Phi,R}\left(t-s\right)}\,\textrm{d}s \quad\forall t\in\left[0,T\right]\end{align*}$$

for (70), where $x(t): = \frac{1}{2}\|u(t)-v(t)\|^2_{L^2}$ is the second term on the left hand side and $\Pi: = \left(\frac{1}{\epsilon}+\frac{1}{2}\right)\|{_{\theta}{\mathbf{F}}}(u)-{_{\theta}{\mathbf{F}}}(v)\|_{\mathcal{U}^*\times H}^2$. We then get

$$\begin{align*} \|u-v\|_{C\left(0,T;L^2\left(\Omega\right)\right)}^2\unicode{x2A7D} \left(\frac{2}{\epsilon}+1+2C_{\Phi,R}T\text{exp}^{2C_{\Phi,R}T}\right)\|{_{\theta}{\mathbf{F}}}\left(u\right)-{_{\theta}{\mathbf{F}}}\left(v\right)\|_{\mathcal{U}^*\times H}^2. \end{align*}$$

In (70), one can moreover bound the term $(\underline{a}/8-\epsilon)\|u-v\|_{L^2(0,T;H^1(\Omega))}$ on the left hand side by the entire right hand side, thus eventually by $\|{_{\theta}{\mathbf{F}}}(u)-{_{\theta}{\mathbf{F}}}(v)\|_{\mathcal{U}*\times H}$. Altogether, we obtain

$$\begin{align} \|u-v\|_{C\left(L^2\right)\cap\, \mathcal{U}}\unicode{x2A7D} \widetilde{C}_{\Phi,R}\|{_{\theta}{\mathbf{F}}}\left(u\right)-{_{\theta}{\mathbf{F}}}\left(v\right)\|_{\mathcal{U}^*\times H} \end{align} \tag{ 71 }$$

with some $\widetilde{C}_{\Phi,R}\gt0$. Importantly, as the parameter θ lies in the ball $B_R(\theta^\dagger)$ with finite radius R, the estimate (71) is locally uniform. If only $(\varphi,u_0)$ are the unknown parameters, the radius R is allowed to be arbitrarily large, as ${C}_{\Phi,R}$, thus $\widetilde{C}_{\Phi,R}$, depends only on (c, a). This shows Lipschitz stability of $u\mapsto{_{\theta}{\mathbf{F}}}(u)$, equivalently the desired coercivity. □

Remark 9 (coercivity with $\mathcal{V}$-norm). By evaluating

$$\begin{align*}\dot{u}-\dot{v} = \Delta \left(u-v\right)- c \left(u-v\right) - \Phi\left(u\right)+\Phi\left(v\right)\end{align*}$$

in $L^2(0,T;H^{-1}(\Omega))$-norm, one can from the preceding result immediately deduce the stronger estimate

$$\begin{align*} \|u-v\|_{\mathcal{V}}\unicode{x2A7D} \widetilde{C}_{\Phi,R}\|{_{\theta}{\mathbf{F}}}\left(u\right)-{_{\theta}{\mathbf{F}}}\left(v\right)\|_{\mathcal{U}^*\times H} \end{align*}$$

for Lipschitz nonlinearity Φ, case (a) in (66). However, as pointed out in remark (3)(iii), a $\mathcal{U}$-norm coercivity is sufficient for the bi-level algorithm.

Remark 10 (Tangential cone condition). In view of corollary 3, the tangential cone condition for the lower-level problem (in u) is clearly satisfied for $\Phi(u) = 0$. When $\Phi(u)$ is as in (66), as coercivity/Lipschitz stability property is already establish in proposition 6, one only needs to assume Lipschitz continuity of its derivative to conclude the tangential cone condition. Thus, convergence of the lower-level iteration is highly achievable.

For the upper-level problem, verification of the tangential condition (in θ) is another technical question. We refer to [30] for a study of such condition in several linear and nonlinear parabolic benchmark problems, to [21] for abstract linear elliptic problems, to [23] for quantitative elastography, to [24] for magnetic resonance elastography, to [36] for the electrical impedance tomography, to [11] for full wave from inversion, and recently, to [1, 46] in the context of machine learning.

This section is dedicated to discussion regarding a real-life application: magnetic particle imaging (MPI). The material of the following discussion is extracted from the two papers [31, 42]. The work in [31] puts forward the mathematical analysis of the parameter identification problem in MPI; subsequently, [42] proceeds with Landweber, Landweber-Kaczmarz algorithms and numerical experiments. Moreover, an ad hoc bi-level algorithm, in which lower-level iteration is used as an approximation method for the nonlinear PDE (73), was successfully implemented without full analysis.

5.2.1. State of the art.

5.2.2. Parameter identification in MPI.

The external magnetic field induces the temporal change $\partial\mathbf{m}/\partial t$ of the particles, which by Faraday's law induces a measurable electric voltage

$$\begin{align} v\left(t\right) = \int_0^T \int_{\Omega} -\mu_0 \widetilde{a}_l\left(t-\tau\right) c\left(x\right) \, \mathbf{p}^{\mathrm{R}}\left(x\right) \cdot \frac{\partial}{\partial \tau} \mathbf{m}\left(x,\tau\right) \, \mathrm{d} x \, \mathrm{d} = : \mathcal{K}\mathbf{m}. \end{align} \tag{ 72 }$$

Here, µ₀ is magnetic permeability in vacuum, $\widetilde{a}_l$ is a band-pass filter, c is the known particle concentration and $\mathbf{p}^{\mathrm{R}}$ is the 3D receive coil sensitivity.

Inspired by [15, 38], we use the Landau–Lifshitz–Gilbert (LLG) equation to model the evolution of the magnetization moment vector m in a microscopic scale

$$\begin{align} \alpha_1 m_{\mathrm{S}}^2 \frac{\textrm{d}}{\textrm{d}t}\mathbf{m} - \alpha_2 \mathbf{m} \times \frac{\textrm{d}}{\textrm{d}t} \mathbf{m}-m_{\mathrm{S}}^2\Delta \mathbf{m} & = \lvert \nabla \mathbf{m} \rvert^2 \mathbf{m} + m_{\mathrm{S}}^2\mathbf{h} - \langle \mathbf{m}, \mathbf{h} \rangle \mathbf{m} \ && \text{in} \ \Omega \times \left[0,T\right], \nonumber\\ 0 & = \partial_{\nu} \mathbf{m} && \text{on} \ \partial\Omega \times \left[0,T\right],\\ \mathbf{m}_0 & = \mathbf{m}\left(t = 0\right), \ \lvert \mathbf{m}_0 \rvert = m_{\mathrm{S}} && \text{in} \ \Omega\nonumber \end{align} \tag{ 73 }$$

with the relaxation parameter $\alpha\in\mathbb{R}^2$, external magnetic field h and the saturation magnetization constant m_S . We highlight that as (73) includes the relaxation effect (via α and the cross product) and the particle interaction (via $\Delta \mathbf{m}$ and $|\nabla \mathbf{m} \rvert^2 \mathbf{m}$), our approach is realistic to a greater extent in comparison to simplified models, e.g. the Langevin function.

Combining (72) and (73), we investigate the inverse problem of identifying α along side the magnetization m

$$\begin{align} \text{Find } \alpha = \left(\alpha_1,\alpha_2\right):\quad \mathcal{K}\mathbf{m}\left(\alpha\right) = v\qquad\text{s.t. }\quad \mathbf{m}\left(\alpha\right) \text{solves LLG }\left(73\right). \end{align} \tag{ 74 }$$

The LLG equation (73) is highly nonlinear, with vectorial solution m; numerical solution of (73) is a greatly challenging task [2, 3, 5, 8] (see also [42, remark 1] for a survey). This initiated our idea of a bi-level approach, in which the lower-level problem is used to approximate m.

We invite the reader to review [31, 42] for various numerical results. To illustrate, figure 1 shows the state approximation $\textbf{m}(t,x)$ in the lower-level and the reconstruction as well as the internal operation of the bi-level scheme . The reconstruction with Landweber step size µ = 1 ran 250 upper-level iterations and a total of 11 095 lower-level iterations over 90 000 seconds. Note that this can be sped up by using the larger step size µ = 10. The final state error is below 1%, and the parameter error is below 2% [42, table 13]. In essence, the proposed bi-level scheme demonstrates its robustness through several test cases with experimental parameters, true physical parameters, noise free and noisy data. We are currently carrying out an extended numerical comparison of the three approaches—reduced, all-at-once and bi-level.

Zoom In Zoom Out Reset image size