Sobolev spaces with non-Muckenhoupt weights, fractional elliptic operators, and applications
SSOBOLEV SPACES WITH NON-MUCKENHOUPT WEIGHTS,FRACTIONAL ELLIPTIC OPERATORS, AND APPLICATIONS ∗ HARBIR ANTIL † AND
CARLOS N. RAUTENBERG ‡ Abstract.
We propose a new variational model in weighted Sobolev spaces with non-standardweights and applications to image processing. We show that these weights are, in general, not ofMuckenhoupt type and therefore the classical analysis tools may not apply. For special cases ofthe weights, the resulting variational problem is known to be equivalent to the fractional Poissonproblem. The trace space for the weighted Sobolev space is identified to be embedded in a weighted L space. We propose a finite element scheme to solve the Euler-Lagrange equations, and for theimage denoising application we propose an algorithm to identify the unknown weights. The approachis illustrated on several test problems and it yields better results when compared to the existing totalvariation techniques. Key words.
Variable weights, non Muckenhoupt weights, new trace theorem, fractional Lapla-cian with variable exponent, image denoising.
AMS subject classifications.
1. Introduction.
In this paper we consider weighted Sobolev spaces where theweight w is not necessarely of Muckenhoupt type, and an associated variational modelwith concrete applications to image processing that shows advantageous features ofsuch weights. The particular weights w that we consider in this paper are closelyrelated to fractional Sobolev spaces of differentiability order s ∈ (0 ,
1) and the frac-tional spectral Laplacian ( − ∆) s , and are further motivated by considering a spatiallydependent order x (cid:55)→ s ( x ).Weighted Sobolev spaces have been a topic of intensive study for around 60 yearswith a variety of focuses; we refer the readers to the monographs [17, 15, 25] andreferences within, and further [28, 11, 27, 16, 20] for an introduction to the subject.In particular, a source of interest for such spaces is that they represent the correctsolution space for several degenerate elliptic partial differential equations ([17]) andproblems in potential theory ([25]). Recently, the topic has received a new impulsemainly associated with the extension result of Caffarelli-Silvestre; [5] (Stinga-Torrea;[24]) for fractional elliptic partial differential operators. In this setting, for example,the solution to ( − ∆) s u = f in the fractional Sobolev space H s (Ω) endowed with zeroboundary conditions, can be equivalently obtained as (the restriction to Ω of) thesolution of a PDE with a non-fractional elliptic operator in a weighted Sobolev spacewith weight w ( x, y ) = y − s in the extended domain C = { ( x, y ) ∈ Ω × (0 , + ∞ ) } .The type of weights w that have been more thoroughly studied correspond totwo main classes: 1) weights belonging to some Muckenhoupt class; see [11, 27] and2) composition of functions (mainly power functions) with the distance function to aparticular set; see [15, 20]. In these cases, several important questions like density of ∗ HA has been partially supported by NSF grant DMS-1521590. CNR has been supported via theframework of
Matheon by the Einstein Foundation Berlin within the ECMath project SE15/SE19and acknowledges the support of the DFG through the DFG-SPP 1962: Priority Programme “Non-smooth and Complementarity-based Distributed Parameter Systems: Simulation and HierarchicalOptimization” within Project 11 † Department of Mathematical Sciences, George Mason University, Fairfax, VA 22030, USA. [email protected] , ‡ Weierstrass Institute for Applied Analysis and Stochastics, Mohrenstr. 39, 10117 Berlin, Ger-many. [email protected]
AND Department of Mathematics, Humboldt-University of Berlin,Unter den Linden 6, 10099 Berlin, Germany. [email protected] a r X i v : . [ m a t h . O C ] M a r H. Antil and C.N. Rautenberg smooth functions and characterization of traces have been, at least partially, answered;see [28, 11, 15, 20]. In this paper we consider weights of the type w ( x, y ) = y − s ( x ) on C = { ( 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 el-liptic operators. The spike of research interest on problems with this type of operatorsis mainly due to the results from Caffarelli and collaborators; see [24, 5, 6]. Followingthe renowned results in [5], a large number of contributions on modeling, numericalmethods, applications, regularity results, different boundary conditions, and controlproblems, among others have been considered.Recently, for image denoising a model was proposed in [1]. The variational prob-lem can be formulated in our setting as(1.1) min u (cid:107) ( − ∆) s v (cid:107) L (Ω) + ζ (cid:107) v − f (cid:107) L (Ω) , where ( − ∆) s denotes the fractional power of the Laplacian with zero Neumann bound-ary conditions and s ∈ (0 , ζ > 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 toproblem (1.1). However, it is desirable, from the reconstruction point of view, thatthe regularity of the solution to (1.1) is low in places in Ω where edges or discontinu-ities are present in u true , and that is high in places where u true is smooth or containshomogeneous features. Hence, it is of interest to consider (1.1) where s : Ω → [0 ,
1] isnot a constant. The first main roadblock in letting s to be spatially dependent is thefact that there is no obvious way to define ( − ∆) s ( x ) . In fact we will not attempt todo so in this paper, we leave this as an open question. Instead of (1.1) we considerthe following variational problem:(1.2) min u ˆ C y − s ( x ) (cid:16) θ | u | + |∇ u | (cid:17) d x d y + µ ˆ Ω s ( x ) | u ( x, − f | d x, where C = { ( x, y ) ∈ Ω × (0 , + ∞ ) } , 0 < θ (cid:28) µ > u ( x, u at Ω × { } . For the solution u of this optimization problem,we expect that u ( · , (cid:39) u true . This problem is the focus of the present paper and thefull motivational link between the problems (1.1) and (1.2), is given in the followingsection.We next summarize the main contributions of this paper and discuss how thepaper is organized in what follows.In section 1.1, we provide a rigorous motivation for the study of problem (1.2)by showing that (1.2) represents a generalization of (1.1), after using the Caffarelli-Silvestre extension, in the case of s non-constant. For a constant s ∈ (0 ,
1) thedefinition 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 ourproblem (1.2) on an unbounded domain C = Ω × (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 , w ( x, y ) = y − s ( x ) on C . Thisclass of Sobolev spaces is studied in section 2 and we establish also two fundamentalresults. The first one is that since we allow s ( x ) = 0, for some x ∈ Ω, the weight y − s ( x ) is not, in general, in the Muckenhoupt class A (cf. Proposition 1). Due tothe lack of this A property we cannot use some of the existing machinery on weighted on-Muckenhoupt weights, fractional elliptic operators, and applications s -weighted L space. In addition, we giveexamples of functions with discontinuous traces in Example 1.In section 3, we establish properties for the variational problem of interest (1.2),and a procedure for the selection of the function s . We finalize the paper with section 4,where we develop a discretization scheme, and provide numerical tests that positivelycompared with other methods for image reconstruction. In particular, in section 4.1we introduce a finite element method for (1.6) on a bounded domain C τ := C × (0 , τ ).Since our targeted application is image denosing so we first state an Algorithm todetermine s in Algorithm 1. We apply this approach to several prototypical examplesin section 5. In all cases we obtain better results when compared to the Total Variation(TV) approach. Our point of departure for motivationto study the proposed variational model (and associated function space properties)is Figure 1 where the left panel shows a noisy image, the middle panel shows the
Fig. 1 . Example 1. Left panel: noisy image. Middle panel: reconstruction using Total Variation(TV) approach [8]. Right panel: reconstruction using the approach introduced in this paper. reconstruction using the Total Variation (TV) model, explained in what follows, (weuse piecewise linear finite element discretization for the TV model and we stop thealgorithm when the relative error between the consecutive iterates is smaller than1 e -8.), and the rightmost panel shows the reconstruction using the new frameworkintroduced in this paper and briefly explained in this section. The Rudin-Osher-Fatemi (ROF) model, also called the TV model, was introduced in [23] and is givenby(1.3) min u ∈ BV(Ω) ˆ Ω | Du | + ζ (cid:107) u − f (cid:107) L (Ω) , where Ω ⊂ R N with N ≥ f ∈ L (Ω) isgiven by f = u true + ξ , where ξ is the “noise”. Moreover, ζ > ´ Ω | Du | denotes the total variation seminorm of u on Ω ; see [4]. Theparameter ζ > u tv to (1.3) is to u true ; see [26]. The fundamental question in such applications is: Is it
H. Antil and C.N. Rautenberg possible to capture the sharp transitions (edges) across the interfaces while removingthe undesirable noise from remainder of the domain.
Such a question is not limited toimage denoising but is fundamental to many applications in science and engineering.For instance multiphase flows (diffuse interface models), fracture mechanics, imagesegmentation etc. Even though the ROF model has been extremely successful inpractice, two main drawbacks are observed in reconstructions: 1) loss of contrast isalways present 2) some corners tend to be rounded (cf. Figure 1, middle panel).The problem (1.1) proposed in [1] is obtained if ´ Ω | D · | is replaced by (cid:107) ( − ∆) s ·(cid:107) L (Ω) where ( − ∆) s denotes the fractional powers of the Laplacian with zero Neu-mann boundary conditions, for s ∈ (0 , H s (Ω), thefractional Sobolev space given by H s (Ω) = (cid:40) v ∈ L (Ω) : ˆ Ω ˆ Ω | v ( x ) − v ( y ) | | x − y | n +2 s d x d y < + ∞ (cid:41) . The minimization problem (1.1) is appealing, in contrast to (1.3), as the Euler-Lagrange equations are linear in u , i.e., if u solves (1.1), then equivalently it solves(1.4) (cid:104) ( − ∆) s u, v (cid:105) + ζ ( u − f, v ) = 0 , for all v ∈ H s (Ω). Both (1.1) and hence (1.4) have unique solutions in H s (Ω). Thiscan be easily deduced from the definition of Neumann fractional Laplacian and theequivalence between the spectral and H s (Ω)-norms (cf. [2, Remark 2.5] and [6]).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 (1.4), then there exists U ∈ H ( C ; y − s ) where H ( C ; y − s ) is theweighted H Sobolev space on C := Ω × (0 , ∞ ) with weight w ( x, y ) = y − s , ( x, y ) ∈ C (see the following section for more details), such that tr Ω U = u , where tr Ω operatoris the restriction of the trace map to Ω, and U solves ˆ C y − s ∇U · ∇W d x d y + d s ζ ˆ Ω (tr Ω U − f ) tr Ω W d x = 0 , (1.5)for all W ∈ H ( C ; y − s ), where d s = 2 − s Γ(1 − s )Γ( s ) . Therefore, for s a constant suchthat s ( x ) = s ∈ (0 ,
1) for all x ∈ Ω, θ = 0 and ζ := µs d − s , if U solves (1.2), thenit equivalently solves (1.5). This establishes the connection of (1.2) with (1.4), andhence with (1.1).Further, if s ( · ) is not a constant and s : Ω → [0 ,
1] (with some additional assump-tions made explicit in the following section), and if
U ∈ H ( C ; y − s ( x ) ) solves (1.2),then ˆ C y − s ( x ) ( ∇U · ∇W + θ UW ) d x d y + µ ˆ Ω s ( x ) (tr Ω U − f ) tr Ω W d x = 0 , (1.6)for all W ∈ H ( C ; y − s ( x ) ); this is rigorously done in section 3 where the key ingredientis the characterization of tr Ω of H ( C ; y − s ( x ) ). It is clear that for θ (cid:28)
1, the aboveproblem (whenever well-posed) represents a generalization of (1.5) and hence of (1.1).Now we are in shape to explain Figure 1. First, note that for s (cid:39)
0, we expectlow regularity of solutions, and for s (cid:39) on-Muckenhoupt weights, fractional elliptic operators, and applications < s < s <
1, we observe H s (Ω) ⊂ H s (Ω). Then, after a roughidentification of edges, we choose s small on edges, and large on other regions; theexact procedure is provided in section 4.2. The right image in Figure 1, is obtainedvia this procedure.
2. Weighted Sobolev space with variable s and their traces. Let Ω ⊂ R N with N ≥ ∂ Ω. Wedenote by C := Ω × (0 , ∞ ) a semi-infinite cylinder with boundary ∂ L C := ∂ Ω × [0 , ∞ ).With s : Ω → R measurable and s ( x ) ∈ [0 ,
1] for almost all x ∈ Ω, we define δ ( x ) := 1 − s ( x ) ∈ [ − , , and w ( x, y ) := y δ ( x ) , so that(2.1) w, w − ∈ L ( C ) . We denote L ( C ; y δ ( x ) ) := L ( C ; w ), where L ( C ; y δ ( x ) ) := (cid:26) v : C → R : v is measurable and ˆ C y δ ( x ) | v | d x d y < + ∞ (cid:27) , which is a well-defined weighted Lebesgue space.For u smooth, define the (extended real-valued) functional (cid:107) · (cid:107) H as(2.2) (cid:107) u (cid:107) H := (cid:18) ˆ C y δ ( x ) (cid:16) | u | + |∇ u | (cid:17) d x d y (cid:19) . Note that if (cid:107) u (cid:107) H < + ∞ , the value of s controls how singular the function u is inthe neighborhood of Ω × { } . Define A α ( s ) := { x ∈ Ω : s ( x ) = α, a.e. } , supposethat for α = 1 its Lebesgue measure satisfies |A ( s ) | >
0, and that u is a constant in A ( s ) × (0 , h ), then (cid:107) u (cid:107) H ≥ | u | |A ( s ) | ˆ h y − d y, so that u is zero in Ω × (0 , h ). On the other hand, u is allowed to have a more singularbehaviour on A ( s ).Consider the set C Ω = (Ω × { } ) ∪ C , that is, C Ω is the open cylinder C together with the Ω cap, and denote by H , H , and W , or equivalently by H ( C ; y δ ( x ) ), H ( C ; y δ ( x ) ), and W ( C ; y δ ( x ) ), the spaces H ( C ; w ), H ( C ; w ) and W ( C ; w ) which are defined as H ( C ; w ) is the closure of the set K := { w ∈ C ∞ ( C ) : (cid:107) w (cid:107) H < + ∞} with respect tothe (cid:107) · (cid:107) H norm. H ( C ; w ) is the closure of the set K := { w ∈ C ∞ c ( C Ω ) : (cid:107) w (cid:107) H < + ∞} with respectto the (cid:107) · (cid:107) H norm. W ( C ; w ) is the space of maps u ∈ L ( C ; y δ ( x ) ) ∩ L ( C ) with distributional gradients ∇ u that satisfy |∇ u | ∈ L ( C ; y δ ( x ) ) ∩ L ( C ). H. Antil and C.N. Rautenberg
A few words are in order concerning the definition of H given the slight differencefrom the existing literature: Usually, the definition of H is given over the completion(of finite energy functions) of C ∞ ( C ), in contrast to C ∞ ( C ), with respect to the H norm; see [27]. Provided that ∂ C is regular enough, C is bounded, and w, w − arebounded above and below ( (cid:15) > ∂ C , thenthe closure in the definition of H can be equivalently taken with C ∞ ( C ) or C ∞ ( 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 ⊂ K , which follows from C ∞ c ( C Ω ) ⊂ C ∞ ( C ), we have that H ⊂ H . The notation of “ H ” comes from the following: If v ∈ K , then v | ∂ Ω ×{ } =0, so that restrictions of elements of H to Ω have “zero boundary conditions”. Inparticular, in the case s ∈ (0 ,
1) is constant, the latter remark can be taken out ofquotation marks and considered in the sense of the trace; see [7]. The condition (2.1)implies (see [27, 16]) that the space W is a proper weighted Sobolev space endowedwith the H -norm. While it holds that H ( ˜ C ; ˜ w ) ⊂ W ( ˜ C ; ˜ w ), for general weights ˜ w and domains ˜ C , it does not necessarily hold that H ( ˜ C ; ˜ w ) = W ( ˜ C ; ˜ w ). For ˜ w = 1,the equality holds provided conditions on ˜ C are given, e.g., if ˜ C is bounded and has aLipschitz boundary. A sufficient condition for H ( ˜ C ; ˜ w ) = W ( ˜ C ; ˜ w ), when ˜ w and ˜ C arebenign enough, is that ˜ w belongs to the A ( ˜ C ) Muckenhoupt class, that is(2.3) (cid:32) Q ˆ Q ˜ w d x d y (cid:33) (cid:32) Q ˆ Q ˜ w − d x d y (cid:33) ≤ M, for some M > Q ⊂ ˜ C , see [25, 10]. Further, note that (2.1),does not imply (2.3). Additionally, if 0 < (cid:15) ≤ s ≤ − (cid:15) a.e., then ( x, y ) (cid:55)→ y − s is aMuckenhoupt A ( C ) weight. However, we are particularly interested in cases where s satisfies(2.4) ess inf x ∈ Ω s ( x ) = 0 , as this will allow almost perfect approximation of functions with jump discontinuities.Remarkably enough condition (2.4) prevents w to be a Muckenhoupt weight in generalas we show next. Proposition Suppose that for x ∈ Ω , s ( x ) (cid:39) | x − x | q (locally) on x forsome q > , that is, there exists positive constants R , M , m such that (2.5) m | x − x | q ≤ s ( x ) ≤ M | x − x | q , for all x such that | x − x | ≤ R . Then, w / ∈ A ( C ) where w ( x, y ) := y − s ( x ) .Proof. First step. Suppose first that x = 0 ∈ Ω , s ( x ) = | x | q for x ∈ B R (0) , andlet Q R = B R (0) × (0 , y ) for < R ≤ R < and < y < . Using cylindrical coordinates we observe ˆ Q R w d x d y ≥ C ˆ R ˆ y r N − y − r q d y d r ≥ C ˆ R r N − y − r q d r ≥ C R N y , on-Muckenhoupt weights, fractional elliptic operators, and applications C , C , and C are positive and independent of R and y . Similarly, ˆ Q R w − d x d y ≥ C ˆ R ˆ y r N − y − (1 − r q ) d y d r ≥ C ˆ R r N − − q y r q d r ≥ C R N − q y R q , where also C , and C are positive and independent of R and y . Since | Q R | (cid:39) R N y ,we have (cid:32) | Q R | ˆ Q R w d x d y (cid:33) (cid:32) | Q R | ˆ Q R w − d x d y (cid:33) ≥ C y R q R q , for some C > R . Taking R ↓
0, we have that the L.H.S. expressionis unbounded.
Second step: Consider x (cid:54) = 0 , s ( x ) (cid:54) = | x | q , and show that the above inequalityholds similarly for a family of cubes R (cid:55)→ ˜ Q R .If x (cid:54) = 0 and s ( x ) = | x | q , then a linear change of coordinates can be performedand the proof above still holds, so the choice of x = 0 is without loss of generality.Additionally, let ˜ Q R ⊃ Q R be defined as ˜ Q R = S R (0) × (0 , y ) where S R (0)is the smallest square that contains B R (0). Hence, the quotient | S R (0) | / | B R (0) | isindependent of R and for c := | ˜ Q R | / | Q R | we observe (cid:32) | ˜ Q R | ˆ ˜ Q R w d x d y (cid:33) (cid:32) | ˜ Q R | ˆ ˜ Q R w − d x d y (cid:33) ≥ c (cid:32) | Q R | ˆ Q R w d x d y (cid:33) (cid:32) | Q R | ˆ Q R w − d x d y (cid:33) . Hence, by taking R ↓
0, we have that ( x, y ) (cid:55)→ y − | x | q is not in A ( C ).Finally, for s ( x ) (cid:54) = | x | q but when (2.5) holds true, we similarly have (cid:32) | Q | ˆ Q w d x d y (cid:33) (cid:32) | Q | ˆ Q w − d x d y (cid:33) ≥ ˜ c (cid:32) | Q | ˆ Q ˜ w d x d y (cid:33) (cid:32) | Q | ˆ Q ˜ w − d x d y (cid:33) , for some ˜ c >
0, where ˜ w ( x, y ) := y − | x | q and w ( x, y ) := y − s ( x ) , i.e., w / ∈ A ( C ). Remark
2. Although the above result is quite general with respect of the rate ofdecrease of s , at this point, we do not know if there are continuous functions s withzeros within the domain Ω, for which ( x, y ) (cid:55)→ y − s ( x ) is an element of A ( C ).The study of the trace space of H is of utmost importance in what follows. Inparticular, we are interested in the restriction of the trace operator to Ω. The paper[7] studies the trace of H when s is independent of x and away from 0. However, thatapproach which is based on [18] cannot be applied here. We next prove a trace char-acterization result, the proof is inspired by [20] where the authors consider boundeddomains and weights of distance to the boundary type. Theorem Let s be such that its zero set A ( s ) := { x ∈ Ω : s ( x ) = 0 , a.e. } haszero measure. Then, there exist a unique bounded linear trace operator tr Ω : H, H → L (Ω; (1 − δ ( x )) ) = L (Ω; s ( x ) ) H. Antil and C.N. Rautenberg such that tr Ω ( u ) = u | Ω ×{ } , for all u ∈ H ∩ C ∞ ( C ) , and hence also for all u ∈ H ∩ C ∞ c ( C Ω ) .Proof. Consider u ∈ H such that u ∈ C ∞ ( C ). Let ( x, y ) ∈ C and x / ∈ A ( s ), thenit follows that u ( x, y ) = u ( x,
0) + ˆ yD N +1 u ( x, ty ) d t, where D N +1 corresponds to the partial derivative with respect to the N +1 coordinate.Then, multiplication by y δ ( x ) / and integrating from 0 to σ > y leadsto | u ( x, | ˆ σ y δ ( x ) / d y ≤ ˆ σ | u ( x, y ) | y δ ( x ) / d y + ˆ σ ˆ y | D N +1 u ( x, ty ) | y δ ( x ) / d y d t. Note that δ ( x ) / ≤ / < σ ≤
1, then y / ≤ y δ ( x ) / for all y ∈ (0 , σ ). Therefore,23 σ / ≤ ˆ σ y δ ( x ) / d y, and hence,(2.6) | u ( x, | ≤ I + I , for I := 32 σ − / ˆ σ | u ( x, y ) | y δ ( x ) / d y, I = 32 σ − / ˆ ˆ σ y | D N +1 u ( x, ty ) | y δ ( x ) / d y d t, and where we have used Tonelli’s theorem to switch the order of integration in theexpression that orginates I .H¨older’s inequality implies that, for I , we have I ≤ σ − (cid:18) ˆ σ d y (cid:19) (cid:18) ˆ σ | u ( x, y ) | y δ ( x ) d y (cid:19) ≤ σ − (cid:32) ˆ R | u ( x, y ) | y δ ( x ) d y (cid:33) , for any R ≥ σ .Further, for I consider the change of variable z = ty and again by H¨older’sinequality we observe: I ≤ σ − ˆ ˆ σt | D N +1 u ( x, z ) | z δ ( x )2 t − − δ ( x )2 d z d t ≤ σ − (cid:32) ˆ (cid:16) t − δ ( x )4 (cid:17) d t (cid:33) ˆ (cid:32) ˆ σt | D N +1 u ( x, z ) | z δ ( x )2 t δ ( x )4 − − δ ( x )2 d z (cid:33) d t ≤ σ − (cid:18) − δ ( x ) (cid:19) ˆ t − δ ( x )2 (cid:32) ˆ σt | D N +1 u ( x, z ) | z δ ( x )2 d z (cid:33) d t ≤ (cid:18) − δ ( x ) (cid:19) ˆ t − δ ( x )2 (cid:32) ˆ σt | D N +1 u ( x, z ) | z δ ( x ) d z (cid:33) d t . on-Muckenhoupt weights, fractional elliptic operators, and applications R ≥ σ , we have I ≤ (cid:18) − δ ( x ) (cid:19) / (cid:32) ˆ t − δ ( x )2 d t (cid:33) / (cid:32) ˆ R | D N +1 u ( x, z ) | z δ ( x ) d z (cid:33) / ≤ − δ ( x ) (cid:32) ˆ R | D N +1 u ( x, z ) | z δ ( x ) d z (cid:33) / . Therefore, from (2.6), we obtain that for any R ≥ σ | u ( x, | (1 − δ ( x )) ≤ σ − (1 − δ ( x )) (cid:32) ˆ R | u ( x, y ) | y δ ( x ) d y (cid:33) + 3 (cid:32) ˆ R | D N +1 u ( x, z ) | z δ ( x ) d z (cid:33) ≤ σ − (cid:32) ˆ R | u ( x, y ) | y δ ( x ) d y (cid:33) + 3 (cid:32) ˆ R | D N +1 u ( x, z ) | z δ ( x ) d z (cid:33) ≤ σ − (cid:32) ˆ R | u ( x, y ) | y δ ( x ) d y + ˆ R | D N +1 u ( x, z ) | z δ ( x ) d z (cid:33) . Integration over Ω with respect to x of the above squared inequality leads to ˆ Ω | u ( x, | (1 − δ ( x )) d x ≤ M σ (cid:32) ˆ Ω ˆ R | u ( x, y ) | y δ ( x ) d y + ˆ Ω ˆ R | D N +1 u ( x, z ) | z δ ( x ) d z (cid:33) ≤ M σ (cid:32) ˆ Ω ˆ R | u ( x, y ) | y δ ( x ) d y + ˆ Ω ˆ R |∇ u ( x, z ) | z δ ( x ) d z (cid:33) , where M σ = 6 σ − , R ≥ σ , with σ ∈ (0 ,
1] and for an arbitrary u ∈ H such that u ∈ C ∞ ( C ).Since R ≥ σ is arbitrary, we have that for an arbitrary u ∈ H ∩ C ∞ ( C ), we have (cid:107) u ( · , (cid:107) L (Ω;(1 − δ ( x )) ) ≤ M σ (cid:107) u (cid:107) H . The operator tr Ω is the above extension to H : Since H is the closure of C ∞ ( C ) withrespect to the energy norm, we can take σ = 1 with M σ = 6, and obtain (cid:107) tr Ω u (cid:107) L (Ω;(1 − δ ( x )) ) ≤ (cid:107) u (cid:107) H , which completes the proof, since (1 − δ ( x )) = 4 s ( x ) . The result for H is analogous.The behavior of s controls the local regularity on space tr Ω H , in particular, tr Ω H contains piecewise constants in any dimension for appropriate choices of s as we seein the following example. Example There exists non-constant s : Ω → R maps such that tr Ω H containspiecewise smooth functions. For the sake of simplicity, let Ω = B / (0) , define s ( x ) =0 H. Antil and C.N. Rautenberg ˆ s ( x ) , where x = ( x , x , . . . , x N ) , and such that < δ ≤ ˆ s ( x ) ≤ κ < / for all | x | < (cid:15) , and < δ ≤ ˆ s ( x ) for all x . Let w ( x, y ) := , x < −√ y ; √ y ( x + √ y ) , √ y ≤ x ≤ −√ y ; , √ y < x .so that |∇ w ( x, y ) | = , x < −√ y ; y + x y , √ y ≤ x ≤ −√ y ; , √ y < x .Define u ( x, y ) := η ( y ) w ( x, y ) for ( x, y ) ∈ C , with η ∈ C ∞ c ( R ) , such that η ( y ) = 1 if y ∈ [0 , R ∗ ] for some R ∗ , and therefore, ˆ B (cid:15) (0) ˆ y − s |∇ w ( x, y ) | d x d y ≤ ˆ B (cid:15) (0) ˆ y − κ |∇ w ( x, y ) | d x d y < + ∞ . Then, u ∈ W and smoothing argument involving mollification shows that u ∈ H , and tr Ω u = χ { x ≥ } .
3. The Optimization Problem .
Throughout this section we suppose that f ∈ L (Ω; s ( x ) ), s ( x ) ∈ [0 ,
1] for almost all x ∈ Ω, and that A ( s ) := { x ∈ Ω : s ( x ) =0 , a.e. } has zero measure.Given a regularization parameter µ >
0, and parameter θ > u ∈ H ( C ; y − s ( x ) ) J ( u, s ) := ˆ C y − s ( x ) (cid:16) θ | u | + |∇ u | (cid:17) d x d y + µ ˆ Ω s ( x ) | tr Ω u − f | d x. Remark θ ). In case the weights w are of class A , we can set θ = 0, given that the Poincar´e type inequalities are available in this case. However,the weights y − s ( x ) in our case are not in A since we allow s ( x ) = 0 for some x ∈ Ω.Such Poincar´e type inequalities on H ( C ; y − s ( x ) ) are not known yet.The model (3.1) provides a completely new approach to approximate nonsmoothfunctions f with tr Ω u . Such a situation often occurs in inverse problems where given f one is interested in a reconstruction tr Ω u which is close to f . For instance in imagedenoising problems where f represents a noisy image, given by f = u true + ξ where u true is the true target to recover, ´ Ω ξ = 0, and ´ Ω | ξ | = σ for a known σ > Ω u is close to u true on Ω.For a fixed s and given f we next state existence and uniqueness of solution to(3.1). This follows immediately from Theorem 3. Theorem
There exist a unique solution to (3.1) .Proof.
Existence follows from application of direct methods of calculus of varia-tions, the norm definition on H and the Theorem 3. Uniqueness follows directly fromconvexity arguments.The first order necessary and sufficient optimality conditions for the minimization on-Muckenhoupt weights, fractional elliptic operators, and applications ˆ C y − s ( x ) ( ∇ u · ∇ w + θuw ) d x d y + µ ˆ Ω tr Ω u tr Ω w s ( x ) d x = µ ˆ Ω f tr Ω w s ( x ) d x, for all w ∈ H. (3.2)The first step for computational implementation of the problem above requires to solvethe problem in a bounded domain. Notice that all results above are valid if we replacethe unbounded domain C with boundary ∂ L C by a bounded domain C τ = Ω × (0 , τ )with 0 < τ < ∞ and boundary ∂ L C τ = ( ∂ Ω × [0 , τ ]) ∪ Ω × { τ } . In that case theproblem (3.2) becomes: find v ∈ H τ := H ( C τ ; y δ ( x ) ) ˆ C τ y − s ( x ) ( ∇ v · ∇ w + θvw ) d x d y + µ ˆ Ω tr Ω v tr Ω w s ( x ) d x = µ ˆ Ω f tr Ω w s ( x ) d x, for all w ∈ H τ . (3.3)It is known that for a constant s ∈ (0 ,
1) the solution v to (3.3) converges to u solving(3.2) exponentially with respect to τ >
1, see [21]. Such exponential approximationis also expected in the case of a non-constant s . Indeed, as we discussed earlier, aconstant s implies that tr Ω u solves the fractional equation (1.4). A relation between(3.2) and fractional PDE of type (1.4) with s ( x ) instead of s is an open question.We expect that studying (3.2) or equivalently (3.1) is a first step in establishing thedefinition of fractional ( − ∆) s ( x ) for x ∈ Ω.For given s and f , Theorem 5 shows existence of a unique solution to (3.1) andequivalently (3.2). Nevertheless for our model (3.1) to be useful in practice we needto determine the function s . For this matter, consider the following important pointsthat will lead to the selection procedure for s .( i ) Precise edge recovery.
In example 1, a family of s functions is identifiedso that u ∈ H ( C ; y − s ( x ) ), and tr Ω u = χ { x ≥ } . This suggests that for theoptimization problem (3.1) to recover edges, or more precisely for tr Ω u to havediscontinuities 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 suggestto utilize a rough edge detection at first and then force s to be close zero onneighborhoods of these regions.( ii ) Homogeneous/flat regions recovery.
We proceed rather formally first. Sup-pose first that θ = 0, and that problem (3.1) is posed not on H ( C ; y − s ( x ) ),but on the Banach space of equivalence classes generated by the seminorm ´ C y − s ( x ) |∇ u | d x d y . Consider s ( x ) = 1 on a certain region Ω . Supposethat |∇ u | ≥ e > × (0 , y ), then ˆ C y − s ( x ) |∇ u | d x d y ≥ (cid:15) ˆ Ω × (0 ,y ) y − d x d y = + ∞ . Hence, the solution to (3.1), needs to satisfy |∇ u | (cid:39) ×{ } . This suggestthat flat regions would be recovered if s = 1 on that same region. However, if θ >
0, then analogously, solutions to problem (3.1), would be forced to satisfy | u | (cid:39) × { } . Hence, we consider 0 < θ (cid:28) s ( x ) (cid:39)
1, but s ( x ) < u is not forced to beidentically zero but |∇ u | (cid:39) H. Antil and C.N. Rautenberg
We provide an algorithm based on the two points above in what follows. Amore detailed bilevel optimization framework (where both u and s are optimizationvariables) will be considered in a forthcoming publication.
4. Numerical Method.
We now focus on the discretization of the truncatedproblem (3.3) and the selection of an appropriate s function. (3.3) . From hereon we will assume that Ω is polygo-nal/polyhedral Lipschitz. We recall that the results of previous sections remains validif we replace the unbounded domain C with boundary ∂ L C by C τ = Ω × (0 , τ ) with τ > ∂ L C = ( ∂ Ω × [0 , τ ]) ∪ Ω × { τ } . The Euler-Lagrange equationsfor the resulting problem on C τ are given by (3.3).We begin by introducing a discretization for (3.3). We will follow the notationfrom [3]. For a constant s such a discretization for v was first considered in [21]. Let T Ω = { E } be a conforming and quasi-uniform triangulation of Ω, where E ∈ R N isan element that is isoparametrically equivalent to either to the unit cube or to theunit simplex in R N . We assume T Ω ∝ M N . Thus, the element size h T Ω fulfills h T Ω ∝ M − .Furthermore, let I τ = { I k } K − k =0 , where I k = [ y k , y k +1 ], is anisotropic mesh in [0 , τ ]in the sense that [0 , τ ] = (cid:83) K − k =0 I k . For a constant s ∈ (0 ,
1) we define the anisotropicmesh in y -direction as(4.1) y k = (cid:18) kK (cid:19) γ , k = 0 , . . . , K, γ > 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 thesecan be used to compensate the singular effects [19, 21]. In all our implementationswe will choose a fixed constant s in (4.1).We construct the triangulations T τ of the cylinder C τ as tensor product triangu-lations by using T Ω and I τ . Let T denotes the collection of such anisotropic meshes T τ . For each T τ ∈ T we define the finite element space V ( T τ ) as V ( T τ ) := { V ∈ C ( C τ ) : V | T ∈ P ( E ) ⊕ P ( I ) ∀ T = E × I ∈ T τ } . In case ∂ Ω has Dirichlet boundary conditions we define our finite element space as V ( T τ ) = V ( T τ ) ∩ { V : V | ∂ C τ = 0 } , i.e., functions with zero boundary conditions. Incase E is a simplex then P ( E ) = P ( E ), the set of polynomials of degree at most 1.If E is a cube then P ( E ) equals Q ( E ), the set of polynomials of degree at most 1in each variable. In our numerical illustrations we shall work with simplices.We define the finite element space for s as S ( T Ω ) := { S ∈ L ∞ (Ω) : S | E ∈ P ( K ) for all E ∈ T Ω } which is a space of piecewise constant functions defined on T Ω . The discrete versionof (3.3) is then given by: Find V ∈ V ( T τ ) ˆ C τ y − S ( x ) ( ∇ V · ∇ W + θV W ) d x d y + µ ˆ Ω tr Ω V tr Ω W S ( x ) d x = µ ˆ Ω f tr Ω W S ( x ) d x, for all W ∈ V ( T τ ) , (4.2) on-Muckenhoupt weights, fractional elliptic operators, and applications S ∈ S ( T Ω ). We compute the stiffness and mass matrices in (4.2) exactly. Thecorresponding forcing boundary term is computed by a quadrature formula whichis exact for polynomials up to degree 5. For a given S the resulting discrete linearsystem (4.2) is solved using the Preconditioned Conjugate Gradient (PCG) methodwith a block diagonal preconditioner. In view of what was stated about θ in section 3,we have set θ = 10 e -10 in all our examples. We let τ = 1 + ( T Ω ) , this choiceis motivated by the fact that for a constant s such a τ balances the finite elementapproximation on C τ and the truncation error from C to C τ [21, Remark 5.5]. Thenumber of points in the y -direction is taken to be K = 20. We use a moderatelyanisotropic mesh in the y -direction by setting s = 0 .
32 in (4.1). Our experimentssuggest 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(4.2). We have observed that µ only affects the “contrast” or magnitude of tr Ω V .Nevertheless one can determine µ using the well established techniques such as L-curve method [12]. In our case we choose a fixed µ for each example (no optimizationwas carried to select µ ). On the other hand selecting S which is a function is muchmore delicate. One option is to use a bilevel strategy as in [13, 14] to determine both S and V in an optimization framework. This is a part of our future work. In thispaper we propose a different approach in Algorithm 1.Even though our targeted application is image denoising the approach we presentis general enough to be applicable to a wider range of applications. We first noticethat a typical image is given on a rectangular grid (pixels). Since for (4.2) we areworking on simplices we first need to interpolate the given image from the rectangulargrid to a simplicial mesh. This is a delicate question especially given the fact thattypically we only have access to a noisy image. Before we interpolate the noisy imageonto a simplicial mesh we perform an intermediate step.We solve ROF model (1.3) using [8] and ζ > tv . We select a mild tolerance tol tv inthis 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 Lagrangeinterpolant I T Ω u tv of u tv .We next evaluate I T Ω u tu on the simplicial mesh. In order to reduce the com-putational cost we use Adaptive Finite Element Method (AFEM) to generate theappropriate mesh. In the nutshell, we start with a coarse mesh T Ω = { E } where E is an element in T Ω . For each E ∈ T Ω we then evaluate gradient of I T Ω u tv on E and we denote this gradient by ∇ I T Ω u tv | E . We use this gradient to define an edgeindicator function, we call this edge indicator function as the estimator on E . Basedon a marking strategy we then mark a subset of elements in T Ω . Subsequently weperform 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 theimage and close to 1 away from sharp edges. Intuitively smaller the S the lesser thesmoothness and otherwise.Once we have S then we can immediately solve the linear equation (4.2) for V .
5. Numerical Examples.
In this section we illustrate the proposed schemewith the help of several examples. We consider f given by f = u true + ξ where u true isthe true image, and ξ is noise with the following properties: ´ Ω ξ = 0, and ´ Ω | ξ | = σ for a known σ >
0. In all cases we find S by using Algorithm 1 and then solve the4 H. Antil and C.N. Rautenberg
Algorithm 1
Selection of S Data: f ∈ L (Ω , s ( x ) ), ζ , tol tv , N refine , λ , β , ν Solve total variation minimization problem (1.3) with regularization parameter ζ and tolerance tol tv using for instance [8]. Generate a piecewise linear Lagrangeinterpolant I T Ω u tv of u tv . Construct an Adaptive Finite Element Method (AFEM) based on the followingiterative loop with N refine iterations: Solve → Estimate → Mark → Refine
We describe each of these modules next:a.
Solve : For a given triangulation T Ω = { E } evaluate the elementwise gradi-ent of I T Ω u tv . We denote the gradient on each element E by ∇ I T Ω u tv | E .b. Estimate : Use the edge indicator function E ( E ; λ ) = 1 − (cid:16) λ − |∇ I T Ω u tv | E | (cid:17) − ∀ E ∈ T Ω as an estimator. Here, λ > Mark : Use the D¨orfler marking strategy [9] (bulk chasing criterion) withparameter β ∈ (0 , M ⊂ T Ω fulfilling (cid:88) E ∈M E ( E ; λ ) ≥ β (cid:88) E ∈ T Ω E ( E ; λ ) . d. Refine : Generate a new mesh T (cid:48) Ω = { E (cid:48) } by bisecting all the elementscontained in M using the newest-vertex bisection algorithm [22]. On each element E (cid:48) ∈ T (cid:48) Ω , set S ( E (cid:48) ) = 1 − E ( E (cid:48) ; ν ).resulting linear equation (4.2) for V . We call tr Ω V as the reconstruction. We alsocompare our reconstruction tr Ω V with the reconstruction we obtain by solving (1.3)with tol ∗ tv = 10 e -8. In this case we choose the parameter ζ so that a normalizedweighted sum of Peak Signal to Noise Ratio (PSNR) and the Structural Similarity(SSIM) index is maximized, as in section 4.2 of [14]. That is, we compare our resultsto the ones of the TV model with an optimized parameter ζ >
0. Note that thisparameter is not the one used in step 1 from algorithm 1.
In our first example we con-sider an image with a circle, triangle and a square as shown in Figure 2 (top row left).Notice that in this case the right pointing edge of the triangle does not align withthe grid. We consider two different noise levels. In both cases the noise is normallydistributed with mean 0 and standard deviations 0.10 and 0.15 respectively. Theresults are shown in Figures 2 and 3 for standard deviations 0 .
10 and 0 .
15, respec-tively. In both cases we first compute S by using Algorithm 1 where we have set: ζ = 0 .
2, tol tv = 10 e -4, N refine = 8, λ = 300, β = 0 . ν = 200 in Algorithm 1. Wefurther set µ = 8050 in (4.2). We then solve for V using PCG. We call tr Ω V as ourreconstruction.The top row of Figure 2 shows the original and noisy images (left to right). Themiddle row shows u tv and S obtained using Algorithm 1. In the bottom row wecompare the results using the TV approach (left) and tr Ω V (right) computed using on-Muckenhoupt weights, fractional elliptic operators, and applications Ω V panels, froma different angle. Similar description holds for Figure 3. Figure 4 is again viewingFigure 3 from a different angle.As we can notice TV tends to round up the corners (cf. Figure 1 (middle)).On the other hand as our theory predicted we can truly capture the edges in tr Ω V (cf. Figure 1 (right)). This is further corroborated by Table 1 where we have showna comparison between PSNR and SSIM for these two approaches for two differentstandard deviations. σ PSNR (TV) PSNR (New) SSIM (TV) SSIM (New)0.1 3.7299e+01 4.8147e+01 9.4016e-01 9.5710e-010.15 3.5712e+01 4.0451e+01 9.3439e-01 9.4890e-01
Table 1
Example 1: PSNR and SSIM using two different standard deviations ( σ = 0 . and σ = 0 . )using TV and proposed scheme (New). In our second example we consider 6 stripes withintensities equal to 0.8, 0.7, 0.6, 0.5, 0.4 and 0.35 (from left to right), respectively(cf. Figure 5 (top row left)). As in the previous example we again consider twoadditive noise levels with mean 0 standard deviations 0 . .
15 respectively. Ourresults are shown in Figures 5 and 6 respectively. Here again we first find S byusing Algorithm 1 using the parameters ζ = 0 .
2, tol tv = 10 e -4, N refine = 8, λ = 15, β = 0 . ν = 100. We further set µ = 2900 in (4.2). We then solve for V . Wecall tr Ω V as our reconstruction. A comparison between PSNR and SSIM is shownin Table 2. As we noticed in the previous example we again obtain visually almostperfect reconstruction using our approach. σ PSNR (TV) PSNR (New) SSIM (TV) SSIM (New)0.1 3.6253e+01 2.7917e+01 9.0799e-01 9.4789e-010.15 3.3616e+01 2.7419e+01 8.7028-01 9.3000e-01
Table 2
Example 2: PSNR and SSIM using two different standard deviations ( σ = 0 . and σ = 0 . )using TV and proposed scheme (New). In the first two examples we considered syn-thetic images. In our final example we consider a more realistic situation. We considera prototypical image of the cameraman (cf. Figure 7). As in the previous exampleswe again consider two additive noise levels with standard deviations 0 . . S . Here we have set the underlyingparameters as ζ = 0 .
2, tol tv = 10 e -4, N refine = 8, λ = 0 . β = 0 . ν = 20.We further set µ = 10 in (4.2). We then solve for V and we call tr Ω V as ourreconstruction. A comparison between PSNR and SSIM is shown in Table 3. As wenoticed in the previous examples we again obtain better reconstructions using ourapproach.
6. Conclusion and further directions.
A new variational model associated tothe fractional Laplacian was introduced. In particular, we have identified a weighted6
H. Antil and C.N. Rautenberg
Fig. 2 . Example 1 ( σ = 0 . ). Top row (from left to right): original and noisy images, re-spectively. Middle row: u tv (left) from Step 1 in Algorithm 1 and the corresponding s (right)from Step 3. Bottom row: reconstruction using total variation with optimized ζ > (left) and ourapproach (right), respectively. Sobolev space with respect to the weight w = y − s ( x ) appropriate for the treatmentof the problem. We have shown that in general these weights are not of Muckenhoupttype, and have established that the trace space embeds in an s − weighted Lebesguespace. We have provided a discretization method for the full problem, and an algo-rithm for its resolution that also builds a selection procedure for s . The full schemeis advantageous when it comes to recovery of discontinuous features, details, homo-geneous regions, and also for contrast preservation, in data perturbed by additivenoise.Future research directions are multiple, we enumerate some of them. on-Muckenhoupt weights, fractional elliptic operators, and applications Fig. 3 . Example 1 ( σ = 0 . ). Top row (from left to right): original and noisy images,respectively. Middle row: u tv (left) from Step 1 in Algorithm 1 and the corresponding s (right)from Step 3. Bottom row: reconstruction using total variation with optimized ζ > (left) and ourapproach (right), respectively.
1) The study of the optimization problem (3.1), seems to be the first step for theidentification of a possible definition for ( − ∆) s ( x ) . Such a task does not seem directlyapproachable via spectral or functional calculus points of view.2) Full characterization of the trace space for H and H . We have identified theembedding of tr Ω H into L (Ω; s ( x ) ), but it would be of interest to understand theSobolev regularity of tr Ω H .3) Differentiability and stability properties of the solution to (3.1) with respect to s . For the optimal selection of s , it is required to identify a topology over an admissi-ble set for s , so that solutions to (3.1), are stable with respect to perturbations. This8 H. Antil and C.N. Rautenberg
Fig. 4 . Example 1 ( σ = 0 . ). Left panel: noisy image. Middle panel: reconstruction usingTotal Variation (TV) approach with optimal ζ > . Right panel: reconstruction using the ourapproach. TV tends to smooth out the edges and corners that are not aligned with the grid. On theother hand it is possible to have perfect recovery using our approach. seems like a complex task where the usual obstacles from homogenization and conver-gence of differential operators are present. In addition, difficulty is increased as thesolutions to (3.1) belong to a state space depending on s as well. The differentiabilityissue is further more complex, but such a study will be the first step in establishingstationarity 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 (3.1), admits the following generalizationmin u ∈ H ( C ; y − s ( x ) ) ˆ C y − s ( x ) (cid:16) θ | u | + |∇ u | (cid:17) d x d y + µ ˆ Ω s ( x ) | K (tr Ω u ) − f | d x, where K is a bounded linear operator on L (Ω; s ( x ) ). This would allow to deal withthe problem of finding y such that Ky = f . In this case, the choice of s could besimply associated with regions of the domain Ω where it is more important to recover y more accurately. σ PSNR (TV) PSNR (New) SSIM (TV) SSIM (New)0.1 2.7054e+01 2.9640e+01 8.0475e-01 8.3653e-010.15 2.4376e+01 2.6884e+01 7.4340e-01 7.9335e-01
Table 3
Example 3:PSNR and SSIM using two different standard deviations ( σ = 0 . and σ = 0 . )using TV and proposed scheme (New). Acknowledgement.
We thank Gunay Dogan for point us to the edge errorindicator and Pablo Raul Stinga for sharing the reference [20].
REFERENCES[1] H. Antil and S. Bartels. Spectral Approximation of Fractional PDEs in Image Processing andPhase Field Modeling.
Comput. Methods Appl. Math. , 17(4):661–678, 2017.[2] H. Antil, J. Pfefferer, and S. Rogovs. Fractional operators with inhomogeneous boundaryconditions: Analysis, control, and discretization. arXiv preprint arXiv:1703.05256 , 2017.[3] H. Antil and C.N. Rautenberg. Fractional elliptic quasi-variational inequalities: the-ory and numerics.
Accepted in Interfaces and Free Boundaries. Available inhttps://arxiv.org/pdf/1712.07001.pdf , 2018.[4] H. Attouch, G. Buttazzo, and G. Michaille.
Variational analysis in Sobolev and BV spaces ,volume 17 of
MOS-SIAM Series on Optimization . Society for Industrial and Appliedon-Muckenhoupt weights, fractional elliptic operators, and applications Fig. 5 . Example 2 ( σ = 0 . ). Top row (from left to right): original and noisy images, re-spectively. Middle row: u tv (left) from Step 1 in Algorithm 1 and the corresponding s (right)from Step 3. Bottom row: reconstruction using total variation with optimized ζ > (left) and ourapproach (right), respectively. Mathematics (SIAM), Philadelphia, PA; Mathematical Optimization Society, Philadelphia,PA, second edition, 2014. Applications to PDEs and optimization.[5] L. Caffarelli and L. Silvestre. An extension problem related to the fractional Laplacian.
Comm.Part. Diff. Eqs. , 32(7-9):1245–1260, 2007.[6] L.A. Caffarelli and P.R. Stinga. Fractional elliptic equations, Caccioppoli estimates and regu-larity.
Ann. Inst. H. Poincar´e Anal. Non Lin´eaire , 33(3):767–807, 2016.[7] A. Capella, J. D´avila, L. Dupaigne, and Y. Sire. Regularity of radial extremal solutions forsome non-local semilinear equations.
Comm. Part. Diff. Eqs. , 36(8):1353–1384, 2011.[8] A. Chambolle. An algorithm for total variation minimization and applications.
J. Math.Imaging Vision , 20(1-2):89–97, 2004. Special issue on mathematics and image analysis.[9] W. D¨orfler. A convergent adaptive algorithm for Poisson’s equation.
SIAM J. Numer. Anal. ,33(3):1106–1124, 1996. H. Antil and C.N. Rautenberg
Fig. 6 . Example 2 ( σ = 0 . ). Top row (from left to right): original and noisy images,respectively. Middle row: u tv (left) from Step 1 in Algorithm 1 and the corresponding s (right)from Step 3. Bottom row: reconstruction using total variation with optimized ζ > (left) and ourapproach (right), respectively. [10] J. Duoandikoetxea. Fourier analysis , volume 29 of
Graduate Studies in Mathematics . AmericanMathematical Society, Providence, RI, 2001.[11] V. Gol’dshtein and A. Ukhlov. Weighted Sobolev spaces and embedding theorems.
Trans.Amer. Math. Soc. , 361(7):3829–3850, 2009.[12] P.C. Hansen and D.P. O’Leary. The use of the L -curve in the regularization of discrete ill-posedproblems. SIAM J. Sci. Comput. , 14(6):1487–1503, 1993.[13] M. Hinterm¨uller and C.N. Rautenberg. Optimal selection of the regularization function in aweighted total variation model. Part I: Modelling and theory.
J. Math. Imaging Vision ,59(3):498–514, 2017.[14] M. Hinterm¨uller, C.N. Rautenberg, T. Wu, and A. Langer. Optimal selection of the regular-ization function in a weighted total variation model. Part II: Algorithm, its analysis andnumerical tests.
J. Math. Imaging Vision , 59(3):515–533, 2017.on-Muckenhoupt weights, fractional elliptic operators, and applications Fig. 7 . Example 3 ( σ = 0 . ). Top row (from left to right): original and noisy images, re-spectively. Middle row: u tv (left) from Step 1 in Algorithm 1 and the corresponding s (right)from Step 3. Bottom row: reconstruction using total variation with optimized ζ > (left) and ourapproach (right), respectively. [15] A. Kufner. Weighted Sobolev spaces . A Wiley-Interscience Publication. John Wiley & Sons,Inc., New York, 1985. Translated from the Czech.[16] A. Kufner and B. Opic. How to define reasonably weighted Sobolev spaces.
Comment. Math.Univ. Carolin. , 25(3):537–554, 1984.[17] A. Kufner and A.-M. S¨andig.
Some applications of weighted Sobolev spaces , volume 100 of
Teubner-Texte zur Mathematik [Teubner Texts in Mathematics] . BSB B. G. TeubnerVerlagsgesellschaft, Leipzig, 1987. With German, French and Russian summaries.[18] J.L. Lions. Th´eor`emes de trace et d’interpolation.
Annali della Scuola Normale Superiore diPisa-Classe di Scienze , 13(4):389–403, 1959.[19] D. Meidner, J. Pfefferer, K. Sch¨urholz, and B. Vexler. hp -finite elements for fractional diffusion. arXiv preprint arXiv:1706.04066 , 2017.[20] A. Nekvinda. Characterization of traces of the weighted Sobolev space W ,p (Ω , d (cid:15)M ) on M . H. Antil and C.N. Rautenberg
Fig. 8 . Example 3 ( σ = 0 . ). Top row (from left to right): original and noisy images,respectively. Middle row: u tv (left) from Step 1 in Algorithm 1 and the corresponding s (right)from Step 3. Bottom row: reconstruction using total variation with optimized ζ > (left) and ourapproach (right), respectively.Czechoslovak Math. J. , 43(118)(4):695–711, 1993.[21] R.H. Nochetto, E. Ot´arola, and A.J. Salgado. A PDE approach to fractional diffusion in generaldomains: A priori error analysis. Found. Comput. Math. , 15(3):733–791, 2015.[22] R.H. Nochetto, K.G. Siebert, and A. Veeser. Theory of adaptive finite element methods:an introduction. In
Multiscale, nonlinear and adaptive approximation , pages 409–542.Springer, Berlin, 2009.[23] L.I. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation based noise removal algorithms.
Phys. D , 60(1-4):259–268, 1992. Experimental mathematics: computational issues in non-linear science (Los Alamos, NM, 1991).[24] P.R. Stinga and J.L. Torrea. Extension problem and Harnack’s inequality for some fractionaloperators.
Comm. Part. Diff. Eqs. , 35(11):2092–2122, 2010.[25] B.O. Turesson.
Nonlinear potential theory and weighted Sobolev spaces . Springer, 2000.on-Muckenhoupt weights, fractional elliptic operators, and applications [26] Z. Wang, A.C. Bovik, H.R. Sheikh, and E.P. Simoncelli. Image quality assessment: from errorvisibility to structural similarity. IEEE transactions on image processing , 13(4):600–612,2004.[27] V. V. Zhikov. On weighted Sobolev spaces.
Mat. Sb. , 189(8):27–58, 1998.[28] V. V. Zhikov. On the density of smooth functions in a weighted Sobolev space.