Solving ill-posed control problems by stabilized finite element methods: an alternative to Tikhonov regularization

Tikhonov regularization is one of the most commonly used methods of regularization of ill-posed problems. In the setting of finite element solutions of elliptic partial differential control problems, Tikhonov regularization amounts to adding suitably weighted least squares terms of the control variable, or derivatives thereof, to the Lagrangian determining the optimality system. In this note we show that stabilization methods for discretely ill--posed problems developed in the setting of convection--dominated convection--diffusion problems, can be highly suitable for stabilizing optimal control problems, and that Tikhonov regularization will lead to less accurate discrete solutions. We consider data assimilation problems for Poisson's equation as illustration and derive new error estimates both for the the reconstruction of the solution from measured data and reconstruction of the source term from measured data. These estimates include both the effect of discretization error and error in measurements.


Introduction
In this note we propose an alternative to the classical Tikhonov regularization approach in finite element approximations of optimal control problems governed by elliptic partial differential equations. We shall, following [2], consider problems of the type where J is a cost functional, A is an elliptic differential operator for the state variable u, and B an impact operator for the control variable q. Introducing costate variables λ, this problem can be formulated as finding saddle points to the Lagrangian functional where (·, ·) denotes the L 2 inner product, determined by the system In a finite element setting, the continuous states, controls, and costates (u, q, λ) ∈ V × Q × V are replaced by their discrete counterparts (u h , q h , λ h ) ∈ V h × Q h × V h , where V h and Q h are finite dimensional counterparts of the appropriate Hilbert spaces V and Q, respectively. Typically, the cost functional measures some distance between the discrete state and a known or sampled function u 0 over a subdomain M ⊆ Ω where Ω ⊂ R d , d = 2, 3 is a polyhedral (polygonal) domain of computation. which may not lead to a well posed problem. A classical regularization method due to Tikhonov, see [11], is to add a stabilizing functional n(q, q), where, typically, (6) n(q, q) := α q − q b 2 + β ∇(q − q b ) 2 , with α and β regularization parameters and q b the background state, or first guess state. The role of the background state is to diminuish the nonconsistent character of the Tikhonov regularization and implies additional a priori knowledge on the system beyond the samples u 0 . In this note we will assume that no such additional a priori data are at hand, q b = 0 and that there is no physical justification for the addition of the term n(q, q), or that the size of α and β given by the application are too small to provide sufficient stabilization of the system for computational purposes. The objective of this note is to show that these regularizations can be improved upon in a finite element framework. The approach that we will follow is to eliminate the Tikhonov regularization on the continuous level and instead regularize the discrete formulation hence making the regularization part of the computational method in the form of a weakly consistent stabilization. The terminology stabilization versus regularization is slightly ambiguous but in classical numerical analysis the method of modified equations [10] provides a link between these concepts. 1.1. Model problems. We will discuss two model problems below. First let M ⊂ Ω be a ball B r 1 (x 0 ) with radius r 1 , centered at x 0 and assume that measurements u 0 are available in this ball. We then wish to reconstruct the solution u in B r 2 (x 0 ) under the a priori assumption that B r 2 (x 0 ) ⊂ B r 3 (x 0 ) ⊂ Ω and u ∈ H 1 (Ω), is a weak solution to −∆u = f, in Ω.
This problem can be cast in the form of a constrained minimization problem: subject to (8) − ∆u = f, in Ω.
We will refer to this problem as the data assimilation problem below, but it is also strongly related to boundary control problems. It is known that in this case if for α = 0 u 0 is such that a solution exists, then by unique continuation of harmonic functions this solution is unique. This statement can be quantified in the following three sphere's inequality: Lemma 1. (Three spheres inequality) Assume that u : Ω → R is a weak solution of (8) with f ∈ H −1 (Ω) such that f H −1 (Ω) ≤ ε for some ε > 0. For every r 1 , r 2 , r 3 , r such that 0 < r 1 < r 2 < r 3 < r and for every x 0 ∈ Ω such that dist(x 0 , ∂Ω) > r there holds where B r i , i = 1, 2, 3 are balls centered at x 0 , C > 0 and τ, 0 < τ < 1 only depend on the geometry of Ω, and the ratios r 2 /r 1 and r 3 /r 2 .
Proof. For a proof in the non-homogeneous case considered here see Allessandrini et al. [1,Theorem 1.10].
Remark 1. We do not track the exact form of the geometric constants that appear in the proof of Lemma 1. They are all considered included in the canonical constant C above. For the precise definition of the result with all exact dependencies we refer to [1].
Introducing the Lagrange multiplier λ ∈ H 1 0 (Ω), we have the optimality system: find (u, λ) ∈ H 1 (Ω)×H 1 0 (Ω) such that The difficulty in this case is that the equation (11) is ill-posed, since u ∈ H 1 (Ω) and µ ∈ H 1 0 (Ω). It follows that coercivity fails for u in H 1 (Ω). Clearly this is not the case when α > 0, but then a nonconsistent perturbation is added to the system that can not easily be quantified. Below we will instead add a regularization on the discrete level. Indeed we will penalize the fluctuations of the gradient over element faces and show that the added coercivity on the high frequency content of the solution is sufficient to obtain a priori estimates leading to error estimates through Lemma 1. This part of the analysis draws on earlier ideas for the elliptic Cauchy problem from [4,5]. Below we will assume that u 0 is the unperturbed measurement for which the unique solution exists and consider a numerical method with perturbed data.
Our second example considers the case where the data is available in the whole domain, M ≡ Ω, but the source term is unknown and must be reconstructed. The challenge here being that the application of the Laplacian is unstable. This case will be referred to as source reconstruction below, but is also related to a distributed control problem. We consider the elementary problem: minimize Here, u 0 is known data and Ω is a convex polygonal (polyhedral) subset of R d , d = 2, 3, with outword pointing normal n. We assume that we wish to solve (12)-(13) in the situation where u 0 is a measurement on a system that is of the form (13). This means that if no perturbations are present in the data and measurements are available in every point of Ω, the minimizer for α = 0 is u = u 0 ∈ H 1 0 (Ω)∩H 2 (Ω) =: W and an associated q = −∇ 2 u ∈ L 2 (Ω) exists so that (13) is satisfied. Assume also that the shift theorem u H 2 (Ω) ≤ C q L 2 (Ω) holds. Below we will consider the problem of reconstructing q from u 0 accounting for both the discretization error and the error due to errors in the measured data u 0 .
Introducing the Lagrange multiplier λ ∈ H 1 0 (Ω), we have the optimality system: find (u, q, λ) ∈ H 1 0 (Ω) × L 2 (Ω) × H 1 0 (Ω) such that α Ω q r dΩ + Ω λ r dΩ = 0 ∀r ∈ L 2 (Ω), We note here that the trace of λ is zero on the boundary but that this does not affect q in the infinite dimensional case. However, when we discretize and solve for discrete counterparts (u h , λ h , q h ), a problem arises in that the finite dimensionality of the problem forces the zero boundary condition on λ h onto q h via the L 2 -projection in (15). This has profound implications for the accuracy close to the boundary of the control q h . As a remedy for this we propose to introduce the regularization of q in the discrete setting so that it acts only on fluctuations of the gradient of the source term, with a particular scaling in the mesh-size. Since the kernel of this operator is so big, the stability of the system (14)-(16) is compromised. We therefore also introduce a stabilization of the type suggested above, where the jumps of the gradient of u are penalized as well. This is sufficient to make the system inf-sup stable in a suitable discrete norm as we shall see below.

Derivation of the discrete model
Let {T h } h denote a family of shape regular and quasi uniform tesselations of Ω into nonoverlapping simplices, such that for any two different simplices K, K ∈ T h , K ∩ K consists of either the empty set, a common face/edge or a common vertex. The outward pointing normal of a simplex K will be denoted n K . We denote the set of interior element faces F in T h by F I . To each face we associate a normal n F whose orientation is arbitrary but fixed. We define the standard finite element space of continuous piecewise affine functions on T h where P 1 (K) denotes the set of affine functions on K. We also define V 0 h := V h ∩ H 1 0 (Ω).
2.1. Data assimilation. Consider the discrete formulation: h , wheref := f + δ f andũ 0 := u 0 + δu 0 , with δ f ∈ H −1 (Ω) and δu 0 ∈ L 2 (M) denote measurement errors in the source term and data. The bilinear forms are given by ) denotes the jump of the quantity y h over the face F, with normal n F and finally This may then be written on the compact form, find with bilinear forms given by (19), (20) and (21) above. Below we will distinguish the stabilization parameters of s 1 (·, ·) and s 5 (·, ·) and then denote them by γ 1 and γ 5 respectively. This may then be written on the compact form, find

Preliminary technical results
First we define the semi norms associated with the stabilization operator We recall the following well known inverse and trace inequalities As an immediate consequence of (25) we have the following stabilities for some C si > 0 depending only on the mesh geometry, projections on the respective finite element spaces. The following error estimate is known to hold both for i h and π h , To prove stability of our formulations below we need to show that for any We prove this result in this technical Lemma.
There exists c 1 , c 2 > 0 such that for any . Proof. The discrete Poincaré type inequality (28) may be proved using a compactness argument similar to that of [8]. To keep down the technical detail we here instead use an approach with continuous Poincaré inequalities and discrete interpolation.
where ∆ K := ∪ K :K∩K ∅ . The following Poincaré inequality is well known (see [9,Lemma B.63]). If f : H 1 (Ω) → R is a linear functional that is non-zero for constant functions then For instance, we may take As an immediate consequence we have the bound that is the set of interior triangles of M. It then follows by the stability of I h that Adding and subtracting I h ∇v in the second term on the right hand side of (31) gives where we used (30) in the second inequality. For the third term on the right hand side of (32) we once again use Poincaré's inequality, the stability of I h and the inverse inequality (25) to conclude that Using the fact that ∇v h is constant on each element, so that where in the second inequality we have applied the inverse inequality (25) to each term in the sum, and finally used (30). We thus obtain which combined with (32) gives By which we have proven (28). The lower bound of (29) is immediate by the stability of the L 2 -projection and the inverse and trace inequalities of equation (25). To prove the upper bound we write from which the upper bound follows, since by the triangle inequality followed by the first inequality of (26), with i = 3, To prove (33) we first define the support of a nodal basisfunctions ϕ i by Ω i := {x ∈ Ω : ϕ i (x) > 0}. Then let N I denote the set of indices of basis functions that are in the interior of the domain, that is, for each value i ∈ N I the closure of the support of the associated basis function has empty intersection with the boundary, Ω i ∩ ∂Ω = ∅. For each i ∈ N I we define the macropatch ∆ i := ∪ j:Ω j ∩Ω i ∅ Ω j . This means that ∆ i consist of Ω i and any other patch Ω j sharing two triangles (in 2D) or several tetrahedra (in 3D) with Ω i . Since T h is shape regular we may map the patch ∆ i to a shape regular∆ i such that diam(∆ i ) = 1. We define the linear map F : ∆ i →∆ i . LetF andÑ denote the set of interior faces and interior nodes respectively of∆ i and define the scalar product on∆ i It follows that ifw h denotes the mapped function w h andφ j denotes the mapped basis function ϕ j then We now prove that | · | s,∆ is a norm onw h |∆ i . It is clearly a semi norm so we only need to prove that |w h | s,∆ = 0 impliesw h |∆ i = 0. If |w h | s,∆ = 0 thenw h is an affine function on∆ i . It is straightforward to check that the only affine function that can satisfy (34) is the zero function. To be precise, otherwisẽ w h has to be odd with respect to the center of mass of all the basis functions.
Since w h is affine it will vanish along a line in 2D and on a plane in 3D. This means that the centers of mass of all the nodal basis functions associated to the nodes inÑ must be on the line (in the plane), which is impossible. The constant of the estimate depends on the shape regularity. Defining where F i denotes the set of interior faces of ∆ i , we obtain, by scaling back to the physical geometry, that there exists C > 0 depending only on the local mesh geometry such that, To conclude we observe that since the overlap between different patches ∆ i is bounded uniformly in h there holds which is the desired inequality.

Error analysisdata assimilation
We introduce the triple norm (27) and (25) the following approximation estimate is straightforward to show Observe that the terms in the above norm do not have matching dimensions. Indeed there is a constant of the dimension of an inverse length scale present in the first two terms in the right hand side. In the term over L 2 (M) this is to avoid a too strong penalty on possibly perturbed data and in the second term of the right hand side it comes from the application of the discrete Poincaré inequality (28). Stability in the norm (35) is sufficient to deduce the existence of a discrete solution to the system (17)-(18), however the norm is too weak to be useful for error estimates.
The analysis takes the following form, following the framework of [3]. First we prove inf-sup stability of the form A DA [·, ·] in the norm (35). From this the existence of discrete solution to the linear system follows. Then we show an error estimate in the norm (35) that is independent of the stability of the data assimilation problem and gives convergence rates for the residuals of the approximation. Finally we show that the error satisfies an equation of the type (8), with the right hand side given by the residual. The a priori error estimates on the residual together with the assumed a priori estimate on the exact solution allows us to deduce error bounds through Lemma 1.
h be the solution of (17)-(18) then there exists c s > 0 such that Proof. First observe that Using the Cauchy-Schwarz inequality, an arithmetic-geometric inequality and the inverse inequality (25) in the last term of the right hand side we get . and similarly . The contribution hu h H 1 (Ω) is added to the right hand side by applying (28). Using once again the Cauchy-Schwarz inequality, an arithmetic-geometric inequality and the inverse inequality (25) we see that which, together with the Poincaré inequality for λ h , concludes the proof. Proof. Let ξ h := u h − i h u, with i h u denoting the Lagrange interpolant of u. By the triangle inequality we have For the second term on the right hand side we apply the inf-sup condition of Proposition 1, we obtain the equality Using the Cauchy-Schwarz inequality in the terms of the right hand side we immediately deduce Applying this inequality in (37) obtain the bound and the result follows from approximation estimate (27) and (36).
h be the solution of (17)-(18) and u ∈ H 2 (Ω) the solution to (7)-(8), with α = 0. Then for some 0 < τ < 1 depending on r 1 /r 3 , r 2 /r 3 and the smallest eigenvalue there holds It follows that e is a weak solution to the problem (8) with the right hand side r ∈ H −1 (Ω) and we are in the framework of Lemma 1. We now only need to show the bound r H −1 (Ω) ≤ ε and then apply the inequality (9) with e in the place of u. By definition of the dual norm we have We proceed using the definition of r, Here we used partial integration in the form a(·, ·), a trace inequality and approximation to obtain We now apply Lemma 1 to obtain Setting ε = C r (h f L 2 (Ω) + |u h | s 1 + δ f H −1 (Ω) ) and observing that by Proposition 2 and once again by Proposition 2 we conclude.
Corollary 1. Assume that for h 0 > 0 there holds Proof. This result follows immediately by applying the assumed bound (38) in the error estimate of Theorem 1.
Observe that in particular the above result implies that if the exact data is available in M we can compute the solution in B r 2 (x 0 ) to arbitrary precision. Also observe that Theorem 1 provides both a priori and a posteriori error bounds. This means that perturbations in data can be compared with the computational residual in an a posteriori procedure to drive adaptive algorithms for the computation of the reconstruction.

Error analysissource reconstruction
The error analysis in this case follows a similar outline, however instead of Lemma 1 we may here use a compactness argument to obtain convergence orders.
We will also use the Ritz-projection defined by h . It is well known that if Ω is convex the following estimate holds u − r h u L 2 (Ω) + h ∇(u − r h u) L 2 (Ω) ≤ Ch 2 |u| H 2 (Ω) .
We will first prove an estimate where we assume that q is more regular and show that in this case the stabilization of the velocity is superfluous.
For the estimate on the source term observe that for γ 5 ≥ 0 If on the other hand γ 5 > 0 then we may use where it is immediate to show that q − π h q H −1 (Ω) ≤ Ch 2 ∇q L 2 (Ω) and . For the second term on the right hand side the estimate (44) holds and for the first term we observe that with

We may then use the equation to deduce
and after a Cauchy-Schwarz inequality and an inverse inequality in the second factor, . It follows that h π 0 h η h L 2 (Ω) ≤ Ch ∇(u − u h ) L 2 (Ω) and the above bounds that (43) holds.
We see that if less regularity is assumed for the source term q ∈ L 2 (Ω) then convergence of the gradient can no longer be deduced, however a priori bounds on u h and q h are nevertheless achieved that may be used to prove convergence in the asymptotic limit. In order to obtain an estimate with convergence order in h also in the case where q ∈ L 2 (Ω) we take γ 1 > 0 and prove that this allows us to obtain stronger control of the approximation of the source term. This in its turn allows us to prove convergence using an interpolation argument between H −2 and L 2 as we shall see below.
Proof. For some β > 0 to be fixed, . Observing that Cauchy-Schwarz inequalities, arithmetic-geometric inequalities and the stability inequality (26), with i = 5 leads to Using partial integration, the fact that π 0 h t h | ∂Ω = 0 and a trace inequality (25) we have the bound We may then fix β = min(C −2 si , C −2 t ) to show that for some c > 0 depending only on the mesh geometry By Lemma 2 and the quasi uniformity of the mesh there holds for some c > 0 depending only on the mesh geometry and it follows that for some c > 0 depending only on the mesh geometry To conclude we need to show that To this end we note that For the second term in the right hand side we may write where the constant only depends on the constant of quasiuniformity of the meshes.
Proposition 5. Let (u, q) ∈ W × L 2 (Ω) satisfy (13) and let (u h , q h λ h ) ∈ V S R h be the solution of (22)-(23). Then there holds and Proof. Considering (40) it is sufficient to prove the estimate for the discrete errors ξ h := u h − i h u and η h := q h − π h q. From Proposition 4 we know Using the definition of (13), (22)-(23) and the orthogonality of the L 2projection we may write We see that using a Cauchy-Schwarz inequality, the approximation properties of the Lagrange-interpolant and the regularity of u we have Similarly for term three we use Cauchy-Schwarz inequality and the approximation result (40) to obtain In term IV we use an integration by parts followed by the stability (26) 1 , with i = 5, to obtain Finally for term V we proceed using Cauchy-Schwarz inequality followed by an inverse inequality to obtain Collecting the bounds for the terms I−V above we deduce that the following bound holds which proves the first claim. The second claim follows immediately by the definition of the triple norm, an inverse inequality and the first inequality.
Now consider the problem: find z ∈ H 1 0 (Ω) such that −∆z = q − q h in the sense of distributions.
It follows that q − q h H −1 (Ω) ≤ ∇z L 2 (Ω) = (∇z, ∇z) Considering the last term in the right hand side we see that Taking the square root of the last expression we obtain the desired estimate for q − q h H −1 (Ω) . We now consider the estimate on the error in the gradient of u. Since u, u h ∈ H 1 0 (Ω) there holds For the first and second terms of the right hand side we use the first estimate of Proposition 5 to obtain ). For the third term we observe that by the definition (45), the bound (46) and the regularity assumption on u we may write . Collecting the bounds of the terms I-IV and taking square roots concludes the proof. 5.1. The effect of measurement error. Both the above estimates take into account the measurement errors. We will here discuss the second case in some more detail. Observe that these estimates are nonstandard in the sense that they measure the distance of the approximate solution to any pair of functions (u, q) that satisfy (13). First assume that u 0 ∈ W then u 0 satisfies the constraint and u h will converge to u 0 and q h to −∆u 0 , according to the estimates above. Otherwise under our assumptions u 0 is a measurement on the form u 0 =ũ 0 + δu 0 and in this case the solution that we wish to approximate is u =ũ 0 and q = −∆ũ 0 . Provided δu 0 L 2 (Ω) is a priori known the estimates above give us a bound on how close u h , q h are toũ 0 and −∆ũ 0 . The a perturbation term δu 0 L 2 (Ω) does not vanish under mesh-refinement and that causes the estimates to blow up if u 0 is not in the range of the Laplace operator.
The above results also gives us quantitative information on how the perturbation δu 0 in the measurements affects the computation. For general measurement errors we see that formally refinement improves the solution as long as δu 0 L 2 (Ω) q L 2 (Ω) <<Ch.
Since q is unknown it has to be replaced by q h in practice, which is be reasonable as long as q h L 2 (Ω) does not grow under mesh refinement. If we assume that u is sampled on some coarse scale H and u 0 is taken to be a piecewise affine interpolant using these samples. Then u − u 0 L 2 (Ω) ≤ CH 2 u H 2 (Ω) and the above estimate becomes We see that we have two terms with different behavior as we refine in h. To ensure that the solution improves under mesh refinement we can check that the discretization error dominates the measurement error. This is true as long as h Another a posteriori criterion for when refining the mesh improves the solution is given by: (1) ∇u h L 2 (Ω) and q h L 2 (Ω) non-increasing under refinement (2) s 1 (u h , u h ), s 5 (q h , q h ) decreasing under refinement. 6. Numerical examples 6.1. Data assimilation. We consider problem 7 on the domain Ω = (−1, 1)× (−1, 1) and the data with right-hand side f = −∆u 0 , and M = (−1/4, 1/4) × (−1/4, 1/4). We set γ = 10 −4 . The first mesh on Ω, in a sequence, and the mesh on M is given in Fig. 1. The different domains B r 2 on which convergence is measured is shown in Fig. 2, and the obtained convergence is shown in Fig. 3. Note that the error constant increases with r 2 and that the convergence rate decreases slowly, cf.   Table 1. Errors and convergence rates underlying Figure 2; meshsize is defined as the inverse square root of the number of nodes.
In Fig. 4 we show the interpolant of the exact source and a typical solution obtained using the gradient jump stabilization method. In Fig. 5 we show the effact of (properly scaled) Tikhonov regularization using L 2 , i.e., q , and H 1 , i.e., ∇q , regularizations, respectively. Note that the L 2 regularization gives the wrong boundary conditions in the discrete scheme, whereas H 1 works better while still giving a spurious boundary effect, well known from similar approaches used in fluid mechanics, cf. Burman and Hansbo [7].
The observed convergence using (22)-(23) using only the stabilization term s 5 (q h , w h ) and for a variety of choices of γ is shown in Fig. 6. We note that the convergence q h − q L 2 (Ω) and u h − u H 1 (Ω) is first order in both cases. The convergence of u h is completely unaffected by the choice of γ.
For the non-smooth case we let the solution be two different constant in the radial direction, u = 1 for 0 ≤ r ≤ 1/4 and u = 0 for 3/4 ≤ r, interconnected by a cubic C 1 -polynomial i the radial direction. This means that the source term will have jumps at r = 1/4 and r = 3/4 so that q ∈ H 1/2− (Ω) for any > 0. In Fig. 7 we show the observed rate of convergence, which drops to about O(h 1/2 ) for q h − q L 2 (Ω) but remains O(h) for u h − u H 1 (Ω) . The error constant is now affected by γ also for the convergence of u h . Finally, we show the effect of perturbing the data randomly, with constant amplitude and with amplitude decreasing O(h). In Fig. 8  Consider Ω = (0, 1) × (0, 1) and the right hand side defined as a discontinuous cross shaped function (see Figure 11, left plot) written using boolean binary functions as The data u 0 is reconstructed using P 4 finite elements on the one hand on a mesh that is fitted to the discontinuities of f (120 × 120 structured) resulting in a very accurate solution ( u − u h L 2 (Ω) ≤ C(120) −5 f L 2 (Ω) ) and on the other hand on a mesh that is not fitted to the discontinuities of f (110 × 110 elements). The unfitted data results in spurious high frequency oscillations with small amplitude in the high order finite element solution, as can be seem in Figure 9 (left fitted data and right unfitted data). The L 2 -norm of the difference of the fitted and the unfitted solution is a good measure of the size of the perturbation. It is 1.7 × 10 −4 .
First we fixed γ = 10 −6 after a few steps of a line search algorithm, using the same stabilization parameter for s 1 and s 5 . We solved the problem using 6 unstructured (Delaunay) meshes with 20, 30, 40, 60, 80 and 100 elements on the domain side. The L 2 -error in q is given in Figure 10. Circle markers indicate the result obtained with the stabilized method and square markers the result obtained taking γ = 0 above. In the left plot the data u 0 is given by the accurate computation and in the right plot the perturbed data is used. As can be seen in the plots, for the unperturbed data the stabilized method performs slightly better than the unstabilized method and has approximately h 1 2 order convergence in the L 2 -norm, which is optimal. For the perturbed data on the other hand the situation is dramatically different, whereas the stabilized method almost has the same convergence, on coarse meshes and only stagnates on finer meshes when the effect of the perturbation becomes important, the unstabilized method diverges.
The exact and reconstructed source function for the case of fitted data (on the 80 × 80 mesh) is given in Figure 11 and with unfitted data in Figure 12.              Figure 11. Contout lines of exact and reconstructed source terms. Left: exact source term. Middle: reconstructed source term using stabilization, unperturbed data. Right: unstabilized reconstruction, unperturbed data. Figure 12. Left: reconstructed source term using stabilization, perturbed data. Right: unstabilized reconstruction, perturbed data.