Quasi-solution of linear inverse problems in non-reflexive Banach spaces
Christian Clason Andrej Klassen






Quasi-solution of linear inverse problems in non-reflexive Banach spaces
    Christian Clason
        University Duisburg-Essen, Faculty of Mathematics,
        Thea-Leymann-Str. 9,
        45127 Essen, Germany
        (christian.clason@uni-due.de, andrej.klassen@uni-due.de)
     
    Andrej Klassen[1]

    June 26, 2018
=======================================================================================================================================================================================================================================




    We consider the method of quasi-solutions (also referred to as Ivanov regularization) for the regularization of linear ill-posed problems in non-reflexive Banach spaces. Using the equivalence to a metric projection onto the image of the forward operator, it is possible to show regularization properties and to characterize parameter choice rules that lead to a convergent regularization method, which includes the Morozov discrepancy principle. Convergence rates in a suitably chosen Bregman distance can be obtained as well.
    We also address the numerical computation of quasi-solutions to inverse source problems for partial differential equations in L^∞(Ω) using a semismooth Newton method and a backtracking line search for the parameter choice according to the discrepancy principle.
    Numerical examples illustrate the behavior of quasi-solutions in this setting.




§ INTRODUCTION


We consider the operator equation Au=y for a linear injective operator A:U → Y with dense range, where Y is a uniformly convex Banach space and U the – not necessarily reflexive – dual of a separable normed vector space. Instead of the exact data y∈ R(A) (the range of A), usually only a noisy measurement y^δ∈ Y∖ R(A) with y^δ-y_Y≤δ for some noise level δ is available, and regularization methods have to be applied to obtain a stable (approximate) solution.
Three related strategies that can be found in the literature are:

    
  * Tikhonov regularization <cit.>, consisting in solving for some α>0 the minimization problem
        
    Tmin_u∈ U1/2Au-y_Y^2 + α/2u_U^2;

    
  * Morozov regularization (also known as the method of residuals) <cit.>, consisting in solving the constrained minimization problem
        
    Mmin_u∈ Uu_U   subject to  Au-y_Y≤δ;

    
  * Ivanov regularization (also known as the method quasi-solutions) <cit.>, consisting in solving for some ρ>0 the constrained minimization problem
        
    Imin_u∈ UAu-y_Y   subject to  u_U≤ρ.


It is the last one that is the focus of this work. Let us first briefly comment on the relation of these approaches. Since A is linear, all three minimization problems are convex, and hence their solutions coincide for a suitable choice of α and ρ, respectively; see, e.g., <cit.>. However, this is no longer the case for nonlinear operators, as was shown by a counterexample in <cit.> (which also, to the best of our knowledge, coined the term Ivanov regularization); see also <cit.>. Furthermore, if U is not a Hilbert space, (<ref>) is no longer necessarily easier to solve than (<ref>) or (<ref>). In fact, for U=L^∞(Ω) (which motivates the problem setting stated in the beginning), <cit.> proposed solving (<ref>) by transforming it into the form (<ref>) with an additional penalty α/2ρ^2 for (the in this case free variable) ρ. This motivates directly considering (<ref>) as a regularization strategy.

A further motivation is parameter identification in partial differential equations, where the domain of the (nonlinear) forward mapping is given by pointwise almost everywhere restrictions on the parameter. If such restrictions already yield a regularization of the ill-posed parameter identification problem, unwanted smoothing from the introduction of a Tikhonov penalty can be avoided. 

In this work, we consider as a first step admissible sets of the form

    M_ρ := u ∈ Uu_U≤ρ,

where we are particularly interested in the case U= L^∞(Ω) for a bounded domain Ω⊂^d.
We then take as the quasi-solution of Au=y the minimizer u_ρ∈ M_ρ⊂ U of 

    inf_u ∈ M_ρAu - y_Y,

where the radius ρ≥ 0 plays the role of regularization parameter. Its proper choice is one of the main issues in showing that the method of quasi-solutions can be considered a (convergent) regularization method, and the characterization of such choice rules – especially of a posteriori-type – is one of the main contributions of this work. The main difficulty here is the fact that U is a non-reflexive Banach space, which requires working in the weak* topology. This is addressed by considering (<ref>) as a projection problem in Y, allowing exploiting the strong topology of this space.
We also discuss the numerical solution of (<ref>) for the case that A is the solution operator to an elliptic partial differential equation and U= L^∞(Ω).

Let us remark on related literature. Quasi-solutions were introduced by Ivanov <cit.> as the stabilization of linear ill-posed problems in normed spaces through a compact restriction of the solution domain based on Tikhonov's theorem <cit.>. 
Further works <cit.> weakened the requirement to weakly compact norm balls in uniformly convex Banach spaces; see also <cit.>. More general compact solution domains of monotone or convex and bounded functions were considered in <cit.>.
More recently, convergence rates for quasi-solutions to nonlinear inverse problems with norm balls in Hilbert scales were proved in <cit.>. 
However, the (a posteriori) parameter choice for the radius was not discussed. 
Finally, in parallel to this work, convergence and convergence rates for both Morozov and Ivanov regularization under variational source conditions were derived in a more abstract setting in <cit.>.
In contrast to this, the current work also treats the numerical computation of quasi-solutions in a non-reflexive Banach space, which to the best of our knowledge has not been done so far. 

This paper is organized as follows. <Ref> collects results on metric projections onto image sets under A that will be needed in the following. In <ref>, we show that the method of quasi-solutions is indeed a regularization strategy, characterize parameter choice rules leading to a convergent regularization method, and derive convergence rates in a suitably chosen Bregman distance. <Ref> is concerned with the numerical solution of an inverse source problem for a partial differential equation with U=L^∞(Ω) based on a semi-smooth Newton method and a line search for the a posteriori choice according to the discrepancy principle. Finally, <ref> contains numerical examples illustrating the behavior of quasi-solutions in this setting.




§ QUASI-SOLUTION AS METRIC PROJECTION IN IMAGE SPACE


To avoid having to work with the weak or weak* topology in U, we consider the problem of finding a quasi-solution as a metric projection onto the range of A. While existence of a quasi-solution can also be shown by standard methods, the analysis in image space will be helpful for analyzing parameter choice rules.

For this purpose, we make the following assumption on the forward operator.

    The linear operator A:U→ Y 
    
        
  * is injective,
        
  * has dense range,
        
  * is bounded in the strong topologies of U and Y,
        
  * is weak*-to-strong closed, i.e., for any sequence {u_n}_n∈ with u_n  u in U and Au_n  y in Y we have that Au = y.
    

From the definition of M_ρ and these assumptions on A, we immediately obtain the following useful property.

    For every ρ≥ 0, the image set Q_ρ := A(M_ρ)⊂ Y is non-empty, closed, convex, and bounded.

 
    Non-emptiness, convexity, and boundedness are obvious since M_ρ is a ball and A is linear and bounded. To see the closedness, let {y_n}_n∈ be a sequence in Q_ρ converging strongly to some y∈ Y. By definition of Q_ρ, there thus exists for every n∈ a u_n∈ M_ρ such that Au_n = y_n.
    Since {u_n}_n∈⊂ M_ρ is bounded and U is the dual of a separable vector space, we can apply the Banach–Alaoglu theorem to extract a subsequence, still denoted by {u_n}_n∈, with u_n  u∈ M_ρ (since the unit ball in U is weakly* sequentially closed). By the closedness <ref> (iv) on A, we thus have y = Au ∈ Q_ρ.

Since Y was assumed to be a uniformly convex Banach space, we can apply the following result from <cit.>.

    For every y∈ Y, there exists a unique metric projection y_ρ∈ Q_ρ such that
    
    y_ρ -y_Y = inf_q∈ Q_ρq-y_Y.
 
    Furthermore, the mapping P_Q_ρ:y↦ y_ρ is continuous.

By the definition of Q_ρ, this yields the existence of a u_ρ∈ M_ρ with y_ρ=Au_ρ that attains the minimum in (<ref>). Together with the injectivity of A, we thus immediately obtain the existence of a unique quasi-solution.
 
    For every y∈ Y and ρ≥ 0 there exists a unique solution u_ρ∈ M_ρ to (<ref>).


We can therefore introduce the distance function

    d:[0,∞) × Y → [0,∞)      (ρ,y) ↦Au_ρ - y_Y = min_u ∈ M_ρAu-y_Y.

The following property of the distance function will be crucial in the following. It was first proved in <cit.> for the case that M_ρ is compact, and our extension to weakly* compact and convex subsets follows the line of their proof. (Here we point out that convex and compact sets are in general not weakly* compact.)
 
    For every y∈ Y, the mapping ρ↦ d(ρ,y) is surjective on 
    (0,y_Y). 


    Let y∈ Y and σ∈ (0,y_Y) be given. Since the range of A is dense in Y, there exists a sufficiently large ρ such that Q_ρ∩ B_σ(y) is non-empty. The set of all such ρ>0 is thus non-empty and bounded from below and therefore admits an infimum
    
    ρ̃ := inf{ρ>0: Q_ρ∩ B_σ(y) ≠∅}.

    By the definition of the infimum, there thus exists a sequence {ρ_n}_n∈ℕ⊂ (0,∞) with
    
    S_n := Q_ρ_n∩ B_σ(y) ≠∅

    and ρ_n→ρ̃. Without loss of generality, we can assume that this sequence is monotonically decreasing, which implies that Q_ρ_n+1⊆ Q_ρ_n and hence that S_n+1⊆ S_n. 

    We now show that the intersection S:=⋂_n∈ S_n is non-empty. First, since every S_n⊂ Y is bounded, convex, and closed (as the intersection of such sets) as well as non-empty, we can apply <ref> to obtain for each n∈ an s_n := P_S_n(0). 
    Furthermore, the monotonicity of {ρ_n}_n∈ implies that  s_n_Y ≤ρ_1. Hence, {s_n}_n∈ is bounded and, since uniformly convex Banach spaces are reflexive, thus contains a subsequence weakly converging to some s∈ Y.
    As S is the intersection of convex and closed sets and therefore weakly closed, we obtain that s∈ S.

    We next show that s∈ Q_ρ̃. First, by construction we have that s∈ S_n⊂ Q_ρ_n for every n∈, i.e., there exists a u∈ M_ρ_n such that Au=s. (Since A is injective, u cannot depend on n.) We can thus pass to the limit to obtain that
    
    u_U ≤ρ_n →ρ̃.

    In fact, we even have that u_U = ρ̃, because if the inequality were strict, we could find a ρ̂<ρ̃ such that s∈ Q_ρ̂∩ B_σ(y), in contradiction to the definition of ρ̃. This shows the claim.

    Similarly, since s∈ S_n for all n∈, we have that s∈ B_σ(y). Assume now that s-y_Y<σ, i.e., s is an interior point of B_σ(y). Then there exists for h:=-s an >0 small enough such that s̅:= s+ h=(1-)s∈ B_σ(y).    
    Since A is linear, we have that 0∈ Q_ρ̃ and hence by <ref> also s̅ =  0 + (1-)s ∈ Q_ρ̃. 
    (Note that the assumption 0<σ<y_Y implies that 0∉ B_σ(y) and thus in particular s≠ 0.) 
    This implies that there exists a u̅∈ M_ρ̃ such that Au̅ = s̅. 
    But from the linearity and injectivity of A it then follows that u̅ = (1-)u and hence that
    
    u̅_U = (1-)u_U = (1-)ρ̃< ρ̃,

    again in contradiction to the definition of ρ̃. This implies that s-y_Y=σ.

    It remains to show that s-y_Y = inf_u∈ M_ρ̃Au-y_Y. For this, we argue as in the last step: If there exists a û∈ M_ρ̃ such that Aû-y_Y<σ, we can again find a u̅∈ M_ρ̃ such that Au̅∈ B_σ(y) but u̅_U < ρ̃ to contradict the definition of ρ̃. Hence,
    
    d(ρ̃,y) = s-y_Y = σ

    as claimed.
 

The proof of <ref> shows in particular that unless d(ρ,y)=0 – i.e., y∈ Q_ρ⊂ R(A) – the unique quasi-solution u_ρ∈ M_ρ from <ref> always lies on the boundary of M_ρ, a result that was already shown in <cit.> in the case that U and Y are Hilbert spaces.

    For every ρ≥ 0 and y∈ Y∖ Q_ρ, the solution to (<ref>) satisfies u_ρ_U = ρ.

This result also implies a strict monotonicity of the distance function.

    For every y∈ Y∖ R(A), the mapping ρ↦ d(ρ,y) is strictly monotonically decreasing and continuous.


    If ρ_1 ≤ρ_2, the definition of the distance function and the fact that M_ρ_1⊂ M_ρ_2 immediately imply that
    
    d(ρ_2,y)=inf_u ∈ M_ρ_2Au-y_Y ≤inf_u ∈ M_ρ_1Au-y_Y = d(ρ_1,y).

    Assume now that ρ_1< ρ_2 but d(ρ_1,y)=d(ρ_2,y). By <ref>, we have that u_ρ_1_U=ρ_1 < ρ_2 = u_ρ_2_U and in particular u_ρ_1≠ u_ρ_2. 
    But by injectivity of A, we have that Au_ρ_1≠ Au_ρ_2 as well, and thus the uniform convexity of Y implies that u̅ := 1/2(u_ρ_1+u_ρ_2)∈ M_ρ_2 satisfies
    
    Au̅ -y = 12(Au_ρ_1 + Au_ρ_2)-y_Y < 1/2(Au_ρ_1 -y_Y +Au_ρ_2 -y_Y ) = d(ρ_2,y),

    in contradiction to the definition of the distance function.

    Together with the surjectivity from <ref>, the strict monotonicity implies the continuity.


Finally, we will also require the continuity of the distance function with respect to its second argument.

    For each ρ≥ 0, the mapping y↦ d(ρ,y) is non-expansive (i.e., Lipschitz continuous with constant 1).

 
    Consider y_1,y_2 ∈ Y and corresponding projections q_1 := P_Q_ρy_1 and q_2 := P_Q_ρy_2. The optimality property of the projections then imply that
    
    d(ρ,y_1)    = q_1 - y_1_Y ≤q_2 - y_1_Y ≤q_2 - y_2_Y + y_2 - y_1_Y, 
    
            d(ρ,y_2)    = q_2 - y_2_Y ≤q_1-y_2_Y ≤q_1 - y_1_Y + y_1 - y_2_Y,

    which together yield that 
    
    |d(ρ,y_1) - d(ρ,y_2)| ≤y_1 - y_2_Y.

     




§ REGULARIZATION AND PARAMETER CHOICE


In this section, we show that the method of quasi-solutions is a regularization in the classical sense <cit.> and leads – combined with an appropriate parameter choice rule for ρ – to a convergent regularization method. Since in this case the solution domain M_ρ is not compact and U is not assumed to be reflexive, we can only expect weak* convergence.



 §.§ Stability and convergence


To show that quasi-solutions are a regularization of the ill-posed operator equation Au=y, we first need to show that for each ρ≥ 0, the solution to (<ref>) depends continuously on y∈ Y. We again exploit the equivalence of (<ref>) with the metric projection onto Q_ρ:= A(M_ρ).

    Let ρ≥ 0 and y∈ Y be given. Then for every sequence {y_n}_n∈⊂ Y with y_n→ y, the corresponding sequence of solutions {u_n}_n∈ to (<ref>) for y_n in place of y converges weakly* to u_ρ.

 
    First, <ref> yields for every n∈ a unique metric projection q_n∈ Q_ρ and hence, by definition of Q_ρ, a u_n∈ M_ρ such that Au_n = q_n. Furthermore, by continuity of the metric projection, we have that q_n→ q:= P_Q_ρ y. As in the proof of <ref>, we can now extract a subsequence, still denoted by {u_n}_n∈, such that u_n  u∈ M_ρ with Au=q, i.e., u∈ M_ρ minimizes Au-y_Y over M_ρ. This shows that u=u_ρ is the solution to (<ref>). Since this solution is unique by <ref>, a subsequence–subsequence argument yields weak* convergence of the full sequence.


In the case of Ivanov regularization, convergence as ρ→∞ for fixed y∈ R(A) is obvious since for every ρ≥u^†_U with Au^† = y, the choice u_ρ = u^†∈ M_ρ clearly solves (<ref>). Nevertheless, we state the consequence as a theorem for the sake of completeness.

    Let y∈ R(A) be given and u^†∈ U such that Au^† = y. Then for every sequence {ρ_n}_n∈ with ρ_n→∞, we have u_ρ_n→ u^†.




 §.§ Parameter choice


We now address the convergence of Ivanov regularization in combination with a parameter choice strategies for ρ. As we will see, there are fundamental differences to Tikhonov (and Morozov) regularization, and hence the results of this section are one of the main contributions for this work.

We begin by characterizing parameter choices that lead to a convergent regularization method. In the following, let y∈ R(A) be arbitrary, {δ_n}_n∈ be a non-negative sequence with δ_n → 0, and y^δ_n∈ Y with y^δ_n-y_Y≤δ_n for each n∈. Let further u^†∈ U denote the (unique) solution to Au=y and for given ρ_n≥ 0, let u_ρ_n^δ_n∈ U denote the solution to 

    min_u∈ M_ρ_nAu-y^δ_n_Y.


    If A is not continuously invertible, the sequence {u_ρ_n^δ_n}_n∈ weakly* converges to u^† as n→∞ if and only if the following two conditions hold:
    
        
  * {ρ_n}_n∈ is bounded,
        
  * lim inf_n→∞ρ_n ≥u^†_U.
    

 
    We first show that these conditions are sufficient by proceeding similarly to the proof of <ref>. 
    First, condition (i) implies that {u_ρ_n^δ_n}_n∈ is bounded, and hence we can extract a subsequence u_n u for some u∈ U. By passing to a further subsequence, we can in addition assume that ρ_n →ρ̂≥u^†_U by condition (ii). This implies that
    
    P_Q_ρ_ny-y_Y = d(ρ_n,y) → d(ρ̂,y) = 0,

    where we have used the continuity of the distance function from <ref> and the fact that Au^† = y∈ Q_ρ̂. Furthermore, for every such n, <ref> yields the existence of a unique metric projection P_Q_ρ_n y^δ_n = Au_ρ_n^δ_n which by definition satisfies
    
    P_Q_ρ_ny^δ_n - y^δ_n_Y = inf_q ∈ Q_ρ_nq - y^δ_n_Y ≤P_Q_ρ_ny - y^δ_n_Y.

    Combining the above, we obtain that
    
    P_Q_ρ_ny^δ_n - y_Y    ≤P_Q_ρ_ny^δ_n - y^δ_n_Y + y^δ_n - y_Y
       ≤P_Q_ρ_ny - y^δ_n_Y + δ_n 
       ≤P_Q_ρ_ny - y_Y + 2δ_n.

    Passing to the limit then yields that P_Q_ρ_ny^δ_n→ y and hence by the closedness that Au=y, i.e., u=u^†. Convergence of the full sequence now follows again from a subsequence–subsequence argument.

    To show the necessity, first note that by assumption R(A) is not closed and hence has no interior points (otherwise the algebraic interior of R(A) would not be empty either, and the linearity of A would lead to the contradiction R(A)=Y.) For any n∈, we can thus find a y^δ_n∈ B_δ_n(y) ∖ R(A). <Ref> then implies that the corresponding quasi-solution satisfies u_ρ_n^δ_n_U = ρ_n. Hence, if {ρ_n}_n∈ is unbounded, {u_ρ_n^δ_n}_n∈ is unbounded as well and therefore cannot weakly* converge (since weakly* convergent sequences are bounded by the Banach–Steinhaus Theorem). On the other hand, if u_ρ_n^δ_n u^†, then the weak* lower semi-continuity of the norm immediately implies that
    
    u^†_U ≤lim inf_n→∞u_ρ_n^δ_n_U = lim inf_n→∞ρ_n.

 

<Ref> shows in particular that Ivanov regularization converges under the constant parameter choice ρ_n ≡u^†_U.
However, this is not a valid parameter choice in the strict sense, which may only depend on δ and y^δ, but not on the exact solution u^†. In fact, <ref> implies that there can be no convergent a priori parameter choice ρ = ρ(δ) unless the inverse problem is well-posed.

    Let ρ(δ) be an a priori parameter choice. If u_ρ_n^δ_n u^† for ρ_n = ρ(δ_n) with δ_n→ 0 as n→∞, then A is continuously invertible.


    Since the claim must hold independently of u^†, condition (ii) of <ref> has to be satisfied for any u^†∈ U. This implies a fortiori that {ρ_n}_n∈ is unbounded, in contradiction to condition (i). Hence, A must be continuously invertible.

Of course, under the additional a priori information u^†∈ M_ρ^† for some ρ^†≥ 0, the choice ρ_n≡ρ^† leads to a convergent regularization method (cf. <cit.>).

This leaves a posteriori choice rules ρ=ρ(δ,y^δ). Here we consider the Morozov discrepancy principle <cit.>, which for some τ >1 chooses ρ(δ,y^δ) such that

    δ≤Au_ρ(δ,y^δ)^δ - y^δ_Y ≤τδ.

This choice indeed leads to a convergent method.
    
    Let ρ_n:=ρ(δ_n,y^δ_n) be chosen according to (<ref>) for some τ>1. Then u_ρ_n^δ_n u^† as n→∞.


    We make use of the properties of the distance function (<ref>) from <ref>. First note that for any 0<δ_n<y^δ_n_Y (which holds for n sufficiently large provided y_Y>0), <ref> yields the existence of a ρ_n satisfying (<ref>). 
    It remains to verify the conditions of <ref>.

    We first show that {ρ_n}_n∈ is bounded by u^†_U. Assume that there exists an n∈ such that ρ_n > u^†_U, which implies that u^†∈ M_ρ_n.
    Since u_ρ_n^δ_n_U = ρ_n for δ_n>0 by <ref>, we further have that u_ρ_n^δ_n≠ u^†. Hence it follows from the definition  and uniqueness of the quasi-solution that
    
    Au_ρ(δ,y^δ)^δ - y^δ_Y = min_u∈ M_ρ_nAu-y^δ_n_Y < Au^† -y^δ_n_Y ≤δ_n,

    in contradiction to the choice of ρ_n.

    Assume now that
    
    lim inf_n →∞ρ_n < u^†_U=:ρ^†.

    Since the inequality is strict, there exists – after passing to a subsequence – for any >0 an N ∈ such that ρ_n < ρ̃:=u^†_U - for all n > N. From <ref> we thus obtain that d(ρ_n,y^δ_n) > d(ρ̃,y^δ_n) for all n>N.
    The parameter choice rule (<ref>) and <ref> then imply that
    
    lim_n →∞τδ_n ≥lim_n →∞ d(ρ_n,y^δ_n) ≥lim_n →∞ d(ρ̃,y^δ_n) = d(ρ̃,y) > 0

    since u^†∉ M_ρ̃, in contradiction to δ_n→ 0.

The proof of <ref> in fact shows that for ρ_n:=ρ(δ_n,y^δ_n) chosen according to the discrepancy, we have ρ_n→u^†_U from below.




 §.§ Convergence rates


Since we have only shown weak* convergence of u_ρ^δ to u^† as δ→ 0, we cannot expect convergence rates for the strong error u_ρ^δ-u^†_U. As usual for inverse problems in Banach spaces, we thus consider rates for the error measured by the Bregman distance; see, e.g., <cit.>. We recall that for a convex functional J:U→∪{∞}, the convex subdifferential of J at u∈ U with J(u)<∞ is given by

    ∂ J(u) = ξ∈ U^*ξ,v-u_U ≤ J(v)-J(u)  for all v∈ U.

Note that the subdifferential can be empty unless J(u)<∞ for all u∈ U. For given u∈ U and ξ∈∂ J(u), we can now define the Bregman distance

    D_J^ξ(v;u) = J(v)-J(u) - ξ,v-u_U    for all v∈ U.

By the definitions, we have D_J^ξ(v;u) ≥0 for any v∈ U as well D_J^ξ(u; u)=0, which justifies using the Bregman distance to measure the error. Here, we choose J(u) := u_U and point out that this choice is different from the classical one as suggested in <cit.>; to fit into their framework, we would have to take J(u)=δ_M_ρ(u) with δ_M_ρ(u) = 0 if u∈ M_ρ and ∞ else, which however would make the Bregman distance uninformative for v∉ M_ρ.

Together with a classical source condition, we can then derive the expected convergence rate under the a posteriori choice rule (<ref>).
 
    Assume there exists a w ∈ Y^* with ξ := A^*w∈∂ J(u^†). If ρ=ρ(δ,y^δ) is chosen according to (<ref>), then 
    
    D_J^ξ(u_ρ^δ;u^†) ≤ Cδ.
 


    First note that the proof of <ref> together with <ref> implies that
    
    u_ρ_U=ρ≤u^†_U.

    Hence, we can use the source condition together with the parameter choice rule to estimate
    
    [b]
                D^ξ_J(u_ρ;u^†) 
                    = u_ρ_U - u^†_U- A^*w, u_ρ - u^†_U
       ≤| A^*w , u_ρ - u^†_U| = | w , A(u_ρ - u^†)_Y| 
       ≤w_Y^*Au_ρ - Au^†_Y
       ≤w_Y^*(Au_ρ - y^δ_Y + y^δ - y_Y )
       ≤w_Y^*(τ + 1)δ.
 



    We end this section by remarking on the interpretation of the convergence in Bregman distance in our context. Recall that the subgradient ξ in (<ref>) satisfies by definition.
    
    ξ∈∂·_L^∞(u^†) = ξ∈ L^∞(Ω)^*ξ , u^†_L^∞(Ω) = u^†_L^∞(Ω),
            ξ_L^∞(Ω)^* = 1.

    Assume now that the set C:=x∈Ω|u(x)|=ρ is measurable. 
    In this case, it is straightforward to verify that one possible choice is
    ξ∈  L^1(Ω) ⊂ (L^∞)^*
    defined pointwise almost everywhere by
    
    ξ(x) :=
            -|C|^-1   if x ∈Ωu^†(x) = ρ
    -0    if x ∈Ω|u^†(x)| < ρ
    
                -|C|^-1   if x ∈Ωu^†(x) = -ρ,

    where |C| denotes the Lebesgue measure of C.
    In this case, convergence in Bregman distance entails pointwise convergence u^δ_ρ(x)→ u^†(x) on the active set where |u^†(x)| = ρ. Of course, this choice does not satisfy the source condition ξ∈ R(A^*)⊂ H^1(Ω), so <ref> does not imply a pointwise convergence on this set with rate 𝒪(δ).




§ NUMERICAL SOLUTION


We now address the numerical computation of the quasi-solution problem (<ref>) for the case U=L^∞(Ω) and Y=L^2(Ω) for a bounded Lipschitz domain Ω⊂^d, d∈{1,2,3}, using a combination of a semi-smooth Newton method and a backtracking procedure for the a posteriori parameter choice of ρ. As a model problem, we consider the identification of the source term u∈ L^∞(Ω) from observation of the state y∈ H^1(Ω) given as the weak solution to

    {
                -Δ y + cy     = u   in Ω, 
    ∂_ν y    = 0   on ∂Ω,
            .

for given c> 0. Setting now

    A:L^∞(Ω)→ L^2(Ω),     u↦ y,

it is easy to see that A is linear, injective, continuous, and weakly* closed. It is well-known from PDE-constrained optimization that in this case u_ρ is a solution to (<ref>) (which is known in that context as a bang-bang control problem, see, e.g., <cit.>) if and only if

    u_ρ = _M_ρ(u_ρ - p_ρ),

where 

    _M_ρ:L^2(Ω)→ L^2(Ω),    
        [_M_ρ(v)](x) =: 
        ρ   if    v(x)    >  ρ
    
            v(x)    if   |v(x)| ≤ρ
    
            -ρ   if    v(x)    < -ρ

and p_ρ is the solution to the adjoint problem

    {
                -Δ p + cp    = f   in Ω, 
    ∂_ν p    = 0   on ∂Ω,
            .

for f=Au_ρ-y^δ; see, e.g., <cit.>. (Here we rely on the fact that the adjoint state p_ρ∈ H^1(Ω) to allow a pointwise almost everywhere evaluation.)

We now wish to solve (<ref>) using a locally superlinearly convergent semi-smooth Newton method <cit.>. However, it is known that the pointwise projection is semi-smooth from L^p(Ω)→ L^q(Ω) if and only if p>q, and hence (<ref>) is not semi-smooth with respect to u_ρ. Usually, this is addressed by including additional (small) Tikhonov regularization in (<ref>), which allows canceling u_ρ inside the projection; see, e.g., <cit.> as well as <cit.>. Since the focus of this work is Ivanov regularization, we do not follow this route here and instead consider a finite element discretization of (<ref>).

Let Y_h ⊂ L^2(Ω) be a finite-dimensional space spanned by the usual continuous piecewise linear nodal basis (“hat”) functions based on the vertices {x_j}_j=1^N of a triangulation of Ω. Taking now y_h,u_h,p_h∈ Y_h (which we identify with the vector of their basis coefficients, which coincide with the values at the nodes in this case), noting that functions in Y_h attain their maximum and minimum at the nodes, and introducing the stiffness matrix K_h∈^N× N and the mass matrix M_h∈^N× N, we obtain the discretized conditions

    {
                (K_h+cM_h) y_h    = M_h u_h,
    
                (K_h+cM_h) p_h    = M_h (y_h-y^δ_h),
    
                u_h    = _[-ρ,ρ](u_h-p_h),
            .

where y^δ_h∈ Y_h denotes the L^2 projection of y^δ and the projection is now understood as componentwise. This is now semi-smooth from ^N to ^N, where a Newton derivative at v∈^N direction w∈^N is given componentwise by

    [D_N _[-ρ,ρ](v)w]_i =
        
            w_i       if   |w_i| ≤ρ, 
    
            0       else.
  
Introducing the active and inactive sets as 

    {_+     :=  iu_i - p_i > ρ
    _-     :=  iu_i - p_i <-ρ
        :=  i|u_i - p_i| ≤ρ.

and the corresponding characteristic functions via
[_]_i = 1 if i∈ and 0 else, we can write a Newton step as 

    B_k := 
        [       K_h+cM_h              O           -M_h;           -M_h       K_h+cM_h              O;              O     __+^k∪_-^k _^k-__+^k∪_-^k ][ δ y; δ p; δ u ]
    
        =  -
        [         (K_h+cM_h)y_h^k - M_hu_h^k; (K_h+cM_h)p_h^k - M_h(y_h^k-y^δ_h);           u_h-_[-ρ,ρ](u_h^k-p_h^k) ]
        =:  - b^k

and setting y_h^k+1=y_h^k+δ y, p_h^k+1=p_h^k+δ p, and u_h^k+1=u_h^k+δ u. The iteration is terminated if either a fixed number k_max of iterations is exceeded or if both _+^k=_+^k-1 and _-^k=_-^k-1. In the latter case, (p_h^k,y_h^k,u_h^k) is already a solution to (<ref>); see <cit.>.
To improve the robustness of the procedure, we also include dampening, which requires an additional criterion that b^k_2≤tol for some small tol>0 for successful termination. The full procedure is given in Algorithm <ref>.


In practice, we have to deal with the fact that Newton methods converge only locally. We therefore use a backtracking-type method for computing a regularization parameter ρ satisfying the discrepancy principle (<ref>) as a continuation strategy. Our approach is based on the following two observations: First, since we are dealing with a discretized problem, the operator A_h:=(K_h+cM_h)^-1 is invertible and thus any y^δ_h is in the range of A_h. There thus exists a ρ̂ such that the constraint in the discrete version of (<ref>) is inactive and the semi-smooth Newton method will therefore converge in a single step. Second, the solution for a given ρ_k will be a good starting point for ρ_k+1<ρ_k provided the difference is not too large. We thus proceed in three phases:

    
  I. Starting from some ρ_0≥ 0, we quickly increase ρ_k until the Newton method converged in a fixed number of iterations and the residual is smaller than the noise level.
    
  II. We then quickly decrease ρ_k until the Newton method fails to converge.
    
  III. Starting from the last converged iteration, we adaptively decrease ρ_k further to ensure convergence of the Newton method, until the discrepancy principle (<ref>) is satisfied.

The full procedure is given in Algorithm <ref>. Of course, we also terminate the iteration if the discrepancy principle is already satisfied in phase I or II. Furthermore, to reduce the number of computed quasi-solutions, the increase of ρ_k in phase I can be made adaptive as well.




§ NUMERICAL EXAMPLE


We illustrate the regularization properties of the quasi-solution using a numerical example for the setting considered in <ref>, which is a moderately ill-posed problem. 
Specifically, we take Ω=[-1,1]^2⊂^2 and discretize it using a uniform (Friedrichs–Keller) triangulation with N=128× 128 vertices. The potential coefficient is fixed at c=1. We then choose a piecewise constant true parameter u^† consisting of inclusions of different (positive and negative) heights with ρ^†:=u^†_L^∞(Ω)=4; see <ref>. (For the sake of conciseness, we omit the subscript h in the following.) From this, noisy data is generated by setting

    y^δ = y^† + s/100y^†_∞η/η_2,

where y^† = A_h u^†, s>0 is a chosen relative noise percentage, η∈^N is a random vector with η_i, 1≤ i≤ N, independently normally distributed with mean 0 and standard deviation 1, and η_2 := η^TM_h η is the discrete L^2(Ω)-norm.

We then compute the quasi-solution u^δ_ρ using a Matlab implementation of the procedure discussed in <ref>, where the parameters are chosen as follows. In <ref>, we set q=0.7, i_max=10, k_max=30 and tol=10^-9. In <ref>, we set τ = 1.1 and ρ_0=10.

fig:recon:1fig:recon:0001 show the quasi-solutions for s = 1,0.02,0.01,0.001,0.0001 with corresponding ρ=ρ(δ,y^δ)≈ 2.5,3.896,3.952,3.995,3.999 according to the discrepancy principle. First, we note that in all cases u^δ_ρ_L^∞(Ω) = ρ as expected from <ref>. For s=1, we have in fact that u^δ_ρ(x) = ρ almost everywhere, which leads to a poor reconstruction quality. However, as s and therefore δ:=y^δ-y^†_2 decreases, the quality increases until the quasi-solution is virtually identical to u^† for s=0.0001. While the reconstruction only becomes really acceptable for s<0.01, two properties are of note: First, ignoring the noisy background, the location and shape of the inclusions is recovered well even for large s; in particular, no smoothing is visible in any of the reconstructions. Second, this is especially the case for the inclusions of largest magnitude for which u^†_ρ(x) = ρ^†, for which location and shape are reconstructed perfectly (as can be seen – with a bit of effort – even for s=1 in <ref>). This ties in with <ref>, which indicates improved convergence behavior at these points. 


This is further illustrated by quantitative results for a larger selection of values of the noise percentage s in <ref>, where we give the actual noise level δ, the discrepancy d(ρ,y^δ) achieved by the parameter choice <ref>, and the corresponding reconstruction error of the quasi-solution. For the sake of completeness, we give here both the L^2(Ω)-error u^δ_ρ-u^†_2 and the L^∞(Ω)-error u^δ_ρ-u^†_∞ as well as the error of the duality pairing ξ,u^δ_ρ-u^† in the discrete Bregman distance D^ξ_J(u^δ_ρ;u^†) for ξ∈^N constructed according to <ref>, i.e.,

    ξ_i = 
        -d    if u^†_i = ρ^†,
    -0    if |u^†_i| < ρ^†,
    -d    if u^†_i = -ρ^†,
            for     d=|i|u^†_i| = ρ^†|.

(Note that in the discrete setting, A_h is invertible and hence ξ∈ R(A_h^*).) It can be observed that at least for s≤ 0.01, all these errors behave as 𝒪(δ), corroborating both <ref> and the observed pointwise convergence.




§ CONCLUSION


We have extended the method of quasi-solutions, or Ivanov regularization, to non-reflexive Banach spaces and shown weak* stability and convergence as well as convergence rates in Bregman distances. In particular, we have characterized parameter choice rules that lead to a convergent regularization method. While it turns out that in this setting, a true a priori choice is not feasible, the classical a posteriori choice according to the discrepancy principle is possible and can be exploited as a continuation strategy for the efficient numerical computation of quasi-solutions to inverse source problems using a Newton-type method. Numerical examples illustrate the regularization properties of this approach.

The results in <ref> show that the method of quasi-solutions is indeed a very weak regularization. On the one hand, this is beneficial since it doesn't introduce additional smoothing as Tikhonov regularization with an L^2(Ω) penalty would. On the other hand, the reconstructions show large residual noise unless the data noise is relatively small. This indicates that at least for inverse source problems, Ivanov regularization is likely too weak for practical application in general; an exception would be problems where it is expected that u^†∈{-ρ^†,ρ^†} almost everywhere or where the main interest lies in the correct localization of strong inclusions.

Further study of quasi-solutions in Banach spaces would therefore be justified. Besides the extension to nonlinear parameter identification problems for PDEs mentioned in the introduction, it would be interesting to consider more general choices of M_ρ such as, e.g., M_ρ=u0≤ u(x)≤ρ or M_ρ=uρ^-1≤ u(x) ≤ρ; the former would be relevant for the problem of identifying a potential coefficient in an elliptic PDE (and furnish a link to sparsity regularization), while the latter would be appropriate for identification of a diffusion coefficient. Finally, it is an open question whether heuristic parameter choice rules for ρ are possible, e.g., based on a model function approach for the distance function d(ρ,y^δ).



§ ACKNOWLEDGMENTS


This work was supported by the German Science Foundation (DFG) under grant CL 487/1-1. 
The authors wish to thank Barbara Kaltenbacher for helpful remarks and discussions.



