# A Preconditioned Richardson Regularization for the Data Completion Problem and the Kozlov-Maz'ya-Fomin Method

Duc Thang Du*, Faten Jelassi*†

* LMAC, EA 2222, Université de Technologie de Compiègne Compiègne, France.
duc-thang.du@utc.fr

† LAMSIN, Faculté des Sciences de Bizerte, Jarzouna, 7021 Bizerte, Tunisia.
faten.jelassi@lamsin.rnu.tn

**ABSTRACT.** Using a preconditioned Richardson iterative method as a regularization to the data completion problem is the aim of the contribution. The problem is known to be exponentially ill posed that makes its numerical treatment a hard task. The approach we present relies on the Steklov-Poincaré variational framework introduced in [Inverse Problems, vol. 21, 2005]. The resulting algorithm turns out to be equivalent to the Kozlov-Maz'ya-Fomin method in [Comp. Math. Phys., vol. 31, 1991]. We conduct a comprehensive analysis on the suitable stopping rules that provides some optimal estimates under the General Source Condition on the exact solution. Some numerical examples are finally discussed to highlight the performances of the method.

**RÉSUMÉ.** L'objectif est d'utiliser une méthode itérative de Richardson préconditionnée comme une technique de régularisation pour le problème de complétion de données. Le problème est connu pour être sévèrement mal posé qui rend son traitement numérique ardu. L'approche adoptée est basée sur le cadre variationnel de Steklov-Poincaré introduit dans [Inverse Problems, vol. 21, 2005]. L'algorithme obtenu s'avère être équivalent à celui de Kozlov-Maz'ya-Fomin parû dans [Comp. Math. Phys., vol. 31, 1991]. Nous menons une analyse complète pour le choix du critère d'arrêt, et établissons des estimations optimales sous les Conditions Générale de Source sur la solution exacte. Nous discutons, enfin, quelques exemples numériques qui confortent les pertinence de la méthode.

**KEYWORDS:** Cauchy problem, Regularization, iterative method, Morozov's discrepancy principle.

**MOTS-CLÉS:** Problème de Cauchy, Régularisation, Méthode itérative, Principe de Morozov.
# 1. Introduction and variational formulation

Let $\Omega$ be a bounded domain in $\mathbb{R}^d$ ($d = 2, 3$) with $n$ the unit normal to the boundary $\Gamma = \partial\Omega$, oriented outward. Assume that $\Gamma$ is the union of the disjoint $\Gamma_C$ and $\Gamma_I$ that are disjoint for simplification (see Figure 1).

<figure><img src="image_1.png" /><figcaption>**Figure 1.** The boundary $\Gamma_C$ where measurements are possible and $\Gamma_I$ is unreachable.</figcaption></figure>

Assume given a datum and a flux $(g, \varphi)$ in $H^{1/2}(\Gamma_C) \times H^{-1/2}(\Gamma_C)^{(1)}$. The Cauchy problem we deal with reads as: find $u$ such that

$$(-\Delta)u = 0 \quad \text{in } \Omega, \qquad (1)$$

$$u = g \quad \text{on } \Gamma_C, \qquad (2)$$

$$\partial_n u = \varphi \quad \text{on } \Gamma_C, \qquad (3)$$

$$u = ? \quad \text{on } \Gamma_I. \qquad (4)$$

Hadamard J. demonstrates through an example given in (see [7]) that the data completion problem is ill-posed. Its severe ill-posedness for general geometries is proved by Ben Belgacem F. in [2]. A variational framework of it is proposed in [3], which consists in the duplication of the unknown $u$ into $(u_D, u_N)$ as follows: let $\lambda \in H^{1/2}(\Gamma_I)$, then $u_D = u_D(\lambda, g)$ is solution of

$$(-\Delta)u_D = 0 \quad \text{in } \Omega, \\ u_D = g \quad \text{on } \Gamma_C, \\ u_D = \lambda \quad \text{on } \Gamma_I,$$

while $u_N = u_N(\lambda, \varphi)$ satisfies the problem

$$(-\Delta)u_N = 0 \quad \text{in } \Omega, \\ \partial_n u_N = \varphi \quad \text{on } \Gamma_C, \\ u_N = \lambda \quad \text{on } \Gamma_I.$$

The key idea is to consider the common trace $(u_D)|_{\Gamma_I} = (u_N)|_{\Gamma_I} = \lambda(\in H^{1/2}(\Gamma_I))$ as the main unknown of the problem. Finding $\lambda$ allows to complete the boundary data and to obtain thereby the solution $u$ of the Cauchy problem. The $\lambda$ we look for should satisfy the flux equation,

$$\partial_n u_D(\lambda, g) = \partial_n u_N(\lambda, \varphi) \quad \text{on } \Gamma_I. \qquad (5)$$

<footnote n="1">1. The Sobolev spaces are defined in [9].</footnote>
Indeed, if $\lambda$ solves (5), it may be checked by Holmgren's Uniqueness Theorem that $u_D = u_N = (u)(\in H^1(\Omega))$ is the solution of the Cauchy problem.

In the subsequent, we are interested in the preconditioned Richardson iterative method to approximate problem (5). An outline of the paper is as follows: We recall the Steklov-Poincaré variational formulation. Then we describe the iterative Richardson method and establish the connection with the work of Kozlov V.A., Maz'ya V.G. and Fomin A.V. (see [8], [13]). Afterward, we conduct an *a-priori* and *a-posteriori* analysis of the Richardson algorithm used as a regularization strategy when associated to the Discrepancy principle. We conclude by numerical examples to illustrate the reliability of the iterative regularization.

In order to construct a variational formulation of the Cauchy problem (1)-(4), we use the notation provided in [3] for the solutions, that is

$$ (u_D(\mu), u_N(\mu)) := (u_D(\mu, 0), u_N(\mu, 0)), $$

$$ (\check{u}_D(g), \check{u}_N(\varphi)) := (u_D(0, g), u_N(0, \varphi)). $$

The variational formulation has been provided in [3] and consists in: find $\lambda \in H^{1/2}(\Gamma_I)$ such that (2): for all $\mu \in H^{1/2}(\Gamma_I)$,

$$ \int_{\Omega} \nabla u_D(\lambda) \nabla u_D(\mu) \, dx - \int_{\Omega} \nabla u_N(\lambda) \nabla u_N(\mu) \, dx = -\langle \varphi, u_N(\mu) \rangle_{1/2, \Gamma_C} \\ - \int_{\Omega} \nabla \check{u}_D(g) \nabla u_D(\mu) \, dx. $$

It may be put under the following compact form (with obvious notations)

$$ s(\lambda, \mu) = \ell(\mu), \quad \forall \mu \in H^{1/2}(\Gamma_I), \qquad (6) $$

where $s(\cdot, \cdot) = s_D(\cdot, \cdot) - s_N(\cdot, \cdot)$ and $\ell(\cdot) = \ell_D(\cdot) - \ell_N(\cdot)$. The forms $s_D(\cdot, \cdot)$, $s_N(\cdot, \cdot)$ are continuous, symmetric and elliptic on $H^{1/2}(\Gamma_I)$, and $\ell_D(\cdot)$, $\ell_N(\cdot)$ are continuous linear on $H^{1/2}(\Gamma_I)$ (see [12]). It is proved in [3, Lemma 3.3] that the bilinear form $s(\cdot, \cdot)$ is symmetric, non-negative definite which means that that $s_D(\cdot, \cdot) > s_N(\cdot, \cdot)$. In addition, $s(\cdot, \cdot)$ is compact and its eigenvalues, though non-negative, are clustered around zero which arises serious difficulties in the numerical treatment of the Cauchy problem.

REMARK. —

Each of the bilinear forms $s_D(\cdot, \cdot)$ and $s_N(\cdot, \cdot)$ are related to an inner-product on $H^{1/2}(\Gamma_I)$, and their corresponding norms on $H^{1/2}(\Gamma_I)$ are equivalent to the natural norm $\|\cdot\|_{H^{1/2}(\Gamma_I)}$, that is to say, for $\mu \in H^{1/2}(\Gamma_I)$,

$$ \|\mu\|_{s_D} = (s_D(\mu, \mu))^{1/2} \approx \|\mu\|_{H^{1/2}(\Gamma_I)}. $$

In all the sequel we use the norm $\| \cdot \|_{s_D}$ instead of $\| \cdot \|_{H^{1/2}(\Gamma_I)}$.

Before switching to the regularization issues, we provide a stability property of the linear form $\ell(\cdot)$ with respect to $s(\cdot, \cdot)$ that will play an important role in the analysis we have in mind. Let us denote first(3)

$$ \eta = |\check{u}_D(g) - \check{u}_N(\varphi)|_{H^1(\Omega)}, \qquad (7) $$

<footnote n="2">2. $\langle \cdot, \cdot \rangle_{1/2, \Gamma_C}$ is the duality pairing bilinear of $H^{-1/2}(\Gamma_C)$ and $H^{1/2}(\Gamma_C)$.</footnote>

<footnote n="3">3. The symbol $| \cdot |_{H^1(\Omega)}$ stands for the semi-norm in $H^1(\Omega)$.</footnote>
There holds that (see [5])

$$ \ell(\mu) \leq \eta \sqrt{s(\mu, \mu)}, \quad \forall \mu \in H^{1/2}(\Gamma_I). \qquad (8) $$

## 2. Preconditioned Richardson Algorithm

We aim to build a stable approximate solution of the exact one (in some sense) by the iterative Richardson method. We do not however apply it directly to the equation (6), but use it in combination with a preconditioner. We describe the preconditioner we have in mind after the construction of the operator of iterations. It is denoted by $T \in \mathcal{L}(H^{1/2}(\Gamma_I))$ and is defined as follows: for any $\lambda \in H^{1/2}(\Gamma_I)$, $T\lambda \in H^{1/2}(\Gamma_I)$ is the solution of

$$ s_D(T\lambda, \mu) = s_N(\lambda, \mu), \quad \forall \mu \in H^{1/2}(\Gamma_I). $$

$T$ is symmetric non-negative and contracting on $(H^{1/2}(\Gamma_I), s_D(\cdot, \cdot))$. As a result, the operator $(I - T)$, acting on $H^{1/2}(\Gamma_I)$, satisfies

$$ s_D((I - T)\lambda, \mu) = s(\lambda, \mu), \quad \forall \mu \in H^{1/2}(\Gamma_I). $$

It is symmetric, non-negative and contracting. It is also compact. The data $f$ of the preconditioned system is constructed as follows: find $f \in H^{1/2}(\Gamma_I)$ such that

$$ s_D(f, \mu) = \ell(\mu), \quad \forall \mu \in H^{1/2}(\Gamma_I). \qquad (9) $$

Once, all this done, it is readily checked that the Steklov-Poincaré problem (6) may be rewritten under the following form: find $\lambda \in H^{1/2}(\Gamma_I)$ such that

$$ (I - T)\lambda = f, \quad \text{in } H^{1/2}(\Gamma_I). \qquad (10) $$

We are now in position to perform the Richardson method, to the preconditioned problem (10) to compute a sequence $(\lambda_n)_n \subset H^{1/2}(\Gamma_I)$ satisfying

$$ \lambda_{n+1} - T\lambda_n = f, \quad \text{in } H^{1/2}(\Gamma_I). \qquad (11) $$

### REMARK. —

An equivalent form of (11) is as follows

$$ s_D(\lambda_{n+1}, \mu) = s_N(\lambda_n, \mu) + \ell(\mu), \quad \forall \mu \in H^{1/2}(\Gamma_I). \qquad (12) $$

This is the variational form of the iterative flux equation: find $(\lambda_n)_n \subset H^{1/2}(\Gamma_I)$ such that

$$ \partial_n u_D(\lambda_{n+1}, g) = \partial_n u_N(\lambda_n, \varphi). \qquad (13) $$

The main advantage of the variational formulation of the problem shows up in the discrete level. When the problem is approximated, say by the finite element method, we obtain a stiffness squared matrix that inherits all the properties of $s(\cdot, \cdot)$, symmetry and non-negativity definiteness. The overall regularizing tools developed for singular and ill-conditioned matrices can hence be tested on this matrix.

The general theory of the Richardson algorithm shows that under some sufficient conditions we obtain a converging regularization method. Our aim is to establish that under
some less stringent constraints the convergence is still guaranteed. In addition, when combined with the discrepancy principle the theory predicts that the Richardson method does not necessarily converge. Nevertheless, we are able to prove that it actually converges. The mile stone for the analysis we have in mind is the following stability on $f$ which is straightforwardly issued from (8) and (9).

**Lemma 1.** We have that

$$ (f, \mu)_{s_D} \le \eta \|(I-T)^{1/2} \mu\|_{s_D}, \quad \forall \mu \in H^{1/2}(\Gamma_I). \quad (14) $$

REMARK. —

A particular consequence of (14) is that

$$ \|f\|_{s_D} \le \eta. $$

### 3. Connection with Kozlov-Maz'ya-Fomin's Method

The first step of our analysis is the equivalence between the preconditioned Richardson procedure and the one proposed by Kozlov, Maz'ya and Fomin in 1991 (see [8, 13]), currently pointed at as the (KMF)-method. Actually, they are the same method written in different ways. By this result, one can see why the convergence of the algorithm occurs only when the Cauchy data is exempted from errors, i.e. when $\ell \in \mathcal{R}(s)$ (or equivalently when $f \in \mathcal{R}(I-T)$). We recall the construction of the (KMF)-sequences. Assume $(u_n)$ is known, then $(v_n, u_{n+1})$ are calculated as solutions of the following iterative boundary value problems, the function $v_n$ is the solution of

$$ \begin{aligned} (-\Delta)v_n &= 0, & \text{in } \Omega \\ \partial_n v_n &= \varphi, & \text{on } \Gamma_C, \\ v_n &= u_n, & \text{on } \Gamma_I. \end{aligned} \quad (15) $$

and $u_{n+1}$ satisfies the boundary value problem

$$ \begin{aligned} (-\Delta)u_{n+1} &= 0, & \text{in } \Omega \\ u_{n+1} &= g, & \text{on } \Gamma_C, \\ \partial_n u_{n+1} &= \partial_n v_n, & \text{on } \Gamma_I. \end{aligned} \quad (16) $$

The following equivalence holds.

**Proposition 1.** Let $(\lambda_n)_n \subset H^{1/2}(\Gamma_I)$ be the solution of the preconditioned Richardson iterative procedure (11) and $(v_n, u_n)_n \subset H^1(\Omega) \times H^1(\Omega)$ be provided by the KMF method. Then we have that $u_n = u_D(\lambda_n, g)$ and $v_n = u_N(\lambda_n, \varphi)$.

*Proof.* Let the Cauchy data $(g, \varphi) \in H^{1/2}(\Gamma_C) \times H^{-1/2}(\Gamma_C)$ be given. Set $v_n = u_N(\lambda_n, \varphi)$ and $u_n = u_D(\lambda_n, g)$. On one hand side, the two first lines of (15) are fulfilled and the condition on $\Gamma_I$ is given by

$$ v_n = u_n (= \lambda_n), \quad \text{on } \Gamma_I. $$
On the other hand side, the two first lines of (16) hold true. Furthermore, due to (5) and (13), we are allowed to write that

$$ \partial_{\mathbf{n}} u_{n+1} = \partial_{\mathbf{n}} v_n, \quad \text{on } \Gamma_I. $$

The proof is complete.

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

The convergence of the sequences $(u_n)_n$ and $(v_n)_n$ towards the solution $(u)$ of the Cauchy problem when the data $(g, \varphi)$ are exact comes directly from the Proposition 1 and ([8, 13]).

**Corollary 1.** Assume that $f \in \mathcal{R}(I - T)$, then the sequence $(\lambda_n)_n$ provided by (11) converges towards $\lambda$ as $n \to \infty$. In the contrary, if $f \notin \mathcal{R}(I - T)$, the sequence blows-up.

## 4. Stopping rules

In reality, the exact linear form $\ell(\cdot)$ (and by then $f$) is not available which means that $(g, \varphi)$ are not known exactly. Only perturbed data are accessible $(\tilde{g}, \tilde{\varphi}) = (g + \delta g, \varphi + \delta \varphi)$. The noise level affecting Cauchy data is assumed to be $\epsilon' > 0$, that is,

$$ \|\delta g\|_{H^{1/2}(\Gamma_C)} + \|\delta \varphi\|_{H^{-1/2}(\Gamma_C)} \le \epsilon'. \quad (17) $$

This induces a perturbation on the data $f$ in equation (10) since it suffers also from a deviation. Instead of the exact $f$, we dispose of $\tilde{f} = f + \delta f$. Formula (14) allows to show that

$$ \|\delta f\|_{s_D} \le C\epsilon' \quad (:= \epsilon) $$

We denote by $(\tilde{\lambda}_n)_n$ the perturbed sequence obtained by the Richardson method with $f$ replaced by $\tilde{f}$. We start from $\tilde{\lambda}_0 = \lambda_0(= 0)$ and the induction henceforth reads as

$$ \tilde{\lambda}_{n+1} - T\tilde{\lambda}_n = \tilde{f}. $$

The aim is to discuss now the issue of the stopping rule $n^* = n^*(\epsilon)$ to ensure that the sequence $(\tilde{\lambda}_{n^*})_{\epsilon>0}$ converges toward the exact solution $\lambda$ in $H^{1/2}(\Gamma_I)$ when $\epsilon$ decays. The parameter choice $n^*$ guarantees that the Richardson algorithm results in a convergent regularization. The analysis relies on the bias-variance decomposition (with respect to the $s_D$-norm $\| \cdot \|_{s_D}$)

$$ \|\lambda - \tilde{\lambda}_n\|_{s_D} \le \|\lambda - \lambda_n\|_{s_D} + \|\lambda_n - \tilde{\lambda}_n\|_{s_D}. \quad (18) $$

The *bias error* is the error caused by the iterative method, while the *variance error* describes the effect of the erroneous measurements. The choice of $n^*$ depends on the smoothness of the exact solution $\lambda$ and can be made *a-priori* if that information is available. More likely, it may be achieved thanks to an *a-posteriori* artifice, the Discrepancy Principle of Morozov, that does not need any other information than the approximated solution. Henceforth we set $e_n = (\lambda_n - \lambda)$ and $\tilde{e}_n = (\tilde{\lambda}_n - \lambda_n)$.
## 4.1. An a-priori choice

To fix a-priori the stopping rule of the algorithm we conduct the convergence, starting by the variance error.

**Lemma 2.** There holds that

$$ \|\tilde{\lambda}_n - \lambda_n\|_{s_D} \le (\sqrt{n})\epsilon. \qquad (19) $$

If the iterations are stopped at $n^* = n^*(\epsilon)$ such that

$$ \lim_{\epsilon \to 0} n^* = +\infty, \quad \lim_{\epsilon \to 0} (\sqrt{n^*})\epsilon = 0. \qquad (20) $$

The variance error decays to zero.

*Proof.* We have that

$$ \tilde{\lambda}_{n+1} - T\tilde{\lambda}_n = \tilde{f}, \quad \text{and} \quad \lambda_{n+1} - T\lambda_n = f. $$

Subtracting one equation from the other, we obtain that

$$ \tilde{e}_{n+1} = T\tilde{e}_n + (\delta f) = \sum_{k=0}^{n} T^k (\delta f) = (I - T^{n+1})(I - T)^{-1}(\delta f). $$

In view of the symmetry of $T$ and $(I - T)$, we have

$$ \|\tilde{e}_n\|_{s_D}^2 = ((\delta f), (I - T^n)(I - T)^{-1}(\tilde{e}_n))_{s_D}. $$

Now, owing to the stability (14), we derive that

$$ \begin{aligned} \|\tilde{e}_n\|_{s_D}^2 &\le (\delta\eta)\|(I - T)^{1/2}(I - T^n)(I - T)^{-1}(\tilde{e}_n)\|_{s_D} \\ &\le (\delta\eta)\|(I - T^n)(I - T)^{-1/2}\|_{s_D} \|\tilde{e}_n\|_{s_D}. \end{aligned} $$

We obtain therefore

$$ \|\tilde{e}_n\|_{s_D} \le (\delta\eta)\|(I - T^n)(I - T)^{-1/2}\|_{s_D} \le (\delta\eta)\sqrt{n}. $$

Next, by (17), we have

$$ (\delta\eta) = |\tilde{u}_D(\delta g) - \tilde{u}_N(\delta\varphi)|_{H^1(\Omega)} \le \epsilon, $$

and hence

$$ \|\tilde{e}_n\|_{s_D} \le (\sqrt{n})\epsilon. $$

The proof is complete with $n = n^*(\epsilon)$ chosen as in (20). □

## REMARK. —

In the general case, the convergence of the variance error is ensured under the sufficient condition $(n\epsilon)$ goes to zero. Our result is better and remind the behavior of the Landweber method (see [6]). The combination of bias-variance decomposition, Proposition 1 and Lemma 2 make out of the preconditioned Richardson regularization, with the a-priori choice of $n^*$, a convergent strategy.

Now, the convergence of the bias error being arbitrary slow (see [6]) we may be tempted
to consider the effect of some regularity on the exact solution. We introduce here the General Source Condition (GSC), currently employed for ill-posed problems (see, e.g., [6]) that is $\lambda \in \mathcal{R}(I - T)$:

$$ \text{there exists } \chi \in H \text{ such that } \lambda = (I - T)\chi. \qquad (21) $$

We set $E = \|\chi\|_{s_D}$. The following lemma gives the convergence of the bias error.

**Lemma 3.** Assume that the solution $\lambda$ of problem (12) satisfies the (GSC) assumption. Then the following bound holds

$$ \|\lambda_n - \lambda\|_{s_D} \le \frac{E}{2n}. $$

*Proof.* Without loss of generality, we set $\lambda_0 = 0$. We have that,

$$ e_n = -T^n \lambda = -T^n (I - T) \chi. $$

We obtain therefore,

$$ \|e_n\|_{s_D} = \|T^n(I - T)\chi\|_{s_D} \le E\|T^n(I - T)\|_{s_D} \le \frac{E}{2n}. $$

The proof is complete.

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

The following theorem provides an optimal stopping rule ($n^* = n^*(\epsilon)$) when the (GSC) assumption is fulfilled.

**Theorem 1.** Let the solution $\lambda$ of problem (12) satisfy the (GSC) assumption. Then, choosing $n^*(\epsilon) = E^{2/3}\epsilon^{-2/3}$ yields the bound

$$ \|\tilde{\lambda}_n - \lambda\|_{s_D} \le \frac{3}{2} E \left(\frac{\epsilon}{E}\right)^{2/3}. \qquad (22) $$

*Proof.* The above results give

$$ \|\tilde{\lambda}_n - \lambda\|_{s_D} \le \frac{E}{2n} + \sqrt{n\epsilon}. $$

Choosing $n^* = n^*(\epsilon)$ so that

$$ n^*(\epsilon) = E^{1/3}\epsilon^{-1/3}, $$

yields the expected estimate. The proof is complete.

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

**REMARK.** —

Theorem 1 gives a stopping rule of iteration $n = O(\epsilon^{-2/3})$, to ensure the convergence. An alternative consists in using an a-posteriori rule; the one we select here is the Discrepancy Principle of Morozov.
## 4.2. The Discrepancy Principle

Generally, an a-priori choice of the stopping parameter is impossible, because of the difficulty to access an accurate information on the smoothness of $\lambda$. We may therefore match the Discrepancy Principle of Morozov (see [10]). The criterion we adopt here is based on the minimum value of the Kohn-Vogelius functional, defined as follows (see [11, 3]):

$$ \tilde{K}(\mu) = \frac{1}{2} |u_D(\mu, \tilde{g}) - u_N(\mu, \tilde{\varphi})|_{H^1(\Omega)}^2, \quad \forall \mu \in H^{1/2}(\Gamma_I). $$

That function is called $K(\cdot)$ when $(\tilde{g}, \tilde{\varphi})$ are replaced by $(g, \varphi)$ and $K^\delta(\cdot)$ if they are replaced by $(\delta g, \delta \varphi)$. Before the description of the strategy we intend to use, we need to state the monotonicity of the sequence $(K(\lambda_n))_n$. Some formulas will help in the analysis. We recall first that, in [5], is stated that

$$ K(\mu) = \frac{1}{2} s(\mu, \mu) - \ell(\mu) + \eta^2, \quad \forall \mu \in H^{1/2}(\Gamma_I). \quad (23) $$

Identity (23) still holds for $\tilde{K}(\cdot)$ and $K^\delta(\cdot)$ where $\eta$, given in (7), is changed into $\tilde{\eta}$ and $(\delta\tilde{\eta})$ with obvious notations. Moreover, the following holds (see [4])

$$ 2K(\mu) = s(\mu - \lambda, \mu - \lambda), \quad \forall \mu \in H^{1/2}(\Gamma_I). \quad (24) $$

**Lemma 4.** The sequences $(K(\lambda_n))_n$, $(\tilde{K}(\lambda_n))_n$ and $(K^\delta(\lambda_n))_n$ are decreasing.

*Proof.* Owing to (12) and the Lax-Milgram theorem, we derive that

$$ \frac{1}{2} s_D(\lambda_{n+1}, \lambda_{n+1}) - s_N(\lambda_n, \lambda_{n+1}) - \ell(\lambda_{n+1}) \le \frac{1}{2} s_D(\lambda_n, \lambda_n) - s_N(\lambda_n, \lambda_n) - \ell(\lambda_n). $$

As a result, we obtain that

$$ K(\lambda_{n+1}) + \frac{1}{2} s_N(\lambda_{n+1}, \lambda_{n+1}) + s_N(\lambda_n, \lambda_{n+1}) \le K(\lambda_n) - \frac{1}{2} s_N(\lambda_n, \lambda_n), $$

and therefore

$$ K(\lambda_{n+1}) + \frac{1}{2} s_N(\lambda_{n+1} - \lambda_n, \lambda_{n+1} - \lambda_n) \le K(\lambda_n). $$

The non-negativity of $s_N(\cdot, \cdot)$ completes the proof. $\square$

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

Let us now remark that the exact solution $\lambda \in H^{1/2}(\Gamma_I)$ of the entire problem satisfies the bound

$$ (2\tilde{K}(\lambda))^{1/2} = |\tilde{u}_D(\delta g) - \tilde{u}_N(\delta \varphi)|_{H^1(\Omega)} = (\delta \eta) \le \epsilon. $$

We fix $\sigma > 1$. The Discrepancy Principle suggests to interrupt the iterations at the first time $n^* = n^*(\epsilon)$ where the following inequality holds true

$$ (2\tilde{K}(\tilde{\lambda}_{n^*}))^{1/2} \le \sigma\epsilon. \quad (25) $$

Such a procedure yields a convergence rate that can be compared to the a-priori choice. We need the following preparatory Lemmas.
**Lemma 5.** Assume that the solution $\lambda$ of problem (12) satisfies the (GSC) assumption. Under the stopping iteration $n^* = n(\epsilon)$ provided in (25), the variance bound holds

$$ \|\tilde{\lambda}_{n^*} - \lambda_{n^*}\|_{s_D} \leq CE \left(\frac{\epsilon}{E}\right)^{2/3}. $$

*Proof.* Due to the Discrepancy Principle (25) we have

$$ \begin{aligned} \sigma\epsilon &< |u_D(\tilde{\lambda}_{(n^*-1)}, \tilde{g}) - u_N(\tilde{\lambda}_{(n^*-1)}, \tilde{\varphi})|_{H^1(\Omega)} \\ &< |u_D(\lambda_{(n^*-1)}, g) - u_N(\lambda_{(n^*-1)}, \varphi)|_{H^1(\Omega)} + |u_D(\tilde{e}_{(n^*-1)}, \delta g) - u_N(\tilde{e}_{(n^*-1)}, \delta \varphi)|_{H^1(\Omega)} \\ &= (2K(\lambda_{(n^*-1)}))^{1/2} + (K^\delta(\tilde{e}_{(n^*-1)}))^{1/2}. \end{aligned} $$

That $(K^\delta(\tilde{e}_n))_n$ is decreasing yields

$$ \begin{aligned} \sigma\epsilon &< (2K(\lambda_{(n^*-1)}))^{1/2} + (K^\delta(0))^{1/2} \\ &\leq (2K(\lambda_{(n^*-1)}))^{1/2} + (\delta\eta) \leq (2K(\lambda_{(n^*-1)}))^{1/2} + \epsilon. \end{aligned} $$

We deduce therefore that

$$ (\sigma - 1)\epsilon < (2K(\lambda_{(n^*-1)}))^{1/2}. \qquad (26) $$

Now, on account of (24) we write that

$$ \begin{aligned} (K(\lambda_n))^{1/2} &= (s(\lambda_n - \lambda, \lambda_n - \lambda))^{1/2} = (((I-T)e_n, e_n)_{s_D})^{1/2} \\ &= \|(I-T)^{1/2}e_n\|_{s_D} = \|(I-T)^{3/2}T^n(\chi)\|_{s_D} \leq \frac{E}{(n+1)^{3/2}}. \end{aligned} $$

Back to (26), we come up with the bound

$$ (\sigma - 1)\epsilon \leq \frac{E}{(n^*)^{3/2}}, $$

from which we derive that

$$ n^* \leq C \left(\frac{E}{\epsilon}\right)^{2/3}. $$

Calling back the estimate (19), we establish that

$$ \|\tilde{\lambda}_{n^*} - \lambda_{n^*}\|_{s_D} \leq (\sqrt{n^*})\epsilon \leq CE \left(\frac{\epsilon}{E}\right)^{2/3} $$

This completes the proof.

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

**Lemma 6.** Under the stopping iteration $n^* = n^*(\epsilon)$ chosen by the Discrepancy Principle, there holds that

$$ \|\lambda - \lambda_{n^*}\|_{s_D} \leq DE \left(\frac{\epsilon}{E}\right)^{2/3}. $$

*Proof.* Let us denote $\chi_n \in H^{1/2}(\Gamma_I)$ such that $\lambda_n = (I - T)\chi_n$. The existence of $\chi_n$ may be obtained by induction with $\|\chi_{n^*} - \chi\|_{s_D} \leq E$. Invoking the interpolation inequality (see [6, 4]), we have by (24)

$$ \|e_{n^*}\|_{s_D} \leq (s(e_{n^*}, e_{n^*}))^{1/3} \|\chi_{n^*} - \chi\|_{s_D}^{1/3} = (2K(\lambda_{n^*}))^{1/3} \|\chi_{n^*} - \chi\|_{s_D}^{1/3}. $$
On the other hand side, by the Discrepancy Principle (25), we obtain that

$$ (2K(\lambda_{n^*}))^{1/2} \le |u_D(\tilde{\lambda}_{n^*}, \tilde{g}) - u_N(\tilde{\lambda}_{n^*}, \tilde{\varphi})|_{H^1(\Omega)} + |u_D(\tilde{e}_{n^*}, \delta g) - u_N(\tilde{e}_{n^*}, \delta \varphi)|_{H^1(\Omega)} \\ \le (2\tilde{K}(\tilde{\lambda}_{n^*}))^{1/2} + (2K^\delta(\tilde{e}_{n^*}))^{1/2} \le \sigma\epsilon + (\delta\eta) \le (\sigma + 1)\epsilon. $$

Given that $\|\chi_{n^*} - \chi\|_{s_D} \le E$, replacing in the above formula completes the proof. $\square$

**REMARK.** —

The constants C and D are known explicitly and are given by

$$ C = (\sigma - 1)^{-2/3}, \quad D = (\sigma + 1)^{2/3}. $$

Putting together Lemmas 5 and 6 provides the following global convergence result.

**Theorem 2.** Assume that the solution $\lambda$ of problem (12) satisfies the (GSC) assumption. Then, choosing $n^* = n^*(\epsilon)$, thanks to the Morozov principle (25), the following optimal bound holds

$$ \|\tilde{\lambda}_{n^*} - \lambda\|_{s_D} \le CE \left(\frac{\epsilon}{E}\right)^{2/3}. \qquad (27) $$

**REMARK.** —

The order of the a-posteriori estimate is $\epsilon^{2/3}$. That is to say the preconditioned Richardson iterative method combined to the Discrepancy Principle of Morozov (25) provides a super-convergent regularization strategy.

**REMARK.** —

The General Source Condition $\lambda \in \mathcal{R}(I-T)$ can be weakened to $\lambda \in \mathcal{R}((I-T)^\alpha)$, where $\alpha \in (0, 1)$ or to the logarithmic type, i.e $\lambda \in \mathcal{R}((\log(I-T)^{-1})^\alpha)$. The convergence rates we are able to exhibit remain of optimal order.

**REMARK.** —

Both Theorems 1 and 2 show that the Richardson iterative method behaves like the Landweber algorithm (see [6]). This suggests that whatever the Cauchy boundary conditions $(g, \varphi)$ are in $H^{1/2}(\Gamma_C) \times H^{-1/2}(\Gamma_C)$, the resulting data $f$ obtained by (9) belongs to $\mathcal{R}(\sqrt{I-T})$, which make us believe that the variational formulation data completion problem is actually the Euler-Lagrange of some least-squares problems.

# 5. Numerical Examples

The computations are realized by a finite element method. The meshes we use are all triangular(⁴) and uniform, the finite elements are linear(⁵). The advantage of the Steklov-Poincaré equation (10) is that it yields a symmetric and non-negative definite stiffness matrix (see [5]). In our experiences, the inversion of the algebraic equations is realized by the preconditioned Richardson algorithm. Nevertheless, different methods may tested to solve that equation such as Krylov subspaces type methods.

The purpose of the first experience is to assess the solution obtained by means of the Discrepancy Principle combined to the iterative preconditioned Richardson method. The

<footnote n="4">4. The meshes are generated by the EMC2 public software; it can be downloaded from the INRIA Web-site,</footnote>

<footnote n="5">5. The primary procedures, in Fortran 77, are provided on the O. Pironneau home-page,</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 2.** The squared domain with two holes.</figcaption></figure>

domain $\Omega$ is the unity square with two holes in it as indicated in Figure 2. The aim is the reconstruction of the solution of the data completion problem.

The Cauchy boundary $\Gamma_C$ is the external boundary of the square and the two holes are the incomplete boundary $\Gamma_I$. Let $w$ be defined by

$$w(x, y) = \cosh(\frac{\pi}{2}x) \cos(\frac{\pi}{2}y) + \sin(\frac{\pi}{2}x) \sinh(\frac{\pi}{2}y).$$

The data on $\Gamma_C$ are $(g, \varphi) = (w, \partial_n w)$ so that the exact solution of the Cauchy problem (1)-(4) is $u = w$. A random noise is generated on the data $(g, \varphi)$ by the Fortran random function ( ). In this computation the magnitude of the noise is 10% with respect to

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Computed solutions, $u_D$ to the left and $u_N$ to the right for 10% -noise. The gap between $u$ and $u_N$ is measured to 0.061 with respect to the maximum norm.</figcaption></figure>

the max norm. Notice that the $\epsilon$ is so that

$$\epsilon = (\delta\eta) = |\check{u}_D(\delta g_C) - \check{u}_N(\delta\varphi_C)|_{H^1(\Omega)}, \quad (28)$$

necessary in the Discrepancy Principle is not at hand. The deviation (28) is not available to us and needs to be estimated. It is dependent on multiple factors and especially of the domain $\Omega$ and has to be re-evaluated for each experiment through some preparatory calculations. The evaluation of the incidence of the noise level (0.1) on the value of $\epsilon$ is achieved by conducting several computations where there is no signal on the data
<figure><img src="image_9.png" /><figcaption>**Figure 4.** Curves of exact and computed solutions $\lambda$ and $\tilde{\lambda}_{n*}$ along both holes contours with 10% noise. Quadrangular's contour to the left and circular's contour to the right.</figcaption></figure>

$(g, \varphi) = (0, 0)$ which suffers however from different noises the magnitude of which does not exceed 0.1. Next, sorting the obtained results and taking the mean-value of the different $(\delta\eta)$ allows to approximate $\epsilon$. The first row of Table 1 provides $\epsilon$ for different $(L^\infty)$-magnitudes of the noise, it seems to grow linearly with respect to the maximum norm of the pollution.

<table><caption>**Table 1.** The third row gives the stopping iteration by the discrepancy principle. The last row indicates the accuracy of the computed solution. The notation $(u-u_N)_r$ is for $\max |u-u_N|/\max |u|$.</caption><tr><td>max (Noise)</td><td>0.01</td><td>0.05</td><td>0.025</td><td>0.075</td><td>0.1</td></tr><tr><td>ε</td><td>0.1466</td><td>0.3741</td><td>0.7518</td><td>1.003</td><td>1.544</td></tr><tr><td>n*</td><td>48</td><td>17</td><td>8</td><td>6</td><td>5</td></tr><tr><td>(u - u<sub>N</sub>)<sub>r</sub></td><td>0.038</td><td>0.045</td><td>0.05</td><td>0.056</td><td>0.061</td></tr></table>

Currently, the parameter $\sigma$ is fixed very close to unity for the discrepancy equation (25). Although $\sigma = 1$ is not tolerated by the analysis, we observe that doing so yields very satisfactory results (the selection $\sigma = 1.01$ results also in good results). In fact, those in Table 1 are obtained by $\sigma = 1$. Figure 3 depicts the approximated solutions $u_D(\tilde{\lambda}_{n*}, \tilde{g})$ and $u_N(\tilde{\lambda}_{n*}, \tilde{\varphi})$ when the data are polluted by a 10%-noise. For this experiment the iterations are stopped at $n^* = 5$ and the computed solutions seem reliable. Observe that the noise on the Neumann conditions (the solution $u_N$ is regular at the boundary $\Gamma_C$) is less visible than the one affecting the Dirichlet conditions (the oscillations on $u_D$ at the boundary $\Gamma_C$ may be observed). Notice that when the algorithm is not interrupted by the discrepancy principle, the accuracy of the computations is slightly improved in a first stage, before blowing-up. We also represent in Figure 4 the curves of the exact solution $u$ so as the approximated one $u_N(\tilde{\lambda}_{n*}, \tilde{\varphi})$ on the two components of $\Gamma_I$, the quadrangular and the circular part. The (local) relative gap between the two curves is evaluated independently for each component. It equals 0.065 for the quadrangle (curves to the left in Figure 4) and is a little higher for the circular contour 0.12 (curves to the right in Figure 4). These may be considered as acceptable results for computations with data influenced by 10%-noise.

The second experiment involved in completing the data where the solution $u$ suffers from the proximity of a singularity. The domain here is $\Omega = ] -0.5, 0.5[ \times ]0.05, 1[$. The
<figure><img src="image_10.png" /><figcaption>Figure 5. The approximated solutions for the nearly singular test, $u_D$ to the left and $u_N$ to the right.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 6. The solution $\lambda$ and the computed one $\tilde{\lambda}_{n*}$. The front is smoothed for $\tilde{\lambda}_{n*}$.</figcaption></figure>

Cauchy boundary $\Gamma_C$ is the union of the vertical sides $\{-0.5, 0.5\} \times ]0.05, 1[$ and the upper boundary $]-0.5, 0.5[ \times \{1\}$, and the lower side is the incomplete boundary $\Gamma_I = [-0.5, 0.5] \times \{0.05\}$. The goal is the reconstruction by the preconditioned Richardson method of the solution

$$v(x, y) = \sqrt{x+r}, \quad r = \sqrt{x^2 + y^2}$$

That solution corresponds to the singularity generated by horizontal crack with the tip located at the origin. That singular point is narrowly close to $\Gamma_I$, the solution may present a stiff front at the vicinity of the origin and the reconstruction of $\lambda = u_{|\Gamma_I}$ is expected to be harder than for the first test. In Figure 5 and Figure 6 we observe a smoothing of the computed solution at the middle of $\Gamma_I$. The magnitude of the noise is 5% which corresponds to $\epsilon = 0.3332$. Here again we fix $\sigma$ to unity for the Discrepancy Principle stopping rule. The solver requires $n^* = 81$ iterations to achieve the convergence. The relative maximum error for the global solution on the whole computational domain is 0.073 which is pretty satisfactory.
## 6. Conclusion

The purpose here is to apply the preconditioned Richardson iterative method to the variational formulation, introduced in [3], for the Cauchy problem. When the Cauchy data are not compatible, no solution is available, thus the approximate solution can not converge. With a suitable choice of the iteration stopping index $n$ depending on the noise level $\epsilon$, we can guarantee the convergence of computed solution $\tilde{\lambda}_n$ toward the exact one $\lambda$. We address here the two ways to achieve the parameter selection. An optimal *a-priori* criterion has been discussed and studied. Afterward, the *a-posteriori* rule of the discrepancy principle of Morozov, based on the Kohn-Vogelius functional, provides a monotonic method which result in a convergent strategy. As checked out in Section 3, the regularizing algorithm we obtain is nothing else than a different form of the Kozlov-Maz'ya-Fomin method (see [8]). The advantages are twofold. It is possible to conduct a numerical analysis on a standard equation with symmetric and non-negative operator. Similarly, in the computational ground, users are left to handle symmetric algebraic systems to which they may apply a solver of their choice (see [14]).

**Acknowledgments.** We are deeply grateful to Ben Belgacem F. for his careful reading of the paper. We also thank very much anonymous reviewers for their valuable comments. This work was partially supported by la région de Picardie, (FRANCE) under the program Appui à l'émergence and le Ministère de l'Enseignement Supérieur, de la Recherche Scientifique et de la Technologie (TUNISIA) under the LR99ES-20 program.

## 7. References

[1] AZAÍEZ M., BEN ABDA A., BEN ABDALLAH J., Revisiting the Dirichlet-Neumann Solver for the Cauchy System. *Intern. J. Appl. Math. and Mech.* vol. 1, 2005.

[2] BEN BELGACEM F., "Why is the Cauchy's problem severely ill-posed?", *Inverse Problems*, vol. 23, num. 2, 2007.

[3] BEN BELGACEM F., EL FEKIH H., "On Cauchy's problem: I. A variational Steklov-Poincaré theory", *Inverse Problems*, vol. 21, num. 6, 2005.

[4] BEN BELGACEM F., EL FEKIH H., JELASSI F., "The Lavrentiev regularization of the data completion problem", *Inverse Problems*, vol. 24, num. 4, 2008.

[5] AZAÍEZ M., BEN BELGACEM F., EL FEKIH H. "On Cauchy's problem: II. Completion, regularization and approximation", *Inverse Problems*, 22 (doi:10.1088/0266-5611/22/4/012):1307-1336, 2006.

[6] ENGL H.W., HANKE M., NEUBAUER A., "Regularization of Inverse Problems", 2000.

[7] HADAMARD J., "Lectures on Cauchy's problem in linear partial differential equation.", New York, Dover, 1953.

[8] KOZLOV V.A., MAZ'YA V.G., FOMIN A.V., "An iterative method for solving the Cauchy problem for elliptic equations", *Comp. Math. Phys.*, vol. 31, num. 1, 1991.

[9] LIONS J-L., MAGENES E., "Problème aux limites nonhomogènes et applications", *Travaux et Recherches Mathématiques*, vol. 1, num. 17, 1968.

[10] MOROZOV V.A., "On the solution of functional equations by the method of regularization", *Doklady*, vol. 167, num. 3, 1966.

[11] KOHN R., VOGELIUS M. "Determining conductivity by boundary measurements II. Interior results", *Comm. Pure and Appl. Math.*, XXXVIII:643-667, 1985.
[12] QUARTERONI A., VALLI A. “Domain Decomposition Methods in Partial Differential Equations”, *Oxford University Press*, 2004.

[13] ENGL H. W., LEITAŌ A. “A Mann iterative regularization method for elliptic Cauchy problems”, *Numer. Funct. Anal. Optimiz.*, 22 (7&8: 861 - 884) 2001.

[14] BEN ABDALLAH J., “A Conjugate Gradient Type Method for the Steklov-Poincaré Formulation of the Cauchy-Poisson Problem”, *Intern. J. Appl. Math. and Mech.* vol. 3, 27-40, 2007.