Total generalized variation

Total Generalized Variation Kristian Bredies Karl Kunisch Thomas Pock May 14, 2010 Abstract The novel concept of total generalized variation of a function u is introduced and some of its essential properties are proved. Differ- ently from the bounded variation semi-norm, the new concept involves higher order derivatives of u. Numerical examples illustrate the high quality of this functional as a regularization term for mathematical imaging problems. In particular this functional selectively regularizes on different regularity levels and, as a side effect, does not lead to a staircasing effect. Keywords: Bounded variation, total generalized variation, tensor fields, regularization, image denoising. AMS Subject Classification: 49J52, 49N45, 68U10. 1 Introduction Most mathematical formulations of inverse problems and in particular of mathematical imaging problems are cast in the form min u F(u) +R(u), (1.1) where F represents the data fidelity and R the regularization term. If G denotes the forward modeling operator then the most common fidelity term is of the form F(u) = 1 2 ‖G(u)− z‖2, (1.2) where z stands for the possibly error-prone data and ‖ · ‖ denotes an ap- propriately chosen Hilbertian norm. Similarly the most frequently chosen regularization term is given by R(u) = α 2 |u|2, (1.3) where α is the regularization parameter and | · | again denotes a Hilbertian norm or semi-norm. It is now becoming well-accepted that the mathematical 1 and computational simplicity of the norm-of-squares terms must be put into perspective with some serious shortcomings. If the errors in the data contain outliers or if the error is of impulsive type the fidelity terms suggested by methods from robust statistics should be preferred over (1.2). Similarly (1.3) is frequently not an appropriate choice. In fact, the regularization term penalizes a certain property of u, which is quantified in the choice of the norm | · |, and the natural proportionality by which this should enter into R would be 1-homogeneous rather than quadratic. In the present paper the focus will be on the choice of R. One of the early proposal for a refined choice was given in [ROF]. It uses the bounded variation semi-norm R(u) = α ∫ Ω |Du|, (1.4) where u is defined on the bounded domain Ω ⊂ IRd. This choice is highly effective when compared to e.g. R(u) = α2 ∫ Ω |∇u|2 dx if the data z to be reconstructed are piecewise constant, since it is more apt to preserve corners and edges. The bounded variation semi-norm, however, also has some shortcomings, most notably the staircasing phenomenon. To briefly explain this effect, we assume that G = I so that (1.1) describes the imaging denoising problem. If the true image contains not only flat, but also slanted regions, then the image reconstructed on the basis of bounded variation semi-norm tends to be piecewise constant (staircasing). This staircasing effect was rigorously established in [Ni, CCN, R], for example. For diverse other aspects on the topic of constructing appropriate regularization or filter functionals in image reconstruction we refer to [SGGHL] and the references given there. In this paper we propose and analyze the regularization term of the form TGVkα(u) = sup {∫ Ω udivk v dx ∣∣∣ v ∈ Ckc (Ω, Symk(IRd)), ‖divl v‖∞ ≤ αl, l = 0, . . . , k − 1 } , (1.5) where Symk(IRd) denotes the space of symmetric tensors of order k with arguments in IRd, and αl are fixed positive parameters. For the definition of the remaining quantities, we ask for the readers’ patience until Section 2. Suffice it to say at this moment that for k = 1, α0 = 1 the semi-norm TGVkα coincides with the bounded variation semi-norm. We refer to TGV k α as total generalized bounded variation of order k with weight α ∈ IRk. From the definition of TGVkα it is immediately clear that it involves (generalized) derivatives of u of order i = 1, . . . , k, and that the kernel of TGVkα is the set of polynomials of degree less than k. Intuitively the total generalized bounded variation further automatically balances the first to the k-th derivatives of u among themselves. It will be shown that TGVkα shares some properties 2 of TV: It is also rotationally invariant and for k = 2, the total generalized variation of the indicator function of a smooth set Ω′ ⊂⊂ Ω equals α1 PerΩ′ = α1 TV(χΩ′), where PerΩ′ denotes the perimeter of Ω′. If differs, however, for functions which are not piecewise constant. As a further preview we point out that in dimension 1, with Ω = ]0, L[, k = 2, α0, α1 > 0 we have for u(x) = 2∑ i=1 pi(x)χΩi , with pi(x) = aix+ bi, for each a1, a2, b1, b2 ∈ IR and Ω1 = ]0, c[, Ω2 = ]c, 1[ that TGV2α(u) = α1|p2(c)− p1(c)|+ α0|p′1(c)− p′2(c)| provided that the jump point c is not near the boundary, i.e., c ∈ [α0/α1, L− α0/α1]. In particular, TGV2α(u) does not penalize the derivative of order l = 0, 1 unless it jumps at c. As already mentioned our motivation for studying TGV2α(u) is based on the fact that it involves and balances higher-order derivatives of u. As a consequence it reduces the staircasing effect of the bounded variation func- tional. This will be demonstrated in our numerical experiments. The use of higher-order derivatives with the aim of reducing staircasing is not new. In [CL] the inf-convolution functional min u1+u2=u ∫ Ω |∇u1|+ α|∇(∇u2)| dx was proposed and proved to be practically efficient, eliminating the stair- casing effect, for denoising problems with images which contain various grey levels as well as edges and corners. This idea was followed upon in a modified form in [CEP] where the regularization term is of the form min u1+u2=u ∫ Ω |∇u1|+ α|∆u2| dx, i.e. the second derivative is replaced by the Laplacian and a dual method for its numerical realization is derived. A different functional was proposed and tested in [CMM]. It is given by∫ Ω |∇u|+ αΦ(|∇u|)(L(u))2 dx, where Φ is a real-valued function that reflects the presence of edges in the sense that its value approaches 0 when the gradient |∇(u)| is large, and L(u) is an elliptic operator. For this choice of regularization functional the ab- sence of the staircasing effect was verified in [DFLM]. In [PS] regularization terms of the form R(u) = ∫ Ω |D∇l−1(u)|, 3 where ∫ Ω |D∇l−1(u)| denotes the total variation of the (l − 1)-th derivative of u ∈ W l−1,1, are considered, and special structure of minimizers of the resulting problems (1.1) are investigated. Higher-order regularization func- tionals in the discrete setting which are related to our computations for the second-order case were further proposed and tested in [SS]. As we shall see in Section 3, even for the case k = 2 the proposed functional TGVkα does not agree with those regularization functionals which were considered earlier in the literature. In particular, although convex, it is structurally different from any infimal-convolution functional, especially the approaches of [CL, CEP]. Let us give a briefly outline of the following sections. Section 2 contains a compact treatise of tensor and, in particular, symmetric tensor analysis in a manner that is useful for the variational analysis context. The defi- nition of total generalized variation norms and some of its basic properties are given in Section 3. Based on Fenchel duality an equivalent description of TGVkα(u) is derived for u sufficiently regular. Moreover the relationship to an appropriately defined k-fold inf-convolution is obtained. A subsection is devoted to the special case k = 2. The description of the numerical pro- cedure that was employed as well as carefully selected numerical denoising experiments are contained in Section 4. The Appendix contains some basic results involving symmetric k-tensor fields. 2 Spaces of symmetric tensor fields This section is mainly devoted to the introduction of the notions we are going to utilize in the main parts of this article. For many of the considerations which are going to follow, the concept of symmetric tensor fields plays a central role. Therefore, we give a rather extensive introduction to make this paper more self-contained and also for the convenience for those readers, who are familiar with tensor analysis. Most of the results, although scattered in several chapters, can also be found in the books [BG, S]. We mainly restrict ourselves to a general introduction of symmetric ten- sor fields. Some more specific results can be found in the appendix. Throughout the paper, d ≥ 1 denotes the dimension which is typically 2 or 3, in applications. Let T k(IRd) = {ξ : IRd × · · · × IRd︸ ︷︷ ︸ k-times → IR ∣∣ ξ k-linear} Symk(IRd) = {ξ : IRd × · · · × IRd︸ ︷︷ ︸ k-times → IR ∣∣ ξ k-linear and symmetric} be the vector space of k-tensors and symmetric k-tensors, respectively (ac- tually, these are spaces of (0, k)-tensors, but since we only deal with co- variant vectors, we omit the 0). Here ξ ∈ T k(IRd) is called symmetric, 4 if ξ(a1, . . . , ak) = ξ(api(1), . . . , api(k)) for all pi ∈ Sk, where Sk denotes the permutation group of {1, . . . , k}. The case k = 0 corresponds to scalar values, for k = 1, Sym(IRd, IR) = IRd while for k = 2, Sym2(IRd) = Sd×d, i.e. the space corresponds to sym- metric matrices. Note three basic operations for k-tensors. For ξ ∈ T k(IRd) and η ∈ T l(IRd) the tensor product (ξ ⊗ η)(a1, . . . , ak+l) = ξ(a1, . . . , ak)η(ak+1, . . . , ak+l) yields an element of T k+l(IRd). We define the trace tr(ξ) ∈ T k−2(IRd) of ξ ∈ T k(IRd) with k ≥ 2 by tr(ξ)(a1, . . . , ak−2) = d∑ i=1 ξ(ei, a1, . . . , ak−2, ei) where ei denotes the i-th standard basis vector. This operation can be iterated, for example trl(ξ⊗η) for ξ ∈ Symk+l(IRd) and η ∈ Syml(IRd) yields a symmetric k-tensor. Every k-tensor ξ ∈ T k(IRd) can be symmetrized by (9ξ)(a1, . . . , ak) = 1 k! ∑ pi∈Sk ξ(api(1), . . . , api(k)). The symmetrization is a projection, i.e. 92ξ = 9ξ. The spaces T k(IRd) and, consequently, Symk(IRd) will be equipped with the scalar product ξ · η = trk(ξ¯ ⊗ η) = ∑ p∈{1,...,d}k ξ(ep1 , . . . , epk)η(ep1 , . . . , epk), ξ¯(a1, . . . , ak) = ξ(ak, . . . , a1), leading canonically to the norm |ξ| = √ξ · ξ. Again, this is the absolute value for k = 0, for k = 1, this corresponds to the Euclidean norm in IRd and in case k = 2, we can identify ξ ∈ Sym2(IRd) with ξ = ξ11 · · · ξ1d... . . . ... ξ1d · · · ξdd  , |ξ| = ( d∑ i=1 ξ2ii + 2 ∑ i 0. Then, the total generalized variation of order k with weight α for u ∈ L1loc(Ω) is defined as the value of the functional TGVkα(u) = sup {∫ Ω udivk v dx ∣∣∣ v ∈ Ckc (Ω, Symk(IRd)), ‖divl v‖∞ ≤ αl, l = 0, . . . , k − 1 } (3.1) with taking the value ∞ if the respective set is unbounded from above. The space BGVkα(Ω) = {u ∈ L1(Ω) ∣∣ TGVkα(u) <∞} , ‖u‖BGVkα = ‖u‖1 + TGVkα(u) is called the space of functions of bounded generalized variation of order k with weight α. Remark 3.2. For k = 1 and α > 0, we see that TGV1α(u) = α sup {∫ Ω udiv v dx ∣∣∣ v ∈ C1c (Ω, Sym1(IRd)), ‖v‖∞ ≤ 1} = αTV(u). Thus one can indeed speak of a generalization of the total variation. In the following, we will derive some basic properties of the total gener- alized variation. Proposition 3.3. The following statements hold: 1. TGVkα is a semi-norm on the normed space BGV k α(Ω), 10 2. for u ∈ L1loc(Ω), TGVkα(u) = 0 if and only if u is a polynomial of degree less than k, 3. for fixed k and positive weights α, α˜ ∈ IRk, the semi-norms TGVkα and TGVkα˜ are equivalent, 4. TGVkα is rotationally invariant, i.e. for each orthonormal matrix O ∈ IRd×d and u ∈ BGVkα(Ω) we have that u˜ ∈ BGVkα(OTΩ) with TGVkα(u˜) = TGVkα(u), where u˜(x) = u(Ox), 5. for r > 0 and u ∈ BGVkα(Ω), we have, defining u˜(x) = u(rx), that u˜ ∈ BGVkα(r−1Ω) with TGVkα(u˜) = r −d TGVkα˜(u) , α˜l = αlr k−l. Proof. Let us begin proving the first statement. Note that TGV can be interpreted as the dual semi-norm in which the set Kkα(Ω) = {divk v ∣∣ v ∈ Ckc (Ω, Symk(IRd)), ‖divl v‖∞ ≤ αl, l = 0, . . . , k − 1} (3.2) is taking the role of the “predual unit ball”: TGV(u) = sup w∈Kkα ∫ Ω uw dx. (3.3) It is easy to see that Kkα(Ω) is balanced and convex. The former implies that TGVkα is positively one-homogeneous