Sobolev spaces with non-Muckenhoupt weights, fractional elliptic operators, and applications
Harbir Antil Carlos N. Rautenberg


Automated tuning of inter-dot tunnel couplings in quantum dot arrays
    L. M. K. Vandersypen
    December 30, 2023
====================================================================





We propose a new variational model 
in weighted Sobolev spaces with non-standard weights and  applications to image processing. We show that these weights 
are, in general, not of Muckenhoupt type and therefore the classical analysis tools may not apply. For special cases of the weights,
the resulting variational problem is known to be equivalent to the fractional Poisson 
problem. The trace space for the weighted Sobolev space is identified to be embedded in a weighted L^2 
space.  We propose a finite element scheme to solve the Euler-Lagrange 
equations, and for the image denoising application we propose an algorithm to identify the 
unknown weights. The approach is illustrated on several test problems
and it yields better results when compared to the existing total variation techniques. 





 Variable weights, non Muckenhoupt weights, new trace theorem, fractional Laplacian with variable exponent, image denoising.




 35S15, 
 26A33, 
 65R20, 
 65N30  








§ INTRODUCTION




In this paper we consider weighted Sobolev spaces where the weight w is not necessarely of Muckenhoupt type, and an associated variational model with concrete applications to image processing that shows advantageous features of such weights. The particular weights w that we consider in this paper are closely related to fractional Sobolev spaces of differentiability order s∈ (0,1) and the fractional spectral Laplacian (-Δ)^s, and are further motivated by considering a spatially dependent order x↦ s(x). 

Weighted Sobolev spaces have been a topic of intensive study for around 60 years with a variety of focuses; we refer the readers to the monographs <cit.> and references within, and further <cit.> for an introduction to the subject. In particular, a source of interest for such spaces is that they represent the correct solution space for several degenerate elliptic partial differential equations (<cit.>) and problems in potential theory (<cit.>). Recently, the topic has received a new impulse mainly associated with the extension result of Caffarelli-Silvestre; <cit.> (Stinga-Torrea; <cit.>) for fractional elliptic partial differential operators. In this setting, for example, the solution to (-Δ)^su=f in the fractional Sobolev space H^s(Ω)  endowed with zero boundary conditions, can be equivalently obtained as (the restriction to Ω of) the solution of a PDE with a non-fractional elliptic operator in a weighted Sobolev space with weight w(x,y)=y^1-2s in the extended domain 𝒞={(x,y)∈Ω× (0,+∞)}.


The type of weights w that have been more thoroughly studied correspond to two main classes: 1) weights belonging to some Muckenhoupt class; see <cit.> and 2) composition of functions (mainly power functions) with the distance function to a particular set; see <cit.>. In these cases, several important questions like density of smooth functions and characterization of traces have been, at least partially, answered; see <cit.>. In this paper we consider weights of the type w(x,y)=y^1-2s(x) on 𝒞={(x,y)∈Ω× (0,+∞)}, which are, in general, neither of 1) or 2) type.

The variational problem of interest in this paper is closely related to fractional elliptic operators. The spike of research interest on problems with this type of operators is mainly due to the results from Caffarelli and collaborators; see <cit.>. Following the renowned  results in <cit.>, a large number of contributions on modeling, numerical methods, applications, regularity results, different boundary conditions, and control problems, among others have been considered. 

Recently, for image denoising a model was 
proposed in <cit.>. The variational problem can be formulated in our setting as 

    min_u1/2(-Δ)^s/2 v_L^2(Ω)^2 
         + ζ/2v-f^2_L^2(Ω)  ,

where (-Δ)^s denotes the fractional power of the Laplacian with zero Neumann boundary
conditions and  s ∈ (0,1); see <cit.>. Additionally, ζ>0 is a given constant, and f=u_true+noise where u_true is the desired target of the optimization procedure. The choice of s has a direct consequence on the global regularity of the solution to problem (<ref>). However, it is desirable, from the reconstruction point of view, that the regularity of the solution to (<ref>) is low in places in Ω where edges or discontinuities are present in u_true, and that is high in places where  u_true is smooth or contains homogeneous features. Hence, it is of interest to consider (<ref>) where s:Ω→ [0,1] is not a constant.  The first main roadblock in letting s to be spatially dependent is the fact that 
there is no obvious way to define (-Δ)^s(x). In fact we will not attempt to 
do so in this paper, we leave this as an open question. 
Instead of (<ref>) we consider the following variational problem:

    min_u1/2∫_y^1-2s(x)(θ |u|^2 + |∇ u|^2 )  x y+ μ/2∫_Ω s(x)^2 |u(x,0) - f| ^2   x,

where 𝒞={(x,y)∈Ω× (0,+∞)}, 0<θ≪ 1,   μ>0 is a given constant, and u(x,0) stands for the trace of u at Ω×{0}.  For the solution u of this optimization problem, we expect that u(·,0)≃ u_true. This problem is the focus of the present paper and the full motivational link between the problems (<ref>) and (<ref>), is given in the following section.

We next summarize the main contributions of this paper and discuss how the paper is organized in what follows. 

In section <ref>, we provide a rigorous motivation for the study of problem (<ref>) by showing that  (<ref>) represents a generalization of (<ref>), after using the Caffarelli-Silvestre extension, in the case of s non-constant. For a constant s ∈ (0,1) the definition of fractional Laplacian in terms of the
Caffarelli-Silvestre or the Stinga-Torrea extension is by now well-known. However such a result remains open 
when s(x) ∈ [0,1] for x∈Ω. Towards this end we emphasize that 
we have formulated our problem (<ref>) on an unbounded domain =Ω× (0,+∞), and not directly over Ω.

In order to handle varying 
s in space, for almost all (f.a.a) x ∈Ω we assume s(x) ∈ [0,1], and consider weighted Sobolev with weights w(x,y)=y^1-2s(x) on 𝒞. This class of Sobolev spaces is studied in section <ref> and  we establish also two fundamental results. The first one is that since we allow 
s(x) =0, for some x∈Ω, the weight y^1-2s(x) is not, in general, in the Muckenhoupt class A_2
(cf. Proposition <ref>). Due to the lack of this A_2 property we cannot use some of the 
existing machinery on weighted Sobolev spaces. In particular, the density of smooth functions is not always given in weighted spaces in the absence of the properties inherited by the Muckenhoupt class.  The second important result of section <ref> characterizes the trace space for this class of weighted Sobolev spaces, this is given in Theorem <ref>, and it is established that the trace space embeds continuously on the s-weighted L^2 space. In addition, we  give examples of functions with discontinuous traces in Example <ref>. 

In section <ref>, we establish properties for the variational problem of interest (<ref>), and a procedure for the selection of the function s. We finalize the paper with section <ref>, where we develop a discretization scheme, and provide numerical tests that positively compared with other methods  for image reconstruction. In particular, in section <ref> we introduce a finite element method for (<ref>) on a bounded domain _τ := × (0,τ). 
Since our targeted application is image denosing so we first state an 
Algorithm to determine s in Algorithm <ref>. 
We apply this approach to several prototypical examples in section <ref>.
In all cases we obtain better results when compared to the Total Variation (TV) approach.






 §.§ Motivation and some notation



Our point of departure for motivation to study the proposed variational model (and associated function space properties) is Figure <ref> 



where the left panel 
shows a noisy image, the middle panel shows the reconstruction using the Total Variation (TV) model, explained in what follows,
(we use piecewise linear finite element discretization for the TV model and we stop
the algorithm when the relative error between the consecutive iterates is smaller than 1e-8.), and the rightmost panel shows the reconstruction using 
the new framework introduced in this paper and briefly explained in this section. The Rudin-Osher-Fatemi (ROF) model, also called the  
TV model, was introduced  in <cit.> and is given by


    min_u ∈ BV(Ω)∫_Ω|D u| + ζ/2u-f^2_L^2(Ω) ,


where Ω⊂ℝ^N with N ≥ 1 
is a open bounded Lipschitz domain, and f ∈ L^2(Ω) is given by f=u_true+ξ, where ξ is the “noise”. Moreover, 
ζ > 0 is the regularization parameter, and  
 ∫_Ω|D u| denotes the total variation seminorm of u on Ω 
; see <cit.>.  The parameter ζ > 0 in Figure <ref> is chosen so that it maximizes a weighted sum of the Peak Signal to 
Noise Ratio (PSNR) and the Structural Similarity index (SSIM); as in <cit.>. These two indexes measure with different criteria how close the solution u_ tv to (<ref>) is to u_true; see <cit.>.


The fundamental question in such applications is:  Is it possible to capture the sharp 
transitions (edges) across the interfaces while removing the undesirable noise from remainder
of the domain. 


Such a question is not limited to image denoising but is fundamental to many applications
in science and engineering. For instance multiphase flows (diffuse interface models),
fracture mechanics, image segmentation etc. 

Even though the ROF model has been extremely successful in practice, two main drawbacks are observed in reconstructions: 1) loss of contrast is always present 2) some corners tend to be rounded  (cf. Figure <ref>, 
middle panel). 


The problem (<ref>) proposed in <cit.> is obtained if  ∫_Ω|D · | is replaced by 1/2(-Δ)^s/2·_L^2(Ω)^2  
where (-Δ)^s denotes the fractional powers of the Laplacian with zero Neumann boundary
conditions, for s∈ (0,1). Solutions to (<ref>) are given in H^s(Ω), the fractional Sobolev space given by

    H^s(Ω)={ v∈ L^2(Ω): ∫_Ω∫_Ω|v(x)-v(y)|^2/|x-y|^n+2s x y<+∞}.

The minimization problem 
(<ref>) is appealing, in contrast to (<ref>),  as the Euler-Lagrange equations are linear in u, i.e., if u solves (<ref>), then equivalently it solves


    ⟨ (-Δ)^s u, v⟩+ ζ (u-f, v)  = 0,


for all v∈ H^s(Ω). Both (<ref>) and hence (<ref>) have unique solutions 
in H^s(Ω). This can be easily deduced from the definition of Neumann fractional Laplacian 
and the equivalence between the spectral and H^s(Ω)-norms (cf. <cit.> and <cit.>). 
The H^s(Ω) space further provides flexibility in terms of the regularity by choosing s appropriately. 

The Caffarelli-Silvestre (or Stinga-Torrea) extension technique establishes that if u∈ H^s(Ω) solves (<ref>), then there exists 𝒰∈ H(;y^1-2s) where H(;y^1-2s) is the weighted H^1 Sobolev space on := Ω× (0,∞) with weight w(x,y)=y^1-2s, (x,y)∈𝒞 (see the following section for more details), such that 𝒰=u, where  operator is the restriction of the trace map to Ω, and 𝒰 solves

    ∫_ y^1-2s∇𝒰·∇𝒲    x  y
      +  d_s ζ∫_Ω  (𝒰-f)    𝒲  x=0,

for all 𝒲∈ H(;y^1-2s), where  d_s=2^1-2sΓ(1-s)/Γ(s). Therefore, for s a constant such that  s(x)=s∈ (0,1) for all x∈Ω, θ=0 and ζ := μ s^2 d_s^-1, if 𝒰 solves (<ref>), then it equivalently solves (<ref>). This establishes the connection of  (<ref>) with (<ref>), and hence with (<ref>).

Further, if s(·) is not a constant and s:Ω→[0,1] (with some additional assumptions made explicit in the following section), and if 𝒰∈ H(;y^1-2s(x)) solves (<ref>), then 

    ∫_ y^1-2s(x)(∇𝒰·∇𝒲 + θ𝒰𝒲)     x  y
      +  μ∫_Ω s(x)^2 (𝒰-f)    𝒲  x=0 ,

for all 𝒲∈ H(;y^1-2s(x)); this is rigorously done in section <ref> where the key ingredient is the characterization of  of H(;y^1-2s(x)). It is clear that for θ≪1, the above problem (whenever well-posed) represents a generalization of (<ref>) and hence of (<ref>).

Now we are in shape to explain Figure <ref>. First, note that for s≃ 0, we expect low regularity of solutions, and for s≃ 1 we expect higher regularity of solutions: In fact, for constant 0<s_1< s_2<1, we observe H^s_2(Ω)⊂ H^s_1(Ω).  Then, after a rough identification of edges, we choose s small on edges, and large on other regions; the exact procedure is provided in section <ref>. The right image in Figure <ref>, is obtained via this procedure. 






§ WEIGHTED SOBOLEV SPACE WITH  VARIABLE S AND THEIR TRACES



Let Ω⊂ℝ^N with N ≥ 1 be a non-empty, open, bounded domain 
with Lipschitz boundary ∂Ω. We denote by := Ω× (0,∞) a semi-infinite 
cylinder with boundary ∂_L := ∂Ω× [0,∞). With s:Ω→ℝ measurable and s(x) ∈ [0,1] for almost all x ∈Ω, we define 

    δ(x) := 1-2s(x) ∈ [-1,1],      and       w(x,y):=y^δ(x),

so that

    w, w^-1∈ L^1_loc().

We denote L^2(;y^δ(x)):=L^2(;w), where

    L^2(;y^δ(x)):={v:→ℝ: v  is measurable and ∫_ y^δ(x) |v|^2  x y<+∞},

which is a well-defined weighted Lebesgue space.

For u smooth, define the (extended real-valued) functional ·_H  as

    u_H 
      := ( ∫_ y^δ(x)(|u|^2 + |∇ u|^2 )  x y )^1/2.

Note that if u_H <+∞, the value of s controls how singular the function u is in the neighborhood of Ω×{0}. Define 𝒜_α(s):={x∈Ω: s(x)=α,   a.e.}, suppose that for α=1 its Lebesgue measure satisfies |𝒜_1(s)|>0, and that  u is a constant in 𝒜_1(s)× (0,h), then

    u^2_H≥ |u|^2|𝒜_1(s)|∫_0^h y^-1 y,

so that u is zero in Ω× (0,h). On the other hand, u is allowed to have a more singular behaviour on 𝒜_0(s). 


Consider the set 

    _Ω=(Ω×{0})∪,

that is, _Ω is the open cylinder  together with the Ω cap,
and denote by H, H_0, and W, or equivalently by  H(;y^δ(x)), H_0(;y^δ(x)), and W(;y^δ(x)), the spaces 
H(;w), H_0(;w) and W(;w) which are defined as


  

   H(;w) is the closure of the set K:={w∈ C^∞(): w_H<+∞} with respect to the ·_H norm.

  

   H_0(;w)  is the closure of the set  K_0:={w∈ C_c^∞(_Ω): w_H<+∞} with respect to the ·_H  norm. 

  

   W(;w) is the space of maps u∈ L^2(;y^δ(x))∩ L_loc^1() with distributional gradients ∇ u that satisfy |∇ u| ∈ L^2(;y^δ(x))∩ L_loc^1().


  


A few words are in order concerning the definition of H given the slight difference from the existing literature: Usually, the definition of H is given over the completion (of finite energy functions) of  C^∞(), in contrast to C^∞(), with respect to the H norm; see  <cit.>. Provided that ∂ is regular enough,  is bounded, and w,w^-1 are bounded above and below (ϵ>0 away from zero) on a neighborhood of ∂, then the closure in the definition of H can be equivalently taken with C^∞() or C^∞(), indistinctly. Since in our case the two last conditions are not satisfied, we consider H the way it is defined here. 


First note that since K_0⊂ K, which follows from C_c^∞(_Ω)⊂ C^∞(), we have that H_0⊂ H. The notation of “H_0” comes from the following: If v∈ K_0, then v|_∂Ω×{0} =0, so that restrictions of elements of H_0 to Ω have “zero boundary conditions”. In particular, in the case s∈ (0,1) is constant, the latter remark can be taken out of quotation marks and considered in the sense of the trace; see <cit.>. The  condition (<ref>) implies (see <cit.>) that  the space W
is a proper weighted Sobolev space endowed with the H-norm.  While it holds that H(;w̃)⊂ W(;w̃), for general weights w̃ and domains , it does not necessarily hold that H(;w̃)=W(;w̃). For w̃=1, the equality holds provided conditions on  are given, e.g., if  is bounded and has a Lipschitz boundary. A sufficient condition for H(;w̃)=W(;w̃), when w̃ and  are benign enough, is that w̃ belongs to the A_2()  Muckenhoupt class, that is

    (1/Q∫_Qw̃ x y)(1/Q∫_Qw̃^-1 x y)≤ M,

for some M>0 and all open cubes Q⊂, see <cit.>.  Further, note that  (<ref>), does not imply (<ref>). Additionally, if 0<ϵ≤ s≤ 1-ϵ a.e., then (x,y)↦ y^1-2s is a 
Muckenhoupt A_2() weight.   However, we are particularly interested in cases where s satisfies 

    ess inf_x∈Ω s(x) =0 ,

as this will allow almost perfect approximation of functions with jump discontinuities.
Remarkably enough condition (<ref>) prevents w to be a Muckenhoupt weight in 
general as we show next.



Suppose that for x_0∈Ω, s(x)≃ |x-x_0|^q (locally) on x_0 for some q>0, that is, there exists positive constants R_0, M_0, m_0 such that

    m_0 |x-x_0|^q≤ s(x)≤  M_0|x-x_0|^q,

for all x such that |x-x_0|≤ R_0. Then, w∉ A_2() where w(x,y):=y^1-2s(x).


 First step.
Suppose first that x_0=0∈Ω, s(x)=|x|^q for x∈ B_R_0(0), and let Q_R=B_R(0)× (0,y_0) for 0<R≤ R_0<1 and 0<y_0<1. 

Using cylindrical coordinates we observe

    ∫_Q_Rw   x y   ≥  C_1 ∫_0^R∫_0^y_0r^N-1y^1-2 r^q y  r≥ C_2  ∫_0^R r^N-1y_0^2-2 r^q r 
       ≥ C_3 R^N y_0^2,

where C_1, C_2, and C_3 are positive and independent of R and y_0.  Similarly, 

    ∫_Q_Rw^-1 x y   ≥ C_1 ∫_0^R∫_0^y_0r^N-1y^-(1-2 r^q) y  r≥ C_4  ∫_0^R r^N-1-qy_0^2 r^q r 
       ≥ C_5 R^N-q y_0^2R^q,

where also  C_4, and C_5 are positive and independent of R and y_0. Since |Q_R|≃ R^N y_0, we have

    (1/|Q_R|∫_Q_Rw   x y)(1/|Q_R|∫_Q_Rw^-1 x y)≥ C_6  y_0^2R^q/R^q,

for some C_6>0 independent of R. Taking R↓ 0, we have that the L.H.S. expression is unbounded.


Second step: Consider x_0≠ 0, s(x)≠ |x|^q, and show that the above inequality holds similarly for a  family of cubes R↦Q̃_R.

If x_0≠ 0 and s(x)=|x|^q, then a linear change of coordinates can be performed and the proof above still holds, so the choice of x_0=0 is without loss of generality. 

Additionally, let Q̃_R⊃ Q_R be defined as Q̃_R=S_R(0)× (0,y_0) where S_R(0) is the smallest square that contains B_R(0). Hence,  the quotient |S_R(0)|/|B_R(0)| is independent of R and for c:=|Q̃_R|/|Q_R| we observe

    (1/|Q̃_R|∫_Q̃_Rw   x y)(1/|Q̃_R|∫_Q̃_Rw^-1 x y)
               ≥1/c^2(1/|Q_R|∫_Q_Rw   x y)(1/|Q_R|∫_Q_Rw^-1 x y).

Hence, by taking R↓ 0, we have that  (x,y)↦ y^1-2|x|^q is not in A_2(). 

Finally, for s(x)≠ |x|^q but when (<ref>) holds true, we similarly have 

    (1/|Q|∫_Qw   x y)(1/|Q|∫_Qw^-1 x y)
               ≥c̃(1/|Q|∫_Qw̃ x y)(1/|Q|∫_Qw̃^-1 x y),

for some c̃>0, where w̃(x,y):=y^1-2|x|^q and w(x,y):=y^1-2s(x), i.e., w∉ A_2(). 
 

Although the above result is quite general with respect of the rate of decrease of s, at this point, we do not know if there are continuous functions s with zeros within the domain Ω, for which (x,y)↦ y^1-2s(x) is an element of A_2().


 The study of the trace space of H is of utmost importance in what follows. In particular, we are interested in the restriction of the trace operator to Ω. The paper <cit.> studies the trace of H when s 
is independent of x and away from 0. However, that approach which is based on 
<cit.> cannot be applied here. We next prove a trace characterization result, the proof is inspired by <cit.> where the authors 
consider bounded domains and weights of distance to the boundary type.


Let s be such that its zero set 𝒜_0(s):={x∈Ω: s(x)=0,   a.e.}  has zero measure. Then, there exist a unique bounded linear trace operator 

    : H, H_0 → L^2(Ω; (1-δ(x))^2)=L^2(Ω; s(x)^2)

such that (u) = u|_Ω×{0}, for all u∈ H∩ C^∞(), and hence also for all u ∈ H_0∩ C_c^∞(_Ω).




Consider u∈  H such that u∈ C^∞(). Let (x,y)∈ and x∉𝒜_0(s), then it  follows that 
    u(x,y)=u(x,0)+∫_0^1  y D_N+1 u(x,ty)  t,

where D_N+1 corresponds to the partial derivative with respect to the N+1 coordinate. Then, multiplication by y^δ(x)/2 and integrating from 0 to σ>0 with respect to y leads to

    |u(x,0)|∫ _0^σ y^δ(x)/2 y≤∫ _0^σ  |u(x,y)|y^δ(x)/2 y+ ∫ _0^σ∫_0^1   y |D_N+1 u(x,ty)| y^δ(x)/2 y t.

Note that δ(x)/2≤ 1/2 and suppose without loss of generality 0<σ≤ 1, then y^1/2≤ y^δ(x)/2 for all y∈ (0,σ). Therefore,

    2/3σ^3/2≤∫ _0^σ y^δ(x)/2 y,

and hence,

    |u(x,0)|≤ I_1+I_2,

for 

    I_1:=3/2σ^-3/2∫ _0^σ  |u(x,y)|y^δ(x)/2 y ,    
    I_2=3/2σ^-3/2∫_0^1∫ _0^σ   y |D_N+1 u(x,ty)| y^δ(x)/2 y t ,

and where we have used Tonelli's theorem to switch the order of integration in the expression that orginates I_2. 


Hölder's inequality implies that, for I_1, we have

    I_1   ≤3/2σ^-3/2 (∫ _0^σ y)^1/2(∫ _0^σ  |u(x,y)|^2y^δ(x) y)^1/2≤3/2σ^-1(∫ _0^R  |u(x,y)|^2y^δ(x) y)^1/2,

for any R≥σ. 

Further, for I_2 consider the change of variable z=ty and again by Hölder's inequality we observe:

    I_2≤3/2σ^-1/2∫_0^1∫ _0^σ t  |D_N+1 u(x,z)| z^δ(x)/2 t^-1-δ(x)/2 z t
       ≤3σ^-1/2/2 (∫_0^1 (t^-1+δ(x)/4)^2 t )^1/2 (∫_0^1 ( ∫ _0^σ t  |D_N+1 u(x,z)| z^δ(x)/2 t^1+δ(x)/4-1-δ(x)/2 z)^2 t )^1/2
       ≤3/2σ^-1/2 (2/1-δ(x) )^1/2 (∫_0^1 t^-3+δ(x)/2( ∫ _0^σ t  |D_N+1 u(x,z)| z^δ(x)/2 z)^2  t )^1/2
       ≤3/2 (2/1-δ(x) )^1/2 (∫_0^1 t^-1+δ(x)/2( ∫ _0^σ t  |D_N+1 u(x,z)|^2 z^δ(x) z)  t   )^1/2.

Then, for any R≥σ, we have

    I_2   ≤3/2 (2/1-δ(x) )^1/2 (∫_0^1 t^-1+δ(x)/2 t   )^1/2( ∫ _0^R  |D_N+1 u(x,z)|^2 z^δ(x) z    )^1/2
       ≤3/1-δ(x)( ∫ _0^R  |D_N+1 u(x,z)|^2 z^δ(x) z )^1/2.



Therefore, from (<ref>), we obtain that for any R≥σ

    |u(x,0)|(1-δ(x))
         ≤3/2σ^-1 (1-δ(x)) (∫ _0^R  |u(x,y)|^2y^δ(x) y)^1/2+3 ( ∫ _0^R  |D_N+1 u(x,z)|^2 z^δ(x) z )^1/2
         ≤ 3σ^-1(∫ _0^R  |u(x,y)|^2y^δ(x) y)^1/2+3 ( ∫ _0^R  |D_N+1 u(x,z)|^2 z^δ(x) z )^1/2
         ≤ 6σ^-1(∫ _0^R |u(x,y)|^2y^δ(x) y + ∫ _0^R  |D_N+1 u(x,z)|^2 z^δ(x) z )^1/2.

Integration over Ω with respect to x of the above squared inequality leads to

    ∫_Ω |u(x,0)|^2(1-δ(x))^2 x
               ≤ M_σ^2 (∫_Ω∫ _0^R |u(x,y)|^2y^δ(x) y + ∫_Ω∫ _0^R  |D_N+1 u(x,z)|^2 z^δ(x) z )
               ≤ M_σ^2 (∫_Ω∫ _0^R |u(x,y)|^2y^δ(x) y + ∫_Ω∫ _0^R  |∇ u(x,z)|^2 z^δ(x) z ),

where M_σ= 6σ^-1, R≥σ, with σ∈ (0,1] and for an arbitrary u∈  H such that u∈ C^∞(). 


Since R≥σ is arbitrary, we have that for an arbitrary u∈  H∩ C^∞(), we have

    u(·,0)_L^2(Ω; (1-δ(x))^2)≤ M_σu_H.

The operator  is the above extension to H: Since H is the closure of  
C^∞()
with respect to the energy norm, we can take σ=1 with  M_σ=6, and obtain

    u_L^2(Ω; (1-δ(x))^2)≤ 6 u_H,

which completes the proof, since   (1-δ(x))^2=4s(x)^2. The result for H_0 is analogous. 


  



The behavior of s controls the local regularity on space H, in particular, H contains piecewise constants in any dimension for appropriate choices of s as we see in the following example. 

There exists non-constant s:Ω→ℝ maps  such that H contains piecewise smooth functions. For the sake of simplicity, let Ω=B_1/2(0), define s(x)=ŝ(x_1), where x=(x_1,x_2,…,x_N), and such that 0<δ≤ŝ(x_1)≤κ<1/4 for all |x_1|<ϵ, and 0<δ≤ŝ(x_1) for all x_1.  Let

    w(x,y):={[                 0,       x_1< -√(y) ;; 1/2√(y)(x_1+√(y)),  √(y)≤ x_1≤-√(y) ;;                 1,          √(y)<x_1. ].

so that 

    |∇ w(x,y)|^2={[                  0,        x_1< -√(y) ;; 1/4y+x_1^2/16 y^3 ,   √(y)≤ x_1≤-√(y) ;;                  0,           √(y)<x_1. ].

Define u(x,y):=η(y)w(x,y) for (x,y)∈, with η∈ C^∞_c(ℝ), such that η(y)=1 if y∈ [0,R^*] for some R^*, and therefore,

    ∫_B_ϵ(0)∫_0^1y^1-2s|∇ w(x,y)|^2 x  y ≤∫_B_ϵ(0)∫_0^1y^1-2κ|∇ w(x,y)|^2 x  y<+∞.

Then, u∈ W and smoothing argument involving mollification shows that u∈ H, and u = χ_{x_1≥ 0}. 







§ THE OPTIMIZATION PROBLEM 


 Throughout this section we suppose that f∈ L^2(Ω; s(x)^2), s(x) ∈ [0,1] 
for almost all x ∈Ω, and that  𝒜_0(s):={x∈Ω: s(x)=0,   a.e.}  has zero measure. 

Given a regularization parameter μ > 0, and parameter 
θ > 0 we consider the following variational problem 

    min_u∈ H(𝒞; y^1-2s(x)) J(u,s):=∫_y^1-2s(x)(θ |u|^2 + |∇ u|^2 )  x y+ μ/2∫_Ω s(x)^2 | u - f| ^2   x.


 




In case the weights w are of class A_2, we can set 
 θ = 0, given that  the Poincaré type 
 inequalities are available in this case.  However, the weights y^1-2s(x) in our case are not in 
 A_2 since we allow s(x) =0 for some x∈Ω. Such Poincaré type inequalities on H(;y^1-2s(x)) are not known yet.
 


The model (<ref>) provides a completely new approach to approximate nonsmooth
functions f with u. Such a situation often occurs in inverse 
problems where given f one is interested in a reconstruction u which is 
close to f. For instance in 
image denoising problems where f represents a noisy image, given by f=u_true+ξ where u_true is the true target to recover, ∫_Ωξ=0, and ∫_Ω|ξ|^2=σ^2 for a known  σ>0. The first term in (<ref>) is the regularization and the second term  is the so-called data fidelity, together they  ensure that the reconstruction u is close 
to u_true on Ω. 


For a fixed s and given f we next state existence and uniqueness of solution to (<ref>). This follows immediately from Theorem <ref>.


There exist a unique solution to (<ref>). 


Existence follows from application of direct methods of calculus of variations, the norm definition on H and the Theorem <ref>. Uniqueness follows directly from convexity arguments.

The first order necessary and sufficient optimality conditions for the minimization 
problem (<ref>) are given by: 

    ∫_ y^1-2s(x)(∇ u ·∇ w + θ u w)  x  y
      +    μ∫_Ω u     w  s(x)^2 x  
       = μ∫_Ω f  w  s(x)^2 x,    w∈ H .

The first step for computational implementation of the problem above requires to solve the problem in a bounded domain. Notice that all results above are valid if we replace the unbounded domain  with boundary ∂_L 
by a bounded domain _τ = Ω× (0,τ) with 0< τ < ∞ and boundary ∂_L_τ = (∂Ω× [0,τ]) ∪Ω×{τ}. In that case the problem (<ref>) becomes: find v ∈ H^τ := H(_τ;y^δ(x)) 


    ∫__τ y^1-2s(x)(∇ v ·∇ w + θ v w)  x  y
      +    μ∫_Ω v     w  s(x)^2 x  
       = μ∫_Ω f  w  s(x)^2 x,    w∈ H^τ .

It is known that for a constant s ∈ (0,1) the solution v to (<ref>) converges to u solving (<ref>) exponentially with respect to τ > 1, see 
<cit.>. Such exponential approximation is also expected in the case of a non-constant s. Indeed, as we discussed earlier, a constant s implies that 
u solves the fractional equation (<ref>). A relation between (<ref>)
and fractional PDE of type (<ref>) with s(x) instead of s is an open 
question.  We expect that studying (<ref>) or equivalently (<ref>) is a first step in establishing the definition of fractional 
(-Δ)^s(x) for x ∈Ω. 


For given s and f, Theorem <ref> shows existence of a unique solution to  (<ref>) and equivalently (<ref>). Nevertheless for our
model (<ref>) to be useful in practice we need to determine the function s. For this matter, consider the following important points that will lead to the selection procedure 
for s.


 
  * 
  Precise edge recovery. In example <ref>, a family of s functions    is identified so that u∈ H(𝒞; y^1-2s(x)), and u = χ_{x_1≥ 0}. This suggests that for the optimization problem (<ref>) to recover edges, or more precisely for u to have discontinuities in the same places where u_true has them  (for f=u_true+noise) the function s needs to be close to or equal to zero in these regions. This suggest to utilize a rough edge detection at first and then force s to be close zero on neighborhoods of these regions.

 
  * 
  Homogeneous/flat regions recovery. We proceed rather formally first. Suppose first  that θ=0, and  that  problem (<ref>) is posed not on H(𝒞; y^1-2s(x)), but on the Banach space of equivalence classes generated by the seminorm ∫_y^1-2s(x)  |∇ u|^2  x y. Consider s(x)=1 on a certain region Ω_0. Suppose that |∇ u|≥ e>0 on Ω_0× (0,y_0), then

    ∫_y^1-2s(x)  |∇ u|^2  x y≥ϵ∫_Ω_0× (0,y_0)y^-1 x y =+∞.

Hence, the solution to  (<ref>), needs to satisfy |∇ u|≃ 0 near Ω×{0}. This suggest that flat regions would be recovered if s=1 on that same region. However, if θ>0, then analogously, solutions to  problem (<ref>), would be forced to satisfy |u|≃ 0 near Ω×{0}. Hence, we consider 0<θ≪ 1 and s(x)≃ 1, but s(x)<1, on regions where homogeneous features are present so that u is not forced to be identically zero but |∇ u|≃ 0 is enforced there.  



We provide an algorithm based on the two points above in what follows. A more detailed 
bilevel optimization framework (where both u and s are optimization variables) 
will be considered in a forthcoming publication. 





§ NUMERICAL METHOD



We now focus on the discretization of the truncated problem (<ref>) and the selection of an appropriate s function.




 §.§ Discretization of (<ref>)



From hereon we will assume that Ω is polygonal/polyhedral Lipschitz. 
We recall that the results of previous sections remains valid if we replace the 
unbounded domain  with boundary ∂_L by _τ = Ω× (0,τ) 
with τ > 0 and boundary ∂_L = (∂Ω× [0,τ]) ∪Ω×{τ}. The Euler-Lagrange equations for the resulting problem on _τ are given
by (<ref>).

We begin by introducing a discretization for (<ref>). We will follow the notation from <cit.>. For a constant s such a discretization for v was  first 
considered in <cit.>. Let 𝒯_Ω={E} be a conforming and 
quasi-uniform triangulation of Ω, where E∈ℝ^N is an element that is 
isoparametrically equivalent to either to the unit cube or to the unit simplex in ℝ^N. 
We assume #𝒯_Ω∝ M^N. Thus, the element size h_𝒯_Ω 
fulfills h_𝒯_Ω∝ M^-1. 

Furthermore, let ℐ_τ={I_k}_k=0^K-1, where I_k = [y_k,y_k+1], 
is anisotropic mesh in [0,τ] in the sense that [0,τ]=⋃_k=0^K-1 I_k. For a constant s ∈ (0,1) we define the anisotropic mesh in y-direction as 

    y_k = (k/K)^γ , 
        k = 0, …,K,   γ > 1/s .

This choice is motivated by the singular behavior of the solution towards the boundary Ω for a constant s. In that case anisotropically refined meshes are preferable as these can be used to compensate the singular effects <cit.>. 
In all our implementations we will choose a fixed constant s in (<ref>).
 
We construct the triangulations 𝒯_τ of the cylinder 𝒞_τ as tensor product triangulations by using 𝒯_Ω and ℐ_τ. Let 𝕋 denotes the collection of such anisotropic meshes 𝒯_τ. For each 𝒯_τ∈𝕋 we define the finite element space 𝕍(𝒯_τ) as 

    𝕍(𝒯_τ):={V∈ C^0(𝒞_τ):V|_T∈𝒫_1(E)⊕ℙ_1(I) ∀ T=E× I∈𝒯_τ}.

In case ∂Ω has Dirichlet boundary conditions we define our finite element space as 𝕍_0(𝒯_τ) = 𝕍(𝒯_τ) ∩{V : V|_∂_τ = 0}, i.e., functions with zero boundary conditions. 
In case E is a simplex then 𝒫_1(E)=ℙ_1(E), the set of polynomials of degree at most 1. If E is a cube then 𝒫_1(E) equals ℚ_1(E), the set of polynomials of degree at most 1 in each variable. In our numerical illustrations we shall work with simplices. 

We define the finite element space for s as   

    𝕊(𝒯_Ω) := {S ∈ L^∞(Ω) : S|_E ∈ℙ_0(K)
         E ∈𝒯_Ω}

which is a space of piecewise constant functions defined on 𝒯_Ω. The discrete 
version of (<ref>) is then given by: 
Find V ∈𝕍(𝒯_τ) 


    ∫__τ y^1-2S(x)   (∇ V ·∇ W + θ V W)  x  y
      + μ∫_Ω V     W  S(x)^2 x  
       = μ∫_Ω f  W  S(x)^2 x,    W∈𝕍(𝒯_τ) ,


for S ∈𝕊(𝒯_Ω). We compute the stiffness and mass matrices 
in (<ref>) exactly. The corresponding forcing boundary term is computed by a quadrature formula which is exact for polynomials up to degree 5. For a given S the resulting discrete linear system (<ref>) is solved 
using the Preconditioned Conjugate Gradient (PCG) method with a block diagonal preconditioner. 





 §.§ Parameter selection



In view of what was stated about θ in section <ref>, we have set θ = 10e-10 in all our examples. 
We let 
τ = 1+1/3(#𝒯_Ω) ,
this choice is motivated by the fact that for a constant s such a τ 
balances the finite element approximation on _τ and the truncation error 
from  to _τ <cit.>. 
The number of points in the y-direction is taken to be K=20. We use
a moderately anisotropic mesh in the y-direction by setting s = 0.32 in (<ref>).
Our experiments suggest that the results are stable under reducing θ or s
or increasing K further. 

It then remains to specify the constant μ and the function S in order to realize 
(<ref>). We have observed that μ only affects the 
“contrast" or magnitude of V. Nevertheless one can determine μ using 
the well established techniques such as L-curve method <cit.>. In our case 
we choose a fixed μ for each example (no optimization was carried to select μ). 
On the other hand selecting S which is a function is much more delicate. One option is to use 
a bilevel strategy as in <cit.> to 
determine both S and V in an optimization framework. This is a part of our future work. In 
this paper we propose a different approach in Algorithm <ref>.

Even though our targeted application is image denoising the approach we present is general 
enough to be applicable to a wider range of applications. We first notice that a typical 
image is given on a rectangular grid (pixels). Since for (<ref>) we are working on simplices we first need to interpolate the given image from the rectangular grid to a 
simplicial mesh. This is a delicate question especially given the fact that typically we only 
have access to a noisy image. Before we interpolate the noisy image onto a simplicial mesh
we perform an intermediate step. 

We solve ROF model (<ref>) using <cit.> and ζ>0 chosen to a fixed value for all examples. We stop the algorithm 
when the relative difference between two consecutive iterates is smaller 
than the given tolerance tol_tv. We select a mild tolerance tol_tv 
in this step as we want to preserve the sharp features but still remove a certain noise.
We call the resulting solution as u_ tv. We then generate a piecewise linear Lagrange interpolant I_𝒯_Ω u_ tv of u_ tv. 
 

We next evaluate 
I_𝒯_Ω u_ tu on the simplicial mesh. In order to reduce the computational
cost we use Adaptive Finite Element Method (AFEM) to generate the appropriate mesh. In the
nutshell, we start with a coarse mesh 𝒯_Ω = {E} where E is an element in
𝒯_Ω. For each E∈𝒯_Ω we then evaluate gradient of 
I_𝒯_Ω u_ tv on E and we denote this gradient by 
∇ I_𝒯_Ω u_ tv|_E. We use this gradient to define an edge 
indicator function, we call this edge indicator function as the estimator on E. Based on 
a marking strategy we then mark a subset of elements in 𝒯_Ω. Subsequently
we perform the mesh refinements. We execute this loop N_ refine times. 

Finally we set S so that it is close to 0 at the sharp edges (large gradient) 
in the image and close to 1 away from sharp edges. Intuitively smaller the S the lesser
the smoothness and otherwise. 




Once we have S then we can immediately solve the linear equation (<ref>) 
for V. 






§ NUMERICAL EXAMPLES



In this section we illustrate the proposed scheme with the help of several examples. We consider f given by f=u_true+ξ where u_true is the true image, and ξ is noise with the following properties: ∫_Ωξ=0, and ∫_Ω|ξ|^2=σ^2 for a known  σ>0.
In all cases we find S by using Algorithm <ref> and then solve the 
resulting linear equation (<ref>) for V. We call V as
the reconstruction. We also compare our reconstruction V with the 
reconstruction we obtain by solving (<ref>) with tol^*_ tv = 10e-8. 
In this case we choose the parameter ζ so that a normalized weighted sum of Peak Signal to 
Noise Ratio (PSNR) and the Structural Similarity (SSIM) index is maximized, as in  
section 4.2 of <cit.>. That is, we compare our results to the ones of the TV model with an optimized parameter ζ>0. Note that this parameter is not the one used in step <ref> from algorithm <ref>.





 §.§ Example 1: circle, triangle and square



In our first example we consider an image with a circle, triangle and a square as shown in
Figure <ref> (top row left). Notice that in this case the right pointing edge of
the triangle does not align with the grid. 
We consider two different noise levels. In both cases the noise is normally distributed with mean 0 and standard deviations 
0.10 and 0.15 respectively. The results are shown in Figures <ref> and 
<ref> for standard deviations 0.10 and 0.15, respectively. In both cases
we first compute S by using Algorithm <ref> where we have set:
ζ = 0.2, 
tol_ tv = 10e-4, 
N_ refine = 8,
λ = 300, 
β = 0.99,
ν = 200
in Algorithm <ref>. We further set μ = 8050 in (<ref>).
We then solve for V using PCG. We call V as our reconstruction.  

The top row of Figure <ref> shows the original and noisy images (left to 
right). The middle row shows u_ tv and S obtained using Algorithm <ref>. 
In the bottom row we compare the results using the TV approach (left) and V (right) 
computed using
our approach by solving (<ref>). Notice that Figure <ref> is simply obtained by viewing Figure <ref>, 
in particular, the noisy, reconstruction using TV and V panels, from a different angle.
Similar description holds for Figure <ref>. Figure <ref> 
is again viewing Figure <ref> from a different angle.

As we can notice TV tends to round up the corners (cf. Figure <ref> (middle)). On the other hand as our theory predicted we can truly capture the edges in 
V (cf. Figure <ref> (right)). This is further corroborated by Table <ref> where we have shown a comparison between PSNR and SSIM for these two 
approaches for two different standard deviations. 
















 §.§ Example 2: stripes



In our second example we consider 6 stripes with intensities equal to 0.8, 0.7, 0.6, 0.5, 0.4 
and 0.35 (from left to right), respectively (cf. Figure <ref> (top row left)). As in the previous example
we again consider two additive noise levels with mean 0 standard deviations 0.1 and 0.15 respectively. Our results are shown in Figures <ref> and <ref> respectively. Here again we first find S by using Algorithm <ref> using the parameters
ζ = 0.2, 
tol_ tv = 10e-4, 
N_ refine = 8,
λ = 15, 
β = 0.99,
ν = 100. 
We further set μ = 2900 in (<ref>).
We then solve for V. We call V as our reconstruction. 
A comparison between PSNR and SSIM is shown in Table <ref>.
As we noticed in the previous example we again obtain visually almost
perfect reconstruction using our approach. 













 §.§ Example 3: cameraman



In the first two examples we considered synthetic images. In our
final example we consider a more realistic situation. We consider
a prototypical image of the cameraman (cf. Figure <ref>). As in the previous examples 
we again consider two additive noise levels with standard deviations 0.1 and 0.15 respectively. Our results are shown in Figures <ref> and <ref>, respectively

At first we apply the Algorithm <ref> to find S. Here we have set the underlying parameters as 
ζ = 0.2, 
tol_ tv = 10e-4, 
N_ refine = 8,
λ = 0.7, 
β = 0.99,
ν = 20.
We further set μ = 10^4 in (<ref>).
We then solve for V and we call V as our reconstruction. 
A comparison between PSNR and SSIM is shown in Table <ref>.
As we noticed in the previous examples we again obtain better reconstructions using our approach.




§ CONCLUSION AND FURTHER DIRECTIONS

A new variational model associated to the fractional Laplacian was introduced. In particular, we have identified a weighted Sobolev space with respect to the weight w=y^1-2s(x) appropriate for the treatment of the problem. We have shown that in general these weights are not of Muckenhoupt type, and have established that the trace space embeds in an s-weighted Lebesgue space. We have provided a discretization method for the full problem, and an algorithm for its resolution that also builds a selection procedure for s. The full scheme is advantageous when it comes to recovery of discontinuous features, details, homogeneous regions, and also for contrast preservation, in data perturbed by additive noise.

Future research directions are multiple, we enumerate some of them. 

1) The study of the optimization problem (<ref>), seems to be the first step for the identification of a possible definition for (-Δ)^s(x). Such a task does not seem directly approachable via spectral or functional calculus points of view.

2) Full characterization of the trace space for H and H_0. We have identified the embedding of H into L^2(Ω;s(x)^2), but it would be of interest to understand the Sobolev regularity of H.

3) Differentiability and stability properties of the solution to (<ref>) with respect to s. For the optimal selection of s,  
it is required to identify a topology over an admissible set for s, so that solutions to (<ref>), are stable with respect to perturbations. This seems like a complex task where the usual obstacles from homogenization and convergence of differential operators are present. In addition, difficulty is increased as the solutions to (<ref>) belong to a state space depending on s as well. The differentiability issue is further more complex, but such a study will be the first step in establishing stationarity systems useful for implementation in the optimal selection of s.

4) The extension of the presented methods to a general class of inverse problems. Problem (<ref>), admits the following generalization

    min_u∈ H(𝒞; y^1-2s(x))∫_y^1-2s(x)(θ |u|^2 + |∇ u|^2 )  x y+ μ/2∫_Ω s(x)^2 |K( u) - f| ^2   x,

where K is a bounded linear operator on L^2(Ω; s(x)^2). This would allow to deal with the problem of finding y such that Ky=f. In this case, the choice of s could be simply associated with regions of the domain Ω where it is more important to recover y more accurately.
















§ ACKNOWLEDGEMENT



We thank Gunay Dogan for point us to the edge error indicator and Pablo Raul Stinga
for sharing the reference <cit.>.



plain
