Iterative Refinement of A Modified Lavrentiev Regularization Method for De-convolution of the Discrete Helmholtz Type Differential Filter
Nathaniel Mays Ming Zhong
math.NA

IRMLR]Iterative Refinement of A Modified Lavrentiev Regularization Method for De-convolution of the Discrete Helmholtz Type Differential Filter


Verona, WI
nate.mays@gmail.com


Department of Applied Mathematics & Statistics, Johns Hopkins University, Baltimore, MD 21218, USA
mzhong5@jhu.edu








We propose and analyze an iterative refinement of a modified Lavrentiev regularization method for deconvolution of the discrete Helmholtz-type differential filter.  The modification for the Lavrentiev regularization method exploits the properties of the Helmholtz filter, and we prove that the modification reduces the error bound between the original solution and the approximated solution.  Furthermore, we derive an optimal stopping condition on the number of iterations necessary for the regularization.  We provide numerical examples demonstrating the benefits of this iterative modified Lavrentiev regularization over a family of Tikhonov regularization methods.


[
    Ming Zhong^*
    December 30, 2023
=====================




§ INTRODUCTION


Since its introduction in 1963, Large Eddy Simulation (LES),  a model in between the direct simulation of the Navier-Stokes equations (DNS) and Reynolds-averaged Navier-Stokes equations (RANS), has been widely applied in meteorology, astrophysics, aerospace, mechanical, chemical and environmental engineering <cit.>.  With a certain low-pass filter, LES is able to reduce the spatial scales (and sometimes even the temporal scales),  therefore reducing the computational load of doing DNS <cit.>.  Most LES filters use convolution, defined as follows,

    u̅(x; δ) = ∫_x' ∈Ω G(x - x'; δ)u(x')  dx'

where Ω is the spatial domain, and δ is the filtering radius (also known as the length scale cutoff) <cit.>.  However the LES filter we are interested in uses a Helmholtz-type differential equation.  In this paper, we will address the issue on how to deconvolve this Helmholtz filter so that we can resolve the smaller length scales below the spatial cutoff, δ.  We are given a Helmholtz filter G and a filtered solution u̅, where Gu = u̅ with u being the desired solution, then we will use an iterative regularization algorithm to find an approximation ũ, such that ũ≈ u, and the error, u - ũ, is well contained.  The mathematical definition of the filtering process is given as follows,  
[Noise Free Model Problem]  
Let X and Y be Hilbert spaces.  Given a linear filter operator G: X → Y and a filtered signal u̅∈Range(G).  The noise free model problem is to find u ∈ X which satisfies

    Gu = u̅.



[Noisy Model Problem] 
Let X and Y be Hilbert spaces.  Given a linear filter operator G: X → Y and a filtered signal u̅∈Range(G) and noise ϵ∈ Y.
The noisy model problem is to find u ∈ X satisfying

    Gu = u̅ + ϵ.


We will mainly consider the case where the linear filter G is a Helmholtz filter, hence solving either Problem <ref> or Problem <ref> becomes a deconvolution problem.  This kind of deconvolution problem is an important inverse problem <cit.>.  This problem occurs in many applications including parameter identification <cit.>, the deconvolution problem of image processing <cit.>, and the closure problem in turbulence modeling <cit.>.  The deconvolution problem gets more complicated when noise is present in the filtering process.
It is known that if G is compact and Range(G) is infinite dimensional, then Problems <ref> and <ref> are ill-posed <cit.>.  Tikhonov-Lavrentiev regularization, a regularization method described further in Algorithm <ref>, which introduces a regularization parameter α, is a well-established method used to solve Problems <ref> and <ref> <cit.>.  However, Tikhonov-Lavrentiev regularization is a general method which does not exploit the properties of the Helmholtz-type differential filter (with filtering radius δ).  Therefore we introduce a new regularization method with a modification to iterative Tikhonov-Lavrentiev regularization, which in Theorem <ref> shows that it is able to exploit the properties of the Helmholtz filter G and improves the error bounds for Problems <ref>. That is, a small algorithmic modification leads to a large improvement in the error bounds.


In section <ref>, we introduce the necessary notation and inequalities which are used in the proofs of our theorems.  Section <ref> describes the Modified Iterative Tikhonov-Lavrentiev Regularization (Mitlar) algorithm in details and shows the error bounds in both the continuous case and discrete case.  Section <ref> provides an optimal stopping condition for the total number of iterations to counter the presence of noise in the the filtering process.  Finally, in section <ref>, we show a number of numerical examples to verify the convergence rate, the optimal stopping condition and also compare the performance of our new algorithm to a family of existing Tikhonov-Lavrentiev regularization methods: original Tikhonov-Lavrentiev regularization, iterative Tikhonov-Lavrentiev regularization, and modified Tikhonov-Lavrentiev regularization.



§ PRELIMINARIES AND NOTATION


Throughout this paper, we use the standard notation for Lebesgue and Sobolev spaces and their norms.  Also, Ω will be a regular, bounded, polyhedral domain in ℝ^n.  We define the following space

    X = H^1_0(Ω)^d = { v ∈ L^2(Ω)^d: ∇ v ∈ L^2 (Ω)^d× d and  v = 0  on ∂Ω}.

The norm · (when subscript is not present) will also denote the L^2(Ω) norm unless otherwise specified in a proof.  Similarly the inner product (·, ·) will denote the L^2(Ω) inner product.
We will use the notation X^h ⊂ X to denote a finite dimensional subset of X.  An example of X^h is the set of continuous polynomials of degree k.  We also assume that we have homogeneous boundary data throughout.
We use the following approximation inequalities, see <cit.>,

    inf_v ∈ X^h u - v _L^2(Ω)   ≤ C h^k+1 u _H^k+1(Ω),    u ∈ H^k+1(Ω)^n, 
    inf_v ∈ X^h u - v _H^1(Ω)   ≤ C h^k u _H^k(Ω),    u ∈ H^k+1(Ω)^n.

Other well known inequalities used herein include:

  
  * Cauchy-Schwartz inequality: | ( f,g ) | ≤fg,   ∀ f,g ∈ L^2(Ω).
  
  * Young's inequality: ab ≤ϵ/p a^p + ϵ^-q/p/qb^q, where 1 < p,q < ∞, 1/p + 1/q = 1, ϵ > 0, and a,b ≥ 0.
  
  * Poincare-Friedrich's inequality: v ≤ C_PF∇ v ,   ∀ v ∈ X.
  
  * Triangle inequality: a + b ≤a  + b.




 §.§ The differential filter


The differential filter (also known as the Helmholtz-type differential filter) is used in multiple large eddy simulation models <cit.>.  This filter is equivalent to the Pao filter used in image processing <cit.>.  

The differential filter G is defined as Gu=u̅, where u and u̅ satisfy

    -δ^2 Δu̅ + u̅ = u    in Ω.








The differential filter is equivalent (see <cit.>) to the following variational formulation.  Find u∈ H^1_0(Ω) satisfying

    δ^ 2 ( ∇u, ∇ v ) + (u, v) = (u,v),  ∀ v ∈ X.









  Let X^h be a finite dimensional subspace of X.  We define G^h:L^2(Ω)^d → X^h where
u^h = G^h u which is the unique solution in X^h to 

    δ^2 u^h, v^h + u^h,v^h =
    u,v^h,  ∀ v^h ∈ X^h.



Lemmas <ref>, <ref>, <ref>, and <ref> are quoted from <cit.> for completeness.


If u ∈ L^2(Ω)^d, the following stability estimate for problem (<ref>) holds:

    δ^2 ∇u^2 + 1/2u^2 ≤1/2 u ^2.




The operator G: L^2(Ω)^d → X is self-adjoint.



 If ∇ u ∈  L^2(Ω)^d and u satisfies (<ref>), then
 
    δ^2/2∇( u - u ) ^2 +  u - u^2 ≤δ^2/2∇ u ^2.

 If, additionally Δ u ∈ L^2(Ω)^d, then
  
    δ^2 ∇( u - u ) ^2 + 1/2 u - u^2 ≤δ^4 Δ u ^2.




The operator G^h:L^2(Ω)^d → X^h is self-adjoint and positive semi-definite on L^2(Ω) and positive definite on X^h.




 §.§ Tikhonov regularization


A major tool which is often used to solve inverse problems is the Tikhonov regularization <cit.>.  In order to solve an ill-posed linear system Gu = u̅, a general Tikhonov regularization will try to find the minimizer of the following problem, 

    u_T = u{Gu - u̅^2 + Γ u^2}

where Γ is a suitably chosen linear operator, called Tikhonov operator (or Tikhonov matrix when G is a matrix).  The Tikhonov minimizer of (<ref>), u_T, has a closed form expression, 

    u_T = (G^*G + Γ^*Γ)^-1G^*u

where G^* and Γ^* are Hermitian transposes of G and Γ respectively.  In most cases, Γ is picked as a multiple of the identity operator, that is, Γ = α I.  In this case, (<ref>) simplifies down to, 

    u_T = (G^*G + α^2 I)^-1G^*u

If the operator G is monotone[In the case of G being the Helmholtz filter operator, G is self-adjoint and positive semi-definite, hence monotone.], then instead of solving Gu = u̅ with the the perturbed normal equation from Tikhonov regularization, Tikhonov-Lavrentiev regularization can be used.  This is known in the literature as the method of Lavrentiev Regularization <cit.> or the method of Singular Perturbation <cit.>,
 
  Choose a regularization parameter α > 0.  Solve for u_0 satisfying

    ( G + α I ) u_0    = u,   in Ω.


This regularization method can be improved by an iterative method.  

Choose a regularization parameter α>0 and fix the number of updates J ≥ 1.  The iterative Tikhonov-Lavrentiev approximations u_j (0 ≤ j ≤ J) are found by solving

    ( G + α I ) u_0    = u,   in Ω,
    
      ( G + α I ) ( u_j - u_j-1 )    = u - G u_j-1,   in Ω.


Given a source condition, it is shown in <cit.> that Tikhonov-Lavrentiev and iterative Tikhonov-Lavrentiev regularization converge to u as α→ 0 (in the noise free case) and as ϵ→ 0 and α→ 0 (in the noisy case).

  Suppose that G is non-negative definite.  Fix α>0.  Let e_j = u - u_j for all j = 0, …, J.  Suppose, for some β≥ 0 that u∈ Range(G^β) and the noise is bounded ϵ≤ϵ_0 < ∞.  Then, there exists a constant C(J) < ∞ such that, for any 0 ≤ J ≤β,

    e_J ≤(J+1)ϵ_0/α + α^J+1 C(J).

Moreover, if α = α(ϵ_0) = C ϵ_0^1/(J+2) we have that e_J ≤ C ϵ_0^1 - 1/(J+2).

The error bounds for Tikhonov-Lavrentiev regularization are similar to that of Tikhonov regularization, see <cit.>.



§ MODIFICATION TO ITERATIVE TIKHONOV-LAVRENTIEV REGULARIZATION


Algorithm <ref> describes the modification which we add to iterative Tikhonov-Lavrentiev Regularization to develop our new algorithm, the Modified Iterative Tikhonov-Lavrentiev Regularization (Mitlar).  We analyze the error in the continuous case by separating it into the following components: the regularization error in the Mitlar algorithm and the amount of noise amplification due to our regularization.  We then discretize Mitlar in Algorithm <ref>. We analyze the error in the discretized case by separating it into the following components: the regularization error in the continuous Mitlar algorithm, the discretization error in the solution, and the discretized noise amplification due to the discrete Mitlar algorithm.

We define the following regularization operators D_α and D_α, j for convenience of notation.

  For α > 0 define the modified Tikhonov-Lavrentiev operator D_α to be

    D_α  = [(1-α)G+α I]^-1.

For j > 0, define the j^th modified iterative Tikhonov-Lavrentiev operator D_α,j by

    D_α,ju = u_j,

where u_j is obtained via Algorithm <ref>.


Assume G is the differential filter defined in (<ref>). Algorithm <ref> is equivalent  to the following variational formulation.  Given u∈ L^2(Ω), then u_J = D_J u̅ is the unique solution to the following equations

    αδ^2 ∇ u_0,∇ v + u_0,v   =  δ^2∇u̅,∇ v + u̅, v,   ∀ v ∈ X,  and
    αδ^2 ∇ u_j,∇ v + u_j,v   =  δ^2∇u̅,∇ v + u̅, v + αδ^2 ∇ u_j-1,∇ v,  ∀ v ∈ X.


Theorem <ref> shows that this modification to Tikhonov-Lavrentiev regularization provides a higher order deconvolution error compared to iterated Tikhonov-Lavrentiev regularization.  The following lemmas and propositions are needed for the proof of the error bound in Theorem <ref>.
 
For 0 ≤α≤ 1, the function f(x) = ((1-α)x + α)^-1 maps the interval (0,1] to [1,1/α), and the function g(x) = x((1-α)x + α)^-1 maps the interval (0,1] to (0,1].


The term (1-α)x + α is a convex combination of x and 1, so 

    α < (1-α)x + α≤ 1,  and
       1 ≤1/(1-α)x + α < 1/α.

For the bounds on g(x), consider

    g'(x) = α ((1-α)x + α)^-2.

So g'(x) has no critical points in the interval (0,1).  Therefore g(x) attains its extrema on the boundary of [0,1].  Note that  g(1) = 1, and g(x) = 0 if and only if x = 0.  Therefore g: (0,1] → (0,1].


  For 0 ≤α≤ 1, the operators D_α, D_α G, and I -D_α G  are bounded.  In particular, they satisfy

    D_α≤1/α,    D_α G ≤ 1,  and  I - D_α G ≤ 1.



The method of proof is similar to that employed in <cit.>.  
The differential filter operator G has a spectrum that lies in (0,1]. 
Therefore by Lemma <ref>, the spectrum of D_α = ((1- α)G + α I)^-1 lies between [1, 1/α).  
Also by Lemma <ref>, the spectrum of D_α G = ((1- α)G + α I)^-1 G lies between (0,1].  Similarly, the spectrum of I - D_α G lies between [0,1).

 
  The error equation e_J = u - u_J is given by

    e_J = (-αδ^2)^J+1 (D_α G)^J+1(Δ^J+1u),

and the error is bounded

    e_J ≤ (αδ^2)^J+1Δ^J+1 u.



  For 0 < j ≤ J, we start with (<ref>) and an identity for the original solution u,

    [(1-α)G+α I](u_j - u_j-1)     = u̅-Gu_j-1  and
    
      [(1-α)G+α I](u - u)     = u̅ - Gu.

Subtracting these equations and rearranging gives

    e_j = α D_α (I - G) e_j-1 = α D_α G (G^-1 - I)
    e_j-1.

For j = 0, we use (<ref>) and the true solution 

    [(1-α)G+α I]u_0   =u̅  and
    
      [(1-α)G+α I]u    =(1-α)u̅ + α u.

Subtraction gives

    e_0 = α D_α (I - G) u = α D_α G (G^-1 - I) u.

Thus, we arrive at,

    e_J = α^J + 1(D_αG)^J + 1(G^-1 - I)^J + 1u

Using the fact that (G^-1 - I) u = -δ^2Δ u, we have shown equation (<ref>).  The norm of the error is bounded by taking the norm of the error equation (<ref>) and using the bound on D_α G in (<ref>) to obtain

    e_J   = (-αδ^2)^J+1 (D_α G)^J+1(Δ^J+1u)
       ≤ (αδ^2)^J+1(D_α G)^J+1Δ^J+1u
       ≤ (αδ^2)^J+1Δ^J+1u.


 
The j^th step of the Mitlar algorithm, u_j, is given by

    u_j:= D_α,ju 
         =D_α∑_i=0^j(α D_α(I-G))^i u.



  Starting with (<ref>), solve for u_j with the equations 

    I - D_α G = α D_α (I - G)   and   u_0 = D_αu.


    u_j    = u_j-1 + D_α(u - G u_j-1)
       = D_αu + (I - D_α G) u_j-1
       = D_αu + α D_α (I - G) u_j-1
       ⋮
       = D_α∑_i=0^j(α D_α(I-G))^i u.

as claimed.

Noise amplification is one of the fundamental difficulties in solving ill-posed inverse problems <cit.>.  The noise amplification is studied in Problem <ref> where u̅ has additive noise ϵ.  The Mitlar algorithm applied to this problem gives an improvement over iterated Tikhonov regularization in the noise free portion of the error as shown in Proposition <ref>.  The bound on the error in the noisy data is no worse as shown in Theorem <ref>.

  Under the conditions of Algorithm <ref> and (<ref>) and
if there exists some ϵ_0 such that ϵ <
ϵ_0, then the error in the j^th step of the Mitlar algorithm is

    e_J = (-αδ^2)^J+1 (D_α G)^J+1(Δ^J+1u) 
          + D_α∑_j=0^J(I - D_α G)^j ϵ.

The error is bounded,

    e_J ≤ (αδ^2)^J+1Δ^J+1 u +
    (J+1)ϵ_0/α.



Using the definition the J^th modified iterative Tikhonov-Lavrentiev operator D_α, J and Gu = u̅  + ϵ, we have

    u_J    = D_α,Ju
       = D_α,J (Gu - ϵ) 
       = D_α,JGu - D_α,Jϵ

We then divide the error, e_J = u - u_J, into two parts,

    e_J = (u - D_α,JGu) + D_α,Jϵ

Together with Proposition <ref> and <ref>, we have

    e_J = (-αδ^2)^J + 1(D_αG)^J + 1(Δ^J + 1u) + D_α∑_j = 0^J(I - D_αG)^j ϵ

as claimed.  To get a bound on the norm of the error, start with the error equation and take the norm and use the inequalities in (<ref>).

    e_J    =  (-αδ^2)^J+1 (D_α G)^J+1(Δ^J+1u) 
          + D_α∑_j=0^J(I - D_α G)^j ϵ
       ≤ (αδ^2)^J+1Δ^J+1u 
          + D_α∑_j=0^J(I - D_α G)^jϵ
       ≤ (αδ^2)^J+1Δ^J+1u 
          + 1/α∑_j=0^Jϵ_0
       ≤ (αδ^2)^J+1Δ^J+1 u +
    (J+1)ϵ_0/α.


We see that this is an improvement over iterative Tikhonov-Lavrentiev regularization because of its double asymptotic behavior in α and δ of the error bound.  Each update step in the method contributes an extra factor of αδ^2 , whereas each update step of iterative Tikhonov-Lavrentiev contributes only an extra factor of α.  However, we also notice that each iterative step adds an term, ϵ_0/α, to the error bound, thus possibly increasing the error.  How to balance the total number of iterations and the noise becomes significant in reducing the error bound.  And we will discuss such relationship in section <ref>.



 §.§ Discrete Mitlar Applied to the differential operator


The results of the previous section are now extended to the discrete form of the Mitlar algorithm.  The modified iterated Tikhonov-Lavrentiev regularization operator applied to the differential filter is defined in Definition <ref> variationally on a finite dimensional space. 


Given a filter radius δ > 0 of the differential filter operator G, and fix a regularization parameter 0 ≤α≤ 1 and stopping number J≥ 0.  
If Δ^j u _L^2(Ω) is bounded for all j ≤ J+1, then the error to the problem in (<ref>) using the discrete Mitlar algorithm is bounded.  In particular,

    u - D_α, J^h G^h u _L^2(Ω)   ≤ (αδ^2)^J+1Δ^J+1 u_L^2(Ω)
          + C (√(α)δ h^k + h^k+1) max_0≤ j ≤ J D_α,j G u _H^k+1(Ω).



We denote ·_L^2(Ω) by ·.  Then we add and subtract the exact deconvolution term, and use the triangle inequality,

    u - D_α, J^h G^h u ≤ u - D_α, J G u  +  D_α, J G u - D_α, J^h G^h u .

The first term of (<ref>) is bounded by (<ref>)

    u - D_α, J G u ≤ (αδ^2)^J+1Δ^J+1 u .

For the second term of (<ref>), start with (<ref>) and take v = v^h, then subtract equation (<ref>). For j=1, …, J, we have

    αδ^2 ∇ (u_j-u_j^h),∇ v^h + u_j-u_j^h,v^h 
      =   αδ^2 ∇ (u_j-1 - u_j-1^h) ,∇ v^h.

The case when j=0 follows similarly or see <cit.>.
We define η_j = u_j - w_j^h and ϕ^h_j = u^h_j - w_j^h for some w_j^h ∈ X^h to be chosen later for each j=1,…,J.  Using these definitions, we write (<ref>) as

    αδ^2 ∇ (η_j-ϕ_j^h),∇ v^h + η_j-ϕ_j^h,v^h 
      =   αδ^2 ∇ (η_j-1 - ϕ_j-1^h) ,∇ v^h.

Take v^h = ϕ_j^h, denote e_j = u_j - u_j^h = η_j - ϕ_j^h, and separate the terms to get

    αδ^2 ∇ϕ_j^h ^2 + ϕ_j^h ^2    = 
          αδ^2 ∇η_j, ∇ϕ_j^h + η_j,ϕ_j^h - αδ^2 ∇ e_j-1, ∇ϕ_j^h
       ≤αδ^2 ∇η_j ^2 + αδ^2/4∇ϕ_j^h ^2 
          + 1/2η_j ^2 + 1/2ϕ_j^h^2
            + αδ^2 ∇ e_j-1^2 + αδ^2/4∇ϕ_j^h^2.

Move the ∇ϕ_j^h^2 and ϕ_j^h^2 terms from the left hand side to the right, and then multiply by 2 to get

    αδ^2 ∇ϕ_j^h ^2 +  ϕ_j^h ^2 
             ≤ 2 αδ^2 ∇η_j ^2  
          +  η_j ^2 
          + 2αδ^2 ∇ e_j-1^2.

Use e_j≤η_j  + ϕ_j^h and ∇ e_j≤∇η_j  + ∇ϕ_j^h to obtain the recursion 

    αδ^2 ∇ e_j^2 + e_j^2 ≤ 3 αδ^2 ∇η_j ^2  
          +  2 η_j ^2 
          + 2 αδ^2 ∇ e_j-1^2.

Thus

    e_J≤ C(J) max_0≤ j ≤ J( √(αδ^2)∇η_j   
          +   η_j )

This inequality holds for all w_j^h ∈ X^h, so take the infimum over X^h and apply the approximation inequalities (<ref>) to obtain

    D_α, J G u - D_α, J^h G^h u ≤   C (√(α)δ h^k + h^k+1) max_0≤ j ≤ J D_α,j G u _H^k+1(Ω).

Combining equations (<ref>) and (<ref>) proves the claim.

Problem <ref> still needs to be addressed.  If our data consists of discrete measurements that contain noise ϵ, making G^hu = u^h + ε, then approximations of the error from that noise are needed.  This problem is addressed by applying the discretized modified iterated Tikhonov-Lavrentiev algorithm to the discretized and noisy data u^h.
First, we prove the boundedness of operators G^h, D^h, and D^h_J.
 
The operators G^h:X^h → X^h, D^h:X^h → X^h, and D^h_J:X^h → X^h are bounded and furthermore they satisfy

    G^h   ≤ 1,
    D^h   ≤1/α, 
    D^h_J   ≤J+1/α, 
    D^h G^h   ≤ 1,  and
    I - D^h G^h   ≤ 1.



For the first, take u ∈ X^h then G^h u≤u^h by Cauchy-Schwartz and Young inequalities to equation (<ref>).  For the second, note that D^h = [(1-α)G^h + α I]^-1 is the convex combination of positive operators, so its spectrum is bounded by 1/α.  For the third, we write out D^h_J = D^h ∑_i=0^J (α D^h(I - G^h))^i.  Then taking u ∈ X^h we obtain,

    D^h_J u   = D^h ∑_j=0^J (α D^h(I - G^h))^j u 
       ≤D^h∑_j=0^J (α D^h(I - G^h))^ju
       ≤1/α∑_j=0^J u
       ≤J+1/αu.

The spectrum of D^h G^h lies in between (0,1] and the spectrum of I - D^h G^h lies in between [0,1) proving the result.


  If the noise ε∈ X^h is bounded ε≤ϵ_0, then the error e_j between the desired solution, u,  and the discretized Mitlar solution applied to noisy data u^h with G^hu =  u^h+ ε is bounded, and

    e_j   :=  u - D_J^h u^h 
       ≤J+1/αϵ_0+ (αδ^2)^J+1Δ^J+1 u
            + C (√(α)δ h^k + h^k+1) max_0≤ j ≤ J D_α,j G u _H^k+1(Ω).



Use the triangle inequality to separate the error into two pieces, the true discretization error and the error associated with noise
  
    u - D_J^h u̅^h    ≤ u - D_J^h(G^hu - ε)
       ≤ u - D_J^hG^hu  +  D_J^h ε.

Using Theorem <ref> to bound the first term,

    u - D_J^h G^hu    ≤ (αδ^2)^J+1Δ^J+1 u
          + C (√(α)δ h^k + h^k+1) max_0≤ j ≤ J D_α,j G u _H^k + 1(Ω).

Lemma <ref> gives a bound on the second term,

    D_J^h ε≤J+1/αε.

Combine these results to prove the claim.




§ DESCENT PROPERTIES OF MODIFIED ITERATED TIKHONOV-LAVRENTIEV REGULARIZATION


For a self-adjoint and positive definite G, solving Problem <ref> is equivalent to solving the following minimization problem

    w = argmin_v ∈ X E_0(v),  where  E_0(v)1/2(G v,v) - (u̅,v).

We analyze when the Mitlar approximations, {u_0, u_1, …, u_J}, will form a minimizing sequence for E_0.

Let G be self-adjoint and positive definite and 0 < α≤1/2. Then the Modified Iterative Tikhonov-Lavrentiev iterates are a minimizing sequence for E_0. In particular, 

    E_0(u_j) - E_0(u_j+1)=([(1/2-α)G+α I](u_j+1-u_j),u_j+1-u_j)≥ 0.

with equality achieved only when u_j = u_j + 1.  Thus

    E_0(u_j+1) < E_0(u_j),  unless  u_j+1=u_j.



Expand using the definition of E_0(·) and cancel terms to prove the identity along with the fact that G is self-adjoint.

    E_0(u_j) - E_0(u_j+1)    = 1/2(G u_j,u_j)-(u̅,u_j) - 1/2(G u_j+1,u_j+1)+(u̅,u_j+1)
       =1/2(G u_j, u_j - u_j+1) + 1/2(G (u_j - u_j+1), u_j) 
            - (u̅,u_j - u_j+1 )
       = 1/2(G (u_j - u_j+1), u_j + u_j+1) 
            - ( [ (1-α)G + α I ] (u_j+1 - u_j) + G  u_j, u_j - u_j+1)
       = 1/2(G (u_j - u_j+1), u_j + u_j+1) - ( G (u_j - u_j+1),  u_j+1 ) 
            + α ( [I - G] (u_j - u_j+1), u_j - u_j+1)
       = ([(1/2-α)G+α I](u_j+1-u_j),u_j+1-u_j).


Equation (<ref>) stays positive as long as 0 < α≤1/2, hence E_0(u_j+1)<E_0(u_j) unless u_j+1=u_j as claimed.  

Equation (<ref>) implies that if u_j = u_j+1, then G u_j = u̅.  However, in the noisy case, such convergence is not desired, since u̅ = Gu - ϵ according to Problem <ref>.  Thus, as j→∞, u_j → Gu - ϵ.  This implies that it is critical to stop after a finite number of update steps.
We seek the desired solution, u, to Problem <ref> when noise is present in the filtering process.  Similarly, finding such u is equivalent to the finding the minimizer of the following minimization problem,

    w = argmin_v ∈ X E_ϵ(v),  where  E_ϵ(v)1/2(G v,v) - (u̅ + ϵ,v).

Then we analyze the sequence of noisy Mitlar approximations u_j's in the noisy functional, E_ϵ(·).  Again we expand the difference, E_ϵ(u_j)-E_ϵ(u_j+1), and the following is obtained,

    E_ϵ(u_j)-E_ϵ(u_j+1) = (ϵ,u_j+1-u_j) + ([(1/2-α)G+α I](u_j+1-u_j),u_j+1-u_j)


Let G be self-adjoint and positive definite. Suppose
an estimate on the noise ϵ≤ϵ_0 is known.  Then the solutions from the Modified
Iterative Tikhonov-Lavrentiev Regularization are a minimizing sequence for the noisy functional E_ϵ as long as

    ϵ_0/u_j+1-u_j≤α≤1/2.



First, the Cauchy-Schwartz inequality implies

    | (ϵ, u_j+1 - u_j) | ≤ϵ u_j+1 - u_j

Then, if (<ref>) holds, then

    0    ≤ϵ u_j+1 - u_j - | (ϵ, u_j+1 - u_j) | 
       ≤ϵ_0 u_j+1 - u_j - | (ϵ, u_j+1 - u_j) | 
       ≤α u_j+1 - u_j^2 -  | (ϵ, u_j+1 - u_j) |
       ≤( [ (1/2 - α) G + α I] ( u_j+1 - u_j ), u_j+1 - u_j) - | (ϵ, u_j+1 - u_j) | 
       ≤( [ (1/2 - α) G + α I] ( u_j+1 - u_j ), u_j+1 - u_j) + (ϵ, u_j+1 - u_j) 
        = E_ϵ(u_j)-E_ϵ(u_j+1).


Theorem <ref> implies that when the size of the updates is larger than twice the noise, the updates move the solutions closer to desired solution, u. As the updates become smaller, u_j begins to deviate from an
approximation of  u unless α is increased.
This result can be extended if more is known about the noise or its statistical distribution.  In particular if there exists a projection operator P where Pϵ⊥ (u_j+1 - u_j), then

    (ϵ,u_j+1-u_j) = (ϵ,(I-P)[u_j+1-u_j]).

In other words, if a component of the Mitlar update is in the range of the projection, then that updated component will reduce the error to the desired solution, u.  This suggests the following small algorithmic modification.




§ NUMERICAL ILLUSTRATIONS


We investigate several applications. In section <ref>, we verify the use of our stopping criterion.  In section <ref>, we compare the family of Tikhonov-Lavrentiev regularizations: Tikhonov-Lavrentiev, iterated Tikhonov-Lavrentiev, modified Tikhonov-Lavrentiev, and Mitlar in the application of deconvolving the differential filter.  Section <ref> verifies and compares the convergence rates of the 4 different Tikhonov-Lavrentiev regularization methods.





























































































































 §.§ An example of the stopping criterion



We implement and verify the optimal stopping criterion (Theorem <ref>)  in MATLAB (version: R2014b) with the following details: first we choose a true solution to be u = sin(π x) + sin(200π x), plotted in Figure <ref>, over the interval [0, 2].

Then we discretize the interval with a step size of h = 2/1000 (hence 1001 sample points) and choose the filtering radius for the differential filter to be δ = 6h.  To implement the Helmholtz filter, we begin with approximating the Laplacian operator with a center differencing scheme,

    Δ u ≈Δ^h u = u(x - Δ x) - 2u(x) + u(x + Δ x)/(Δ x)^2.

and define our discrete operator, namely A^h (the inverse operator to G^h) as

    A^h u̅ = -δ^2Δ^hu̅ + u̅  and   A^hG^h = I.

Our simulated data was obtained by filtering the true solution and adding 1% random noise to the filtered data, that is, u̅ = G^hu - ϵ where ϵ = 0.01 *  G^hu (ϵ is generated in Matlab using the command “randn", and normalized to have L^2 norm of 1).  And for calculating the L^2 norm of a function f over [a, b], we use either the composite Trapezoidal rule or the composite Simpson's rule.  We select the regularization parameter α = 0.1 (α≤ 0.5 to have a decreasing sequence for the noise free energy functional) for Mitlar.  And we use the following guideline for finding the optimal stopping J:

  
  * At the j^th iterate, except the initial iterate, we calculate u_j+1-u_j.
  
  * We then compare α, the regularization parameter, to ϵ_0/u_j+1-u_j.
  
  * According to (<ref>): if α≥ϵ_0/u_j+1-u_j, we proceed to next iterate; when α < ϵ_0/u_j+1-u_j, we stop the iteration.

The actual simulation which we did for this demonstration, on the other hand, will not stop once we find the stopping J; instead the optimal J will be recorded.

Figure <ref> shows the noisy energy functional, E_ϵ, calculated with the Mitlar approximation u_j's when there is noise in the filtering process, i.e., Gu = u̅ + ϵ. The calculated optimal stopping point (via Theorem <ref>) occurs after J = 4 iteration steps and is shown as a green dot.  The figure shows that the algorithm stops right where functional reaches its minimum and starts increasing again, hence avoiding the convergence to the noisy solution.  However, because we do not have precise information of the noise (as it is always the case in real life applications), the algorithm will always try to stop before the energy functional reaches its minimum.




 §.§ Comparison of four deconvolution algorithms


We check the efficiency of Algorithm <ref> by comparing the relative error of a solution for a given parameter α to the relative errors found with Tikhonov-Lavrentiev, iterative Tikhonov-Lavrentiev, and Modified Tikhonov-Lavrentiev using the same α.  We implement the codes in MATLAB (version R2014b) with the following details: we start out with the original data as u = sin(π x) + 0.1sin(100π x), with 1001 sample points taken over the interval [0, 2], see figure <ref>; hence the step size is h = 2/1000.  We set our filtering radius at δ = 0.01.  We let the α vary from 1 to 10^-3 and calculate 1, 2, and 3 steps for the two iterative methods.   For an approximation u_appr. to a desired solution u,  the relative error, Err_rel, is defined as, Err_rel = u - u_appr./u.  The results are shown in Figures <ref>, <ref>, and <ref> respectively.  We see that the error in Mitlar is the lowest for any given α; meanwhile as the number of iterates increases, the accuracy of Mitlar increases.  We would like to point out that, for bigger J (total iteration numbers), the difference between Mitlar and the other 3 regularization methods widens for small α.











 §.§ Verification of convergence rates


We calculate the convergence rates of the family of Tikhonov-Lavrentiev regularization methods to verify the convergence rates predicted in Theorem <ref>.  We take a true solution over the domain [0, 2]× [0, 2] of

    u = sin(π x) sin(π y) + sin(20π x)sin(20π y).

see Figure <ref>.  

We discretize using the square command in FreeFEM++ <cit.> with n intervals in each of x and y coordinates and use piecewise continuous linear polynomials.  We use a filter radius of δ = 0.1(2π/n)^1/4 = 𝒪(h^0.25) and regularization parameter α = 0.1(2π/n)^1/2 = 𝒪(h^0.5).  And the results are presented in the following tables.







§ CONCLUSION


We introduced a novel tool for solving some special types of inverse problems, that is, to deconvolve solutions filtered by a Helmholtz-type differential filter.  We show that the noise free errors in using Mitlar are doubly asymptotic in α and δ, that is 𝒪((αδ^2)^J+1).  However, using the Tikhonov-Lavrentiev regularization or iterated Tikhonov-Lavrentiev regularization only results in the noise free errors depending solely on α, which are 𝒪(α) and 𝒪(α^J+1) respectively.  


We also introduce a tool for calculating when to stop the iterations for our iterative algorithm.  We show that continuing to iterate until the solution converges gives the unwanted and noisy solution.  However, our stopping criterion guarantees that the iteration steps are getting closer to the original solution.  The example chosen to illustrate the stopping criterion showed the exact optimal stopping, which is not always guaranteed due to insufficient knowledge of the noise.  The actual implementation, on the other hand, will always try to stop before reaching the optimal number of iterations.  And if we can incorporated more knowledge about the noise added into the model, then we would be able to get a more accurate bound on the optimal number of iterations.  


Furthermore, we would like to point out that the analysis done in Mitlar can be applied to other filters whose transfer function behaves like the Pao filter for small wave lengths, due to the robustness of the Tikhonov-Lavrentiev regularization method.


siam

