Luca Calatroni
email: luca.calatroni@polytechnique.edu

Kostas Papafitsoros
email: kostas.papafitsoros@wias-berlin.de

Analysis and optimisation of a variational model for mixed Gaussian and Salt & Pepper noise removal

. The data discrepancy term of the model combines L 1 and L 2 terms in an infimal convolution fashion and it is appropriate for the joint removal of Gaussian and Salt & Pepper noise. In this work we perform a finer analysis of the model which emphasises on the balancing effect of the two parameters appearing in the discrepancy term. Namely, we study the asymptotic behaviour of the model for large and small values of these parameters and we compare it to the corresponding variational models with L 1 and L 2 data fidelity. Furthermore, we compute exact solutions for simple data functions taking the total variation as regulariser. Using these theoretical results, we then analytically study a bilevel optimisation strategy for automatically selecting the parameters of the model by means of a training set. Finally, we report some numerical results on the selection of the optimal noise model via such strategy which confirm the validity of our analysis and the use of popular data models in the case of "blind" model selection.

Introduction

Image denoising is a classical problem in imaging which is defined as the task of removing oscillations and interferences from a given image f . Given a regular image domain Ω ⊆ R 2 the problem can be formulated mathematically as the task of retrieving a noise-free image u from its noisy measurement f , the latter being the result of a (possibly non-linear) degradation operator T . In its general form, this problem can be written in the following way find u such that f = T (u).

Here the operator T introduces noise in the image, not only in an additive way, and it is not to be confused with the forward operator in the context of inverse problems. In order to obtain a noise-free image, a classical technique consists of minimising an appropriate energy functional J over a suitable Banach space X where the image functions are assumed to lie. In its general form the problem reads as follows: min u∈X {J (u) := R(u) + λΦ(u, f )} .

(1.1)

Here, R(u) stands for the regularisation term encoding a priori information on the regularity of the solution, Φ(u, f ) for the data-fitting measure that depends on the statistical and physical assumptions in the data and λ > 0 is a scalar parameter whose magnitude balances the regularisation against trust in the data. Since the seminal work of Rudin, Osher, Fatemi [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], a popular choice for R in (1.1) is R(u) = |Du|(Ω), the Total Variation (TV) seminorm [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF], due to its ability of preserving salient structures in the image, i.e., edges, while removing noise at the same time. Here Du represents the distributional derivative of the function u ∈ BV(Ω), the space of functions of bounded variation, and |Du|(Ω) is the total variation of this measure. In the recent years, also higher-order regularisation terms that improve upon TV-induced artefacts -notably, the creation of piecewise constant structures -have been proposed in the literature. Among those, the most prominent is the Total Generalized Variation (TGV) [START_REF] Bredies | Total Generalized Variation[END_REF], see also [START_REF] Papafitsoros | Novel higher order regularisation methods for image reconstruction[END_REF] for a comprehensive review.

In this work we will mostly focus on the standard TV regularisation energy, as our work focuses on the choice of Φ rather than of R. Classical data fidelity terms for denoising images with Gaussian or impulsive Salt & Pepper noise are based on the use of the L 2 and L 1 norm, respectively, i.e.,

Φ L 2 (u, f ) := 1 2 ∫ Ω (f -u) 2 dx or Φ L 1 (u, f ) := ∫ Ω |f -u| dx.
These discrepancies are statistically consistent with the assumptions on the noise since they can be derived as the MAP estimators of the underlying likelihood function [START_REF] Benning | Error estimates for general fidelities[END_REF].

There exists a considerable amount of work in the literature regarding the structures of solutions of variational problems with pure L 1 or L 2 fidelity terms. As a result, the differences between the effects that these terms have on the solutions of the corresponding denoising models are well understood. See for instance [START_REF] Allard | Total Variation Regularization for Image Denoising, I. Geometric Theory[END_REF][START_REF] Allard | Total Variation Regularization for Image Denoising, II. Examples[END_REF][START_REF] Allard | Total Variation Regularization for Image Denoising, III. Examples[END_REF][START_REF] Caselles | The discontinuity set of solutions of the TV denoising problem and some extensions[END_REF][START_REF] Chambolle | Geometric properties of solutions to the total variation denoising problem[END_REF][START_REF] Chan | Aspects of Total Variation regularized L 1 function approximation[END_REF][START_REF] Cristoferi | Exact solutions for the denoising problem of piecewise constant images in dimension one[END_REF][START_REF] Duval | The TV-L1 model: a geometric point of view[END_REF][START_REF] Jalalzai | Some Remarks on the Staircasing Phenomenon in Total Variation-Based Image Denoising[END_REF][START_REF] Meyer | Oscillating patterns in image processing and nonlinear evolution equations: the fifteenth Dean Jacqueline B[END_REF][START_REF] Ring | Structural properties of solutions to total variation regularization problems[END_REF][START_REF] Valkonen | The jump set under geometric regularization. Part 1: Basic technique and first-order denoising[END_REF] for TV regularisation and [START_REF] Bredies | Properties of L 1 -TGV 2 : The one-dimensional case[END_REF][START_REF] Papafitsoros | A study of the one dimensional total generalised variation regularisation problem[END_REF][START_REF] Papafitsoros | Asymptotic behaviour of total generalised variation[END_REF][START_REF] Pöschl | Exact solutions of one-dimensional total generalized variation[END_REF][START_REF] Valkonen | The jump set under geometric regularisation. Part 2: Higher-order approaches[END_REF] for TGV. For example, in the case of TV regularisation, it is known that the use of L 2 fidelity does not introduce new discontinuities in the solution, which is not the case for the L 1 fidelity. Moreover, the L 1 model is capable of exact data recovery, in contrast to the L 2 one, where always some loss of contrast occurs.

Image denoising for noise mixtures

Due to different image acquisition and transmission faults, the given image f may be corrupted by a mixture of noise statistics. This is typical, for instance, whenever the presence of noise is due to electronic faults and/or photon-counting processes (such as in astronomy and microscopy applications) combined with actual damages in the receiving sensors, resulting in a lack of information transmittance in only a few image pixels ("burned pixels"). The modelling of Φ in (1.1) is therefore expected to encode such noise combination. In this work we consider the special case when Gaussian and sparse Salt & Pepper noise are present in the data.

Several authors have considered previously such noise combination. In [START_REF] Hintermüller | Subspace correction methods for a class of nonsmooth and nonadditive convex variational problems with mixed L 1 /L 2 data-fidelity in image processing[END_REF][START_REF] Langer | Automated parameter selection in the L 1 -L 2 -TV model for removing gaussian plus impulse noise[END_REF], for instance, a combined model with a L 1 + L 2 data fidelity and TV regularisation is considered for the joint removal of impulsive and Gaussian noise. Another approach is considered in [START_REF] Cai | Two-phase approach for deblurring images corrupted by impulse plus Gaussian noise[END_REF] where two sequential steps having L 1 and L 2 as data fidelity are performed to remove the impulsive and the Gaussian component of the mixed noise, respectively. Framelet-based approaches combining L 1 and L 2 data fidelities in a discrete setting have also been proposed in [START_REF] Dong | Wavelet-based blind image inpainting[END_REF][START_REF] Yan | Restoration of images corrupted by impulse noise and mixed gaussian impulse noise using blind inpainting[END_REF]. However, despite the observed good practical performance of the models described above, the use of an additive and sequential combination of L 1 and L 2 data fitting terms lacks a rigorous statistical interpretation in terms, for instance, of a MAP estimation.

We remark that the combination of other noise distributions such as, for instance, Gaussian and Poisson is also frequent in astronomy and microscopy applications and have been studied in several works such as, for instance, [START_REF] Benvenuto | The study of an iterative method for the reconstruction of images corrupted by Poisson and Gaussian noise[END_REF][START_REF] Jezierska | A Convex Approach for Image Restoration with Exact Poisson-Gaussian Likelihood[END_REF][START_REF] Lanza | Image restoration with Poisson-Gaussian mixed noise[END_REF]. Such noise mixtures, however, are outside the scope of this work.

Infimal convolution modelling of data discrepancies

Recently in [START_REF] Calatroni | Infimal convolution of data discrepancies for mixed noise removal[END_REF], a non-standard variational model for noise removal of mixtures of Salt & Pepper and Gaussian, and Gaussian and Poisson noise has been studied. The model, which will be referred to as TV-IC model, is based on the minimisation of an energy functional which is the sum of |Du|(Ω) and an infimal convolution of single noise data discrepancy terms. Given two positive parameters λ 1 , λ 2 it reads:

min u∈BV(Ω) |Du|(Ω) + Φ λ 1 ,λ 2 (u, f ),
where the data fidelity Φ λ 1 ,λ 2 (u, f ) is defined as

Φ λ 1 ,λ 2 (u, f ) := inf v λ 1 Φ 1 (v) + λ 2 Φ 2 (v, f -u).
(1.2)

Here, Φ 1 , Φ 2 denote standard data fidelity terms typically used for single noise removal such as the L 1 , L 2 norm and the Kullback-Leibler functional. In the particular case of a mixture of Salt & Pepper and Gaussian noise, (1.2) specifies into

Φ 1 (v) = ∥v∥ L 1 (Ω) and Φ 2 (v, f -u) = 1 2 ∥f -u -v∥ 2 L 2 (Ω) . In this case, the minimisation in (1.2) is done over L 1 (Ω) and Φ λ 1 ,λ 2 (u, f ) reads Φ λ 1 ,λ 2 (u, f ) = min v∈L 1 (Ω) λ 1 ∥v∥ L 1 (Ω) + λ 2 2 ∥f -u -v∥ 2 L 2 (Ω) , f, u ∈ L 1 (Ω). (1.3)
It can be easily checked that the minimisation in (1.3) is indeed well-defined. The reader should not be alerted by the fact that the L 1 functions f, u appear in the L 2 part of (1.3), as the variable v takes care of any non-integrability issue, see Proposition 2.1. In fact, as we are going to remark in the same proposition, the functional Φ λ 1 ,λ 2 can be written equivalently as

Φ λ 1 ,λ 2 (u, f ) = ∫ Ω φ(f (x) -u(x)) dx,
where φ is the well-known Huber-regularisation of the absolute value function. As a consequence, the functional Φ λ 1 ,λ 2 can be simply seen as a Huberised L 1 norm.

In [START_REF] Calatroni | Infimal convolution of data discrepancies for mixed noise removal[END_REF], it was shown that the data discrepancy (1.3) corresponds to the joint MAP estimator for a denoising problem featuring a mixture of Laplace and Gaussian noise distributions. The effectiveness of the model for the removal of such noise mixture as well as the additional property of decomposing the noise into its sparse (Salt & Pepper) and distributed (Gaussian) component was there confirmed with extended numerical examples.

Note that a Huber smoothing of the L 1 norm has previously been considered in order to apply fast second order minimisation algorithms such as semismooth Newton method in [START_REF] Hintermüller | Expected absolute value estimators for a spatially adapted regularization parameter choice rule in L 1 -TV-based image restoration[END_REF]. Also, in the purely discrete setting, smoothed TV-L 1 models have been studied in [START_REF] Nikolova | Exact histogram specification for digital images using a variational approach[END_REF] for exact histogram specification. Similar models (among which also Huber-type) were also considered in [START_REF] Baus | Fully smoothed ℓ 1 -T V models: Bounds for the minimizers and parameter choice[END_REF], where the authors obtained bounds on the infinity norm of the difference between data and solutions.

Our contribution. In this work we examine in depth the similarities and the differences between the TV-IC model and the pure TV-L 1 , TV-L 2 ones. We first provide detailed asymptotic results as λ 1 or λ 2 tend either to infinity or to zero and describe how the solution of the model varies in these cases. Note that these results are proved for a general regularisation term. Secondly, in order to have a better insight on the type of solutions one can expect, we do a fine scale analysis of the one-dimensional TV regularised model by computing exact solutions for simple data functions f . Up to our knowledge, this is the first time that the effect of the Huberised L 1 fidelity term is studied in the continuous setting.

In the second part of the paper, we focus on the optimal selection of the parameters λ 1 , λ 2 appearing in (1.3). In order to do that, we consider a bilevel optimisation approach [START_REF] Calatroni | Bilevel approaches for learning of variational imaging models[END_REF][START_REF] De Los Reyes | Image denoising: learning the noise model via nonsmooth PDE-constrained optimization[END_REF] which in its general formulation reads

min λ 1 ,λ 2 ≥0 F (u λ 1 ,λ 2 ) subject to u λ 1 ,λ 2 ∈ argmin u∈BV(Ω) |Du|(Ω) + Φ λ 1 ,λ 2 (u, f ).
(1.4)

Here F denotes a cost functional which measures how far the solution u λ 1 ,λ 2 is from some ground truth (training) image. The parameters λ 1 , λ 2 selected within this framework are therefore those producing the closest reconstruction u λ 1 ,λ 2 to the training image, see Section 4.1 for more details. We perform a rigorous analysis on the existence of solutions of (1.4) as well as a regularised version of it, proving the Fréchet differentiability of the solution map S : (λ 1 , λ 2 ) → u λ 1 ,λ 2 and the existence of adjoint states. This allows to derive a handy characterisation of the gradient of the reduced form of F which can be used for efficient numerical implementations. Our analysis justifies rigorously the formal Lagrangian approach considered in [14, Section 7].

We conclude our study with some numerical experiments connecting the analysis on the structure of solutions discussed above with the problem of learning the optimal noise model for a given noisy image with unknown noise intensity. Our numerical findings show that in case of pure Salt & Pepper and Gaussian denoising, the bilevel optimisation approach applied to the TV-IC model computes optimal parameters which enforces pure L 1 fidelity. In the case of noise mixture, a combination of L 1 and L 2 data fitting is preferred. Interestingly, in the case of pure Gaussian, it is not the pure L 2 data fitting that is selected but still a combination of L 1 and L 2 , indicating the benefit of the use of L 1 discrepancy even in the case of Gaussian noise.

We emphasise that the two parts of the paper are intrinsically connected, see for instance Propositions 4.2 and 4.3. In the former, by making use of the analytical results of the first part, we show with a help of a counterexample, that in order to show existence of solutions for the bilevel optimisation problem (1.4), it is necessary to enforce an upper bound on the parameters λ 1 , λ 2 . The existence of solutions in this case, is shown in Proposition 4.3 also by making use of the results of the first part of the paper.

Overall, this study motivates further the use of the TV-IC model and in general the use of the infimal convolution based fidelity term, by (i) describing the structure of the expected solutions and (ii) by proposing an automated parameter selection strategy making this model more flexible and applicable to mixed denoising problems.

Analysis of the L 1 -L 2 IC model: characterisation and asymptotics

We start this section by observing that the L 1 -L 2 infimal convolution term can be equivalently formulated as a Huberised L 1 norm. This provides an interesting motivation on its effectiveness in the removal of mixed Salt & Pepper and Gaussian noise. As usual, Ω ⊆ R d denotes an open, bounded, connected domain with Lipschitz boundary.

Proposition 2.1. Let λ 1 , λ 2 > 0, f, u ∈ L 1 (Ω) and consider the data fitting term, Φ λ 1 ,λ 2 (u, f ) defined in (1.3). Then Φ λ 1 ,λ 2 (u, f ) = ∫ Ω φ(f (x) -u(x)) dx, (2.1)
where for t ∈ R

φ(t) = { λ 1 |t| - λ 2 1 2λ 2 , if |t| ≥ λ 1 λ 2 , λ 2 2 |t| 2 , if |t| < λ 1 λ 2 .
(2.2)

Proof. The proof is straightforward, in view of

Φ λ 1 ,λ 2 (u, f ) = min v∈L 1 (Ω) ∫ Ω ( λ 1 |v(x)| + λ 2 2 (f (x) -u(x) -v(x)) 2 ) dx,
and noticing that the minimisation in the definition of Φ λ 1 ,λ 2 can be considered pointwise. Immediate calculations show that the optimal v can be computed explicitly as

v opt (x) = { f (x) -u(x) -λ 1 λ 2 f (x)-u(x) |f (x)-u(x)| , if |f (x) -u(x)| ≥ λ 1 λ 2 , 0, if |f (x) -u(x)| < λ 1 λ 2 .
(2.3)

From the formulation (2.1)-(2.2) one sees that the infimal convolution of L 1 and L 2 data fidelities coincides with a smoothed L 1 norm, see Figure 1. This is in fact well-known in the context of optimisation in Hilbert spaces [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF], where an explicit expression like (2.3) is often used to compute the soft-thresholding operators. Furthermore, let us point out here that Proposition 2.1 above is analogous to a similar result about the Huberised total variation functional, which has also been shown to be equivalent to a corresponding infimal convolution functional involving L 1 and L 2 norms, see [START_REF] Burger | Infimal convolution regularisation functionals of BV and L p spaces. Part I: The finite p case[END_REF] and [START_REF] Hintermüler | An Infeasible Primal-Dual Algorithm for Total Bounded Variation-Based Inf-Convolution-Type Image Restoration[END_REF]. The formulation (2.1)-(2.2) of Φ λ 1 ,λ 2 provides an interesting insight on its interpretation and motivates its effectiveness for mixed noise removal. When |f (x)-u(x)| is large, the noise component is interpreted as Salt & Pepper by the model and then Φ λ 1 ,λ 2 behaves locally as ∥f -u∥ L 1 (Ω) . On the other hand, if |f (x) -u(x)| is small, then the model assumes that the noise is Gaussian and enforces a data fidelity Φ λ 1 ,λ 2 locally equal to ∥f -u∥ 2 L 2 (Ω) , see Figure 2 for a visualisation.

Asymptotic behaviour

In this section, we investigate the asymptotic behaviour of the L 1 -L 2 IC model. Here, we do not need to restrict to the TV regulariser, but we can consider a more general regularisation functional J with the following properties: (i) J : L 1 (Ω) → R ∪ {∞} is positive, proper, convex, lower semicontinuous with respect to the strong convergence in L 1 (Ω).

(ii) There exist constants

C 1 , C 2 > 0 such that C 1 |Du|(Ω) ≤ J(u) ≤ C 2 |Du|(Ω) for every u ∈ BV(Ω).
Classical regularisers such as TV, Huber-TV and TGV of any order, satisfy the above properties. Note that this is also true for a large class of structural TV-type functionals that are commonly used in inverse problems, see [START_REF] Hintermüller | A function space framework for structural total variation regularization with applications in inverse problems[END_REF].

We are interested in the following general problem:

min u∈BV(Ω) v∈L 1 (Ω) J(u) + λ 1 ∥v∥ L 1 (Ω) + λ 2 2 ∥f -u -v∥ 2 L 2 (Ω) , (2.4)
which is a more general version of the TV-IC model for Gaussian and Salt & Pepper noise removal

min u∈BV(Ω) v∈L 1 (Ω) |Du|(Ω) + λ 1 ∥v∥ L 1 (Ω) + λ 2 2 ∥f -u -v∥ 2 L 2 (Ω) . (2.5)
The well-posedness of (2.5) has been studied in [START_REF] Calatroni | Infimal convolution of data discrepancies for mixed noise removal[END_REF]. For the more general model (2.4) existence of minimisers (u * , v * ) ∈ BV(Ω) × L 1 (Ω) follows from a simple application of the direct method of calculus of variations. Note however, that the solution u * is not necessary unique as Φ λ 1 ,λ 2 (•, f ) is not strictly convex (J is not necessarily strictly convex either). The same holds for v * due to formula (2.3) which connects it with u * .

One natural question one may ask is in what degree we can expect to recover the single noise models by sending the parameters λ 1 , λ 2 (or their ratio) to infinity. In the following we answer this question by taking advantage of the formulation (2.1)-(2.2) and using some Γ-convergence arguments. Firstly, we extend a corresponding proposition that was shown in [14, Proposition 5.1], to the general regulariser case, adjusted for our purposes. Proposition 2.2. Let (u * , v * ) ∈ BV(Ω) × L 1 (Ω) be an optimal pair for (2.4). Then, the following assertions hold:

(i) If λ 1 → ∞, f ∈ L 1 (Ω) then v * → 0 in L 1 (Ω). (ii) If λ 2 → ∞, f ∈ L 2 (Ω) then ∥f -u * -v * ∥ L 2 (Ω) → 0. If in addition λ 1 is fixed, then the same result holds with f ∈ L 1 (Ω). (iii) If both λ 1 , λ 2 → ∞ and f ∈ L 2 (Ω) then (i) holds and we have that u * → f in L 1 (Ω). If f ∈ BV(Ω) this convergence is also weakly * in BV(Ω).
Proof. (i) We notice that

λ 1 ∥v * ∥ L 1 (Ω) ≤ J(u) + λ 1 ∥v∥ L 1 (Ω) + λ 2 2 ∥f -u -v∥ 2 L 2 (Ω) , ∀u ∈ BV(Ω), v ∈ L 1 (Ω), which by setting v = f -u, implies λ 1 ∥v * ∥ L 1 (Ω) ≤ J(u) + λ 1 ∥f -u∥ L 1 (Ω) , ∀ u ∈ BV(Ω), v ∈ L 1 (Ω). (2.6) Given ϵ > 0, we can find u ϵ ∈ C ∞ c (Ω) such that ∥f -u ϵ ∥ L 1 (Ω) < ϵ. Thus (2.6) becomes ∥v * ∥ L 1 (Ω) ≤ 1 λ 1 J(u ϵ ) + ϵ ⇒ lim sup λ 1 →∞ ∥v * ∥ L 1 (Ω) ≤ ϵ.
Since ϵ was arbitrary, the result follows.

(ii) In this case we have that for every u ∈ BV(Ω) and v ∈ L 1 (Ω)

λ 2 2 ∥f -u * -v * ∥ 2 L 2 (Ω) ≤ J(u) + λ 1 ∥v∥ L 1 (Ω) + λ 2 2 ∥f -u -v∥ 2 L 2 (Ω) , (2.7)
which by setting v = 0, implies

λ 2 2 ∥f -u * -v * ∥ 2 L 2 (Ω) ≤ J(u) + λ 2 2 ∥f -u∥ 2 L 2 (Ω) , ∀ u ∈ BV(Ω).
Then, proceeding as in step (i) the result follows. Notice that if we assume that λ 1 is fixed (or more generally bounded from above) and by merely assuming f ∈ L 1 (Ω), we can have the same result by setting v = f and u = 0 in (2.7).

(iii

) Notice that if both λ 1 , λ 2 → ∞ and f ∈ L 2 (Ω) then from (i), (ii) we get that ∥v * ∥ L 1 (Ω) → 0 and ∥f -u * -v * ∥ L 1 (Ω) → 0 and hence an application of the triangle inequality implies that u * → f in L 1 (Ω). If in addition f ∈ BV(Ω), then we have that |Du * |(Ω) is uniformly bounded by setting v = 0, u = f in C 1 |Du * |(Ω) ≤ J(u * ) ≤ J(u) + λ 1 ∥v∥ L 1 (Ω) + λ 2 2 ∥f -u -v∥ 2 L 2 (Ω) , ∀ u ∈ BV(Ω), v ∈ L 1 (Ω). From compactness in BV(Ω) and the fact that ∥f -u * ∥ L 1 (Ω) → 0 we infer that u * → f weakly * in BV(Ω).
In what follows, we refine the result above and prove convergence of the minimisers of (2.4) to the minimisers of the single noise models. To do so, we apply Γ-convergence arguments [START_REF] Dal Maso | Introduction to Γ-convergence[END_REF] to the IC term Φ λ 1 ,λ 2 .

Proposition 2.3. Let f ∈ L 1 (Ω) and let us define the functional

F λ 1 ,λ 2 : L 1 (Ω) → R + by F λ 1 ,λ 2 (u) := Φ λ 1 ,λ 2 (u, f ). Then (i) For any fixed λ 1 , F λ 1 ,λ 2 Γ-converges to F 1 (•) := λ 1 ∥f -•∥ L 1 (Ω) as λ 2 → ∞. (ii) For any fixed λ 2 , F λ 1 ,λ 2 Γ-converges to F 2 (•) as λ 1 → ∞, where for every u ∈ L 1 (Ω), F 2 is defined as F 2 (u) := λ 2 2 ∥f -u∥ 2 L 2 (Ω) := { λ 2 2 ∥f -u∥ 2 L 2 (Ω) , if f -u ∈ L 2 (Ω), +∞, if f -u ∈ L 1 (Ω) \ L 2 (Ω).
Proof.

For (i), let ( λ (n) 2
) n∈N be a sequence with λ

(n) 2 → ∞ and set F n := F λ 1 ,λ (n) 2 . We notice that F n converges uniformly to λ 1 ∥f -•∥ L 1 (Ω) . Indeed, for u ∈ L 1 (Ω), we have F n (u) -λ 1 ∥f -u∥ L 1 (Ω) = ∫ Ω φ(f -u) -λ 1 |f -u| dx ≤ ∫ |f -u|≥ λ 1 λ (n) 2 λ 2 1 2λ (n) 2 dx + ∫ |f -u|< λ 1 λ (n) 2 λ (n) 2 2 |f -u| 2 -λ 1 |f -u| dx ≤ ∫ Ω λ 2 1 2λ (n) 2 dx + ∫ Ω λ 2 1 2λ (n) 2 dx + ∫ Ω λ 2 1 λ (n) 2 dx → 0 as n → ∞.
Since the last limit is independent of u, the convergence of the functionals is indeed uniform. Moreover F n is continuous with respect to the L 1 topology, see for instance [START_REF] Calatroni | Infimal convolution of data discrepancies for mixed noise removal[END_REF]. Thus from [21, Proposition 5.2] we immediately get that F n Γ-converges to F 1 as n → ∞. For (ii), we now set

F n := F λ (n) 1 ,λ 2 with λ (n) 1 → ∞ and observe that F n converges pointwise to F 2 . Indeed if f -u ∈ L 1 (Ω) \ L 2 (Ω) then F n (u) ≥ ∫ |f -u|< λ (n) 1 λ 2 λ 2 2 |f -u| 2 dx → ∞ as n → ∞.
On the other hand, if f -u ∈ L 2 (Ω), proceeding as before we have

F n (u) - λ 2 2 ∥f -u∥ 2 L 2 (Ω) ≤ ∫ |f -u|≥ λ (n) 1 λ 2 λ (n) 1 |f -u| - ( λ (n) 1 ) 2 2λ 2 dx + ∫ |f -u|≥ λ (n) 1 λ 2 λ 2 2 |f -u| 2 dx ≤ ∫ |f -u|≥ λ (n) 1 λ 2 3 2 λ 2 |f -u| 2 dx → 0 as n → ∞. If ( λ (n) 1
)

n∈N is increasing, it can be easily verified that the sequence (F n ) n∈N is increasing. Moreover the L 2 norm is lower semicontinuous with respect to the strong L 1 topology. Hence from [21, Remark 5.5], we have that F n Γ-converges to F 2 . In the case where ( λ

(n) 1
) n∈N is non-monotonically going to infinity, we can find an increasing subsequence and then the result follows from the Urysohn property of Γ-convergence, see [START_REF] Dal Maso | Introduction to Γ-convergence[END_REF]Proposition 8.3].

As a corollary, we obtain the following result on the convergence of minimisers.

Corollary 2.4 (Convergence to single noise models).

The following two results hold:

(i) Let f ∈ L 1 (Ω), λ 1 > 0 fixed and ( λ (n) 2 ) n∈N with λ (n) 2 → ∞. If (u n ) n∈N ⊆ BV(Ω) is a sequence of minimisers of (2.4), then every subsequence of (u n ) n∈N has a weak * in BV(Ω) cluster point which is a minimiser of min u∈BV(Ω) J(u) + λ 1 ∥f -u∥ L 1 (Ω) .
(2.8)

Moreover if the solution u * of (2.8) is unique, then u n → u * weakly * in BV(Ω) and
J(u n ) + Φ λ 1 ,λ (n) 2 (u n , f ) → J(u * ) + λ 1 ∥f -u * ∥ L 1 (Ω) .
(2.9)

(ii) Let f ∈ L 2 (Ω), λ 2 > 0 fixed and ( λ (n) 1 ) n∈N with λ (n) 1 → ∞. If (u n ) n∈N ⊆ BV(Ω) is a sequence of minimisers of (2.4), then u n → u * weakly * in BV(Ω), where u * is the unique minimiser of min u∈BV(Ω) J(u) + λ 2 2 ∥f -u∥ 2 L 2 (Ω) .
(2.10)

Moreover J(u n ) + Φ λ (n) 1 ,λ 2 (u n , f ) → J(u * ) + λ 2 2 ∥f -u * ∥ 2 L 2 (Ω) .
Proof. (i) Since the functional J is lower semicontinuous with respect to L 1 , then [START_REF] Dal Maso | Introduction to Γ-convergence[END_REF]Proposition 6.25] we have that the minimising functionals

Φ λ 1 ,λ (n) 2 (•, f ) + J(•) also Γ-converge to the functional in (2.8). Moreover, note that (u n ) n∈N is uniformly bounded in BV(Ω).
Since u n is a minimiser we have in fact that for every u ∈ BV(Ω) and every v ∈ L 1 (Ω), the following three inequalities hold

C 1 |Du n |(Ω) ≤ J(u n ) ≤ J(u) + λ 1 ∥v∥ L 1 (Ω) + λ (n) 2 2 ∥f -u -v∥ 2 L 2 (Ω) , (2.11) λ (n) 2 2 ∥f -u n -v n ∥ 2 L 2 (Ω) ≤ J(u) + λ 1 ∥v∥ L 1 (Ω) + λ (n) 2 2 ∥f -u -v∥ 2 L 2 (Ω) ,
(2.12)

λ 1 ∥v n ∥ L 1 (Ω) ≤ J(u) + λ 1 ∥v∥ L 1 (Ω) + λ (n) 2 2 ∥f -u -v∥ 2 L 2 (Ω) . (2.13)
By setting u = 0, v = f in (2.11) one obtains a uniform bound for the sequence

(|Du n |(Ω)) n∈N .
From (2.12) and from the fact that Ω is bounded, one obtains a uniform bound for

(∥f -u n -v n ∥ L 1 (Ω) ) n∈N . Similarly, from (2.13), a uniform bound on (∥v n ∥ L 1 (Ω) ) n∈N is obtained and this means that (∥u n ∥ L 1 (Ω)
) n∈N is also bounded. Thus every subsequence of u n has a cluster point in BV(Ω) with respect to the weak * topology, which must be a minimiser of (2.8), [START_REF] Dal Maso | Introduction to Γ-convergence[END_REF]Corollary 7.20]. Furthermore if (2.8) has a unique minimiser, then every subsequence of u n has a further subsequence that converges to u weakly * in BV(Ω). Thus, in this case u n converges to u weakly * in BV(Ω) and moreover (2.9) holds by [START_REF] Dal Maso | Introduction to Γ-convergence[END_REF]Corollary 7.20].

(ii) The Γ-convergence of the energies follows as above. By setting, v = u = 0 in (2.11) one obtains a uniform bound on (|Du n |(Ω)) n∈N and similarly as in (i) a bound on L 1 (Ω) and consequently in BV(Ω) is obtained for u. The rest of the proof follows as in (i), bearing in mind that the solution of (2.10) is unique.

We summarise our findings so far in Table 1, where we have combined the results of Proposition 2.2 and Corollary 2.4.

In the case of bounded data and TV regularisation, the results obtained above can be refined. We first recall the following well-known result, see [START_REF] Choksi | A few remarks on variational models for denoising[END_REF]Lemma 3.5]. Proposition 2.5. Let u be a solution of (2.5), with f ∈ L ∞ (Ω). Then the following maximum principle holds:

ess inf x∈Ω f (x) ≤ ess inf x∈Ω u(x) ≤ ess sup x∈Ω u(x) ≤ ess sup x∈Ω f (x).
We can now prove the following result for the TV-IC minimisation problem (2.5). Proposition 2.6. Suppose that f ∈ L ∞ (Ω) and the parameters λ 1 , λ 2 > 0 satisfy the following condition

λ 1 λ 2 ≥ 2∥f ∥ ∞ . (2.14) λ 1 → ∞ f ∈ L 1 (Ω) λ 1 → ∞ λ 2 fixed f ∈ L 2 (Ω) λ 2 → ∞ f ∈ L 2 (Ω) (or f ∈ L 1 (Ω) & λ 1 fixed) λ 2 → ∞ λ 1 fixed f ∈ L 1 (Ω) J-λ 1 L 1 has !sol. u * λ 1 → ∞ λ 2 → ∞ f ∈ L 2 (Ω) v v → 0 in L 1 (Ω) v → 0 in L 1 (Ω) ∥f -u -v∥ L 2 (Ω) → 0 v → f -u * in L 1 (Ω) v → 0 in L 1 (Ω) u cannot say u → solution J-λ2 2 L 2 , w * in BV(Ω) cannot say u → u * , w * in BV(Ω) u → f in L 1 (Ω), (w * in BV(Ω) if f ∈ BV(Ω))
Table 1: Summary of all the asymptotic results concerning the solution pair u, v of (2.4) when one or both parameters λ 1 and λ 2 are let to infinity.

Then, if u is a solution of (2.5), there holds

Φ λ 1 ,λ 2 (u, f ) = λ 2 2 ∥f -u∥ 2 L 2 (Ω) . (2.

15)

As a result, the problem (2.5) is equivalent to a standard TV-L 2 minimisation problem.

Proof. This is a direct consequence of Proposition 2.5 and the formulation (2.1)-(2.2) of Φ λ 1 ,λ 2 (u, f ). Indeed, using a translation argument and Proposition 2.5 one shows directly that for the solution u of (2.5), it holds ∥f -u∥ ∞ ≤ 2∥f ∥ ∞ . Thus if (2.14) holds, (2.1) implies that φ(f -u) = λ 2 2 |f -u| 2 so that (2.15) holds as well. We note here that the adaptation of Proposition 2.5 and, consequently, of Proposition 2.6 to other widely used regularisers is not immediate. For instance, it remains an open problem to show that the solution u of the TGV-L 2 problem with data f ∈ L ∞ (Ω) is also an L ∞ function, see for instance the corresponding discussion in [START_REF] Valkonen | The jump set under geometric regularisation. Part 2: Higher-order approaches[END_REF]. However, in dimension one this fact is true when f ∈ BV(Ω), by taking advantage of the estimate ∥u∥ L ∞ (Ω) ≤ C∥u∥ BV(Ω) , see for instance [START_REF] Papafitsoros | Asymptotic behaviour of total generalised variation[END_REF]Proposition 2].

In view of the Proposition above, one sees that in the case of TV regularisation, the Gaussian noise model can be recovered simply by fixing either λ 2 and setting λ 1 large enough or by fixing λ 1 and setting λ 2 small enough. In Figure 3 we graphically depict this behaviour.

Convergence of the parameters to zero and non-exact recovery of the data

We now study the asymptotic behaviour of the model when the parameters are sent to zero. For this analysis and for the sake of simplicity, we focus on the TV minimisation model (2.5) but the results can be easily extended to the general regulariser case. First, we recall the

λ 2 λ 1 λ 1 = 2∥f ∥ ∞ λ 2 TV-L 2 TV-L 1 Figure 3: If λ 1 λ 2 ≥ 2∥f ∥ ∞ then (2.
u Ω := ∫ Ω u dx, u Ω ∈ argmin c∈R ∫ Ω |u -c| dx.
Remark 2.7. Note that the median value is not necessarily unique. Moreover, if u ∈ L 2 (Ω)

then u Ω = argmin c∈R ∫ Ω |u -c| 2 dx.
We have the following result:

Proposition 2.8. Let f ∈ L 1 (Ω) and ( λ (n) 1 ) n∈N , ( λ (n) 2 ) n∈N two sequences such that λ (n) 1 → 0 and λ (n) 1 λ (n) 2 → 0, as n → ∞. (2.16)
Then, denoting by (u n ) n∈N the sequence of the corresponding solutions of (2.5) we have that

u n → f Ω weakly * in BV(Ω).
By this we mean that every subsequence of (u n ) n∈N has a further subsequence converging to a median of f .

Proof. Let (u n ) n∈N ⊆ BV(Ω) be a sequence of the corresponding solutions for the parameters ( λ

(n) 1 , λ (n) 2
)

. Notice that the sequence (u n ) n∈N is uniformly bounded in BV(Ω). Indeed the bound on TV is obtained again from (2.11), while for the L 1 bound, we first observe that:

- λ 2 1 2λ 2 |Ω| + λ 1 ∥f -u∥ L 1 (Ω) ≤ ∫ Ω φ(f -u) dx, ∀u ∈ L 1 (Ω).
(2.17) Using (2.17), we get

- ( λ (n) 1 ) 2 2λ (n) 2 |Ω| + λ (n) 1 ∥f -u n ∥ L 1 (Ω) ≤ ∫ Ω φ(f -u n ) dx ≤ λ (n) 1 ∥f ∥ L 1 (Ω) ⇒ - λ (n) 1 2λ (n) 2 |Ω| + ∥f -u n ∥ L 1 (Ω) ≤ ∥f ∥ L 1 (Ω) ,
where the bound is obtained using (2.16). Thus every subsequence of (u n ) n∈N has a further (not relabelled) subsequence converging to an element u * ∈ BV(Ω). We will show that u * is a median of f . Notice first that since λ

(n)
1 → 0, then from (2.11) we get that |Du n |(Ω) → 0. Thus, from the lower semicontinuity of total variation and the fact that Ω is connected we get that u * is a constant. By thus setting u = c ∈ R and v = f -c along with (2.17) in

|Du n |(Ω) + ∫ Ω φ(f -u n ) dx ≤ |Du|(Ω) + λ (n) 1 ∥v∥ L 1 (Ω) + λ (n) 2 2 ∥f -u -v∥ 2 L 2 (Ω) , we get - ( λ (n) 1 ) 2 2λ (n) 2 |Ω| + λ (n) 1 ∥f -u n ∥ L 1 (Ω) ≤ λ (n) 1 ∥f -c∥ L 1 (Ω) ⇒ - λ (n) 1 λ (n) 2 |Ω| + ∥f -u n ∥ L 1 (Ω) ≤ ∥f -c∥ L 1 (Ω) ⇒ (by taking limits) ∥f -u * ∥ L 1 (Ω) ≤ ∥f -c∥ L 1 (Ω) .
Hence since u * is constant and c ∈ R was arbitrary, we have that u is a median of f .

Similarly, we have the following result:

Proposition 2.9. Let f ∈ L 2 (Ω) and ( λ (n) 1 ) n∈N , ( λ (n) 2 ) n∈N two sequences such that λ (n) 2 → 0 and λ (n) 2 λ (n) 1 → 0, as n → ∞.
Then for the corresponding solutions (u n ) n∈N of (2.5) we have that

u n → f Ω weakly * in BV(Ω).
Proof. The proof follows the same steps as in Proposition 2.8. First, observe that the sequence of solutions (u n ) n∈N is bounded in BV(Ω). Indeed, from (2.11) we have that |Du n |(Ω) → 0. From (2.13) we further get that (v n ) n∈N is bounded in L 1 and from (2.12) we get that

(f -u n -v n ) n∈N is bounded in L 1 .
Thus, from the triangle inequality we have that (u n ) n∈N is bounded in L 1 . Hence, there exists a subsequence of (u n ) n∈N that converges to a function u * weakly * in BV(Ω) with u * being a constant. It remains to show that u * is the mean value of f . As before, we have for an arbitrary c ∈ R

|Du n |(Ω) + ∫ Ω φ(f -u n ) dx ≤ λ (n) 2 2 ∥f -c∥ 2 L 2 (Ω) ⇒ 1 2 ∫ |f -un|< λ (n) 1 λ (n) 2 |f -u n | 2 dx ≤ 1 2 ∥f -c∥ 2 L 2 (Ω)
⇒ (using Fatou's Lemma)

1 2 ∥f -u * ∥ 2 L 2 (Ω) ≤ 1 2 ∥f -c∥ 2 L 2 (Ω) .
Since u is a constant and c ∈ R was arbitrary the proof is complete.

The following proposition states that with the L 1 -L 2 infimal convolution fidelity model we can never expect exact recovery of the data. This is similar to the pure L 2 model, see also [START_REF] Choksi | A few remarks on variational models for denoising[END_REF]Proposition 4.1].

Proposition 2.10. Let f ∈ L 1 (Ω) and u * to be a solution of the minimisation problem (2.5). Then u * = f if and only if f is a constant.

Proof. One direction is straightforward. Suppose now that f is a solution of (2.5). Note that in this case necessarily we must have f ∈ BV(Ω) ⊆ L d * (Ω), where d * = d/(d -1), see [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF]. It follows that for every 0 < ϵ < 1, the function f ϵ := ϵf is suboptimal. Thus we have

|Df |(Ω) ≤ ϵ|Df |(Ω) + Φ λ 1 ,λ 2 (u, ϵf ) =⇒ 0 ≤ (ϵ -1)|Df |(Ω) + Φ λ 1 ,λ 2 (u, ϵf ).
We continue

0 ≤ (ϵ -1)|Df |(Ω) + Φ λ 1 ,λ 2 (u, ϵf ) 0 ≤ (ϵ -1)|Df |(Ω) + ∫ |f -ϵf |≥ λ 1 λ 2 λ 1 |f -ϵf | - λ 2 1 2λ 2 dx + ∫ |f -ϵf |< λ 1 λ 2 λ 2 2 |f -ϵf | 2 dx+ =⇒ 0 ≥ |Df |(Ω) - ∫ |f |≥ λ 1 λ 2 |1-ϵ| λ 1 |f | dx + ∫ |f |≥ λ 1 λ 2 |1-ϵ| λ 2 1 2λ 2 |1 -ϵ| dx -|1 -ϵ| ∫ |f |< λ 1 λ 2 |1-ϵ| λ 2 2 |f | 2 dx.
(2.18)

Now working with each one of the first three terms in (2.18), using dominated convergence, we have

lim ϵ→1 ∫ |f |≥ λ 1 λ 2 |1-ϵ| |f | dx = 0, (2.19
)

lim ϵ→1 ∫ |f |≥ λ 1 λ 2 |1-ϵ| λ 2 1 2λ 2 |1 -ϵ| dx = λ 1 2 lim t→∞ ∫ |f |≥t t dx ≤ λ 1 2 lim t→∞ ∫ |f |≥t |f | dx = 0, (2.20)
and finally

lim ϵ→1 |1 -ϵ| ∫ |f |< λ 1 λ 2 |1-ϵ| λ 2 2 |f | 2 dx = λ 2 2 lim ϵ→1 |1 -ϵ| ∫ |f |< λ 1 λ 2 |1-ϵ| |f | d d-1 |f | d-2 d-1 ≤ λ 2 2 ( λ 1 λ 2 ) d-2 d-1 lim ϵ→1 |1 -ϵ| |1 -ϵ| d-2 d-1 ∫ |f |< λ 1 λ 2 |1-ϵ| |f | d d-1 dx ≤ λ 2 2 ( λ 1 λ 2 ) d-2 d-1 ∥f ∥ d * L d * (Ω) lim ϵ→1 |1 -ϵ| 1 d-1 = 0. ( 2
but if d = 1 then f ∈ L ∞ (Ω) ⊆ L 2 (Ω)
and the analogous calculation to (2.21) follows more easily.

It is clear from all the analysis above that the TV-IC model, at least when f ∈ L ∞ (Ω), can reproduce the TV-L 2 solutions, but as far as the TV-L 1 solutions are concerned, these are only (guaranteed to be) recovered in the limit λ 2 → ∞, see again Figure 3.

The one-homogeneous analogue

We conclude this section by briefly presenting an alternative form of (1.3), i.e., its onehomogeneous analogue, by which the TV-L 1 solutions can also be recovered for finite parameters. This discussion is motivated by some analogous results in [START_REF] Burger | Infimal convolution regularisation functionals of BV and L p spaces. Part I: The finite p case[END_REF]. We define:

Φ λ 1 ,λ 2 1-hom (u, f ) := min v∈L 1 (Ω) λ 1 ∥v∥ L 1 (Ω) + λ 2 ∥f -u -v∥ L 2 (Ω) ,
which, via a straightforward computation gives

min ( λ 1 , λ 2 2|Ω| 1/2 ) ∥f -u∥ L 1 (Ω) ≤ min v∈L 1 (Ω) λ 1 ∥v∥ L 1 (Ω) + λ 2 2|Ω| 1/2 ∥f -u -v∥ L 1 (Ω) ≤ min v∈L 1 (Ω) λ 1 ∥v∥ L 1 (Ω) + λ 2 2 ∥f -u -v∥ L 2 (Ω) ≤ λ 1 ∥f -u∥ L 1 (Ω) .
Hence, it is clear that

Φ λ 1 ,λ 2 1-hom (u, f ) = λ 1 ∥f -u∥ L 1 (Ω) , if λ 1 λ 2 ≤ 1 2|Ω| 1/2 ,
Thus, under such choice the TV-L 1 model can be recovered. We state here without any proof the relationship between these two versions. Let us define the following sets for f ∈ L 1 (Ω):

S IC = { u * ∈ BV(Ω) : u * = argmin u∈BV(Ω) |Du|(Ω) + Φ λ 1 ,λ 2 (u, f ), for some λ 1 , λ 2 > 0 } , S 1-hom IC = { u * ∈ BV(Ω) : u * = argmin u∈BV(Ω) |Du|(Ω) + Φ λ 1 ,λ 2 1-hom (u, f ), for some λ 1 , λ 2 > 0 } , S L 1 = { u * ∈ BV(Ω) : u * = argmin u∈BV(Ω) |Du|(Ω) + λ 1 ∥f -u∥ L 1 (Ω) , for some λ 1 > 0 } .
Then one can show by using similar techniques as in [START_REF] Burger | Infimal convolution regularisation functionals of BV and L p spaces. Part I: The finite p case[END_REF] that

S 1-hom IC = S IC ∪ S L 1 and that, in general, S L 1 \ S IC ̸ = ∅.

Exact solutions

In order to get more insights about the relationship of the TV-IC model with the pure TV-L 1 and TV-L 2 models, we compute in this section some exact solutions for simple one dimensional data functions f . In particular, we set here Ω = (-2L, 2L) for some L > 0, and we consider as data f the following step function

f (x) = { 0, if x ∈ (-2L, -L) ∪ (L, 2L), h, if x ∈ [-L, L], ( 3.1)
where h > 0. Using similar techniques as in [START_REF] Bredies | Properties of L 1 -TGV 2 : The one-dimensional case[END_REF][START_REF] Burger | Infimal convolution regularisation functionals of BV and L p spaces. Part I: The finite p case[END_REF][START_REF] Papafitsoros | A study of the one dimensional total generalised variation regularisation problem[END_REF][START_REF] Ring | Structural properties of solutions to total variation regularization problems[END_REF], we can easily show using primaldual optimality conditions, that a function u ∈ BV(Ω) is a solution of (2.5) if and only if there exists a function v ∈ H 1 0 (Ω) such that

v ′ = { λ 1 f -u |f -u| , if |f -u| ≥ λ 1 λ 2 , λ 2 (f -u), if |f -u| < λ 1 λ 2 , (3.2) v ∈ Sgn(Du), (3.3)
where

Sgn(Du) = { v ∈ L ∞ (Ω) ∩ L ∞ (Ω, Du) : ∥v∥ ∞ ≤ 1, v = dDu d|Du| , |Du|-a.e.
} .

Here dDu d|Du| denotes the Radon-Nikodým density of Du with respect to |Du|. Compared to the aforementioned references, the only difference here is the right-hand side of (3.2) which is the subdifferential of Φ λ 1 ,λ 2 evaluated at f -u. With the help of the optimality conditions above, we are able to compute analytically all the solutions to the problem (2.5) for the data (3.1), and for all combinations of the parameters λ 1 , λ 2 . Note that similarly to pure L 1 and L 2 cases one can show that no new jump discontinuities are created for the solution u, which will be constant in the areas where {f ̸ = u}. Thus all the solutions will be either constants or piecewise constants with jumps at x = -L and x = L which must also have the same orientations with the jumps of f . We first examine the case u = h 2 , i.e., the mean value of f . In such case we have |f -u| = h 2 everywhere and thus if h 2 ≤ λ 1 λ 2 then v ′ = λ 2 (f -u) everywhere. In order for the condition (3.2) to hold we must also have λ 2 ≤ 2 hL . Note that if h 2 ≤ λ 1 λ 2 and λ 2 < 2 hL , then the solution must be constant otherwise one can check that in every case the condition (3.3) would be violated. One can further show in this case that if u is a constant with h 2 < λ 1 λ 2 , the only possibility is u = h 2 . Suppose now that h 2 ≥ λ 1 λ 2 and also λ 1 < 1 L . Observe that in this case, every constant

function u = c with λ 1 λ 2 ≤ c ≤ h -λ 1 λ 2 satisfies (3.2)-(3.
3). In that case we have |f -u| ≥ λ 1 λ 2 everywhere. Notice again that no other constant function is a solution. By contradiction, that would mean that |f -u| > λ 1 λ 2 and |f -u| < λ 1 λ 2 on (-2L, L) ∪ (L, 2L) and (-L, L) respectively (or vice versa). With the help of (3.2) and the fact that v ∈ H 1 0 (Ω) one would then arrive to a contradiction. Furthermore, one can check that discontinuous solutions cannot occur in this case either.

We concentrate now on the case h 2 ≥ λ 1 λ 2 and λ 1 = 1 L . Note that the constant functions u = c with λ 1 λ 2 ≤ c ≤ h -λ 1 λ 2 are solutions in this case as well. However, one can also verify that the following family of discontinuous functions are also solutions:

u(x) =        c 1 , if x ∈ (-2L, -L), h -d, if x ∈ [-L, L], c 2 , if x ∈ (L, 2L), λ 1 λ 2 ≤ c i < h -d ≤ h - λ 1 λ 2 , i = 1, 2.
It can be checked similarly as before that no other solutions can occur. Finally, we consider the case λ 1 > 1 L and λ 2 ≥ 2 hL . We claim that in that case the unique solution is given by

u(x) = { 1 Lλ 2 , if x ∈ (-2L, -L) ∪ (L, 2L), h -1 Lλ 2 , if x ∈ [-L, L].
One can similarly check that no other solution is possible. We summarise our findings in the following proposition.

|Du|(Ω) + Φ λ 1 ,λ 2 (u, f ),
are given by the following formulae:

(i) If h 2 < λ 1 λ 2 and λ 2 ≤ 2 hL
, then the solution is unique and given by u = h 2 .

(ii) If h 2 ≥ λ 1 λ 2 and λ 1 < 1 L , then there exist infinitely many constant solutions given by

u = c, λ 1 λ 2 ≤ c ≤ h - λ 1 λ 2 . λ 1 λ 2 0 2 hL 1 L λ 1 = h 2 λ 2 Unique solution λ1 λ2 λ1 λ2 λ1 λ2 λ1 λ2

Multiple solutions

Multiple solutions

1 Lλ2 1 Lλ2
Unique solution

(iii) If h

2 ≥ λ 1 λ 2 and λ 1 = 1 L , then there there exist infinitely many solutions given by

u(x) =        c 1 , if x ∈ (-2L, -L), h -d, if x ∈ [-L, L], c 2 , if x ∈ (L, 2L), λ 1 λ 2 ≤ c i ≤ h -d ≤ h - λ 1 λ 2 , i = 1, 2.
(iv) If λ 2 > 2 hL and λ 1 > 1 L , then the solution is unique and given by

u(x) = { 1 Lλ 2 , if x ∈ (-2L, -L) ∪ (L, 2L), h -1 Lλ 2 , if x ∈ [-L, L].
A visualisation of these solutions is depicted in Figure 4. Observe that for large enough ratio λ 1 λ 2 all the TV-L 2 solutions are recovered as Proposition 2.6 predicts -compare also Figures 3 and4. Moreover, observe that as λ 1 and λ 1 λ 2 goes to zero, the solutions indeed converge to a median of f , as shown in Proposition 2.8. Note however, that apart from some medians, the solutions of the TV-L 1 model are not recovered. More precisely, the ones that perfectly fit the data in the whole domain or part of it, i.e., u = h, u = 0 and u = f cannot be obtained here. This is also in accordance to Proposition 2.10.

Automatic selection of parameters

We describe now a bilevel optimisation strategy for the estimation of optimal parameters λ 1 and λ 2 in the TV-IC model (2.5) based on the use of training sets, see [START_REF] Calatroni | Bilevel approaches for learning of variational imaging models[END_REF][START_REF] De Los Reyes | Image denoising: learning the noise model via nonsmooth PDE-constrained optimization[END_REF][START_REF] De Los Reyes | Bilevel parameter learning for higher-order Total Variation regularisation models[END_REF]. This approach has been heuristically considered for the TV-IC model in [14, Section 7] with little theoretical justification. To fill this gap, we prove in this section existence results for the solution of the bilevel minimisation problem and for the corresponding adjoint problem, thus making the derivation the optimality systemin [START_REF] Calatroni | Infimal convolution of data discrepancies for mixed noise removal[END_REF] rigorous.

We point out that in [START_REF] Langer | Automated parameter selection for total variation minimization in image restoration[END_REF][START_REF] Langer | Automated parameter selection in the L 1 -L 2 -TV model for removing gaussian plus impulse noise[END_REF] an adaptive optimisation approach has been proposed for the automatic selection of parameters when a linear combination of L 1 and L 2 data fidelities is considered. However, differently to our setting, in that approach the noise level is assumed to be known.

Bilevel optimisation

Learning approaches have become very popular over the recent years due to their ability of combining data-and model-driven algorithms for the optimal design of imaging models. In particular, bilevel optimisation techniques have been proposed by several different authors in discrete [START_REF] Klatzer | Continuous hyper-parameter learning for support vector machines[END_REF][START_REF] Kunisch | A bilevel optimization approach for parameter learning in variational models[END_REF][START_REF] Ochs | Bilevel optimization with nonsmooth lower level problems[END_REF][START_REF] Yu | Learning Reaction-Diffusion Models for Image Inpainting[END_REF] and functional [START_REF] Calatroni | Bilevel approaches for learning of variational imaging models[END_REF][START_REF] De Los Reyes | Image denoising: learning the noise model via nonsmooth PDE-constrained optimization[END_REF][START_REF] De Los Reyes | The structure of optimal parameters for image restoration problems[END_REF][START_REF] De Los Reyes | Bilevel parameter learning for higher-order Total Variation regularisation models[END_REF][START_REF] Hintermüller | Optimal selection of the regularization function in a weighted total variation model. part I: Modelling and theory[END_REF][START_REF] Hintermüller | Optimal selection of the regularization function in a weighted total variation model. part II: Algorithm, its analysis and numerical tests[END_REF] settings as a tool to estimate the "best" variational image restoration model within a certain class by means of training examples. Typically, such examples consist of images obtained in standard acquisition settings, and thus corrupted by noise with equal (unknown) intensity, paired with their corresponding versions ideally acquired in a very low-noise setting. In medical imaging, for instance, such training set can be provided by means of real and/or simulated phantoms. Note that in order to make the estimation robust, a large training set is often desirable; for that, stochastic optimisation techniques and sampling approaches can be used to reduce the computational costs, see, e.g., [START_REF] Calatroni | Dynamic sampling schemes for optimal noise learning under multiple nonsmooth constraints[END_REF].

The general bilevel optimisation problem can be formulated as:

min λ∈[0,∞) m F (u λ ) (4.1)
subject to:

u λ ∈ argmin u∈X { J(u, λ) := |Du|(Ω) + Φ λ (u, f ) } , ( 4.2)
where λ = (λ 1 , . . . , λ m ) ∈ [0, ∞) m are the parameters to optimise and F ≥ 0 is an appropriate quality measure which is minimised under the constraint that the function u λ is a solution of the denoising model (4.2) in a suitable function space X.

In [START_REF] Calatroni | Bilevel approaches for learning of variational imaging models[END_REF][START_REF] De Los Reyes | Image denoising: learning the noise model via nonsmooth PDE-constrained optimization[END_REF] this approach is used to estimate the optimal parameters λ in the case when single data models Φ i , i = 1, . . . , m are linearly combined, i.e., when

Φ λ (u, f ) = m ∑ i=1 λ i Φ i (u, f ).
There, theoretical results showing existence of minima and the adjoint state for the problem (4.1)-(4.2) are shown and Newton-type methods are proposed for its efficient numerical solution. Similar results and algorithms are further studied in [START_REF] De Los Reyes | The structure of optimal parameters for image restoration problems[END_REF][START_REF] De Los Reyes | Bilevel parameter learning for higher-order Total Variation regularisation models[END_REF] for the estimation of optimal parameters of higher-order regularisers (e.g., TGV) combined with Gaussian fidelity.

In the following we set m = 2 and consider the problem of estimating the optimal parameters λ 1 and λ 2 in (2.5).

General framework:

The non-smooth TV-IC bilevel problem reads:

min λ 1 ,λ 2 ≥0 F (u λ 1 ,λ 2 ) subject to u λ 1 ,λ 2 ∈ argmin u∈BV(Ω) |Du|(Ω) + Φ λ 1 ,λ 2 (u, f ), (4.3)
where Φ λ 1 ,λ 2 is the IC fidelity (1.3). We follow [START_REF] Calatroni | Bilevel approaches for learning of variational imaging models[END_REF][START_REF] De Los Reyes | Image denoising: learning the noise model via nonsmooth PDE-constrained optimization[END_REF][START_REF] De Los Reyes | The structure of optimal parameters for image restoration problems[END_REF][START_REF] De Los Reyes | Bilevel parameter learning for higher-order Total Variation regularisation models[END_REF] and introduce an appropriate smoothing of the TV semi-norm combined with a further quadratic smoothing. This is crucial for the following proofs and for the design of the gradient-based optimisation algorithms we intend to use. For ϵ ≪ 1, we then consider the following regularised version of (4.3):

min λ 1 ,λ 2 ≥0 F (u λ 1 ,λ 2 ) subject to u λ 1 ,λ 2 ∈ argmin u∈H 1 (Ω) ϵ 2 ∥u∥ 2 H 1 (Ω) + ∥∇u∥ γ,L 1 (Ω) + Φ λ 1 ,λ 2 ϵ,γ (u, f ). (4.4)
Here, we denote by ∥∇u∥ γ,L 1 (Ω) = ∫ Ω |∇u| γ dx, with | • | γ being a smooth Huber-type regularisation depending on a parameter γ > 0 whose C 1 -derivative reads:

h γ (z) :=        z |z| if γ|z| -1 ≥ 1 2γ , z |z| (1 -γ 2 (1 -γ|z| + 1 2γ ) 2 ) if γ|z| -1 ∈ (-1 2γ , 1 2γ ), γz if γ|z| -1 ≤ -1 2γ . (4.5)
Note that | • | γ has one degree higher regularity than the classical Huber function φ in (2.2). This higher-order Huber-type smoothing has been previously used in [START_REF] Calatroni | Bilevel approaches for learning of variational imaging models[END_REF][START_REF] De Los Reyes | Bilevel parameter learning for higher-order Total Variation regularisation models[END_REF] for similar bilevel problems since it endows the problem (4.4) with further regularity, compare Theorem 4.6.

We then similarly regularise the IC fidelity term as:

Φ λ 1 ,λ 2 ϵ,γ (u, f ) := min v∈L 2 (Ω) { G λ 1 ,λ 2 ϵ,γ (v, u, f ) := ϵ 2 ∥v∥ 2 L 2 (Ω) + λ 1 ∥v∥ γ,L 1 (Ω) + λ 2 2 ∥f -u -v∥ 2 L 2 (Ω)
} , (4.6) where ∥v∥ γ,L 1 (Ω) is defined analogously as above. For simplicity, from now on, we will assume that f ∈ L 2 (Ω).

Inspired by [START_REF] Calatroni | Bilevel approaches for learning of variational imaging models[END_REF][START_REF] De Los Reyes | Bilevel parameter learning for higher-order Total Variation regularisation models[END_REF], we focus on two main choices of F . Namely, we consider the L 2 cost corresponding to Peak Signal to Noise Ratio (PSNR) optimisation

F L 2 (u λ 1 ,λ 2 ) := ∥u λ 1 ,λ 2 -ũ∥ 2 L 2 (Ω) , for training data ũ ∈ L 2 (Ω), (4.7)
and the Huberised TV cost, which is related to quality measures that are more adjusted to actual human perception, such as the Structural Similarity Index (SSIM):

F L 1 γ D (u λ 1 ,λ 2 ) := ∥D(u λ 1 ,λ 2 -ũ)∥ M,γ , for training data ũ ∈ BV(Ω). (4.8)
Here, ∥Du∥ M,γ := ∫ Ω |∇u| γ dx + |D s u|(Ω) so that the smooth Huber-type regularisation is applied on the absolutely continuous part of |Du|. For the abstract formulation of (4.4) in terms of a general F , we refer the reader to [START_REF] De Los Reyes | Bilevel parameter learning for higher-order Total Variation regularisation models[END_REF].

Remark 4.1. Note that if λ 1 = 0 and/or λ 2 = 0, then Φ λ 1 ,λ 2 (u, f ) = 0 for every u ∈ L 1 (Ω). In that case every constant function c is a minimiser of the lower level problem of (4.3).

Thus we trivially have

inf λ 1 ,λ 2 ≥0 F (u λ 1 ,λ 2 ) < ∞ where F is taken to be either F L 2 or F L 1 γ D .

Well-posedness of the TV-IC bilevel problem

We now discuss the well-posedness of the bilevel problems (4.3) and (4.4). For this type of problems, it is a common practice to impose an extra box constraint on the parameters λ 1 , λ 2 in order to ensure existence of solutions (see [START_REF] Calatroni | Bilevel approaches for learning of variational imaging models[END_REF]), although generalisations to unbounded intervals are also possible [START_REF] De Los Reyes | Bilevel parameter learning for higher-order Total Variation regularisation models[END_REF]. The following proposition says that for problem (4.3) (i.e., with no upper bound constraints on the parameter domain) existence of solutions may fail. Proof. Take f to be any non-constant function in BV(Ω) ∩ L 2 (Ω) with the property that there exists λ * such that

f = argmin u∈BV(Ω) |Du|(Ω) + λ 1 ∥f -u∥ L 1 (Ω)
for all λ 1 ≥ λ * . We note that there are a plethora of such functions, in particular this holds for any one-dimensional function in BV(Ω) ∩ L ∞ (Ω), see also [START_REF] Chan | Aspects of Total Variation regularized L 1 function approximation[END_REF]. Now set ũ = f . Then from Proposition 2.10, we have that

F L 2 (u λ 1 ,λ 2 ) > 0 for every λ 1 , λ 2 ≥ 0. Now fix λ 1 ≥ λ * and let λ (n) 2
→ ∞ as n → ∞. Then according to Corollary 2.4 we have that

u λ 1 ,λ (n) 2 → f weakly * in BV(Ω)
as n → ∞. The function f can in fact be chosen such that this convergence is even uniform, see for instance the example in Proposition 3.1. Then it is clear that

F L 2 (u λ 1 ,λ (n) 2 ) → 0, which means that inf λ 1 ,λ 2 ≥0 F L 2 (u λ 1 ,λ 2 ) = 0.
Since F (u λ 1 ,λ 2 ) > 0 for every λ 1 , λ 2 ≥ 0, we have that for such choices of f and ũ, the bilevel problem (4.3) does not have solutions.

Similarly, we can use the example of Proposition 3.1 to show that the bilevel problem (4.3) with F L 1 γ D cost may not have a solution either, in the case when λ 1 and λ 2 are unbounded.

On the other hand, in the following Proposition we show that the existence of solutions of (4.3) is always guaranteed whenever box constraints on λ 1 , λ 2 are considered.

Proposition 4.3 (Well-posedness of (4.3) with box constraints). The bilevel problem (4.3)

with the extra box constraints 0 ≤ λ i ≤ L i , L i > 0 for i = 1, 2 admits an optimal solution ( λ1 , λ2 ) for both choices of cost functionals (4.7)-(4.8).

Proof. Let (λ

(n) 1 , λ (n) 2 ) ∈ C := {(λ 1 , λ 2 ) : 0 ≤ λ 1 ≤ L 1 , 0 ≤ λ 2 ≤ L 2 } be a minimising sequence for (4.

3). Let us denote by u

n := u λ (n) 1 ,λ (n) 2
the corresponding solution to the lower level problem corresponding to the parameter pair (λ

(n) 1 , λ (n) 2 ).
First suppose that after some index n 0 at least one of the terms (λ

(n) 1 ) n∈N and (λ (n) 2 ) n∈N is zero. This means that u n = c n are constants for n ≥ n 0 . Now if F = F L 2 , due to the coercivity of this functional we have that (c n ) n∈N is bounded, so (u n ) n∈N is bounded in BV(Ω). If F = F L 1 γ then it is obvious that F (u n ) = ∥Dũ∥ M,γ for n ≥ 0.
Thus in this case every constant function trivially solves the bilevel problem (4.3).

We can then assume that λ

(n) 1 , λ (n) 2
> 0 for every n ∈ N. We claim again that the sequence u n is bounded in BV(Ω). Indeed, we have that for every n ∈ N we can bound the TV term as

|Du n |(Ω) ≤ |Du n |(Ω) + Φ λ (n) 1 ,λ (n) 2 (u n , f ) ≤ Φ λ (n) 1 ,λ (n) 2 (0, f ) ≤ L 2 2 ∥f ∥ 2 L 2 (Ω) .
To bound u n in L 1 (Ω), we separate the two cases depending on whether the sequence

λ (n) 1 λ (n) 2 is bounded or not.
If this sequence is bounded by some K > 0, we observe that u n is also a miminiser of

min u∈BV(Ω) Φ 1, λ (n) 2 λ (n) 1 (u, f ) + 1 λ (n) 1 |Du|(Ω), which implies that Φ 1, λ (n) 2 λ (n) 1 (u n , f ) ≤ Φ 1, λ (n) 2 λ (n)
1 (0, f ). Then we have the following successive bounds

∫ Ω |f -u n | dx = ∫ |f -un|< λ (n) 1 λ (n) 2 |f -u n | dx + ∫ |f -un|≥ λ (n) 1 λ (n) 2 λ (n) 1 2λ (n) 2 dx + ∫ |f -un|≥ λ (n) 1 λ (n) 2 |f -u n | - λ (n) 1 2λ (n) 2 dx ≤ K|Ω| + 1 2 ∫ Ω |f -u n | dx + Φ 1, λ (n) 2 λ (n) 1 (u n , f ),
whence we get:

1 2 ∫ Ω |f -u n | dx ≤ K|Ω| + Φ 1, λ (n) 2 λ (n) 1 (0, f )
≤ K|Ω| + ∫ |f |≥ λ (n) 1 λ (n) 2 |f | dx + λ (n) 2 2λ (n) 1 ∫ |f |< λ (n) 1 λ (n) 2 |f | 2 dx ≤ K|Ω| + ∫ Ω |f | dx + K 2 |Ω|.
Suppose now that the sequence

λ (n) 1 λ (n) 2
is unbounded. This means that there exists a (nonrelabelled) subsequence

λ (n) 2 λ (n) 1 → 0. Since λ (n) 1
is bounded, this further implies that λ (n) 2 → 0. Then from Proposition 2.9 we get that u n → f Ω weakly * in BV(Ω) and that, in particular, u n is bounded in L 1 (Ω). So in both cases we have that u n is bounded in BV(Ω).

Having this combined with the boundedness of the sequence (λ

2 ) n∈N implies that we can then extract a further non-relabelled subsequence (λ

(n) 1 , λ (n) 2 , u n ) n∈N converging weakly * in R × R × BV(Ω) to a limit point ( λ1 , λ2 , û).
In particular, this entails that u n → û strongly in L 1 (Ω). Similarly as in the proof of Proposition 2.3 one can now show that the sequence J(•, λ

(n) 1 , λ (n)
2 ) Γ-converges to J(•, λ1 , λ2 ), with respect to the strong topology in L 1 (Ω), which means that û is a minimiser of the lower level problem with parameters ( λ1 , λ2 ). Now, using [START_REF] Dal Maso | Introduction to Γ-convergence[END_REF]Corollary 7.20] we finally have:

Φ λ1 , λ2 (û, f ) + |Dû|(Ω) = lim n→∞ Φ λ (n) 1 ,λ (n) 2 (u n , f ) + |Du n |(Ω),
which completes the proof.

Remark 4.4. Similarly, one can prove the existence of solutions to the regularised bilevel problem (4.4). Note that in [START_REF] De Los Reyes | The structure of optimal parameters for image restoration problems[END_REF] similar techniques are used to prove analogous results for general regularisers and data fidelities. There, the authors proved also the outersemicontinuity property of the solution map S which guarantees that the minimisers of the regularised problem converge towards the minimisers of the one where ϵ = 0.

The quadratic H 1 and L 2 regularisations in (4.4) and (4.6) are required to ensure the differentiability of the solution map, as we are going to highlight in the following. Note that in this section we make use of the formulation (4.4) where the two variables u λ 1 ,λ 2 , v λ 1 ,λ 2 are treated jointly so that the lower-level problem actually reads:

(u λ 1 ,λ 2 , v λ 1 ,λ 2 ) ∈ argmin u∈H 1 (Ω) v∈L 2 (Ω) ϵ 2 ( ∥u∥ 2 H 1 (Ω) + ∥v∥ 2 L 2 (Ω) ) + ∥∇u∥ γ,L 1 (Ω) (4.9)
+ λ 1 ∥v∥ γ,L 1 (Ω) + λ 2 2 ∥f -u -v∥ 2 L 2 (Ω) .
An alternative analysis could exploit the characterisation of the IC data fidelity term given by Proposition 2.1 and consider the minimisation over u only. In such case, we believe that only an H 1 -regularisation on u would be enough for the following proofs. Here, however, we stick with the joint approach to compare our results with the ones derived formally in [START_REF] Calatroni | Infimal convolution of data discrepancies for mixed noise removal[END_REF]. Note that whenever u λ 1 ,λ 2 is given, one can compute the corresponding v λ 1 ,λ 2 by simply solving the optimisation problem (4.6). In particular, the following proposition makes explicit a property of v λ 1 ,λ 2 which will be needed in the following. We refer the reader to [14, Remark 2.1] for a similar characterisation in the case non-regularised case.

Proposition 4.5. Let u, f in L 2 (Ω), 0 ≤ λ i ≤ L i with L i > 0, i = 1, 2 and ϵ > 0. Let v λ 1 ,λ 2 the minimiser of the functional G λ 1 ,λ 2 ϵ,γ (•, u, f ) defined in (4.6
). There holds:

v λ 1 ,λ 2 = v λ 1 ,λ 2 (u) = argmin v∈L 2 (Ω) G λ 1 ,λ 2 ϵ,γ (v, u, f ) = prox λ 1 ϵ+λ 2 ∥•∥ γ,L 1 (Ω) ( λ 2 (f -u) ϵ + λ 2 ) , (4.10)
where for any z ∈ L 2 (Ω), prox τ g (z) denotes the proximal-mapping operator in L 2 (Ω) of the function g with parameter τ . In particular,

v λ 1 ,λ 2 : z → prox λ 1 ϵ+λ 2 ∥•∥ γ,L 1 (Ω)
(z) is a firmly non-expansive operator and it is therefore 1-Lipschitz continuous, i.e.:

∥v λ 1 ,λ 2 (u 1 ) -v λ 1 ,λ 2 (u 2 )∥ L 2 (Ω) ≤ ∥u 1 -u 2 ∥ L 2 (Ω) , for all u 1 , u 2 ∈ L 2 (Ω).
Proof. Straightforward calculations give:

v λ 1 ,λ 2 = argmin v∈L 2 (Ω) ϵ 2 ∥v∥ 2 L 2 (Ω) + λ 1 ∥v∥ γ,L 1 (Ω) + λ 2 2 ∥f -u -v∥ 2 L 2 (Ω) = argmin v∈L 2 (Ω) ∥v∥ γ,L 1 (Ω) + ϵ + λ 2 2λ 1 ∥v∥ 2 L 2 (Ω) + 1 2 λ 1 λ 2 ( 1 + ϵ λ 2 ) ∥f -u∥ 2 L 2 (Ω) - λ 2 λ 1 ∫ Ω (f -u)v dx = argmin v∈L 2 (Ω) ∥v∥ γ,L 1 (Ω) + 1 2 λ 1 ϵ+λ 2 λ 2 (f -u) ϵ + λ 2 -v 2 L 2 (Ω) = prox λ 1 ϵ+λ 2 ∥•∥ γ,L 1 (Ω) ( λ 2 (f -u) ϵ + λ 2
) .

The firm non-expansiveness property (4.10) follows then directly from [6, Proposition 12.27].

Optimality system

We now study in more detail the bilevel problem (4.4) with (4.9) as lower level problem and prove the existence of Lagrange multipliers by deriving the optimality system characterising its stationary points. As a by product, we find an easy formula to compute the gradient of the cost functional in terms of its adjoint state, which simplifies the design of the gradient-based algorithm which we are going to use to solve (4.4) in an efficient way. We start defining the Hilbert space H := H 1 (Ω)×L 2 (Ω) endowed with the scalar product (z, w) H := (z 1 , w 1 ) H 1 (Ω) + (z 2 , w 2 ) L 2 (Ω) for all z, w ∈ H. We further denote by S : R 2 → H the solution map S : (λ 1 , λ 2 ) → (u λ 1 ,λ 2 , v λ 1 ,λ 2 ) which assigns to the optimal parameters (λ 1 , λ 2 ) the corresponding solution pair of (4.9). To avoid heavy notations, we will omit the explicit dependence of the pair (u, v) on the parameters (λ 1 , λ 2 ) unless explicitly needed. Note that we take the whole space R 2 as differentiability set although from an imaging point of view the use of negative parameters clearly does not make sense. However, the following differentiability result holds in such case as well.

Recalling the definition of h γ in (4.5), we have that in correspondence with an optimal pair y := (u, v) ∈ H, we have that the following variational equality holds true for all test functions Ψ := (ψ 1 , ψ 2 ) ∈ H:

ϵ (y, Ψ) H + ∫ Ω h γ (∇u)∇ψ 1 dx + λ 1 ∫ Ω h γ (v)ψ 2 dx + λ 2 ∫ Ω (u + v -f )(ψ 1 + ψ 2 ) dx = 0.
ϵ (z, Ψ) H + ∫ Ω h ′ γ (∇u)∇z 1 ∇ψ 1 dx + λ 1 ∫ Ω h ′ γ (v)z 2 ψ 2 dx + θ 1 ∫ Ω h γ (v)ψ 2 dx + λ 2 ∫ Ω (z 1 + z 2 )(ψ 1 + ψ 2 ) dx + θ 2 ∫ Ω (u + v -f )(ψ 1 + ψ 2 ) dx = 0, (4.12)
for all Ψ = (ψ 1 , ψ 2 ) ∈ H.

Proof. Thanks to the ellipticity of the scalar product in H and the monotonicity of h γ , existence and uniqueness of z of (4.12) are guaranteed by Lax-Milgram theorem. Now, we want to show that z is the Fréchet derivative of S. To do that, given θ = (θ 1 , θ 2 ) ∈ R 2 let us define y + = (u + , v + ) := S(λ 1 + θ 1 , λ 2 + θ 2 ). We aim to show that ξ = (ξ 1 , ξ 2 ) := y + -y -z ∈ H satisfies ∥ξ∥ H = o(|θ|).

Writing (4.11) for y and y + , and (4.12) for z and combining them together, we get that for every Ψ ∈ H there holds:

ϵ (ξ, Ψ) H + ∫ Ω ( h γ (∇u + ) -h γ (∇u) ) ∇ψ 1 dx - ∫ Ω h ′ γ (∇u)∇z 1 ∇ψ 1 dx + λ 1 ∫ Ω ( h γ (v + ) -h γ (v) ) ψ 2 dx + θ 1 ∫ Ω ( h γ (v + ) -h γ (v) ) ψ 2 dx -λ 1 ∫ Ω h ′ γ (v)z 2 ψ 2 dx + λ 2 ∫ Ω (ξ 1 + ξ 2 )(ψ 1 + ψ 2 ) dx + θ 2 ∫ Ω ( (u + -u) + (v + -v) ) (ψ 1 + ψ 2 ) dx = 0.
We now add and subtract the terms:

∫ Ω h ′ γ (∇u) ( ∇(u + -u) ) ∇ψ 1 dx, and λ 1 ∫ Ω h ′ γ (v) ( v + -v ) ψ 2 dx,
thus getting:

ϵ (ξ, Ψ) H + ∫ Ω h ′ γ (∇u)∇ξ 1 ∇ψ 1 dx + λ 1 ∫ Ω h ′ γ (v)ξ 2 ψ 2 dx + λ 2 ∫ Ω (ξ 1 + ξ 2 )(ψ 1 + ψ 2 ) dx = - ∫ Ω ( h γ (∇u + ) -h γ (∇u) -h ′ γ (∇u)(∇(u + -u)) ) ∇ψ 1 dx -λ 1 ∫ Ω ( h γ (v + ) -h γ (v) -h ′ γ (v)(v + -v) ) ∇ψ 2 dx -θ 1 ∫ Ω ( h γ (v + ) -h γ (v) ) ψ 2 dx -θ 2 ∫ Ω ( (u + -u) + (v + -v) ) (ψ 1 + ψ 2 ) dx, for all Ψ ∈ H.
We now choose Ψ = ξ. Under this choice, by monotonicity of h ′ γ we have that the last three terms on the left-hand side become non-negative. We further deduce:

∥ξ∥ H ≤ C ( ∥h γ (∇u + ) -h γ (∇u) -h ′ γ (∇u)(∇(u + -u))∥ L 2 (Ω) + |λ 1 | ∥h γ (v + ) -h γ (v) -h ′ γ (v)(v + -v)∥ L 2 (Ω) + |θ 1 | ∥h γ (v + ) -h γ (v)∥ L 2 (Ω) + |θ 2 | ∥u + -u∥ L 2 (Ω) + |θ 2 | ∥v + -v∥ L 2 (Ω)
)

.

where C is a generic positive and finite constant which may change from line to line. By the differentiability and Lipschitz continuity of h γ and h ′ γ , we have:

∥ξ∥ H ≤ C ( o ( ∥∇u + -∇u∥ L 2 (Ω) ) + |λ 1 | o ( ∥v + -v∥ L 2 (Ω) ) (4.13)
+ |θ 2 | ∥u + -u∥ L 2 (Ω) + (|θ 1 | + |θ 2 |) ∥v + -v∥ L 2 (Ω)
) .

We now focus on the terms depending on the difference between

v + = v λ 1 +θ 1 ,λ 2 +θ 2 (u + ) and v = v λ 1 ,λ 2 (u)
. By triangle inequality we have:

∥v + -v∥ L 2 (Ω) ≤ ∥v λ 1 +θ 1 ,λ 2 +θ 2 (u + ) -v λ 1 +θ 1 ,λ 2 +θ 2 (u)∥ L 2 (Ω) + ∥v λ 1 +θ 1 ,λ 2 +θ 2 (u) -v λ 1 ,λ 2 (u)∥ L 2 (Ω) ≤ ∥u + -u∥ L 2 (Ω) + |θ|,
where the last inequality follows from Proposition 4.5 and from the continuity property of the proximal mapping with respect to its parameter which can be easily checked in our case recalling that |h γ (•)| ≤ 1. We then deduce in (4.13) that:

∥ξ∥ H ≤ C ( o(|θ|) + o(∥u + -u∥ H 1 (Ω) ) ) (4.14)
Thus, to conclude it only remains to show that o(∥u

+ -u∥ H 1 (Ω) ) = o(|θ|).
To do that, we use standard Sobolev embeddings and the regularity result of Gröger for second-order systems [29, Theorem 1] and get that for p > 2:

∥u + -u∥ H 1 (Ω) ≤ ∥u + -u∥ W 1,p (Ω) ≤ C|θ| ( ∥div(h γ (∇u))∥ W -1,p (Ω) + ∥h γ (∇u)∥ W -1,p (Ω) ) ≤ C|θ|∥h γ (∇u)∥ L ∞ (Ω) ≤ C|θ|, since |h γ (•)| ≤ 1.
Combining with (4.14) this finishes the proof.

Remark 4.7. The regularity result by Gröger is a classical argument for the proof of Fréchet differentiability in similar bilevel problems (see, e.g., [START_REF] Calatroni | Bilevel approaches for learning of variational imaging models[END_REF]). Note, however, that the original result in [START_REF] Gröger | A W 1,p -estimate for solutions to mixed boundary value problems for second order elliptic differential equations[END_REF] was proved for C 2 -regular domains, while its extension to convex Lipschitz domains (such as image domains) has been proved by Dauge in [START_REF] Dauge | Neumann and mixed problems on curvilinear polyhedra[END_REF].

We now prove the existence and uniqueness of the adjoint state of the problem (4.4).

ϵ(Π, W ) H + ∫ Ω h ′ γ (∇u)∇w 1 ∇p 1 dx + λ 1 ∫ Ω h ′ γ (v)w 2 p 2 dx + λ 2 ∫ Ω (p 1 + p 2 )w 1 dx + λ 2 ∫ Ω (p 1 + p 2 )w 2 dx = - ∫ Ω F ′ (u)w 1 dx, (4.15)
for any W := (w 1 , w 2 ) ∈ H.

Proof. For w ∈ H, let us consider the following bilinear form on H × H:

a(Π, W ) := ϵ(Π, W ) H + ∫ Ω h ′ γ (∇u)∇w 1 ∇p 1 dx + λ 1 ∫ Ω h ′ γ (v)w 2 p 2 dx + λ 2 ∫ Ω (p 1 + p 2 )w 1 dx + λ 2 ∫ Ω (p 1 + p 2 )w 2 dx.
The form a(•, •) is trivially symmetric and coercive, since by taking W = Π, we get:

a(Π, Π) = ϵ∥Π∥ 2 H + ∫ Ω h ′ γ (∇u)∇p 1 ∇p 1 + λ 1 ∫ Ω h ′ γ (v)p 2 p 2 dx + λ ∫ Ω (p 1 + p 2 ) 2 dx ≥ ϵ∥Π∥ 2 H .
by monotonicity of h ′ γ . By Lax-Milgram theorem, we infer that there exists a unique solution of (4.15). Note that by taking w 2 = 0 in (4.15), we get the optimality condition for p 1 , i.e.:

ϵ(p 1 , w 1 ) H 1 (Ω) + ∫ Ω h ′ γ (∇u)∇w 1 ∇p 1 dx + λ 2 ∫ Ω (p 1 + p 2 )w 1 dx = - ∫ Ω F ′ (u)w 1 dx,
for any w 1 ∈ H 1 (Ω). Similarly, for w 1 = 0, we get the optimality condition for p 2 :

ϵ ∫ Ω w 2 p 2 dx + λ 1 ∫ Ω h ′ γ (v)w 2 p 2 dx + λ 2 ∫ Ω (p 1 + p 2 )w 2 dx = 0, for any w 2 ∈ L 2 (Ω).
Finally, we now combine the results above to derive the optimality system of the bilevel problem (4.4). We recall that by y = (u, v) ∈ H we denote the solution pair. Theorem 4.9 (Optimality system). Let ( λ1 , λ2 ) ∈ R + × R + an optimal solution of the problem (4.4). Then, there exists a Lagrange multiplier Π := (p 1 , p 2 ) ∈ H and µ 1 , µ 2 ∈ R + such that the following system holds:

ϵ (y, Ψ) H + ∫ Ω h γ (∇u)∇ψ 1 dx + λ 1 ∫ Ω h γ (v)ψ 2 dx + λ 2 ∫ Ω (u + v -f )(ψ 1 + ψ 2 ) dx = 0, for all Ψ = (ψ 1 , ψ 2 ) ∈ H, ϵ(Π, W ) H + ∫ Ω h ′ γ (∇u)∇w 1 ∇p 1 dx + λ 1 ∫ Ω h ′ γ (v)w 2 p 2 dx (4.16) + λ 2 ∫ Ω (p 1 + p 2 )(w 1 + w 2 ) dx = - ∫ Ω F ′ (u)w 1 dx, for all W = (w 1 , w 2 ) ∈ H, µ 1 := ∫ Ω h ′ γ (v)p 2 dx, µ 2 := ∫ Ω (f -v -u)(p 1 + p 2 ) dx, µ 1 ≥ 0, µ 2 ≥ 0, µ 1 • λ1 = µ 2 • λ2 = 0.
Proof. We can write the bilevel problem (4.4) in a reduced form as:

min λ 1 ,λ 2 ≥0 F(λ 1 , λ 2 ) := F (u λ 1 ,λ 2 ).
Using [57, Theorem 3.1], we deduce the existence of multipliers µ 1 , µ 2 ∈ R + such that:

µ 1 = ∇ λ 1 F( λ1 , λ2 ), µ 2 = ∇ λ 2 F( λ1 , λ2 ), µ 1 ≥ 0, µ 2 ≥ 0, µ 1 • λ1 = µ 2 • λ2 = 0.
Computing the gradient of F by using the chain rule we get:

∇F(λ 1 , λ 2 )[θ 1 , θ 2 ] = ( F ′ (u λ 1 ,λ 2 ), S ′ (λ 1 , λ 2 )[θ 1 , θ 2 ] ) L 2 (Ω) = ∫ Ω F ′ (u λ 1 ,λ 2 )z dx
where z ∈ H is the linearised state provided by Theorem 4.6. Theorem 4.8 ensures that there exist a Lagrange multiplier Π := (p 1 , p 2 ) satisfying the adjoint equation (4.15), which entails:

∫ Ω F ′ (u λ 1 ,λ 2 )z dx = -ϵ(Π, z) H - ∫ Ω h ′ γ (∇u)∇z 1 ∇p 1 dx -λ 1 ∫ Ω h ′ γ (v)z 2 p 2 dx -λ 2 ∫ Ω (p 1 + p 2 )(z 1 + z 2 ) dx = θ 1 ∫ Ω h γ (v)p 2 dx + θ 2 ∫ Ω (u + v -f )(p 1 + p 2 ) dx,
which completes the proof.

Remark 4.10. Theorem 4.9 provides also the following handy formula for the computation of the gradient of the bilevel problem (4.4) with data fidelity (4.6) in a reduced form:

∇F(λ 1 , λ 2 )[θ 1 , θ 2 ] = θ 1 ∫ Ω h γ (v)p 2 dx + θ 2 ∫ Ω (u + v -f )(p 1 + p 2 ) dx. ( 4
.17)

Remark 4.11. Compared to the optimality system derived in [START_REF] Calatroni | Infimal convolution of data discrepancies for mixed noise removal[END_REF]Section 7] via Lagrangian formalism, we note that the optimality system (4.16) presents an additional quadratic ϵregularisation on v. This is needed in the proof of Theorem 4.6 to get a uniform, finite estimate for ξ 2 not depending on λ 2 which may be zero.

Numerical experiments

In this section we report some numerical results on the computation of the optimal parameters ( λ1 , λ2 ) of the TV-IC model by means of the bilevel optimisation strategy described in the previous section. For given training images with a mixture of Gaussian and Salt & Pepper noise with various intensities, we describe in Algorithm 1 the main steps to compute ( λ1 , λ2 ) by solving the optimality system (4.16) via a second-order BFGS optimisation approach. A general review of second-order numerical methods for PDE-constrained optimisation models can be found in the book [START_REF] De Los Reyes | Numerical PDE-constrained optimization[END_REF], while more details on the numerical realisation of similar bilevel models can be found in [START_REF] Calatroni | Bilevel approaches for learning of variational imaging models[END_REF]Section 8.3] and in [26, Section 4.1].

In the following numerical computations:

-We consider test images of size N × N ≡ 256 × 256 pixels. The differential operators are discretised using finite difference schemes with mesh step size h = 1/N . Standard forward/backward differences are used for the discretisation of the divergence/gradient operator, respectively.

-For illustrative purposes, we report the results obtained using a training set consisting of one training pair (ũ, f ) only. We recall that f represents the noisy version of ũ corrupted by a mixture of Gaussian and Salt & Pepper noise of various intensities which we will specify in each case. The efficient extension to multiple constraints can be done similarly as in [START_REF] Calatroni | Bilevel approaches for learning of variational imaging models[END_REF][START_REF] Calatroni | Dynamic sampling schemes for optimal noise learning under multiple nonsmooth constraints[END_REF].

-The lower-level regularised TV-IC problem (4.4) is solved by means of the SemiSmooth Newton (SSN) algorithm described in [START_REF] Calatroni | Infimal convolution of data discrepancies for mixed noise removal[END_REF] with a warm start. The regularisation parameters are chosen as ϵ = 10 -10 and γ = 10 3 . The algorithm is stopped if either the difference between two consecutive iterates is below tol = 10 -6 or if the maximum number of iterations maxiter = 35 is reached.

-For the outer BFGS iterations, an Armijo line-search with parameter η = 10 -4 is employed together with a curvature verification. The Armijo rule:

F(λ k + α k d k ) -F(λ k ) ≤ ηα k ∇F(λ k ) ⊤ d k is checked at any iteration k ≥ 2 Here λ k = (λ 1 k , λ k 2 )
stands for the parameter pair updated along the iterations, d k for the quasi-Newton descent direction and α k for the line step. The expression of ∇F is given in Remark 4.10. The outer algorithm is stopped when the maximum between the norm of the gradient of the cost functional and the difference of two subsequent iterates is smaller than tol 1 = 10 -6 . model. It is then interesting to notice that in the case θ = 1 the bilevel optimisation strategy does not enforce a TV-L 2 model (i.e., a large λ1 ), but rather a combination of TV-L 1 and TV-L 2 . This might be an indication that in practice the TV-L 1 model works well also in the case of pure Gaussian noise removal. This is reflected in our experiment where the estimated optimal data model turns out to be indeed a combination of the two discrepancies. Note that the two parameters λ 1 and λ 2 scale differently, with λ1 only slightly varying across the different simulations. We report in Table 2 the numerical values of ( λ1 , λ2 ) and the corresponding PSNR values of the noisy images fθ and of the optimal reconstructions uλ 1 , λ2 .

Conclusions and outlook

In this paper we have presented a fine analysis of the TV-IC denoising model originally proposed in [START_REF] Calatroni | Infimal convolution of data discrepancies for mixed noise removal[END_REF] for mixed Gaussian and Salt & Pepper noise removal. Our study started with a characterisation of the IC data-discrepancy as a Huber-regularised L 1 discrepancy. We then studied in detail the asymptotic behaviour of the solutions of the TV-IC model using Γ-convergence arguments, showing that the solutions of the single TV-L 1 and TV-L 2 denoising models can be retrieved in the limit as the parameters tend to infinity. We gained more insights on the model by calculating some exact solutions for simple one-dimensional data functions. Using these theoretical results we then formulated and rigorously analysed a bilevel optimisation approach in function spaces for the estimation of the optimal parameters of the model. With the use of a counterexample motivated by our theoretical work, we showed that box constraints on the parameters are indeed necessary to obtain existence.

In the spirit of [START_REF] Calatroni | Bilevel approaches for learning of variational imaging models[END_REF][START_REF] De Los Reyes | Bilevel parameter learning for higher-order Total Variation regularisation models[END_REF], after a suitable regularisation of the non-smooth problem, we then proved the differentiability of the solution map via standard arguments as well as the existence of the adjoint states by which we can derive the corresponding optimality system in a compact form. Thanks to the handy characterisation of the gradient of the bilevel cost functional in terms of its adjoint states, efficient numerical schemes can be easily implemented. In particular, in Section 5 we considered the second-order BFGS Algorithm 1 to numerically compute the optimal parameters of the TV-IC model for noise mixtures with different noise levels. The numerical results show good agreement with the analytical study of the first sections, making this strategy appealing for blind image denoising applications, where the intensity of each noise component is unknown.

Further research could address the validation of the TV-IC model over a set of images in order to estimate the preferred noise model with respect to the denoising application at hand (i.e., the noise intensity, the structure of the image etc.), similarly as done in [START_REF] De Los Reyes | Bilevel parameter learning for higher-order Total Variation regularisation models[END_REF]. The IC discrepancy could further be combined with higher-order regularisers (such as TGV [START_REF] Bredies | Total Generalized Variation[END_REF]) for more visually pleasing reconstructions. Finally, a similar study could be done for more general noise mixtures such as Gaussian & Poisson noise, for which the IC discrepancy has been shown in [START_REF] Calatroni | Infimal convolution of data discrepancies for mixed noise removal[END_REF] to be a statistically consistent model. [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF]. LC acknowledges the joint ANR/FWF Project "Efficient Algorithms for Nonsmooth Optimization in Imaging" (EANOI) FWF n. I1148 / ANR-12-IS01-0003, the Fondation Mathématique Jacques Hadamard (FMJH) and the RISE EU project NoMADS. KP acknowledges the support of the Einstein Foundation Berlin within the ECMath project CH12. Both authors would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the program Methods and Effective Algorithms for Imaging and Vision" when work on this paper was undertaken. This work was supported by EPSRC grant number EP/K032208/1.

  Fixed λ 1 = 1 and varying λ 2 . Fixed λ 2 = 1 and varying λ 1 .

Figure 1 :

 1 Figure 1: Example plots of the Huber function φ.

  |f (x) -u(x)| (d) |f (x) -u(x)| along the line profile

Figure 2 :

 2 Figure 2: Interepretation of the L 1 -L 2 IC term based on the formulation (2.1)-(2.2). First row: original image, noisy version and detail with line profile (red). Second row: difference |f -u| along the line profile. The fidelity functional Φ λ 1 ,λ 2 locally acts like ∥f -u∥ L 1 at points with high values of |f -u| (blue asterisks), thus assuming these points to be corrupted by Salt & Pepper noise. On the other hand, Φ λ 1 ,λ 2 acts locally as ∥f -u∥ 2 L 2 at points with low values of |f -u| (red asterisks), identifying these points as corrupted by Gaussian noise. Parameters: Gaussian variance σ 2 = 0.005, density of pixels corrupted by Salt & Pepper noise d = 5%.

  Figure3: If λ 1 λ 2 ≥ 2∥f ∥ ∞ then (2.5) is equivalent to the TV-L 2 problem, see Proposition 2.6. By fixing λ 1 and sending λ 2 to infinity the solution u converges to a solution of an TV-L 1 in the sense of Corollary 2.4.

Proposition 3 . 1 .

 31 Let Ω = (-2L, 2L) and f ∈ BV(Ω) being the jump function given by (3.1). Then the solutions u to the TV-IC minimisation problem min u∈BV(Ω)

Figure 4 :

 4 Figure 4: Visualisation of all the possible solutions to the TV-IC problem (2.5) for data (3.1) and for all the possible combinations of the parameters λ 1 and λ 2 , see Proposition 3.1.

Proposition 4 . 2 .

 42 There exist data f ∈ L 2 (Ω) and training data ũ ∈ L 2 (Ω) such that the non-smooth bilevel problem (4.3) does not have a solution for the cost function F L 2 .

( 4 . 11 )

 411 We now prove the main differentiability result.

Theorem 4 . 6 (

 46 Fréchet differentiability of the solution map). The solution operator S : R 2 → H which assigns to each parameter pair (λ 1 , λ 2 ) the element y := (u, v) = S(λ 1 , λ 2 ), solution of the TV-IC denoising problem (4.4)-(4.6) is Fréchet differentiable. In particular, for any θ = (θ 1 , θ 2 ) ∈ R 2 its Fréchet derivative is the unique solution z := S ′ (λ 1 , λ 2 )[θ 1 , θ 2 ] ∈ H of the following linearised equation:

Theorem 4 . 8 (

 48 Adjoint equation). Let (u, v) ∈ H.There exists a unique solution Π := (p 1 , p 2 ) ∈ H to the adjoint PDE:

  f (c) uλ 1 , λ2 , F L 2 cost. (d) uλ 1 , λ2 , F L 1 γ D cost.

Figure 5 :

 5 Figure 5: Optimal TV-IC denoising results for initial guess (λ 0 1 , λ 0 2 ) = (1, 10) w.r.t. F L 2 and F L 1 γ D costs (4.7)-(4.8). Noisy image corrupted with Gaussian noise with variance σ 2 = 0.01 and percentage of missing pixels d = 10%. 5b Noisy image f : PSNR=14.40 dB, SSIM=0.17. 5c Optimal uλ 1 , λ2 w.r.t. F L 2 cost: PSNR=28.35 dB, SSIM=0.81. 5d Optimal uλ 1 , λ2 w.r.t. F L 1 γ D cost: PSNR=27.91 dB, SSIM=0.83.

Figure 6 :

 6 Figure 6: Optimal parameters ( λ1 , λ2 ) computed by solving the bilevel system (4.16) with L 2 cost and training pair (ũ, fθ ) depending on a parameter θ ∈ [0, 1] which controls the amount of Gaussian and Salt & Pepper noise in the data, see Figure 7. For θ = 0 only Salt & Pepper noise is present and the model computes optimal parameters enforcing TV-L 1 denoising model. For 0 < θ < 1 the data discrepancies are weighted by smaller parameters λ1 and λ2 . For each θ the bilevel Algorithm 1 is initialised with (λ 0 1 , λ 0 2 ) = (1, 1).

  fθ , θ = 0 (S & P noise). (c) fθ , θ = 0.25. (d) fθ , θ = 0.5. (e) fθ , θ = 0.75. (f) fθ , θ = 1 (Gaussian noise).

Figure 7 :

 7 Figure 7: Training images used for the computation of the optimal parameters ( λ1 , λ2 ) whose values are reported in Figure 6; The noise-free image ũ is corrupted with a mixture of Gaussian and Salt & Pepper noise of different intensities, ranging from pure impulsive noise when θ = 0 to pure Gaussian when θ = 1.

  .21) By combining (2.19), (2.20), (2.21) with (2.18) we get that |Df |(Ω) = 0, and since Ω is connected, f is a constant function. Note that the calculations above assumed that d > 1

Table 2

 2

		PSNR fθ	λ1	λ2	PSNR uλ 1 , λ2
	0	15.41 dB 1.95 123.39	29.71 dB
	0.25 16.19 dB 2.15 61.04	24.66 dB
	0.5	17.49 dB 2.31 39.94	25.35 dB
	0.75 19.28 dB 2.51 45.90	24.01 dB
	1	22.74 dB 2.47 57.35	24.25 dB

Acknowledgments

The authors thank Prof. Juan Carlos De Los Reyes for his useful comments on Section

-The adjoint equation is solved by means of standard sparse linear solvers.

-Our validation images are taken from the public domain.

Algorithm 1 Bilevel optimisation algorithm for computing optimal λ 1 and λ 2 in (4.4)

Output: Optimal parameters λ1 and λ2 . Initialise: λ 0 1 , λ 0 2 , n = 1. while not converging do SSN algorithm to compute (u n , v n ) by solving (4.11) with parameters (λ In Figure 5 we report a numerical experiment confirming the effectiveness of the bilevel optimisation approach on images corrupted by a mixture of Gaussian and Salt & Pepper noise with Gaussian variance σ 2 = 0.01 and percentage of missing pixels d = 10%. We report the result obtained with respect to both the L 2 cost functional (4.7) and the Huberised L 1 gradient cost (4.8). As observed in [START_REF] De Los Reyes | Bilevel parameter learning for higher-order Total Variation regularisation models[END_REF], we remark that minimising with respect to the L 2 cost is indeed equivalent to PSNR optimisation, while the optimisation with respect to the Huberised L 1 gradient cost produces better visual results which is similar to optimising the SSIM.

In order to validate numerically the theoretical insights given by the analysis performed in Section 2, we report in the following Figure 6 a plot of the optimal parameters ( λ1 , λ2 ) computed by solving the bilevel system (4.16) with L 2 cost (4.7) and in correspondence of a training pair (ũ, fθ ) where the noisy image is corrupted by a mixture of Gaussian and Salt & Pepper noise with varying intensity. Namely, for θ ∈ [0, 1], we corrupt the training image ũ in 7a with Gaussian noise with distribution N (0, θσ 2 ), σ 2 = 0.005 and Salt & Pepper noise with a percentage of corrupted pixels equal to d = (1 -θ)10%. Consequently, the noisy image f θ is corrupted by pure impulsive noise for θ = 0 (Figure 7b), by pure Gaussian noise for θ = 1 (Figure 7f) and by a mixture of the two for θ ∈ (0, 1) (Figures 7c-7e). As suggested by the theory (see, in particular, Proposition 2.3 and Corollary 2.4), we observe that when θ = 0 the bilevel strategy selects a large optimal parameter λ2 , enforcing a TV-L 1 denoising model which is well known to be optimal for this type of noise (see, e.g., [START_REF] Duval | The TV-L1 model: a geometric point of view[END_REF][START_REF] Nikolova | A variational approach to remove outliers and impulse noise[END_REF]). Furthermore, by Corollary 2.4 we also have that the lower-level solution of the mixed noise problem approximates (in the sense of Γ-convergence) a solution of the corresponding TV-L 1