Bilevel Parameter Learning for Higher-Order Total Variation Regularisation Models

Journal of Mathematical Imaging and Vision, Jun 2016

We consider a bilevel optimisation approach for parameter learning in higher-order total variation image reconstruction models. Apart from the least squares cost functional, naturally used in bilevel learning, we propose and analyse an alternative cost based on a Huber-regularised TV seminorm. Differentiability properties of the solution operator are verified and a first-order optimality system is derived. Based on the adjoint information, a combined quasi-Newton/semismooth Newton algorithm is proposed for the numerical solution of the bilevel problems. Numerical experiments are carried out to show the suitability of our approach and the improved performance of the new cost functional. Thanks to the bilevel optimisation framework, also a detailed comparison between \(\text {TGV}^2\) and \(\text {ICTV}\) is carried out, showing the advantages and shortcomings of both regularisers, depending on the structure of the processed images and their noise level.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

Bilevel Parameter Learning for Higher-Order Total Variation Regularisation Models

Bilevel Parameter Learning for Higher-Order Total Variation Regularisation Models J. C. De los Reyes 0 1 C.-B. Schönlieb 0 1 T. Valkonen 0 1 0 Department of Mathematical Sciences, University of Liverpool , Liverpool , UK 1 Department of Applied Mathematics and Theoretical Physics, University of Cambridge , Cambridge , UK We consider a bilevel optimisation approach for parameter learning in higher-order total variation image reconstruction models. Apart from the least squares cost functional, naturally used in bilevel learning, we propose and analyse an alternative cost based on a Huber-regularised TV seminorm. Differentiability properties of the solution operator are verified and a first-order optimality system is derived. Based on the adjoint information, a combined quasiNewton/semismooth Newton algorithm is proposed for the numerical solution of the bilevel problems. Numerical experiments are carried out to show the suitability of our approach and the improved performance of the new cost functional. Thanks to the bilevel optimisation framework, also a detailed comparison between TGV2 and ICTV is carried out, showing the advantages and shortcomings of both regularisers, depending on the structure of the processed images and their noise level. Bilevel optimisation; Total variation regularisers; Image quality measures 1 Introduction In this paper, we propose a bilevel optimisation approach for parameter learning in higher-order total variation regu Research Center on Mathematical Modelling (MODEMAT), Escuela Politécnica Nacional, Quito, Ecuador α R(u) + d(K (u), f ), (1.1) where R is a regularising energy that models a-priori knowledge about the image u, d(·, ·) is a suitable distance function that models the relation of the data f to the unknown u, and α > 0 is a parameter that balances our trust in the forward model against the need of regularisation. The parameter α, in particular, depends on the amount of ill-posedness in the operator K and the amount (amplitude) of the noise present in f . A key issue in imaging inverse problems is the correct choice of α, image priors (regularisation functionals R), fidelity terms d and (if applicable) the choice of what to measure (the linear or non-linear operator K ). Depending on this choice, different reconstruction results are obtained. While functional modelling (1.1) constitutes a mathematically rigorous and physical way of setting up the reconstruction of an image—providing reconstruction guarantees in terms of error and stability estimates—it is limited with respect to its adaptivity for real data. On the other hand, data-based modelling of reconstruction approaches is set up to produce results which are optimal with respect to the given data. However, in general, it neither offers insights into the structural properties of the model nor provides comprehensible reconstruction guarantees. Indeed, we believe that for the development of reliable, comprehensible and at the same time effective models (1.1), it is essential to aim for a unified approach that seeks tailor-made regularisation and data models by combining model- and data-based approaches. To do so, we focus on a bilevel optimisation strategy for finding an optimal setup of variational regularisation models (1.1). That is, for a given training pair of noisy and original clean images ( f, f0), respectively, we consider a learning problem of the form min F (u∗) = cost (u∗, f0) subject to α u∗ ∈ arg min {α R(u) + d(K (u), f )} , (1.2) u where F is a generic cost functional that measures the fitness of u∗ to the training image f0. The argument of the minimisation problem will depend on the specific setup (i.e. the degrees of freedom) in the constraint problem (1.1). In particular, we propose a bilevel optimisation approach for learning optimal parameters in higher-order total variation regularisation models for image reconstruction in which the arguments of the optimisation constitute parameters in front of the first- and higher-order regularisation terms. Rather than working on the discrete problem, as is done in standard parameter learning and model optimisation methods, we optimise the regularisation models in infinitedimensional function space. The resulting problems are difficult to treat due to the non-smooth structure of the lower level problem, which makes it impossible to verify standard constraint qualification conditions for Karush–Kuhn–Tucker (KKT) systems. Therefore, in order to obtain characterising first-order necessary optimality conditions, alternative analytical approaches have emerged, in particular regularisation techniques [ 4, 20, 28 ]. We consider such an approach here and study the related regularised problem in depth. In particular, we prove the Fréchet differentiability of the regularised solution operator, which enables to obtain an optimality condition for the problem under consideration and an adjoint state for the efficient numerical solution of the problem. The bilevel problems under consideration are related to the emerging field of generalised mathematical programmes with equilibrium constraints (MPEC) in function space. Let us remark that even for finite-dimensional problems, there are few recent references dealing with stationarity conditions and solution algorithms for this type of problems (see, e.g. [ 18, 30, 33, 34, 38 ]). Let us give an account to the state of the art of bilevel optimisation for model learning. In machine learning, bilevel optimisation is well established. It is a semi-supervised learning method that optimally adapts itself to a given dataset of measurements and desirable solutions. In [ 15, 23, 43 ], for instance, the authors consider bilevel optimisation for finite-dimensional Markov random field models. In inverse problems, the optimal inversion and experimental acquisition setup is discussed in the context of optimal model design in works by Haber, Horesh and Tenorio [ 25, 26 ], as well as Ghattas et al. [ 3, 9 ]. Recently, parameter learning in the context of functional variational regularisation models (1.1) also entered the image processing community with works by the authors [ 10, 22 ], Kunisch, Pock and co-workers [ 14, 33 ], Chung et al. [ 16 ] and Hintermüller et al. [ 30 ]. Apart from the work of the authors [ 10, 22 ], all approaches so far are formulated and optimised in the discrete setting. Our subsequent modelling, analysis and optimisation will be carried out in function space rather than on a discretisation of (1.1). While digitally acquired image data are of course discrete, the aim of high-resolution image reconstruction and processing is always to compute an image that is close to the real (analogue, infinite dimensional) world. Hence, it makes sense to seek images which have certain properties in an infinite dimensional function space. That is, we aim for a processing method that accentuates and preserves qualitative properties in images independent of the resolution of the image itself [45]. Moreover, optimisation methods conceived in function space potentially result in numerical iterative schemes which are resolution and mesh independent upon discretisation [ 29 ]. Higher-order total variation regularisation has been introduced as an extension of the standard total variation regulariser in image processing. As the Total Variation (TV) [ 41 ] and many more contributions in the image processing community have proven, a non-smooth first-order regularisation procedure results in a non-linear smoothing of the image, smoothing more in homogeneous areas of the image domain and preserving characteristic structures such as edges. In particular, the TV regulariser is tuned towards the preservation of edges and performs very well if the reconstructed image is piecewise constant. The drawback of such a regularisation procedure becomes apparent as soon as images or signals (in 1D) are considered which do not only consist of constant regions and jumps but also possess more complicated, higher-order structures, e.g. piecewise linear parts. The artefact introduced by TV regularisation in this case is called staircasing [ 40 ]. One possibility to counteract such artefacts is the introduction of higher-order derivatives in the image regularisation. Chambolle and Lions [ 11 ], for instance, propose a higher-order method by means of an infimal convolution of the TV and the TV of the image gradient called Infimal Convolution Total Variation (ICTV) model. Other approaches to combine first- and second-order regularisation originate, for instance, from Chan et al. [ 12 ] who consider total variation minimisation together with weighted versions of the Laplacian, the Euler-elastica functional [ 13,37 ], which combines total variation regularisation with curvature penalisation, and many more [ 35,39 ] just to name a few. Recently, Bredies et al. have proposed Total Generalized Variation (TGV) [ 5 ] as a higher-order variant of TV regularisation. In this work, we mainly concentrate on two second-order total variation models: the recently proposed TGV [ 5 ] and the ICTV model of Chambolle and Lions [ 11 ]. We focus on second-order TV regularisation only since this is the one which seems to be most relevant in imaging applications [ 6, 31 ]. For ⊂ R2 open and bounded and u ∈ B V ( ), the ICTV regulariser reads ICTVα,β (u) := min v∈W 1,1( ), ∇v∈BV ( ) TGV2α,β (u) := w∈mB Din( ) α Du − w stands for the total variation of u in , BD( ) := {w ∈ L1( ; Rn) | E w M( ;Rn×n) < ∞} is the space of vector fields of bounded deformation on , E denotes the symmetrised gradient and Sym2(R2) denotes the space of symmetric tensors of order 2 with arguments in R2. The parameters α, β are fixed positive parameters and will constitute the arguments in the special learning problem á la (1.2) we (1.4) (1.5) consider in this paper. The main difference between (1.3) and (1.4) is that we do not generally have that w = ∇v for any function v. That results in some qualitative differences of ICTV and TGV regularisation, compare for instance [ 1 ]. Substituting α R(u) in (1.1) by TGV2α,β (u) or ICTVα,β (u) gives the TGV image reconstruction model and the ICTV image reconstruction model, respectively. In this paper, we only consider the case K = I d identity and d(u, f ) = u − f 2L2( ) in (1.1) which corresponds to an image denoising model for removing Gaussian noise. With our choice of regulariser, the former scalar α in (1.1) has been replaced by a vector (α, β) of two parameters in (1.3) and (1.4). The choice of the entries in this vector now do not only determine the overall strength of the regularisation (depending on the properties of K and the noise level), but those parameters also balance between the different orders of regularity of the function u, and their choice is indeed crucial for the image reconstruction result. Large β will give regularised solutions that are close to TV regularised reconstructions, compare Fig. 1. Large α will result in TV2 type solutions, that is solutions that are regularised with TV of the gradient [ 27,39 ], compare Fig. 2. With our approach described in the next section, we propose a learning approach for choosing those parameters optimally, in particular optimally for particular types of images. For the existence analysis of an optimal solution as well as for the derivation of an optimality system for the corresponding learning problem (1.2), we will consider a smoothed version of the constraint problem (1.1)—which is the one in fact used in the numerics. That is, we replace R(u)—being TV, TGV or ICTV in this paper—by a Huber-regularised version and add an H 1 regularisation with a small weight to (1.1). In this setting and under the special assumption of box constraints on α and β, we provide a simple existence proof for an optimal solution. A more general existence result that holds also for the original non-smooth problem and does not require box constraints is derived in [ 19 ], and we refer the reader to this paper for a more sophisticated analysis on the structure of solutions. (a) Too low α, low β. (b) Too low α, optimal β. (c) Too high α, high β. Good match to noisy data optimal T V 2-like behaviour Bad TV2-like behaviour A main challenge in the setup of such a learning approach is to decide what is the best way to measure fitness (optimality) of the model. In our setting this amounts to choosing an appropriate distance F in (1.2) that measures the fitness of reconstructed images to the ‘perfect’, noise-free images in an appropriate training set. We have to formalise what we mean by an optimal reconstruction model. Classically, the difference between the original, noise-free image f0 and its regularised version uα,β is computed with an L22 cost functional FL2 (uα,β ) := 2 uα,β − f0 2L2( ), which is closely related to the PSNR quality measure. Apart from this, we propose in this paper an alternative cost functional based on a Huberised total variation cost FL1η∇ (uα,β ) := |D(uα,β − f0)|γ d x , where the Huber regularisation | · |γ will be defined later on in Definition 2.1. We will see that the choice of this cost functional is indeed crucial for the qualitative properties of the reconstructed image. The proposed bilevel approach has an important indirect consequence: It establishes a basis for the comparison of the different total variation regularisers employed in image denoising tasks. In the last part of this paper, we exhaustively compare the performance of TV, TGV2 and ICTV for various image datasets. The parameters are chosen optimally, according to the proposed bilevel approach, and different quality measures (like PSNR and SSIM) are considered for the comparison. The obtained results are enlightening about when to use each one of the considered regularisers. In particular, ICTV appears to behave better for images with arbitrary structure and moderate noise levels, whereas TGV2 behaves better for images with large smooth areas. Outline of the paper In Sect. 2, we state the bilevel learning problem for the two higher-order total variation regularisation models, TGV and ICTV, and prove existence of an optimal parameter pair α, β. The bilevel optimisation problem is analysed in Sect. 3, where existence of Lagrange multipliers is proved and an optimality system, as well as a gradient formula, is derived. Based on the optimality condition, a BFGS algorithm for the bilevel learning problem is devised in Sect. 5.1. For the numerical solution of each denoising problem, an infeasible semismooth Newton method is considered. Finally, we discuss the performance of the parameter learning method by means of several examples for the denoising of natural photographs in Sect. 5. Therein, we also present a statistical analysis on how TV, ICTV and TGV regularisation compare in terms of returned image quality, carried out on 200 images from the Berkeley segmentation dataset BSDS300. 2 Problem Statement and Existence Analysis We strive to develop a parameter learning method for higherorder total variation regularisation models that maximises the fit of the reconstructed images to training images simulated for an application at hand. For a given noisy image f ∈ L2( ), ⊂ R2 open and bounded, we consider 1 muin Rα,β (u) + 2 u − f 2L2( ) . where, α, β ∈ R. We focus on TGV2, (2.1) Rα,β (u) = TGV2α,β (u) := w∈mB Din( ) α (Du − w) for u ∈ B V ( ). For these models, we want to determine the optimal choice of α, β, given a particular type of images and a fixed noise level. More precisely, we consider a training pair ( f, f0), where f is a noisy image corrupted by normally distributed noise with a fixed variation, and the image f0 represents the ground truth or an image that approximates the ground truth within a desirable tolerance. Then, we determine the optimal choice of α, β by solving the following problem: min (α,β)∈R2 F (uα,β ) s.t. α, β ≥ 0, where F equals the L22 cost (1.6) or the Huberised TV cost (1.7) and uα,β for a given f solves a regularised version of the minimisation problem (2.1) that will be specified in the next section, compare problem (2.3b). This regularisation of the problem is a technical requirement for solving the bilevel problem that will be discussed in the sequel. In contrast to learning α, β in (2.1) in finite dimensional parameter spaces (as is the case in machine learning), we consider optimisation techniques in infinite dimensional function spaces. 2.1 Formal Statement Let ⊂ Rn be an open bounded domain with Lipschitz boundary. This will be our image domain. Usually = (0, w) × (0, h) for w and h the width and height of a twodimensional image, although no such assumptions are made in this work. Our data f and f0 are assumed to lie in L2( ). In our learning problem, we look for parameters (α, β) that for some cost functional F : H 1( ) → R solve the problem 1 J γ,μ(u; α, β) := 2 u − f 2L2( ) + Rαγ,,βμ(u). for u ∈ H 1( ), Here J γ,μ(·; α, β) is the regularised denoising functional that amends the regularisation term in (2.1) by a Huberregularised version of it with parameter γ > 0, and an elliptic regularisation term with parameter μ > 0. In the case of TGV2, the modified regularisation term Rαγ,,βμ(u) then reads, F (uα,β ) min (α,β)∈R2 subject to α, β ≥ 0, where uα,β ∈ arg min J γ,μ(u; α, β) u∈H1( ) (2.3a) (2.3b) (2.3c) TGV2α,,γβ,μ(u) := min w∈H1( ) + μ + 2 α |Du − w|γ dx β |E w|γ dx u 2H1( ) + w 2H1( ) , and in the case of ICTV, we have (2.2) ICTVγα,,βμ(u) := α |Du − ∇v|γ dx min v∈W 1,1( ) ∇v∈BV ( ,Rn)∩H1( ) β |D∇v|γ dx + μ + 2 u 2H1( ) + ∇v 2H1( ) . Here, H1( ) = H 1( ; Rn) and the Huber regularisation |·|γ is defined as follows. Definition 2.1 Given γ ∈ (0, ∞], we define for the norm · 2 on Rm , the Huber regularisation (2.4) |g|γ = g 2 − 21γ , For the cost functional F , given noise-free data f0 ∈ L2( ) and a regularised solution u ∈ H 1( ), we consider in particular the L2 cost 1 FL2 (u) = 2 2 f0 − u 2L2( ;Rd ), as well as the Huberised total variation cost FL1η∇ (u) = |D( f0 − u)|γ d x with noise-free data f0 ∈ BV( ). Remark 2.1 Please note that in our formulation of the bilevel problem (2.3), we only impose a non-negativity constraint on the parameters α and β, i.e. we do not strictly bound them away from zero. There are two reasons for that. First, for the existence analysis of the smoothed problem, the case α = β = 0 is not critical since compactness can be secured by the H 1 term in the functional, compare Sect. 2.2. Second, in [ 19 ], we indeed prove that even for the non-smooth problem (as μ → 0), under appropriate assumptions on the given data, the optimal α, β are guaranteed to be strictly positive. in particular this holds for u = 0. Hence, 1 1 2 un − f 2L2( ) + Rαγn,μ,βn (un) ≤ 2 f 2L2( ). Exemplarily, we consider here the case for the TGV reguαn,βn = TGV2α,nγ,β,μn . The proof for the ICTV lariser, that is Rγ,μ regulariser can be done in a similar fashion. Inequality (2.5) in particular gives (2.5) 2 un H1( ) + 2 1 wn H1( ) ≤ μ f L2( ), where wn is the optimal w for un. This gives that (un, wn ) is uniformly bounded in H 1( ) × H1( ) and that there exists a subsequence {(αn, βn , un , wn)} which converges weakly in R2 × H 1( )×H1( ) to a limit point (αˆ , βˆ, uˆ, wˆ ). Moreover, un → uˆ strongly in L p( ) and wn → wˆ in L p( ; Rn ). Using the continuity of the L2 fidelity term with respect to strong convergence in L2, and the weak lower semicontinuity of the H 1 term with respect to weak convergence in H 1 and of the Huber-regularised functional even with respect to weak∗ convergence in M (cf. Lemma 2.1), we get αˆ |Duˆ − wˆ |γ dx + βˆ |E w|γ dx 2.2 Existence of an Optimal Solution The existence of an optimal solution for the learning problem (2.3) is a special case of the class of bilevel problems considered in [ 19 ], where the existence of optimal parameters in (0, +∞]2N is proven. For convenience of the reader, we provide a simplified proof for the case where additional box constraints on the parameters are imposed. We start with an auxiliary lower semicontinuity result for the Huberregularised functionals. Lemma 2.1 Let u, v ∈ L p( ), 1 ≤ p < ∞. Then, the functional u → |u − v|γ d x , where | · |γ is the Huber regularisation in Definition 2.1, is lower semicontinuous with respect to weak* convergence in M( ; Rd ) Proof Recall that for g ∈ Rm , the Huber-regularised norm may be written in dual form as γ |g|γ = sup q, g − 2 q 22 : q 2 ≤ 1 . Therefore, we find that G(u) := |u − v|γ d x = sup u(x ) · ϕ(x ) dx − γ 2 ϕ(x ) 22 dx : ϕ ∈ Cc∞( ), ϕ(x ) 2 ≤ 1 for every x ∈ . The functional G is of the form G(u) = sup{ u, ϕ −G∗(ϕ)}, i ∞ where G∗ is the convex conjugate of G. Now, let {u }i=1 converge to u weakly* in M( ; Rd ). Taking a supremising sequence {ϕ j } ∞j=1 for this functional at any point u, we easily see lower semicontinuity by considering the sequences { ui , ϕ j − G∗(ϕ j )}i∞=1 for each j . Our main existence result is the following. Theorem 2.1 We consider the learning problem (2.3) for TGV2 and ICTV regularisation, optimising over parameters (α, β) such that 0 ≤ α ≤ α¯ , 0 ≤ β ≤ β¯. Here (α¯ , β¯) < ∞ is an arbitrary but fixed vector in R2 that defines a box constraint on the parameter space. There exists an optimal solution (αˆ , βˆ) ∈ R2 for this problem for both choices of cost functionals, F = FL2 and F = FL1η∇ . 2 Proof Let (αn, βn ) ⊂ R2 be a minimising sequence. Due to the box constraints we have that the sequence (αn, βn ) is bounded in R2. Moreover, we get for the corresponding sequences of states un := u(αn,βn) that J γ,μ(un; αn , βn ) ≤ J γ,μ(u; αn , βn ), ∀u ∈ H 1( ), where in the last step we have used the boundedness of the sequence Rαγn,μ,βn (un) from (2.5) and the convergence of (αn, βn ) in R2. This shows that the limit point uˆ is an optimal solution for (αˆ , βˆ). Moreover, due to the weak lower semicontinuity of the cost functional F and the fact that the set {(α, β) : 0 ≤ α ≤ α¯ , 0 ≤ β ≤ β¯} is closed, we have that (αˆ , βˆ, uˆ) is optimal for (2.3). Remark 2.2 • Using the existence result in [ 19 ], in principle we could allow infinite values for α and β. This would include both TV2 and TV as possible optimal regularisers in our learning problem. 1 2 uˆ − f 2L2( ) + uˆ 2H1( ) + wˆ 2H1( ) 1 ≤ limninf 2 un − f 2L2( ) μ + 2 + μ + 2 + μ + 2 αˆ |Dun − wn|γ dx + βˆ |E wn|γ dx 2 un H1( ) + 2 wn H1( ) 1 ≤ limninf 2 un − f 2L2( ) + βn |E wn|γ dx 2 un H1( ) + 2 wn H1( ) , αn |Dun − wn|γ dx • In [ 19 ], in the case of the L2 cost and assuming that Rαγ,β ( f ) > Rαγ,β ( f0), we moreover show that the parameters (α, β) are strictly larger than 0. In the case of the Huberised TV cost, this is proven in a discretised setting. Please see [ 19 ] for details. • The existence of solutions with μ = 0, that is without elliptic regularisation, is also proven in [ 19 ]. Note that here, we focus on the μ > 0 case since the elliptic regularity is required for proving the existence of Lagrange multipliers in the next section. Remark 2.3 In [ 19 ], it was shown that the solution map of our bilevel problem is outer semicontinuous. This implies, in particular, that the minimisers of the regularised bilevel problems converge towards the minimiser of the original one. 3 Lagrange Multipliers In this section, we prove the existence of Lagrange multipliers for the learning problem (2.3) and derive an optimality system that characterises stationary points. Moreover, a gradient formula for the reduced cost functional is obtained, which plays an important role in the development of fast solution algorithms for the learning problems (see Sect. 5.1). In what follows, all proofs are presented for the TGV2 regularisation case, that is Rα,β = TGV2α,,γβ . However, posγ sible modifications to cope with the ICTV model will also be commented. Moreover, we consider along this section a smoother variant of the Huber regularisation, given by Uγ Aγ Bγ Cγ 1 ⎧⎪ |g| + γ2 Lγ − 2 + γ 2 + γ 3 + 3γ 4 3 + 4γ 2 ⎪ |g|γ = ⎨ Aγ |g| + B2γ |g|2 + C3γ |g|3 + Dγ This modified Huber function is required in order to get differentiability of the solution operator, a matter which is investigated next. ⎪⎪⎩ γ2 |g|2 with 1 Uγ = γ γ Aγ = 1 − 2 2γ + 1 2 2γ , γ γ 3 Bγ = 2 (2γ + 1), Cγ = − 2 , Dγ = − γ33 L3γ − Aγ Lγ . 1 1 + 2γ 1 , Lγ = γ 3.1 Differentiability of the Solution Operator We recall that the TGV2 denoising problem can be rewritten as y = (u, w) = arg min BV ( )×B D( ) 2 1 |u − f |2 + α|Du − w|γ + β|E w|γ . Using an elliptic regularisation, we then get y = arg min H1( )×H1( ) 2 1 1 a(y, y) + 2 |u − f |2 + α|Du − w|γ + β|E w|γ , where a(y, y) = μ u 2H1 + w 2H1 . A necessary and sufficient optimality condition for the latter is then given by the following variational equation: a(y, ) + αhγ (Du − w)(Dφ − ϕ) dx + Theorem 3.1 The solution operator S : R2 → Y , which assigns to each pair (α, β) ∈ R2 the corresponding solution to the denoising problem (3.1), is Fréchet differentiable and its derivative is characterised by the unique solution z = S (α, β)[θ1, θ2] ∈ Y of the following linearised equation: a(z, ) + − hγ (Du − w) (Dφ − ϕ) dx − − β hγ (Ew+) − hγ (Ew) − hγ (Ew)Eδw : Eϕ dx θ2 hγ (Ewα+θ ) − hγ (Ew) : Eϕ dx, for all ∈ Y. (3.3) Testing with get that = ξ and using the monotonicity of hγ (·), we ξ Y ≤ C |α| hγ (Du+ − w+) − hγ (Du − w) αhγ (Du − w)(Dz1 − z2)(Dφ − ϕ) dx θ2 hγ (E w)E ϕ dx βhγ (E w)E z2 E ϕ dx z1φ dx = 0, for all ∈ Y. Proof Thanks to the ellipticity of a(·, ·) and the monotonicity of hγ , the existence of a unique solution to the linearised equation follows from the Lax-Milgram theorem. Let ξ := y+ − y − z, where y = S(α, β) and y+ = S(α + θ1, β + θ2). Our aim is to prove that ξ Y = o(|θ |). Combining the equations for y+, y and z we get that a(ξ, ) + (α + θ1) hγ (Du+ − w+)(Dφ − ϕ) dx α hγ (Du − w)(Dφ − ϕ) dx where ξ := (ξ1, ξ2) ∈ H 1( ) × H1( ). Adding and subtracting the terms αhγ (Du − w)(Dδu − δw)(Dφ − ϕ) dx and βhγ (E w)E δw : E ϕ dx , where δu := uα+θ − u and δw := wα+θ − w, we obtain that a(ξ, ) + αhγ (Du − w)(Dξ1 − ξ2)(Dφ − ϕ) + = − βhγ (Ew)Eξ2 : Eϕ dx + 2 ξ1φ dx α hγ (Du+ − w+) − hγ (Du − w) − hγ (Du − w)(Dδu − δw) (Dφ − ϕ) − θ1 hγ (Du+ − w+) + + + + − − − + − − hγ (Du − w)(Dδu − δw) L2 + |θ1| hγ (Du+ − w+) − hγ (Du − w) L2 + |β| hγ (E w+) − hγ (E w) − hγ (E w)E δw L2 + |θ2| hγ (E wα+θ ) − hγ (E w) L2 , for some generic constant C > 0. Considering the differentiability and Lipschitz continuity of hγ (·), it then follows that ξ Y ≤ C |α| o( y+ − y 1,p) + |θ1| yα+θ − y Y + |β| o( w+ − w 1,p) + |θ2| wα+θ − w H1( ) , (3.4) where · 1, p stands for the norm in the space W1, p( ). From regularity results for second-order systems (see [24, Theorem 1, Remark 14]), it follows that y+ − y 1, p ≤ L|θ |, ξ Y = o(|θ |). ≤ L|θ | Div hγ (Du − w) −1, p + hγ (Du − w) −1, p + Div hγ (E w) −1, p ≤ L|θ | 2 hγ (Du − w) L∞ + hγ (E w) L∞ since |hγ (·)| ≤ 1. Inserting the latter in estimate (3.4), we finally get that Remark 3.1 The extra regularity result for second-order systems used in the last proof and due to Gröger [24, Thm. 1, Rem. 14] relies on the properties of the domain . The result was originally proved for C 2 domains. However, the regularity of the domain (in the sense of Gröger) may also be verified for convex Lipschitz bounded domains [ 17 ], which is precisely our image domain case. Remark 3.2 The Fréchet differentiability proof makes use of the quasilinear structure of the TGV2 variational form, making it difficult to extend to the ICTV model without further regularisation terms. For the latter, however, a Gâteaux differentiability result may be obtained using the same proof technique as in [ 22 ]. 3.2 The Adjoint Equation Next, we use the Lagrangian formalism for deriving the adjoint equations for both the TGV2 and ICTV learning problems. The existence of a solution to the adjoint equation follows from the Lax-Milgram theorem. Defining the Lagrangian associated to the TGV2 learning problem by L(u, w, α, β, p1, p2) = F (u) + μ(u, p1)H1 + μ(w, p2)H1 αhγ (Du − w)(Dp1 − p2) βhγ (E w)E p2 + (u − f ) p1, and taking the derivative with respect to the state variable (u, w), we get the necessary optimality condition L(u,w)(u, w, α, β, p1, p2)[(δu , δw)] = F (u)δu + μ( p1, δu )H1 + μ( p2, δw)H1 + + αhγ (Du − w)(Dδu − δw)(Dp1 − p2) βhγ (E w)E δw E p2 + p1δu = 0. If δw = 0, then Theorem 3.2 Let (u, w) ∈ H 1( ) × H1( ). There exists a unique solution = ( p1, p2) ∈ Y = H 1( ) × H1( ) to the adjoint system μ( , δy )Y + αhγ (Du − w)(Dδu − δw)(Dp1 − p2) + βhγ (E w)E δw E p2 + p1δv = −F (u)δu , for all δy ∈ Y. (3.6) (3.7) p1δu = −∇u F (u)δu , for all δu ∈ H 1( ), (3.5) and the correspondent Lagrangian functional L is given by The corresponding solution is called adjoint state associated to (v, w). Proof We have to show that the left-hand side of equation (3.7) constitutes a bilinear, continuous and coercive form on Y × Y . Linearity and continuity follows immediately. For the coercivity, let us take δy = . Since hγ is a monotone function, the terms αhγ (Du − w)(Dp1 − p2)(Dp1 − p2) and βhγ (E w)E p2 E p2 become positive, yielding μ 2 Y + αhγ (Du − w)(Dp1 − p2)(Dp1 − p2) + βhγ (E w)E p2 E p2 + 2 p1 ≥ μ 2Y . Thus, coercivity holds and, using Lax-Milgram theorem, we conclude that there exists a unique solution to the adjoint system (3.7). Remark 3.3 For the ICTV model, it is possible to proceed formally with the Lagrangian approach. We recall that a necessary and sufficient optimality condition for the ICTV functional is given by μ(u, φ)H1 +μ(∇v, ∇ϕ)H1 + αhγ (Du−∇v)(Dφ−∇ϕ) + βhγ (D∇v)D∇ϕ + (u − f )φ = 0, for all (φ, ϕ) ∈ H 1( ) × H1( ) (3.8) L(u, v, α, β, p1, p2) = F (u) + μ(u, p1)H1 +μ(∇v, ∇ p2)H1 + αhγ (Du − ∇v)(Dp1 − ∇ p2) + βhγ (D∇v)D∇ p2 + (u − f ) p1. Deriving the Lagrangian with respect to the state variables (u, v) and setting it equal to zero yields L(u,v)(u, v, α, β, p1, p2)[(δu , δv)] = F (u)δu + μ( p1, δu )H1 + μ(∇ p2, ∇δv)H1 + + αhγ (Du − ∇v)(Dδu − ∇δv)(Dp1 − ∇ p2) βhγ (D∇v)D∇δv D∇ p2 + p1δu = 0. By taking successively δv = 0 and δu = 0, the following adjoint system is obtained Using the differentiability of the solution operator and the well-posedness of the adjoint equation, we derive next an optimality system for the characterisation of local minima of the bilevel learning problem. Besides the optimality condition itself, a gradient formula arises as byproduct, which is of importance in the design of solution algorithms for the learning problems. Theorem 3.3 Let (α¯ , β¯) ∈ R2 be a local optimal solution + for problem (2.3). Then there exist Lagrange multipliers ∈ Y := H 1( )×H1( ) and λ1, λ2 ∈ R such that the following system holds a(y, ) + α hγ (Du − w)(Dφ − ϕ) dx = (φ, ϕ) ∈ Y, + β for all λ1 = λ2 = λ1 ≥ 0, λ1 · α¯ = λ2 · β¯ = 0. = (φ, ϕ) ∈ Y, hγ (Du − w)(Dp1 − p2), hγ (E w), E p2 λ2 ≥ 0, a( , ) + α hγ (Du − w)(Dp1 − p2)(Dφ − ϕ) dx hγ (E w) E p2 E ϕ dx + p1φ d x = −Fu (u)[φ], min (α,β)∈C F (α, β), Proof Consider the reduced cost functional F (α, β) = F (u(α, β)). The bilevel optimisation problem can then be formulated as where F : R2 → R and C corresponds to the positive orthant in R2. From [47, Thm. 3.1], there exist multipliers λ1, λ2 ∈ R + β hγ (E w)E ϕ dx + (u − f )φ dx = 0, for all Altogether we proved the result. (3.9a) (3.9b) (3.10c) (3.10d) (3.10e) (3.10f) such that λ1 = ∇αF (α¯ , β¯), λ2 = ∇β F (α¯ , β¯), which, taking into account the linearised equation, yields In this section, we propose a second-order quasi-Newton method for the solution of the learning problem with scalar regularisation parameters. The algorithm is based on a BFGS update, preserving the positivity of the iterates through the line search strategy and updating the matrix cyclically depending on the satisfaction of the curvature condition. For the solution of the lower level problem, a semismooth Newton method with a properly modified Jacobi matrix is considered. Moreover, warm initialisation strategies have to be taken into account in order to get convergence for the TGV2 problem. 4.1 BFGS Algorithm Thanks to the gradient characterisation obtained in Theorem 3.3, we next devise a BFGS algorithm to solve the bilevel learning problems with higher-order regularisers. We employ a few technical tricks to ensure convergence of the classical method. In particular, we limit the step length to get at most a fraction closer to the boundary. As shown in [ 19 ], the solution is in the interior for the regularisation and cost functionals we are interested in. Moreover, the good behaviour of the BFGS method depends upon the BFGS matrix staying positive definite. This would be ensured by the Wolfe conditions, but because of our step length limitation, the curvature condition is not necessarily satisfied. (The Wolfe conditions are guaranteed to be satisfied for some step length σ , if our domain is unbounded, but the range, where the step satisfies the criterion, may be beyond our maximum step length and is not necessarily satisfied closer to the current point.) Instead, we skip the BFGS update if the curvature is negative. Overall, our learning algorithm may be written as follows: Algorithm 4.1 (BFGS for denoising parameter learning) Pick Armijo line search constant c, and target residual ρ. Pick initial iterate (α0, β0). Solve the denoising problem (2.3b) for (α, β) = (α0, β0), yielding u0. Initialise B1 = I . Set i := 0, and iterate the following steps: (1) Solve the adjoint equation (3.10b) for i , and calculate ∇F (αi , βi ) from (3.11). (2) If i ≥ 2, do the following: (a) Set s := (αi , βi )−(αi−1, βi−1), and r := ∇F (αi , βi ) − ∇F (αi−1, βi−1). (b) Perform the BFGS update Bi := Bi−1, Bi−1 − (Bi−t1TsB)(iB−i1−s1s)T rr T + sT r sT r ≤ 0, sT r > 0. (3) Compute δα,β from Bi δα,β = gi . (4) Initialise σ := min{1, σmax/2}, where σmax := max{σ ≥ 0 | (αi , βi ) + σ δα,β > 0}. Repeat the following: (a) Let (ασ , βσ ) := (αi , βi ) + σ δα,β , and solve the denoising problem (2.3b) for (α, β) = (ασ , βσ ), yielding uσ . (b) If the residual (ασ , βσ ) − (αi , βi ) / (ασ , βσ ) < ρ, do the following: μ μ + + N (i) If minσ F (ασ , βσ ) < F (αi , βi ) over all σ tried, choose σ ∗ the minimiser, set (αi+1, βi+1) := (ασ ∗ , βσ ∗ ), ui+1 := uσ ∗ , and continue from Step 5. (ii) Otherwise end the algorithm with solution (α∗, β∗) := (αi , βi ). (c) Otherwise, if Armijo condition F (ασ , βσ ) ≤ F (αi , βi )+σ c∇F (αi , βi )T δα,β holds, set (αi+1, βi+1) := (ασ , βσ ), ui+1 := uσ , and continue from Step 5. (d) In all other cases, set σ := σ/2 and continue from Step 4a. (5) If the residual (αi+1, βi+1)−(αi , βi ) / (αi+1, βi+1) < ρ, end the algorithm with (α∗, β∗) := (αi+1, βi+1). Otherwise continue from Step 1 with i := i + 1. Step (4) ensures that the iterates remain feasible, without making use of a projection step. 4.2 An Infeasible Semismooth Newton Method In this section, we consider semismooth Newton methods for solving the TGV2 and the ICTV denoising problems. Semismooth Newton methods feature a local superlinear convergence rate and have been previously successfully applied to image processing problems (see, e.g. [ 21,29,32 ]). The primal-dual algorithm we use here is an extension of the method proposed in [ 29 ] to the case of higher-order regularisers. In variational form, the TGV2 denoising problem can be written as (Du · Dφ + vφ) + αhγ (Du − w)Dφ (u − f )φ = 0, ∀φ ∈ H 1( ) (E w : E ϕ + wϕ) − αhγ (Du − w)Dϕ βhγ (E w) E ϕ = 0, ∀ϕ ∈ H1( ) or, in general abstract primal-dual form, as i=1 max{1/γ , |[ A j y](x )|2}q j (x ) − α j [ A j y](x ) = 0 a.e. in , j = 1, . . . , N . (4.1b) (4.1a) where L ∈ L(H 1( ; Rm ), H 1( ; Rm ) ) is a second-order linear elliptic operator, A j , j = 1, . . . , N , are linear operators acting on y and q j (x ), j = 1, . . . , N , correspond to the dual multipliers. Let us set m j (y) := max{1/γ , |[ A j y](x )|2}. Let us also define the diagonal application D(y) : L2( ; Rm ) → L2( ; Rm ) by [D(y)q](x ) = y(x )q(x ), (x ∈ ) We may derive ∇y [D(m j (y))q j ] being defined by ∇y [D(m j (y)) p j ] = A∗j D(q j )N( A j y) where N(z) := 0, |zz((xx))|2 , |z(x )|2 ≥ 1/γ . |z(x )|2 < 1/γ Then (4.1a), (4.1b) may be written as N L y + (SSN-1)–(SSN-3) converges globally and locally superlinearly near a point where the subdifferentials of the operator on (y, q1, . . . qN ) corresponding (4.1) are non-singular. Further dampening as in [ 29 ] guarantees local superlinear convergence at any point. We do not present the proof, as going into the discretisation and dampening details would expand this work considerably. Remark 4.1 The system (SSN-1) can be further simplified, which is crucial to obtain acceptable performance with TGV2. Indeed, observe that B is invertible, so we may solve δu from Bδy = R1 − A∗j δq j . N j=1 Thus, we may simplify δy out of (SSN-1) and only solve for δq1, . . . , δqN using a reduced system matrix. Finally, we calculate δy from (4.2). For the denoising sub-problem (2.3b), we use the method (SSN-1)–(SSN-3) with the reduced system matrix of Remark 4.1. Here, we denote by y in the case of TGV2 the parameters (4.2) The semismooth Newton method solves (SSN-1) at a current iterate (yi , q1i, . . . qNi ). It then updates (yi+1, q1i+1, . . . qNi+1) := (yi + τ δy, q1i + τ δq1, qNi + τ δqN ), for a suitable step length τ , allowing qi+1 to become infeasible in the process. That is, it may hold that |qij+1(x )|2 > α j , which may lead to non-descent directions. In order to globalise the method, one projects qij+1 := P(qij+1; α j ), where P(q, α)(x ) (SSN-3) := sgn(q(x )) min{α, |q(x )|}, in the building of the Jacobi matrix. Following [ 29,42 ], it can be shown that a discrete version of the method (SSN-2) For the calculation of the step length τ , we use Armijo line search with parameter c = 1e−4. We end the SSN iterations when τ δY i max{1, Y i } ≤ 1e−5, where δY i := (δyi , δq1i, . . . , δqNi ), and Y i := (yi , q1i, . . . , qNi ). 4.3 Warm Initialisation In our numerical experimentation, we generally found Algorithm 4.1 to perform well for learning the regularisation parameter for TV denoising as was done in [ 22 ]. For learning the two (or even more) regularisation parameters for TGV2 denoising, we found that a warm initialisation is needed to obtain convergence. More specifically, we use TV as an aid for discovering both the initial iterate (α0, β0) as well as the initial BFGS matrix B1. This is outlined in the following algorithm: Algorithm 4.2 (BFGS initialisation for TGV2 parameter learning) Pick a heuristic factor δ0 > 0. Then do the following: Table 1 Quantified results for the parrot image ( = 256 = image width/height in pixels) Initial (α,β) (αT∗V/ , αT∗V) 4c 4d 4e 4f 4g 4h (1) Solve the corresponding problem for TV using Algorithm 4.1. This yields optimal TV denoising parameter αT∗V, as well as the BFGS estimate BTV for ∇2F (αT∗V). (2) Run Algorithm 4.1 for TGV2 with initialisation (α0, β0) := (αT∗Vδ0, αT∗V), and initial BFGS matrix B1 := diag(BTVδ0, BTV). With = (0, 1)2, we pick δ0 = 1/ , where the original discrete image has × pixels. This corresponds to the heuristic [ 2,44 ] that if ≈ 128 or 256, and the discrete image is mapped into the corresponding domain = (0, )2 directly (corresponding to spatial step size of one in the discrete gradient operator), then β ∈ (α, 1.5α) tends to be a good choice. We will later verify this through the use of our algorithms. Now, if f ∈ BV((0, )2) is rescaled to BV((0, 1)2), i.e. f (x ) := f (x / ), then with u(x ) := u(x / ) and w(x ) := w(x / )/ , we have the theoretical equivalence 1 2 f − u 2L2((0, )2) + α Du − w M((0, )2;R2) = n2 +β E w M((0, )2;R2×2) 1 2 f − u 2L2((0,1)2) + nα Du − w +n2β E w M((0,1)2;R2×2) . M((0,1)2;R2) (4.3) (4.4) = | |−1/2 between rescaled This introduces the factor 1/ α, β. 5 Experiments In this section, we present some numerical experiments to verify the theoretical properties of the bilevel learning problems and the efficiency of the proposed solution algorithms. In particular, we exhaustively compare the performance of the new proposed cost functional with respect to well-known quality measures, showing a better behaviour of the new cost for the chosen tested images. The performance of the proposed BFGS algorithm, combined with the semismooth Newton method for the lower level problem, is also examined. Moreover, on basis of the learning setting proposed, a thorough comparison between TGV2 and ICTV is carried out. The use of higher-order regularisers in image denoising is rather recent, and the question on whether TGV2 or ICTV performs better has been around. We target that question and, on basis of the bilevel learning approach, we are able to give some partial answers. 5.1 Gaussian Denoising We tested Algorithm 4.1 for TV and Algorithm 4.2 for TGV2 Gaussian denoising parameter learning on various images. Here we report the results for two images, the parrot image in Fig. 4a, and the geometric image in Fig. 5. We applied synthetic noise to the original images, such that the PSNR of the parrot image are 24.7, and the PSNR of the geometric image is 24.8. In order to learn the regularisation parameter α for TV, we picked initial α0 = 0.1/ . For TGV2, initialisation by TV was used as in Algorithm 4.1. We chose the other parameters of Algorithm 4.1 as c = 1e−4, ρ = 1e−5, θ = 1e−8, and = 10. For the SSN denoising method, the parameters γ = 100 and μ = 1e−10 were chosen. We have included results for both the L2-squared cost functional L22 and the Huberised total variation cost functional L1η∇. The learning results are reported in Table 1 for the parrot images, and Table 2 for the geometric image. The denoising results with the discovered parameters are shown in Figs 4 and 5. We report the resulting optimal parameter values, the cost functional value, PSNR, SSIM [ 46 ], as well as the number of iterations taken by the outer BFGS method. Our first observation is that all approaches successfully learn a denoising parameter that gives a good-quality denoised image. Secondly, we observe that the gradient cost functional L1η∇ performs visually and in terms of SSIM significantly better for TGV2 parameter learning than the cost functional L22. In terms of PSNR, the roles are reversed, as should be, since the L22 is equivalent to PSNR. This again confirms that PSNR is a poor-quality measure for images. For TV, there is no significant difference between different cost functionals in terms of visual quality, although the PSNR and SSIM differ. Table 2 Quantified results for the synthetic image ( = 256 = image width/height in pixels) 0.5 Initial α 5c 5d 5e 5f 5g 5h We also observe that the optimal TGV2 parameters (α∗, β∗) generally satisfy β∗/α∗ ∈ (0.75, 1.5)/ . This confirms the earlier observed heuristic that if ≈ 128, 256 then β ∈ (1, 1.5)α tends to be a good choice. As we can observe from Figs. 4 and 5, this optimal TGV2 parameter choice also avoids the staircasing effect that can be observed with TV in the results. In Fig. 3, we have plotted by the red star the discovered regularisation parameter (α∗, β∗) reported in Fig. 4. Studying the location of the red star, we may conclude that Algorithms 4.1 and 4.2 manage to find a nearly optimal parameter in very few BFGS iterations. 5.2 Statistical Testing To obtain a statistically significant outlook to the performance of different regularisers and cost functionals, we made use of the Berkeley segmentation dataset BSDS300 [ 36 ], displayed in Fig. 6. We resized each image to 128 pixels on its shortest edge and take the 128 × 128 top left square of the image. To this dataset, we applied pixelwise Gaussian noise of variance σ 2 = 2, 10, and 20. We tested the performance of both cost functionals, L 1η∇ and L 22, as well as the TGV2, ICTV, and TV regularisers on this dataset, for all noise levels. In the first instance, reported in Figs. 7, 8, 9 and 10 (noise levels σ 2 = 2, 20 only), and Tables 3, 4 and 5, we applied the proposed bilevel learning model on each image individually, to learn the optimal parameters specifically for that image, and a corresponding noisy image for all of the noise levels separately. For the algorithm, we use the same parametrisation as presented in Sect. 5.1. The figures display the noisy images and indicate by colour coding the best result as judged by the structural similarity measure SSIM [ 46 ], PSNR and the objective function value (L 1η∇ or L 22 cost). These criteria are, respectively, the top, middle and bottom rows of colour-coding squares. Red square indicates that TV performed the best, green square indicates that ICTV performed the best and blue square indicates that TGV2 performed the best—this is naturally for the optimal parameters for the corresponding regulariser and cost functional discovered by our algorithms. In the tables, we report the information in a more concise numerical fashion, indicating the mean, standard deviation and median for all the different criteria (SSIM, PSNR and cost functional value), as well as the number of images for which each regulariser performed the best. We recall that SSIM is normalised to [ 0, 1 ], with higher value better. Moreover, we perform a statistical 95 paired t-test on each of the criteria, and (a) Original image (b) Noisy image (a) Original image (b) Noisy image (c) TGV2 denoising, L1η∇ cost (d) TGV2 denoising, L22 cost (c) TGV2 denoising, L1η∇ cost (d) TGV2 denoising, L22 cost (e) ICTV denoising, L1η∇ cost (f ) ICTV denoising, L22 cost (e) ICTV denoising, L1η∇ cost (f ) ICTV denoising, L22 cost (g) TV denoising, L1η∇ cost (h) TV denoising, L22 cost Fig. 4 Optimal denoising results for initial guess α = (αT∗V/ , αT∗V) for TGV2 and α = 0.1/ for TV a pair of regularisers, to see whether any pair of regularisers can be ordered. If so, this is indicated in the last row of each of the tables. Overall, studying the t-test and other data, the ordering of the regularisers appears to be ICTV > TGV2 > TV. (g) TV denoising, L1η∇ cost (h) TV denoising, L22 cost Fig. 5 Optimal denoising results for initial guess α = (αT∗V/ , αT∗V) for TGV2 and α = 0.2/ for TV This is rather surprising, as in many specific examples, TGV2 has been observed to perform better than ICTV, see Figs. 4 and 5, as well as [ 1, 5 ]. Only when the noise is high, appears TGV2 to come on par with ICTV with the L 1η∇ cost functional in Fig. 9 and Table 5. A more detailed study of the results in Figs. 7, 8, 9 and 10 seems to indicate that TGV2 performs better than ICTV when Fig. 7 Ordering of regularisers with individual learning, L1η∇ cost, and noise variance σ 2 = 2, on the 200 images of the BSDS300 dataset, resized. Best regulariser: red TV, green ICTV, blue TGV2; top SSIM, middle PSNR, bottom objective value Fig. 8 Ordering of regularisers with individual learning, L22 cost, and noise variance σ 2 = 2, on the 200 images of the BSDS300 dataset, resized. Best regulariser: red TV, green ICTV, blue TGV2; top SSIM, middle PSNR, bottom objective value the image contains large smooth areas, but ICTV generally seems to perform better for images with more complicated and varying contents. This observation agrees with the results in Figs. 4 and 5, as well as [ 1, 5 ], where the images are of the former type. One possible reason for the better performance of ICTV could be that TGV2 has more degrees of freedom—in ICTV we essentially constrain w = ∇v for some function v— and therefore overfits to the noisy data, until the noise level becomes so high that overfitting would become too high for Fig. 10 Ordering of regularisers with individual learning, L22 cost, and noise variance σ 2 = 20, on the 200 images of the BSDS300 dataset, resized. Best regulariser: red TV, green ICTV, blue TGV2; top SSIM, middle PSNR, bottom objective value any parameter. To see whether this is true, we also performed batch learning, learning a single set of parameters for all images with the same noise level. That is, we studied the model N min α i=1 + Rγ ,μ(u), α 1 Fi (ui,α ) s.t. ui,α ∈ ua∈rgH m1(in) 2 fi − u 2L2( ) Best Std where α = (α, β), f1, . . . , f N are the N = 200 noisy images with the same noise level, and f0,1, . . . , f0,N the original noise-free images. The results are shown in Figs. 11, 12, 13 and 14 (noise levels σ 2 = 2, 20 only), and Tables 6, 7 and 8. The results are still roughly the same as with individual learning. Again, only with high noise in Table 8, TGV2 does not lose to ICTV. Another interesting observation is that TV starts to be frequently the best regulariser for individual images, although still statistically does worse than either ICTV or TGV2. For the first image of the dataset, ICTV does in all of the Figs. 7, 8, 9, 10, 11, 12, 13 and 14 better than TGV2, while for the second image, the situation is reversed. We have highlighted these two images for the L 1η∇ cost in Figs. 15, 16, 17 and 18, for both noise levels σ = 2 and σ = 20. In the case where ICTV does better, hardly any difference can be observed by the eye, while for second image, TGV2 clearly has less staircasing in the smooth areas of the image, especially with the noise level σ = 20. Based on this study, it therefore seems that ICTV is the most reliable regulariser of the ones tested, when the type of image being processed is unknown, and low SSIM, PSNR or L 1η∇ cost functional value is desired. But as can be observed for individual images, it can within large smooth areas exhibit artefacts that are avoided by the use of TGV2. 5.3 The Choice of Cost Functional The L 22 cost functional naturally obtains better PSNR than L 1η∇, as the two former are equivalent. Comparing the results for the two cost funtionals in Tables 3, 4 and 5, we may however observe that for low noise levels σ 2 = 2, 10, and generally for batch learning, L 1η∇ attains better (higher) SSIM. Since SSIM better captures [ 46 ] the visual quality of Fig. 11 Ordering of regularisers with batch learning, L1η∇ cost, and noise variance σ 2 = 2, on the 200 images of the BSDS300 dataset, resized. Best regulariser: red TV, green ICTV, blue TGV2; top SSIM, middle PSNR, bottom objective value Fig. 12 Ordering of regularisers with batch learning, L22 cost, and noise variance σ 2 = 2, on the 200 images of the BSDS300 dataset, resized. Best regulariser: red TV, green ICTV, blue TGV2; top SSIM, middle PSNR, bottom objective value Fig. 13 Ordering of regularisers with batch learning, L1η∇ cost, and noise variance σ 2 = 20, on the 200 images of the BSDS300 dataset, resized. Best regulariser: red TV, green ICTV, blue TGV2; top SSIM, middle PSNR, bottom objective value images than PSNR, this recommends the use of our novel total variation cost functional L 1η∇. Of course, one might attempt to optimise the SSIM. This is however a non-convex functional, which will pose additional numerical challenges avoided by the convex total variation cost. 6 Conclusion and Outlook In this paper, we propose a bilevel optimisation approach in function space for learning the optimal choice of parameters in higher-order total variation regularisation. We present a Fig. 14 Ordering of regularisers with batch learning, L22, cost, and noise variance σ 2 = 20, on the 200 images of the BSDS300 dataset, resized. Best regulariser: red TV, green ICTV, blue TGV2; top SSIM, middle PSNR, bottom objective value rigorous analysis of this optimisation problem as well as a numerical discussion in the context of image denoising. Analytically, we obtain the existence results for the bilevel optimisation problem and prove the Fréchet differentiability of the solution operator. This leads to the existence of Lagrange multipliers and a first-order optimality system characterising optimal solutions. In particular, the existence of an adjoint state allows to obtain a cost functional gradiMed Fig. 15 Image for which ICTV performs better than TGV2, σ = 2 ent formula which is of importance in the design of efficient solution algorithms. We make use of the bilevel learning approach, and the theoretical findings, to compare the performance—in terms of returned image quality—of TV, ICTV and TGV regularisation. A statistical analysis, carried out on a dataset of 200 images, suggests that ICTV performs slightly better than TGV, and both perform better than TV, in average. For denoising of images with a high noise level, ICTV and TGV score comparably well. For images with large smooth areas, TGV performs better than ICTV. Moreover, we propose a new cost functional for the bilevel learning problem, which exhibits interesting theoretical properties and has a better behaviour with respect to the PSNR related L2 cost used previously in the literature. This study raises the question of other, alternative cost functionals. For instance, one could be tempted to used the SSIM as cost, but its non-convexity might present several analytical and numerical difficulties. The new cost functional, proposed in this paper, turns out to be a good compromise between image quality measure and analytically tractable cost term. Fig. 16 Image for which ICTV performs better than TGV2, σ = 20 Fig. 17 Image for which TGV2 performs better than ICTV, σ = 2 (a) Original (b) Individual L1η∇-TGV2, PSNR=28.28, SSIM=0.74 (c) Batch L1η∇-TGV2, PSNR=28.25, SSIM=0.73 Fig. 18 Image for which TGV2 performs better than ICTV, σ = 20 Acknowledgments This research has been supported by King Abdullah University of Science and Technology (KAUST) Award No. KUKI1-007-43, EPSRC grants Nr. EP/J009539/1 “Sparse & Higher-order Image Restoration” and Nr. EP/M00483X/1 “Efficient computational tools for inverse imaging problems”, Escuela Politécnica Nacional de Quito Award No. PIS 12-14, MATHAmSud project SOCDE “Sparse Optimal Control of Differential Equations” and the Leverhulme Trust project on “Breaking the non-convexity barrier”. While in Quito, T. Valkonen has moreover been supported by SENESCYT (Ecuadorian Ministry of Higher Education, Science, Technology and Innovation) under a Prometeo Fellowship. 1. Benning , M. , Brune , C. , Burger , M. , Müller , J.: Higher-order TV methods-enhancement via Bregman iteration . J. Sci. Comput . 54 ( 2-3 ), 269 - 310 ( 2013 ) 2. Benning , M. , Gladden , L. , Holland , D. , Schönlieb , C.-B., Valkonen , T. : Phase reconstruction from velocity-encoded MRI measurements-a survey of sparsity-promoting variational approaches . J. Magn. Reson . 238 , 26 - 43 ( 2014 ) 3. Biegler , L. , Biros , G. , Ghattas , O. , Heinkenschloss , M. , Keyes , D. , Mallick , B. , Tenorio , L., van Bloemen Waanders, B. , Willcox , K. , Marzouk , Y. : Large-Scale Inverse Problems and Quantification of Uncertainty, vol. 712 . Wiley, New York ( 2011 ) 4. Bonnans , J.F. , Tiba , D. : Pontryagin's principle in the control of semilinear elliptic variational inequalities . Appl. Math. Optim . 23 ( 1 ), 299 - 312 ( 1991 ) 5. Bredies , K. , Kunisch , K. , Pock , T. : Total generalized variation . SIAM J. Imaging Sci. 3 , 492 - 526 ( 2011 ) 6. Bredies , K. , Holler , M.: A total variation-based jpeg decompression model . SIAM J. Imaging Sci . 5 ( 1 ), 366 - 393 ( 2012 ) 7. Bredies , K. , Kunisch , K. , Valkonen , T. : Properties of L1 − TGV2: the one-dimensional case . J. Math. Anal. Appl . 398 , 438 - 454 ( 2013 ) 8. Bredies , K. , Valkonen , T. : Inverse problems with second-order total generalized variation constraints . In: Proceedings of the 9th International Conference on Sampling Theory and Applications (SampTA), Singapore ( 2011 ) 9. Bui-Thanh , T. , Willcox , K. , Ghattas , O. : Model reduction for largescale systems with high-dimensional parametric input space . SIAM J. Sci. Comput . 30 ( 6 ), 3270 - 3288 ( 2008 ) 10. Calatroni , L. , De los Reyes, J.C. , Schönlieb , C.-B.: Dynamic sampling schemes for optimal noise learning under multiple nonsmooth constraints . In: Poetzsche, C . (ed.) System Modeling and Optimization , pp. 85 - 95 . Springer Verlag, New York ( 2014 ) 11. Chambolle , A. , Lions , P.-L.: Image recovery via total variation minimization and related problems . Numer. Math. 76 , 167 - 188 ( 1997 ) 12. Chan , T. , Marquina , A. , Mulet , P. : High-order total variation-based image restoration . SIAM J. Sci. Comput . 22 ( 2 ), 503 - 516 ( 2000 ) 13. Chan , T.F. , Kang , S.H. , Shen , J. : Euler's elastica and curvaturebased inpainting . SIAM J. Appl. Math . 63 ( 2 ), 564 - 592 ( 2002 ) 14. Chen , Y. , Pock , T. , Bischof , H.: Learning 1-based analysis and synthesis sparsity priors using bi-level optimization . In: Workshop on Analysis Operator Learning versus Dictionary Learning , NIPS 2012 ( 2012 ) 15. Chen , Y. , Ranftl , R. , Pock , T. : Insights into analysis operator learning: from patch-based sparse models to higher-order mrfs . IEEE Trans. Image Process . ( 2014 ) (to appear) 16. Chung , J. , Español , M.I. , Nguyen , T. : Optimal regularization parameters for general-form tikhonov regularization . arXiv preprint arXiv:1407 . 1911 ( 2014 ) 17. Dauge , M. : Neumann and mixed problems on curvilinear polyhedra . Integr. Equ. Oper. Theory 15 ( 2 ), 227 - 261 ( 1992 ) 18. De los Reyes, J.C., Meyer, C.: Strong stationarity conditions for a class of optimization problems governed by variational inequalities of the second kind . J. Optim. Theory Appl . 168 ( 2 ), 375 - 409 ( 2015 ) 19. De los Reyes, J.C. , Schönlieb , C.-B., Valkonen , T. : The structure of optimal parameters for image restoration problems . J. Math. Anal. Appl . 434 ( 1 ), 464 - 500 ( 2016 ) 20. De los Reyes, J.C.: Optimal control of a class of variational inequalities of the second kind . SIAM J. Control Optim . 49 ( 4 ), 1629 - 1658 ( 2011 ) 21. De los Reyes, J.C. , Hintermüller , M.: A duality based semismooth Newton framework for solving variational inequalities of the second kind . Interfaces Free Bound . 13 ( 4 ), 437 - 462 ( 2011 ) 22. De los Reyes, J.C. , Schönlieb , C.-B.: Image denoising: learning the noise model via nonsmooth PDE-constrained optimization . Inverse Probl. Imaging 7 ( 4 ), 1139 - 1155 ( 2013 ) 23. Domke , J.: Generic methods for optimization-based modeling . In: International Conference on Artificial Intelligence and Statistics , pp. 318 - 326 ( 2012 ) 24. Gröger , K. : A W 1,p-estimate for solutions to mixed boundary value problems for second order elliptic differential equations . Math. Ann. 283 ( 4 ), 679 - 687 ( 1989 ) 25. Haber , E. , Tenorio , L. : Learning regularization functionals-a supervised training approach . Inverse Probl . 19 ( 3 ), 611 ( 2003 ) 26. Haber , E. , Horesh , L. , Tenorio , L. : Numerical methods for the design of large-scale nonlinear discrete ill-posed inverse problems . Inverse Probl . 26 ( 2 ), 025002 ( 2010 ) 27. Hinterberger , W. , Scherzer , O. : Variational methods on the space of functions of bounded hessian for convexification and denoising . Computing 76 ( 1 ), 109 - 133 ( 2006 ) 28. Hintermüller , M. , Laurain , A. , Löbhard , C. , Rautenberg , C.N. , Surowiec , T.M.: Elliptic mathematical programs with equilibrium constraints in function space: Optimality conditions and numerical realization . In: Rannacher, R . (ed.) Trends in PDE Constrained Optimization , pp. 133 - 153 . Springer International Publishing, Berlin ( 2014 ) 29. Hintermüller , M. , Stadler , G. : An infeasible primal-dual algorithm for total bounded variation-based inf-convolution-type image restoration . SIAM J. Sci. Comput . 28 ( 1 ), 1 - 23 ( 2006 ) 30. Hintermüller , M. , Wu , T. : Bilevel optimization for calibrating point spread functions in blind deconvolution . Preprint ( 2014 ) 31. Knoll , F. , Bredies , K. , Pock , T. , Stollberger , R.: Second order total generalized variation (TGV) for MRI . Magn. Reson. Med . 65 ( 2 ), 480 - 491 ( 2011 ) 32. Kunisch , K. , Hintermüller , M. : Total bounded variation regularization as a bilaterally constrained optimization problem . SIAM J. Imaging Sci . 64 ( 4 ), 1311 - 1333 ( 2004 ) 33. Kunisch , K. , Pock , T. : A bilevel optimization approach for parameter learning in variational models . SIAM J. Imaging Sci . 6 ( 2 ), 938 - 983 ( 2013 ) 34. Luo , Z.-Q. , Pang , J.-S. , Ralph , D. : Mathematical Programs with Equilibrium Constraints . Cambridge University Press, Cambridge ( 1996 ) 35. Lysaker , M. , Tai , X.-C.: Iterative image restoration combining total variation minimization and a second-order functional . Int. J. Comput. Vis . 66 ( 1 ), 5 - 18 ( 2006 ) 36. Martin , D. , Fowlkes , C. , Tal , D. , Malik , J.: A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics . In: Proceedings of the 8th International Conference on Computer Vision , vol. 2 , pp. 416 - 423 ( 2001 ). The database is available online at CS/vision/bsds/BSDS300/html/dataset/images.html 37. Masnou , S. , Morel , J.-M.: Level lines based disocclusion . In: 1998 IEEE International Conference on Image Processing (ICIP 98) , pp. 259 - 263 ( 1998 ) 38. Outrata , J.V. : A generalized mathematical program with equilibrium constraints . SIAM J. Control Optim . 38 ( 5 ), 1623 - 1638 ( 2000 ) 39. Papafitsoros , K. , Schönlieb , C.-B.: A combined first and second order variational approach for image reconstruction . J. Math. Imaging Vis . 48 ( 2 ), 308 - 338 ( 2014 ) 40. Ring , W. : Structural properties of solutions to total variation regularization problems . ESAIM 34 , 799 - 810 ( 2000 ) 41. Rudin , L. , Osher , S. , Fatemi , E.: Nonlinear total variation based noise removal algorithms . Phys. D 60 , 259 - 268 ( 1992 ) 42. Sun , D., Han, J .: Newton and quasi-Newton methods for a class of nonsmooth equations and related problems . SIAM J. Optim . 7 ( 2 ), 463 - 480 ( 1997 ) 43. Tappen , M.F. : Utilizing variational optimization to learn Markov random fields . In: 2007 IEEE Conference on Computer Vision and Pattern Recognition (CVPR'07) , pp. 1 - 8 ( 2007 ) 44. Valkonen , T. , Bredies , K. , Knoll , F. : Total generalised variation in diffusion tensor imaging . SIAM J. Imaging Sci . 6 ( 1 ), 487 - 525 ( 2013 ) 45. Viola , F. , Fitzgibbon , A. , Cipolla , R.: A unifying resolutionindependent formulation for early vision . In: 2012 IEEE Conference on Computer Vision and Pattern Recognition (CVPR) , pp. 494 - 501 ( 2012 ) 46. Wang , Z. , Bovik , A.C. , Sheikh , H.R. , Simoncelli , E.P. : Image quality assessment: from error visibility to structural similarity . IEEE Trans. Image Process . 13 ( 4 ), 600 - 612 ( 2004 )

This is a preview of a remote PDF:

J. C. De los Reyes, C.-B. Schönlieb, T. Valkonen. Bilevel Parameter Learning for Higher-Order Total Variation Regularisation Models, Journal of Mathematical Imaging and Vision, 2017, 1-25, DOI: 10.1007/s10851-016-0662-8