Infimal convolution regularisation functionals of BV and $\mathrm{L}^{p}$ spaces. Part I: The finite $p$ case

We study a general class of infimal convolution type regularisation functionals suitable for applications in image processing. These functionals incorporate a combination of the total variation ($\mathrm{TV}$) seminorm and $\mathrm{L}^{p}$ norms. A unified well-posedness analysis is presented and a detailed study of the one dimensional model is performed, by computing exact solutions for the corresponding denoising problem and the case $p=2$. Furthermore, the dependency of the regularisation properties of this infimal convolution approach to the choice of $p$ is studied. It turns out that in the case $p=2$ this regulariser is equivalent to Huber-type variant of total variation regularisation. We provide numerical examples for image decomposition as well as for image denoising. We show that our model is capable of eliminating the staircasing effect, a well-known disadvantage of total variation regularisation. Moreover as $p$ increases we obtain almost piecewise affine reconstructions, leading also to a better preservation of hat-like structures.


Introduction
In this paper we introduce a family of novel TV-L p infimal convolution type functionals with applications in image processing: α Du − w M + β w L p (Ω) , α, β > 0 and p > 1.
Here · M denotes the Radon norm of a measure. The functional (1.1) is suitable to be used as a regulariser in the context of variational non-smooth regularisation in imaging applications. We study the properties of (1.1), its regularising mechanism for different values of p and apply it successfully to image denoising.
1.1. Context. After the introduction of the total variation (TV) for image reconstruction purposes [27], the use of non-smooth regularisers has become increasingly popular during the last decades (cf. [8]). They are typically used in the context of variational regularisation, where the reconstructed image is obtained as a solution of a minimisation problem of the type: (1.2) min u 1 s f − T u s L s (Ω) + Ψ(u). The regulariser is denoted here by Ψ. We assume that the data f , defined on a domain Ω ⊂ R 2 , have been corrupted through a bounded, linear operator T and additive (random) noise. Different values of s can be considered for the first term of (1.2), the fidelity term. For example, models incorporating a L 2 fidelity term (resp. L 1 ) have be shown to be efficient for the restoration of images corrupted by Gaussian noise (resp. impulse noise). Of course, other types of noise can also be considered and in those cases the form of the fidelity term Emails: martin.burger@wwu.de, kp366@cam.ac.uk, ep374@cam.ac.uk, cbs31@cam.ac.uk is adjusted accordingly. Typically, one or more parameters within Ψ balance the strength of regularisation against the fidelity term in the minimisation (1.2).
The advantage of using non-smooth regularisers is that the regularised images have sharp edges (discontinuities). For instance, it is a well-known fact that TV regularisation promotes piecewise constant reconstructions, thus preserving discontinuities. However, this also leads to blocky-like artifacts in the reconstructed image, an effect known as staircasing. Recall at this point that for two dimensional images u ∈ L 1 (Ω), the definition of the total variation functional reads: The total variation uses only first-order derivative information in the regularisation process. This can be seen from that fact that for TV(u) < ∞ the distributional derivative Du is a finite Radon measure and TV(u) = Du M . Moreover if u ∈ W 1,1 (Ω) then TV(u) = Ω |∇u| dx, i.e., the total variation is the L 1 norm of the gradient of u. Higher-order extensions of the total variation functional are widely explored in the literature e.g. [4,5,9,11,12,[20][21][22]25]. The incorporation of second-order derivatives is shown to reduce or even eliminate the staircasing effect. The most successful regulariser of this kind is the second order total generalised variation (TGV) introduced by Bredies et al. [5]. Its definition reads Here α, β are positive parameters and BD(Ω) is the space of functions of bounded deformation, i.e., the space of all L 1 (Ω) functions w, whose symmetrised distributional derivative Ew is a finite Radon measure. This is a less regular space than the usual space of functions of bounded variation BV(Ω) for which the full gradient Du is required to be a finite Radon measure. Note that if the variable w in the definition (1.4) is forced to be the gradient of another function then we obtain the classical infimal convolution regulariser of Chambolle-Lions [9]. In that sense TGV can be seen as a particular instance of infimal convolution, optimally balancing first and second-order information.
In the discrete formulation of TGV (as well as for TV) the Radon norm is interpreted as an L 1 norm. The motivation for the current and the follow-up paper is to explore the capabilities of L p norms within first-order regularisation functionals designed for image processing purposes. The use of L p norms for p > 1 has been exploited in different contexts -infinity and p-Laplacian (cf. e.g. [14] and [19] respectively).

Our contribution.
Comparing the definition (1.1) with the definition of TGV in (1.4), we see that the Radon norm of the symmetrised gradient of w has been substituted by the L p norm of w, thus reducing the order of regularisation. Up to our knowledge, this is the first paper that provides a thorough analysis of TV-L p infimal convolution models (1.1) in this generality. We show that the minimisation in (1.1) is well-defined and that TVL p α,β (u) < ∞ if and only if TV(u) < ∞. Hence TVL p α,β regularised images belong to BV(Ω) as desired. In order to get more insight in the regularising mechanism of the TVL p α,β functional we provide a detailed and rigorous analysis of its one dimensional version of the corresponding L 2 fidelity denoising problem (1.5) min For the denoising problem (1.5) with p = 2 we also compute exact solutions for simple one dimensional data. We show that the obtained solutions are piecewise smooth, in contrast to TV (piecewise constant) and TGV (piecewise affine) solutions. Moreover, we show that for p = 2, the 2-homogeneous analogue of the functional (1.1) (1.6) F (u) = min is equivalent to a variant of Huber TV [17], with the functional (1.6) having a close connection with (1.1) itself. Huber total variation is a smooth approximation of total variation and even though it has been widely used in the imaging and inverse problems community, it has not been analysed adequately. Hence, as a by-product of our analysis, we compute exact solutions of the one dimensional Huber TV denoising problem.
We proceed with exhaustive numerical experiments focusing on (1.5). Our analysis is confirmed by the fact that the analytical results coincide with the numerical ones. Furthermore, we observe that even though a first-order regularisation functional is used, we are capable of eliminating the staircasing effect, similarly to Huber TV. By Bregmanising our method [23], we are also able to enhance the contrast of the reconstructed images, obtaining results very similar in quality to the TGV ones. We observe numerically that high values of p promote almost affine structures similar to second-order regularisation methods. We shed more light of this behaviour in the follow-up paper [7] where we study in depth the case p = ∞. Let us finally note that we also consider a modified version of the functional (1.1) where w is restricted to be a gradient of another function leading to the more classical infimal convolution setting. Even though, this modified model is not so successful in staircasing reduction, it is effective in decomposing an image into piecewise constant and smooth parts.
1.3. Organisation of the paper. After the introduction we proceed with the introduction of our model in Section 2. We prove the well-posedness of (1.1), we provide an equivalent definition and we prove its Lipschitz equivalence with the TV seminorm. We finish this section with a well-posedness result of the corresponding TVL p α,β regularisation problem using standard tools.
In Section 3 we establish a link between the TVL p α,β functional and its p-homogeneous analogue (using the p-th power of · L p (Ω) ). The p-homogeneous functional (for p = 2) is further shown to be equivalent to Huber total variation.
We study the corresponding one dimensional model in Section 4 focusing on the L 2 fidelity denoising case. More specifically, after deriving the optimality conditions using Fenchel-Rockafellar duality in Section 4.1, we explore the structure of solutions in Section 4.2. In Section 4.3 we compute exact solutions for the case p = 2, considering a simple step function as data.
In Section 5 we present a variant of our model suitable for image decomposition purposes, i.e., geometric decomposition into piecewise constant and smooth structures. Section 6 focuses on numerical experiments. Confirmation of the obtained one dimensional analytical results is done in Section 6.2, while two dimensional denoising experiments are performed in Section 6.3 using the split Bregman method. There, we show that our approach can lead to elimination of the staircasing effect and we also show that by using a Bregmanised version we can also enhance the contrast, achieving results very close to TGV, a method considered state of the art in the context of variational regularisation. We finish the section with some image decomposition examples and we summarise our results in Section 7.
In the appendix, we remind the reader of some basic facts from the theory of Radon measures and BV functions.

Total variation and L p regularisation
In this section we introduce the TV-L p functional (1.1) as well as some of its main properties. For α, β > 0 and 1 < p ≤ ∞, we define TVL p α,β : L 1 (Ω) → R as follows: TVL p α,β (u) := min The next proposition asserts that the minimisation in (1.1) is indeed well-defined. We omit the proof, which is based on standard coercivity and weak lower semicontinuity techniques: Proposition 2.1. Let u ∈ BV(Ω) with 1 < p ≤ ∞ and α, β > 0. Then the minimum in the definition (1.1) is attained.
Another useful formulation of the definition (1.1) is the dual formulation: The next Proposition shows that the two expressions coincide indeed.
Proposition 2.2. Let u ∈ BV(Ω) and 1 < p ≤ ∞ then Proof. First notice that in (2.1), we can replace C 1 c (Ω) by C 1 0 (Ω), since C 1 c (Ω) = C 1 0 (Ω) with the closure taken with respect to the uniform norm. We define Then, we can rewrite (2.1) as The Fenchel-Rockafellar duality theory, see [13], allows to establish a relation between the primal problem − inf and its dual min Here F * 1 and F * 2 denote the convex conjugate of F 1 and F 2 respectively. In order to obtain such a connection, we follow [2] where it suffices to show that is a closed vector space. Indeed, we have that λ≥0 λ (domF 2 − domF 1 ) ⊂ X, and for every φ ∈ X, we can write φ = λ(λ −1 φ − 0) with λ −1 φ ∞ ≤ α and 0 ∈ domF 1 . Hence, λ≥0 λ (domF 2 − domF 1 ) = X is a closed vector space and there is no duality gap i.e., and similarly, F * 2 (w) = sup Thus the desired equality is proven. Remark 2.3. The dual formulation of TVL p α,β : L 1 (Ω) → R is useful since one can easily derive that TVL p α,β is lower semicontinuous with respect to the strong L 1 topology since it is a pointwise supremum of continuous functions.
The following lemma shows that the TVL p α,β functional is Lipschitz equivalent to the total variation seminorm.
Having shown the basic properties of the TVL p α,β functional, we can use it as a regulariser for variational imaging problems, by minimising where T : L 2 (Ω) → L 2 (Ω) is a bounded, linear operator and f ∈ L 2 (Ω). We conclude our analysis with existence and uniqueness results for the minimisation problem (2.3).
Proof. The proof is a straightforward application of the direct method of calculus of variations. We simply take advantage of (2.2) and the compactness theorem in BV(Ω) along with the lower semicontinuity property of TVL p α,β . We also refer the reader to the corresponding proofs in [25,28].
Since we are mainly interested in studying the regularising properties of TVL p α,β , from now on we focus on the case where s = 2 and T is the identity function (denoising task) where rigorous analysis can be carried out. We thus define the following problem 3. The p-homogeneous analogue and relation to Huber TV Before we proceed to a detailed analysis of the one dimensional version of (P), in this section we consider its p-homogeneous analogue We show in Proposition 3.2 that there is a strong connection between the models (P) and (P p−hom ). The reason for the introduction of P p−hom is that, in certain cases, it is technically easier to derive exact solutions for (P p−hom ) rather than for (P) straightforwardly, see Section 4.3. Moreover, here we can guarantee the uniqueness of the optimal w * , since and thus w * is unique as a minimiser of a strictly convex functional since 1 < p < ∞. Hence, compared to (P), an optimal solution pair of (P p−hom ) is unique. The next proposition says that, unless f is a constant function then the optimal w in (P p−hom ) cannot be zero but nonetheless converges to 0 as β → ∞. In essence, this means that one cannot obtain TV type solutions with the p-homogeneous model.
(Ω) and let (w * , u * ) be an optimal solution pair of the p-homogeneous problem (P p−hom ). Then w * = 0 if and only if f is a constant function.
For general data f , we have that w * → 0 in L p (Ω) for β → ∞.
Proof. It follows immediately that if f is constant then (0, f ) is the optimal pair for (P p−hom ). Suppose now that (w * , u * ) solve (P p−hom ). It is easy to check that the following also hold: In particular, (3.1) implies that Suppose now that w * = 0. Then (3.2) becomes Furthermore, since (0, u * ) solve (P p−hom ), then for every u ∈ C ∞ c (Ω) and > 0, the pair ( ∇u, u * + u) is suboptimal for (P p−hom ), i.e., By dividing the last inequality by and taking the limit → 0 we have that Ω (f −u * )u dx ≤ 0. By considering the analogous perturbations u * − u , we obtain similarly that Ω (f −u * )u dx ≥ 0 and thus Hence u * = f and since u * solves (3.4) this can only happen when f is a constant function. For the last part of the proposition, (supposing f = 0), simply observe that for every u ∈ BV(Ω) and w ∈ L p (Ω) we have that and setting u = w = 0, we obtain and thus w * p L p (Ω) → 0 when β → ∞.
We can further establish a connection between the 1-homogeneous (P) and the p-homogeneous model (P p−hom ): Proposition 3.2. Let 1 < p < ∞ and f ∈ L 2 (Ω) not a constant. A pair (w * , u * ) is a solution of (P p−hom ) with parameters (α, β p−hom ) if and only if it is also a solution of (P) with parameters (α, Proof. Since f is not a constant by the previous proposition we have that w * = 0. Note that for an arbitrary function u ∈ BV(Ω): This means that w * is an admissible solution for both problems (P) and (P p−hom ), with the corresponding set of parameters (α, β 1−hom ) and (α, β p−hom ) respectively. The fact that the same holds for u * as well, comes from the fact that in both problems we have Finally, it turns out that for p = 2, problem (P p−hom ) is essentially equivalent to the widely used Huber total variation regularisation, [17]. We show that in the next proposition. Then where ϕ : Proof. We have So we focus on the minimisation problem (3.6) min Baring in mind that (as it can easily checked) for c ∈ R d and λ > 0, belongs to L ∞ (Ω) ⊂ L 2 (Ω) and solves (3.6) with optimal value equal to 1 α Ω ϕ(∇u) dx.

The one dimensional case
In order to get more insight into the structure of solutions of the problem (P), in this section we study its one dimensional version. As above, we focus on the finite p case, i.e., 1 < p < ∞. The case p = ∞ leads to several additional complications and will be subject of a forthcoming paper [7]. For this section Ω ⊂ R is an open and bounded interval, i.e., Ω = (a, b). Our analysis follows closely the ones in [6] and [24] where the one dimensional L 1 -TGV and L 2 -TGV problems are studied respectively. 4.1. Optimality conditions. In this section, we derive the optimality conditions for the one dimensional problem (P). We initially start our analysis by defining the predual problem (P * ), proving existence and uniqueness for its solutions. We will employ again the Fenchel-Rockafellar duality theory in order to find a connection between their solution pairs.
We define the predual problem (P * ) as (q Hölder conjugate to p): Existence and uniqueness can be verified by standard arguments: , the predual problem (P * ) admits a unique solution in Observe now that we can also write down the predual problem (P * ) using the following equivalent formulation: We denote the infimum in (P * ) as inf P * . Then, it is immediate that The dual problem of (4.1), see [13], is defined as where K here denotes the adjoint of K. Let (σ, τ ) be elements of H 1 0 (Ω) * × H 1 0 (Ω) * acting as distributions. For the convex conjugate of F 1 , we write However, by standard density arguments we have: Hence, we obtain Therefore, we have proved the following: The dual problem of (P * ) is equivalent to the primal problem (P) in the sense that (w, u) ∈ Y * is a solution of the dual of (P * ) if and only if (w, u) ∈ L p (Ω) × BV(Ω) is a solution of (P).
It remains to verify that we have no duality gap between the two minimisation problems (P) and (P * ). The proof of the following proposition follows the proof of the corresponding proposition in [6]. We slightly modify for our case.
is a closed vector space and thus [2] (4.9) min Proof. Let (φ, ψ) ∈ Y and define ψ 0 (x) = c 1 +c 2 x, where c 1 , c 2 are constants that are uniquely determined by the following conditions (4.10) Since there is no duality gap, we can find a relationship between the solutions of (P * ) and (P) via the optimality conditions, see [13,Prop. 4 and (4.12) Proof. Since there is not duality gap, the optimality conditions read [13, Prop. 4.1 (III)]: for every (φ, ξ) and (w, u) solutions of (P * ) and (P) respectively. Note that in dimension one we have H 1 0 (Ω) ⊆ C 0 (Ω). Hence, for every (σ, τ ) ∈ X * , we have the following: and using (A.3) we can simplify the last expressions with: Indeed, the L p norm is an one homogeneous functional and thus its subdifferential reads: Clearly, for w = 0, the above expression reduces to σ L p (Ω) ≥ z, σ , ∀σ ∈ L p (Ω) which is valid for any z ∈ L q (Ω) with z L q (Ω) ≤ 1, i.e., the unit ball of L q (Ω). If w = 0 then the subdifferential reduces to the Gâteaux derivative of the L p norm, i.e., .
Finally, from (4.14) we have for every Combining all the above results, we obtain the optimality conditions (4.11) and (4.12).
see also [26]. On the other hand when w = 0, the additional condition (4.12) depends on the value of p and as we will see later it allows a certain degree of smoothness in the final solution u.

4.2.
Structure of the solutions. The optimality conditions (4.11) and (4.12) are an important asset since we can determine exactly the structure of the solutions for the problem (P) as this is determined by the regularising parameters α, β and the value of p.
We initially discuss the cases where the solution u of (P) is a solution of a corresponding ROF minimisation problem i.e., w = 0.
then the solution u of (P) coincides with the solution of the ROF minimisation problem (4.17) and w = 0.
Proof. Let (w * , u * ) be a solution pair for (P), then for every (w, u) ∈ L p (Ω) × BV(Ω), Setting w = 0, we get Since w ∈ L p (Ω) with p ∈ (1, ∞], using the inequality (A.1) we have that , and using the condition (4.18) we get From (4.19) and (4.20) we conclude that for every u ∈ BV(Ω), In fact what we have essentially proved above is that when the condition (4.18) holds then TVL p α,β (u) = α Du M , ∀u ∈ BV(Ω). Notice also that when (4.18) holds then we can show that w = 0 is an admissible solution but in general we cannot prove that this solution is unique.
The condition (4.18) is valid for any dimension d ≥ 1. It provides a rough threshold for obtaining ROF type solutions in terms of the regularising parameters α, β and the image domain Ω. However, the condition is not sharp since as we will see in the following sections we can obtain a sharper estimate for specific data f .
The following proposition in the spirit of [6,24] gives more insight into the structure of solutions of (P). The above proposition is formulated rigorously via the use of precise representatives of BV functions, see [1], but for the sake of simplicity we rather not get into the details here. Instead we refer the reader to [6,24] where the analogue propositions are shown for the TGV regularised solutions and whose proofs are similar to the one of Proposition 4.7.
We now consider the case where the solution is constant in Ω, which in fact coincides with the mean valuef of the data f: then the solution of (P) is constant and it is equal tof .
Proof. Clearly, if u is a constant solution of (P), then Du = 0 and from (2.2) we get TVL p α,β (u) = 0. Hence, we have u =f . In order to have u =f , from the optimality conditions (4.11) and (4.12), it suffices to find a function φ ∈ H 1 0 (Ω) such that φ(a) = φ(b) = 0 and Hence, it suffices to choose α and β as in (4.22). Fig. 1. Characterisation of solutions of (P) for any data f : The blue/red areas correspond to the ROF type solutions (w = 0) and the purple area corresponds to the TVL p solutions (w = 0) for 1 < p < ∞. We note that the red/purple areas are potentially larger/smaller as the conditions we have derived are not sharp.
In Figure 1, we summarise our results so far. There, we have partitioned the set {α > 0, β > 0} into different areas that correspond to different types of solutions of the problem (P). The brown area, arising from thresholds (4.22) corresponds to the choices of α and β that produce constant solutions while the blue area corresponds to ROF type solutions, according to threshold (4.18). Therefore, we can determine the area where the non-trivial solutions are obtained i.e., w = 0, see purple region. Note that since the conditions (4.18) and (4.22) are not sharp the red and the purple areas are potentially larger or smaller respectively than it is shown in Figure 1.
The following proposition reveals more information about the structure of solutions in the case w = 0. Proposition 4.9 (TVL p -solutions). Let f ∈ BV(Ω) and suppose that (w, u) ∈ L p (Ω)×BV(Ω) is a solution pair for (P) with p ∈ (1, ∞) and w = 0. Suppose that u > f (or u < f ) on an open interval I ⊂ Ω then the solution u of (P) is obtained by .
Hence, by (4.11) we obtain (4.23) where C = β w p−1 L p (Ω) . 4.3. Exact solutions of (P) for a step function. In what follows we compute explicit solutions of the TVL p denoising model (P) for the case p = 2 for a simple data function. We define the step function in Ω = (−L, L), L > 0 as: We first investigate conditions under which we obtain ROF type solutions, that is w = 0.
4.3.1. ROF type solutions. We are initially interested in solutions that respect the discontinuity at x = 0 and are piecewise constant. From the optimality conditions (4.11)-(4.12), it suffices to find a function v ∈ H 1 0 (Ω) such that (4.25) φ and it is also piecewise affine. It is easy to see that by setting φ(x) = α L (L − |x|), the conditions (4.25) are satisfied and the solution u is piecewise constant. The first condition of (4.12) implies that φ L q (Ω) ≤ β ⇔ β α ≥ ( 2L q+1 ) 1 q and provides a necessary and sufficient condition that need to be fulfilled in order for u to be piecewise constant, that is to say A special case of the ROF-type solution is when u is constant, i.e., when u =f , the mean value of f .
We define φ(x) = h 2 (L − |x|) and in that case we have that φ ∞ ≤ α ⇔ α ≥ hL 2 and q . This implies that Using now (4.26)-(4.27) we can draw the exact regions in the quadrant of {α >, β > 0} that correspond to these two types of solutions, see the left graph in Figure 3 for the special case p = 2. Notice that in these regions w = 0 and the estimates are valid for any p ∈ (1, ∞).

TVL 2 type solutions.
For simplicity reasons, we examine here only the case p = 2 with w = 0 in Ω. However, we refer the reader to Section 6.2 where we compute numerically solutions for p = 2. Using Proposition 4.9, we observe that the solution is given by the following second order differential equation: .
Even though we can tell that the solution of (4.28) has an exponential form, the fact that the constraint on C depends on the solution w, creates a difficult computation in order to recover u analytically. In order to overcome this obstacle, we consider the one dimensional version of the 2-homogeneous analogue of (P) that was introduced in Section 3: (4.29) min u∈BV(Ω) w∈L 2 (Ω) Similarly to Section 4.1, one can derive the optimality conditions for (4.29). A pair (w, u) is a solution of (4.29) if and only if there exists a function φ ∈ H 1 0 (Ω) such that In order to recover analytically the solutions of (P) for p = 2 and determine the purple region in Figure 1  Then, we get u(x) = c 1 e kx +c 2 e −kx with φ(x) = c 1 k e kx − c 2 k e −kx +c 3 for all x ∈ (−L, 0]. Firstly, we examine solutions that are continuous which due to symmetry much have the value h 2 at the x = 0, i.e., u(0) = h 2 . Since φ ∈ H 1 0 (−L, L), we have φ(−L) = 0 and also u (−L) = 0. Finally, we require that φ(0) < α. After some computations, we conclude that where c 1 = c 2 e 2kL , c 2 = h 2(e 2kL +1) and k = 1 √ β . On the other hand, in order to get solutions that preserve the discontinuity at x = 0, we require the following: Then we get  Figure 2. On the other hand in the complementary green region we obtain the solution (4.34). For extreme cases where β → ∞, i.e., k → 0 we obtain tanh(kL) k → L, which means that there is an asymptote of g at α = hL 2 . Although, we know the form of the inverse function of the hyperbolic tangent, we cannot compute analytically the inverse f −1 . However, we can obtain an approximation using a Taylor expansion which leads to where α > 0 and α = hL 2 . Finally, we would like to describe the solution on the limiting case β → ∞. Letting β → ∞ in (4.32), we have that c 1 , c 2 → h 2 and u(x) → h 2 for every x ∈ Ω, which in fact is the mean value obtained from (P). For the discontinuous solutions, we have that c 1 , c 2 → α 2L and Fig. 2. Characterisation of solutions of (4.29) for data f being a step function. The green region corresponds to solutions that preserve the discontinuity at x = 0, (4.34), while the blue region corresponds to continuous solutions, (4.32), both having an exponential form.
i.e., we converge to the solution (4.26). We also get that and c 2 is given either from (4.32) or (4.34). Then, in both cases we have w → 0 as k → 0. Observe that the product of β 2−hom w L 2 (Ω) is bounded as β 2−hom → ∞ for both types of solutions and in fact corresponds to the bounds found in (4.26) and (4.27). Indeed, since The last result is yet another verification of Theorem 3.2 and it shows that there is an one to one correspondence, β 2−hom w L 2 (Ω) ↔ β 1−hom and the purple region of Figure 3 is characterised by the solutions obtained in Figure 2.

An image decomposition approach
In this section, we present another formulation for the problem (P), where we decompose an image into a BV part (piecewise constant) and a part that belongs to W 1,p (Ω) (smooth). Let 1 < p ≤ ∞ and Ω ⊂ R d and consider the following minimisation problem: In this way, we can decompose our image into two geometric components.  smoothness that depends on the value of p. In the one dimensional setting, we can prove that the problems (P) and (5.1) are equivalent.
is a solution of (P).
Even though for d = 1 it is true that every L p function can be written as a gradient, this is not true for higher dimensions. In fact, as we show in the following sections, this constraint is quite restrictive and for example the staircasing effect cannot always be eliminated in the denoising process, see for instance Figure 18.
The existence of minimisers of (5.1) is shown following again the same techniques as in Theorem 2.5. Moreover, due to the strict convexity on the fidelity term of (5.1), one can prove that the sum u + v ∈ BV(Ω) is unique for a solution (u, v) ∈ W 1,p (Ω) × BV(Ω). This result coincides with the uniqueness of (P) problem for u. Finally, if (u 1 , v 1 ), (u 2 , v 2 ) are two minimisers of (5.1), then from the convexity of L(u, v) we have for 0 ≤ λ ≤ 1 Since (u 1 , v 1 ), (u 2 , v 2 ) are both minimisers, the above inequality is in fact an equality. Since If we assume that then we contradict the equality on (5.3). Hence, the Minkowski inequality becomes an equality which is equivalent to the existence of µ > 0 such that ∇v 2 = µ∇v 1 . In other words, we have proved the following proposition that was also shown in [18] in a similar context: Proposition 5.2. Let (u 1 , v 1 ), (u 2 , v 2 ) be two minimisers of (5.1). Then ∃µ > 0 such that ∇v 2 = µ∇v 1 . (5.5)

Numerical Experiments
In this section we present our numerical simulations for the problem (P). We begin with the one dimensional case where we verify numerically the analytical solutions obtained in Section 4.3. We also describe the type of structures that are promoted for different values of p. Finally, we proceed to the two dimensional case where we focus on image denoising tasks and in particular on the elimination of the staircasing effect.
We start by defining the discretised version of problem (P) Here TVL p α,β : R n×m → R is defined as where for x ∈ R n×m , we set x p = ( n,m i,j=1 |x(i, j)| p ) 1 p and for x = (x 1 , x 2 ) ∈ (R n×m ) 2 we define We denote by ∇ = (∇ 1 , ∇ 2 ) the discretised gradient with forward differences and zero Neumann boundary conditions defined as where t denotes the step size. The discrete version of the divergence operator is defined as the adjoint of ∇. That is, for every w = (w 1 , w 2 ) ∈ (R n×m ) 2 and u ∈ R n×m , we have that −divw, u = w, ∇u with We solve the minimisation problem (6.1) in two ways. The first one is by using the CVX optimisation package with MOSEK solver (interior point methods). This method is efficient for small-medium scale optimisation problems and thus it is a suitable choice in order to replicate one dimensional solutions. On the other hand, we prefer to solve large scale two dimensional versions of (6.1) with the split Bregman method [15] which has been widely used for the fast solution of non-smooth minimisation problems.
6.1. Split Bregman for L 2 -TVL p . In this section we describe how we adapt the split Bregman algorithm to our discrete model (6.1). Letting z = ∇u − w, the corresponding unconstrained problem becomes Replacing the constraint, using a Lagrange multiplier λ, we obtain the following unconstrained formulation: The Bregman iteration, see [23], that corresponds to the minimisation (6.6) leads to the following two step algorithm: Since solving (6.7) at once is a difficult task, we employ a splitting technique and minimise alternatingly for u, z and w. This yields the split Bregman iteration for our method: Next, we discuss how we solve each of the subproblems (6.9)-(6.11). The first-order optimality condition of (6.9) results into the following linear system: Here A is a sparse, symmetric, positive definite and strictly diagonal dominant matrix, thus we can easily solve (6.13) with an iterative solver such as conjugate gradients or Gauss-Seidel. However, due to the zero Neumann boundary conditions, the matrix A can be efficiently diagonalised by the two dimensional discrete cosine transform, (6.14) where here W nm is the discrete cosine matrix and D = diag(µ 1 , · · · , µ n * m ) is the diagonal matrix of the eigenvalues of A. In that case, A has a particular structure of a block symmetric Toeplitz-plus-Hankel matrix with Toeplitz-plus-Hankel blocks and one can obtain the solution of (6.9) by three operations involving the two dimensional discrete cosine transform [16] as follows: Firstly, we calculate the eigenvalues of A by multiplying (6.14) with e 1 = (1, 0, · · · , 0) from both sides and using the fact that W nm W nm = W nm W nm = I nm , we get Then, the solution of (6.9) is computed exactly by The solution of the subproblem (6.10) is obtained in a closed form via the following shrinkage operator, see also [15,30]. Indeed, for i = 1, 2 we have Finally, we discuss the solution of the subproblem (6.11). In the spirit of [29], we solve (6.11) by a fixed point iteration scheme. Letting κ = β λ and η = −b k − z k+1 + ∇u k+1 , the first-order optimality condition of (6.11) becomes For given w k , we obtain w k+1 by the following fixed point iteration under the convention that 0/0 = 0. We can also consider solving the p-homogenous analogue (P p−hom ), where for certain values of p, e.g. p = 2, we can solve exactly (6.19), since in that case w k+1 i = η i κ+1 . However, we observe numerically that there is no significant computational difference between these two methods. Let us finally mention that since we do not solve exactly all the subproblems (6.9)-(6.11), we do not have a convergence proof for the split Bregman iteration. However in practice, the algorithm converges to the right solutions after comparing them with the corresponding solutions obtained with the CVX package.
6.2. One dimensional results. For this section, we set m = 1 and thus u ∈ R n×1 , w ∈ (R n×1 ) 2 . Initially, we compare our numerical solutions with the analytical ones, obtained in Section 4.3 for the step function, setting p = 2, h = 100, L = 1 and Ω = [−1, 1]. The domain Ω is discretised into 2000 points. We first examine the cases of where ROF solutions are obtained, i.e., the parameters α and β are selected according to the conditions (4.26) and (4.27), see Figure 4. There we see that the analytical solutions coincide with the numerical ones. Now, we proceed by computing the non-ROF solutions. The numerical solutions are solved using the 2-homogeneous analogue of (4.29), since we have proved that the 1-homogeneous and p-homogeneous problems are equivalent modulo an appropriate rescaling of the parameter β, see Proposition 3.2. In fact, as it is described in Figure 3, in order to obtain solutions from the purple region, it suffices to seek solutions for the 2-homogeneous (4.29). Notice also that these solutions are exactly the solutions obtained solving a Huber TV problem, see Proposition 3.3. The analytical solutions are given in (4.32) and (4.34) and are compared with the numerical ones in Figure 5, where we observe that they coincide. We also verify the equivalence between the 1-homogeneous and 2-homogeneous problems where α is fixed and β is obtained from Proposition 3.2, see Figure 5c.
We continue our experiments for general values of p focusing on the geometric behaviour of the solutions as p increases. In order to compare the solutions for p ∈ (1, ∞), we fix the parameter α and choose appropriate values of β and p. We choose α and β so that they belong to the purple region in Figure 3, i.e., β < ( 2L q+1 ) 1 q α and β < h 2 ( 2L q+1 q+1 ) hold), see Figure 6b. We observe that for p = 4 3 , the solution has a similar behaviour to p = 2, but with a steeper gradient at the discontinuity point. Moreover, the solution becomes almost constant near the boundary of Ω. On the other hand, as we increase p, the slope of the solution near the discontinuity point reduces and it becomes almost linear with a relative small constant part near the boundary. The linear structure of the solutions that appears for large p motivates us to examine the case of a piecewise linear data f defined as (6.20) f Figure 7. We set again Ω = [−1, 1], λ = 1 10 and the data are discretised in 2000 points. As we observe, the reconstruction for p = 15 behaves almost linearly everywhere in Ω except near the boundary. In the follow up paper [7], where the case p = ∞ is examined in detail, the occurrence of this linear structure is justified.
In the last part of this section, we discuss the image decomposition approach presented in Section 5. We treat a more complicated one dimensional noiseless signal with piecewise constant, affine and quadratic components and solve the discretised version of (5.1) using CVX under MOSEK. We verify numerically the equivalence between (5.1) and (P) for p = 2, i.e., (∇v, u+v) corresponds to (w, u) where (v, u) and (w, u) are the solutions of (5.1) and (P) respectively, see Figure 8. We also compare the decomposed parts u, v for two different values of p ( 4 3 and 10). In order to have a reasonable comparison on the corresponding solutions, the parameters α, β are selected such that the residual f − u − v 2 is the same for both values of p. As we observe, the v decomposition with p = 4 3 promotes some flatness on the solution compared to p = 2, compare Figures 8b and 9a. On the other hand for p = 10, the v component promotes again almost affine structures, Figure 9b. Notice, that in both cases the v components are continuous. In fact, this is confirmed analytically for every 1 < p < ∞, since in dimension one W 1,p (Ω) ⊂ C(Ω).  6.3. Two dimensional results. In this section we consider the two dimensional case where u ∈ R n×m , w ∈ (R n×m ) 2 with m > 1 and Ω denotes a rectangular/square image domain. We focus on image denoising tasks and on eliminating the staircasing effect for different values of p. We use here the split Bregman algorithm proposed in Section 6.1.   In Figure 11, we present the best reconstructions results in terms of two quality measures, the Peak Signal to Noise Ratio (PSNR) and the Structural Similarity Index (SSIM), see [31] for the definition of the latter. In each case, the values of α and β are selected appropriately for optimal PSNR and SSIM. Our stopping criterion is the relative residual error becoming less than 10 −6 i.e., Finally, for computational efficiency, we fix λ = 10α when 1 < p < 4 and λ = 1000α when p ≥ 4 (empirical rule). We observe that the best reconstructions in terms of the PSNR have no visual difference among p = 3 2 , 2 and 3 and staircasing is present, Figures 11a, 11b and 11c. This is one more indication that the PSNR -which is based on the squares of the difference between the ground truth and the reconstruction -does not correspond to the optimal visual results. However, the best reconstructions in terms of SSIM are visually better. They exhibit significantly reduced staircasing for p = 3 2 and p = 3 and is essentially absent in the case of p = 2, see Figures 11d, 11e and 11f.
We can also get a total staircasing elimination by setting higher values for the parameters α and β, as we show in Figure 12. There, one observes that on one hand as we increase p,  almost affine structures are promoted -see the middle row profiles in Figure 12 -and on the other hand these choices of α, β produce a serious loss of contrast that however can be easily treated via the Bregman iteration. Contrast enhancement via Bregman iteration was introduced in [23], see also [3] for an application to higher-order models. It involves solving a modified version of the minimisation problem. Setting u 0 = f , for k = 1, 2, . . ., we solve (6.22) Instead of solving (6.1) once for fixed α and β, we solve a sequence of similar problems adding back a noisy residual in each iteration which results to a contrast improvement. For stopping criteria regarding the Bregman iteration we refer to [23]. In Figure 13 we present our best Bregmanised results in terms of SSIM. There, we notice that Bregman iteration leads to a significant contrast improvement, in comparison to the results of Figure 12. In fact, we observe that the Bregmanised TVL 2 (first-order regularisation), can achieve reconstructions that are visually close to the second-order Bregmanised TGV 2 , compare Figures 13e and 13f. The second-order TGV 2 and Bregmanised TGV 2 are solved using the Chambolle-Pock primal-dual method [10].
We continue our experimental analysis with a radially symmetric image, see Figure 14. In Figure 15, we demonstrate that we can achieve staircasing-free reconstructions for p = 3 2 , 2, 3 and 7. In fact, as we increase p, we obtain results that preserve the spike in the centre of the circle, see Figure 15d. This provides us with another motivation to examine the p = ∞ case in [7]. The loss of contrast can be again treated using the Bregman iteration (6.22). The best results of the latter in terms of SSIM are presented in Figure 16, for p = 2, 4 and 7 and they are also compared with the corresponding Bregmanised TGV 2 . We observe that we can obtain reconstructions that are visually close to the TGV 2 ones and in fact notice that for p = 7, the spike on the centre of the circle is better reconstructed compared to TGV 2 , see also the surface plots in Figure 17. We conclude with numerical results for the image decomposition approach of Section 5 which we solve again using the split Bregman algorithm. Recall that in dimension two, the solutions of (5.1) will not necessarily be the same with the ones of (P). In fact, we observe that (5.1) cannot always eliminate the staircasing, see for instance Figure 18. Even though, we can easily eliminate the staircasing both in the square and in the circle by applying TVL p regularisation, Figures 18b and 18d, we cannot obtain equally satisfactory results by solving (5.1). While using the latter we can get rid of the staircasing in the circle, Figure 18c, this is not possible for the square, Figure 18a, where we observe -after extensive experimentationthat no values of α and β lead to a staircasing elimination. This is analogous to the difference between TGV 2 and the TV-TV 2 infimal convolution of Chambolle-Lions [9].
However, as we mentioned before, the strength of the formulation (5.1) lies on its ability to efficiently decompose an image into piecewise constant and smooth parts. We show that in Figure 19, for the image in Figure 18c.

Conclusion
We have introduced a novel first-order, one-homogeneous TV-L p infimal convolution type functional for variational image regularisation. The TVL p functional constitutes a very general class of regularisation functionals exhibiting diverse smoothing properties for different choices of p. In the case p = 2 the well-known Huber TV regulariser is recovered.
We studied the corresponding one dimensional denoising problem focusing on the structure of its solutions. We computed exact solutions of this problem for the case p = 2 for simple one     dimensional data. Hence, as an additional novelty in our paper we presented exact solutions of the one dimensional Huber TV denoising problem.
Numerical experiments for several values of p indicate that our model leads to an elimination of the staircasing effect. We show that we can further enhance our results by increasing the contrast via a Bregman iteration scheme and thus obtaining results of similar quality to those of TGV 2 . Furthermore, as p increases the structure of the solutions changes from piecewise  smooth to piecewise linear and the model, in contrast to TGV 2 , is capable of preserving sharp spikes in the reconstruction. This observation motivates a more detailed study of the TVL p functionals for large p and in particular for the case p = ∞.
This concludes the first part of the study of the TV-L p model for p < ∞. The second part [7], is devoted to the p = ∞ case. There we explore further, both in an analytical and an experimental level, the capability of the TVL ∞ model to promote affine and spike-like structures in the reconstructed image and we discuss several applications.      ∇u denotes the Radon-Nikodym derivative of D a u with respect to L d . When d = 1, ∇u is simply denoted by u . We will also use the following basic inequality regarding inclusions of L p spaces (A.1) h L p 1 (Ω) ≤ |Ω| Unless otherwise stated q denotes the Hölder conjugate of the exponent p, i.e., Regarding the subdifferential of the Radon norm we have that it can be characterised, at least for C 0 functions, as follows [6] (