# Dynamic Sampling Schemes for Optimal Noise Learning Under Multiple Nonsmooth Constraints

Luca Calatroni, Juan De Los Reyes, Carola-Bibiane Schönlieb

► To cite this version:

Luca Calatroni, Juan De Los Reyes, Carola-Bibiane Schönlieb. Dynamic Sampling Schemes for Optimal Noise Learning Under Multiple Nonsmooth Constraints. 26th Conference on System Modeling and Optimization (CSMO), Sep 2013, Klagenfurt, Austria. pp.85-95, 10.1007/978-3-662-45504-3_8 hal-01286460

HAL Id: hal-01286460

https://inria.hal.science/hal-01286460

Submitted on 10 Mar 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamic sampling schemes for optimal noise learning under multiple nonsmooth constraints*

Luca Calatroni<fnref n="1" />, Juan Carlos De Los Reyes<fnref n="2" />, and Carola-Bibiane Schönlieb<fnref n="3" />

**Abstract.** We consider the bilevel optimisation approach proposed in [5] for learning the optimal parameters in a Total Variation (TV) denoising model featuring for multiple noise distributions. In applications, the use of databases (dictionaries) allows an accurate estimation of the parameters, but reflects in high computational costs due to the size of the databases and to the nonsmooth nature of the PDE constraints. To overcome this computational barrier we propose an optimisation algorithm that, by sampling *dynamically* from the set of constraints and using a quasi-Newton method, solves the problem accurately and in an efficient way.

## 1 Introduction

Most images in the real world suffer from noise. In photography noisy images occur when taking a photograph under bad lighting conditions, for instance. Medical imaging applications, such as Magnetic Resonance Imaging (MRI) and Positron Electron Tomography (PET), produce under-sampled and noisy image data. In general, the quality of images obtained from imaging devices in the real world, in the sciences and medicine, is limited by the hardware and the limited time available to measure the image data. Hence, one of the most important tasks in image processing is the reduction of noise in images, called image denoising.

A common challenge in image denoising is the setup of a suitable denoising model. The model depends on the noise distribution and the class of images the denoised solution should belong to. In [5] a bilevel optimisation approach to learn the correct setup for a TV denoising model from a set of noisy and clean test images is proposed. There, optimal parameters $\lambda_i \in \mathbb{R}, i = 1, \dots, d$ are determined by solving the following optimisation problem:

$$ \min_{\lambda_i \ge 0, \ i=1, \dots, d} \frac{1}{2N} \sum_{k=1}^{N} \| \hat{u}_k - u_k \|_{L^2(\Omega)}^2 \quad (1.1) $$

<footnote n="1">¹ Cambridge Centre for Analysis, University of Cambridge, UK</footnote>

<footnote n="2">² Centro De Modelización Matemática, Escuela Politécnica Nacional de Quito, Ecuador</footnote>

<footnote n="3">³ Department of Applied Mathematics and Theoretical Physics, University of Cambridge, UK</footnote>

<footnote>* This work is supported by the King Abdullah University for Science and Technology (KAUST) Award No. KUK-I1-007-43, the EPSRC first grant Nr. EP/J009539/1 and the Cambridge Center for Analysis (CCA).</footnote>
subject to the set of nonsmooth constraints:

$$ \hat{u}_k = \operatorname*{argmin}_{u \in BV(\Omega) \cap A} \left( |Du|(\Omega) + \sum_{i=1}^{d} \lambda_i \int_{\Omega} \phi_i(u, f_k) dx \right), \quad k = 1, \dots, N. \tag{1.2} $$

In (1.1)-(1.2) $\Omega \subset \mathbb{R}^2$ is the image domain, $|Du|(\Omega)$ is the Total Variation (TV) $u$ in $\Omega$ and $BV(\Omega)$ is the space of functions of bounded variation (see [1]). For each $k$, the pair $(u_k, f_k)$ is an element of a set of $N$ pairs of clean and noisy test images, respectively, whereas $\hat{u}_k$ is the TV-denoised version of $f_k$. For $i = 1, \dots, d$ the terms $\phi_i$ represent the different data fidelities, each one modelling one particular type of noise weighted by a parameter $\lambda_i$. Examples of $\phi$ are $\phi(u, f_k) = (u - f_k)^2$ for noise with Gaussian distribution and $\phi(u, f_k) = |u - f_k|$ for the case of impulse noise. The set $A$ is the set of admissible functions such that the data fidelity terms are well defined.

In this paper, we use a simulated database of clean and noisy images. This is not uncommon. Even in real world applications such as MRI, simulated databases are used to tune image retrieval systems, see for instance [7]. Alternatively, we can imagine the retrieval of such a test set for a specific application using phantoms and their noisy acquisitions. Ideally, we would like to consider a very rich database (i.e. $N \gg 1$) in order to get a more robust estimation of the parameters, thus dealing with a very large set of constraints (1.2) that would need to be solved in each iteration of an optimisation algorithm applied to (1.1)-(1.2). The computational solution of such an optimisation problem renders expensive and therefore challenging due to the large-scale nature of the problem (1.1)-(1.2) and due to the nonsmooth nature of each constraint.

In order to deal with such large-scale problems various stochastic optimisation approaches have been presented in literature. They are based on the common idea of solving not all the PDE constraints, but just a sample of them, whose size varies according to the approach one intends to use. In this paper we focus on a stochastic approximation method proposed by Byrd et al. [3] called *dynamic sample size* method. The main idea of this method is to consider an initial, small, training sample of the dictionary to start the algorithm with and *dynamically* increasing its size, if needed, throughout the different steps of the optimisation process. The criterion to decide whether or not the sample size has to be increased is a check on the sample variance estimates on the batch gradient. The desired trade-off between efficiency and accuracy is then obtained by starting with a small sample and gradually increasing its size till reaching the requested level of accuracy. Let us mention that the method of Byrd et al. is one among various stochastic optimisation methods, compare for instance [8,10,2,6,9].

Our work extends the work of [3] in two directions: firstly, in [3] the linearity of the solution map is required which is not fulfilled for our problem (1.1)-(1.2). We are going to show that the strategy of Byrd et al. can be modified for nonlinear solution maps as the one we are considering. Secondly, in [3] the optimization algorithm is of gradient-descent type. Using a BFGS method to solve (1.1)-(1.2) we extend their approach incorporating also *second order* information in form
of an approximation of the Hessian by evaluations of the sample gradient in the iterations of the optimisation algorithm.

**Organisation of the paper.** In the following Section 2 we present the Dynamic Sampling algorithm adapted to the nonlinear framework of problem (1.1)-(1.2), specifying the variance condition on the batch gradient used in our optimisation algorithm. In Section 3 we present the numerical results obtained for the estimation of the optimal parameters in the case of single and mixed noise estimation for the model (1.1)-(1.2) showing significant improvements in efficiency.

**Preliminaries.** We denote the vector of parameters we aim to estimate by $\lambda = (\lambda_1, \dots, \lambda_d) \in \mathbb{R}_{\ge 0}^d$. We also define by $\mathcal{S}$ the solution map that, for each constraint $k = 1, \dots, N$ of (1.2), associates to $\lambda$ and to the noisy image $f_k$ the corresponding Total Variation denoised solution $\hat{u}_k$, that is $\mathcal{S}(\lambda, f_k) = \hat{u}_k$. Let us then define the reduced cost functional $J(\lambda)$ as

$$J(\lambda) := \frac{1}{2N} \sum_{k=1}^{N} \|\mathcal{S}(\lambda, f_k) - u_k\|_{L^2(\Omega)}^2 \quad (1.3)$$

We also define:

$$l(\lambda, f_k) := \|\mathcal{S}(\lambda, f_k) - u_k\|_{L^2(\Omega)}^2, \quad k = 1, \dots, N \quad (1.4)$$

as the loss functions of the functional $J$ defined in (1.3) for each $k = 1, \dots, N$. For every sample $S \subset \{1, \dots, N\}$ of the database, we introduce the batch objective function:

$$J_S(\lambda) := \frac{1}{2|S|} \sum_{k \in S} l(\lambda, f_k) \quad (1.5)$$

## 2 Dynamic sampling schemes for solving (1.1)-(1.2)

To design the optimisation algorithm solving (1.1)-(1.2) we follow the approach used in [5]. There, a quasi-Newton method (namely, the Broyden-Fletcher-Goldfarb-Shanno algorithm BFGS) is considered together with an Armijo backtracking linesearch rule. We combine such algorithm with a modified version of the Dynamic Sampling algorithm presented in [3, Section 3]. In order to compare our algorithm with the Newton-Conjugate Gradient method presented in [3, Section 5], we highlight that in our optimisation algorithm the Hessian matrix is never computed, but approximated efficiently by the BFGS matrix.

Our algorithm starts by selecting from the whole dataset a sample $S$ whose size $|S|$ is small compared to the original size $N$. In the following iterations, if the approximation computed produces an improvement in the cost functional $J$, then the sample size is kept unchanged and the optimisation process continues selecting in the next iteration a new sample of the same size. Otherwise, if the approximation computed is not a good one, a new, larger, sample size is selected and a new sample $S$ of this new size is used to compute the new step. By starting
with small sample sizes it is hoped that in the early stages of the algorithm the solution can be computed efficiently in each iteration. The key point in this procedure is clearly the rule that checks throughout the progression of the algorithm, whether the approximation we are performing is good enough, i.e. the sample size is big enough, or has to be increased. Because of this systematic check on the quality of approximation in each step of the algorithm, such sampling strategy is called *dynamic*.

As in [3], we consider a condition on the batch gradient $\nabla J_S$ which imposes in every stage of the optimisation that the direction $-\nabla J_S$ is a descent direction for $J$ at $\lambda$ if the following condition holds:

$$ \|\nabla J_S(\lambda) - \nabla J(\lambda)\|_{L^2(\Omega)} \le \theta \|\nabla J_S(\lambda)\|_{L^2(\Omega)}, \quad \theta \in [0, 1). \quad (2.6) $$

The computation of $\nabla J$ may be very expensive for applications involving large databases and nonlinear constraints, so we rewrite (2.6) as an estimate of the variance of the random vector $\nabla J_S(\lambda)$. In order to do that, recalling definitions (1.4) and (1.5) we observe that

$$ \nabla J_S(\lambda) = \frac{1}{2|S|} \sum_{k \in S} \nabla l(\lambda, f_k) \quad (2.7) $$

We can compute (2.7) in correspondence to an optimal solution $\hat{\lambda}$ by using [5, Remark 3.4] where a characterisation of $\nabla J$ is given in terms of the adjoint states $p_k$ (see Section 3 for details). By linearity and extending to the multiple-constrained case, we get:

$$ \nabla J_S(\hat{\lambda}) = \sum_{k \in S} \sum_{i=1}^{d} \int_{\Omega} \phi_i'(\hat{u}_k, f_k) p_k dx \quad (2.8) $$

Thanks to this characterisation, we now extend the dynamic sampling algorithm in [3] to the case where the solution map $\hat{S}$ is nonlinear: by taking (2.8) into account and following [3, Section 3] we can rewrite (2.6) as a condition on the variance of the batch gradient that reads as

$$ \frac{\| Var_{k \in S}(\nabla l(\lambda, f_k)) \|_{L^1(\Omega)} N - |S|}{|S| N - 1} \le \theta^2 \| \nabla J_S(\lambda) \|_{L^2(\Omega)}^2. \quad (2.9) $$

For a detailed derivation of (2.9), see [3]. Condition (2.9) is the responsible for possible changes in the sample size in the optimisation algorithm and has to be checked in every iteration. If inequality (2.9) is not satisfied, a larger sample $\hat{S}$ whose size satisfies the descent condition (2.9) needs to be considered. By assuming that the change in the sample size is gradual enough such that, for any given $\lambda$:

$$ \| Var_{k \in \hat{S}}(\nabla l(\lambda, f_k)) \|_{L^1(\Omega)} \approx \| Var_{k \in S}(\nabla l(\lambda, f_k)) \|_{L^1(\Omega)}, $$

$$ \| \nabla J_{\hat{S}}(\lambda) \|_{L^2(\Omega)} \approx \| \nabla J_S(\lambda) \|_{L^2(\Omega)}, $$
we see that condition (2.9) is satisfied whenever we choose $|\hat{S}|$ such that

$$|\hat{S}| \ge \left[ \frac{N - \|Var_{k \in S}(\nabla l(\lambda, f_k))\|_{L^1(\Omega)}}{\|Var_{k \in S}(\nabla l(\lambda, f_k))\|_{L^1(\Omega)} + \theta^2(N-1) \|\nabla J_S(\lambda)\|_{L^2(\Omega)}^2} \right]. \quad (2.10)$$

Conditions (2.9) and (2.10) are the key points in the optimisation algorithm we are going to present: by checking the former, one can control whether the sampling approximation is accurate enough and if this is not the case at any stage of the algorithm, by imposing the latter a new larger sample size is determined.

We remark that these two conditions force the direction $-\nabla J_S$ to be a descent direction. Steepest descent methods are known to be slowly convergent. Algorithms incorporating information coming from the Hessian are generally more efficient. However, normally the computation of the Hessian is very expensive, so Hessian-approximating methods are commonly used. In [3] a Newton-CG method is employed. There, an approximation of the Hessian matrix $\nabla^2 J_S$ is computed only on a subsample $H$ of $S$ such that $|H| \ll |S|$. As the sample $S$ is dynamically changing, the subsample $H$ will change as well (with a fixed, constant ratio) and the computation of the new conjugate gradient direction can be performed efficiently. In this work, in order to compute an approximation of the Hessian we consider the well-known BFGS method which has been extensively used in the last years because of its efficiency and low computational costs.

Before giving a full description of the resulting algorithm solving (1.1)-(1.2), we briefly comment on the linesearch rule that is employed in the update of the BFGS matrix. We choose an Armijo backtracking line search rule with curvature verification: the BFGS matrix is updated only if the curvature condition is satisfied. The Amijo criterion is:

$$J_S(\lambda_k + \alpha_k d_k) - J_S(\lambda_k) \le \alpha_k \eta \nabla J_S(\lambda_k)^\top d_k \quad (2.11)$$

where the value $\eta$ will be specified in Section 3, $d_k$ is the descent direction of the quasi-Newton step, $\alpha_k$ is the length of the quasi-Newton step and $\nabla J_S(\lambda_k)$ is defined in (2.7). The positivity of the parameters is always preserved along the iterations.

We present now the BFGS optimisation with Dynamic sampling for solving (1.1)-(1.2): compared to [3, Algorithm 5.2] we stress once more that the gain in efficiency is obtained thanks to the use of BFGS instead of the Newton-CG sampling method.

## 3 Numerical results

In this section we present the numerical results of the Dynamic Sampling Algorithm 1 applied to compute the numerical solution of (1.1)-(1.2). In our numerical computations we fix the parameter values as follows:

- We consider images of size 150 × 150. We approximate the differential operators by discretising with finite difference schemes with mesh step size
<caption>**Algorithm 1** Dynamic Sampling BFGS for solving (1.1)-(1.2)</caption>

1: Initialize: $\lambda_0$, sample $S_0$ with $|S_0| \ll N$ and model parameter $\theta$, $k = 0$.
2: **while** BFGS not converging, $k \ge 0$
3: &nbsp;&nbsp;&nbsp;&nbsp;sample $|S_k|$ PDE constraints to solve;
4: &nbsp;&nbsp;&nbsp;&nbsp;update of the BFGS matrix;
5: &nbsp;&nbsp;&nbsp;&nbsp;compute direction $d_k$ by BFGS and steplength $\alpha_k$ by Armijo cond. (2.11);
6: &nbsp;&nbsp;&nbsp;&nbsp;define new iterate: $\lambda_{k+1} = \lambda_k + \alpha_k d_k$;
7: &nbsp;&nbsp;&nbsp;&nbsp;**if** condition (2.9) **then**
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;maintain the sample size: $|S_{k+1}| = |S_k|$;
9: &nbsp;&nbsp;&nbsp;&nbsp;**else** augment $S_k$ such that condition (2.10) is verified.
10: **end**

$h = 1/(\text{number of pixels in the } x\text{-direction})$. We use forward difference for the discretisation of the divergence operator and backward differences for the gradient. The Laplace operator is discretised by using the usual five point formula.

- The TV constraints in (1.2) are solved by means of SemiSmooth Newton (SSN) algorithms whose form depends on the $\phi$'s in (1.2) (cf. [5, Section 4]) solving regularised problems which stop if either the difference between two consecutive iterates is small enough or if the maximum number of iterations $\maxiter = 35$ is reached.

- In the Armijo condition (2.11) the value $\eta$ is chosen to be $\gamma = 10^{-4}$.

**Single noise estimation.** As a toy example, we start by considering the case when the noise in the images is normally distributed. In (1.1)-(1.2), this reflects in the estimation of just one parameter $\lambda$ that weights the fidelity term $\phi(u, f_k) = (u - f_k)^2$ in each constraint. Considering the training database $\{(u_k, f_k)\}_{k=1,...,N}$ of clean and noisy images, the problem reduces to:

$$ \min_{\lambda \ge 0} \frac{1}{2N} \sum_{k=1}^{N} \| \hat{u}_k - u_k \|_{L^2(\Omega)}^2 \quad (3.1) $$

where, for each $k$, $\hat{u}_k$ is the solution of the regularised PDE

$$ -\varepsilon \Delta \hat{u}_k - \operatorname{div}(h_{\gamma}(\nabla \hat{u}_k)) + \lambda(\hat{u}_k - f_k) = 0, \quad k = 1, \dots, N. \quad (3.2) $$

In (3.2) $h_\gamma$ arises from a Huber-type regularisation of the subdifferential of $|D\hat{u}_k|$ with parameter $\gamma \gg 1$ and the $\varepsilon$ term is an artificial diffusion term that sets up the problem in the Hilbert space $H_0^1(\Omega)$ (see [5, Section 3] for details).

As shown in [5, Theorem 3.5] the adjoint states $p_k$ can be computed for each constraint as the solution of the following equation

$$ \varepsilon(Dp_k, Dv)_L^2 + (h'_\gamma(D\hat{u}_k)^* Dp_k, Dv)_L^2 \\ + \int_\Omega \lambda p_k v \, dx = -(\hat{u}_k - f_k, v)_L^2, \quad \forall v \in H_0^1(\Omega). \quad (3.3) $$
Recalling also equations (2.7)-(2.8) needed for the computation of the gradient, we can now apply Algorithm 1 to solve (3.1)-(3.2).

For the following numerical tests, the parameters of this model are chosen as follows: $\varepsilon = 10^{-12}$, $\gamma = 100$. The noise in the images has distribution $N(0, 0.05)$. The parameter $\theta$ of the Algorithm 1, is chosen to be $\theta = 0.5$. We will comment on the sensitivity of the method to $\theta$ later on.

Table 1 shows the numerical value of the optimal parameter $\hat{\lambda}$ when varying the size of the dictionary. We measure the efficiency of the algorithms used in terms of the number of nonlinear PDEs solved during the BFGS optimisation and we compare the efficiency of solving (3.1)-(3.2) without and with the Dynamic Sampling strategy. We observe a clear improvement in efficiency when using Dynamic Sampling: the number of PDEs solved in the optimisation process is very much reduced. We note that this corresponds to an increasing number of BFGS iterations which does not appear to be an issue as BFGS iterations are themselves very fast. For the sake of computational efficiency, what really matters is the number of PDEs that need to be solved in each iteration of BFGS. Moreover, thanks to modern parallel computing methods and to the decoupled nature of the constraints in each BFGS iteration, solving such a reduced amount of PDEs makes the computational efforts very reasonable. In fact, we note that the size of the sample is generally maintained very small in comparison to $N$ or just slightly increased. Computing also the relative error between the optimal parameter computed by solving all the PDEs and the one computed with Dynamic Sampling method, we note a good level of accuracy: the difference between the two values remains below 5%.

<table><caption>Table 1: $N$ is the size of the database, $\hat{\lambda}$ is the optimal parameter for (3.1)-(3.2) obtained by solving all the $N$ constraints, whereas $\hat{\lambda}_S$ is the one computed by solving the problem with Algorithm 1. The initial size $S_0$ is chosen to be $|S_0| = 20\%N$. $|S_{end}|$ of the sample at the end of the optimisation algorithm. The efficiency of the algorithms is measured in terms of the PDEs solved. We compare the accuracy of the approximation in terms of the difference $\|\hat{\lambda}_S - \hat{\lambda}\|_1 / \|\lambda_S\|_1$.</caption><thead><tr><th>N</th><th>λ̂</th><th>λ̂<sub>S</sub></th><th>|S<sub>0</sub>|</th><th>|S<sub>end</sub>|</th><th>eff.</th><th>eff. Dyn.S.</th><th>BFGS its.</th><th>BFGS its. Dyn.S.</th><th>diff.</th></tr></thead><tbody><tr><td>10</td><td>3334.5</td><td>3427.7</td><td>2</td><td>3</td><td>140</td><td><b>84</b></td><td>7</td><td>21</td><td>2.7%</td></tr><tr><td>20</td><td>3437.0</td><td>3475.1</td><td>4</td><td>4</td><td>240</td><td><b>120</b></td><td>7</td><td>15</td><td>1.1%</td></tr><tr><td>30</td><td>3436.5</td><td>3478.2</td><td>6</td><td>6</td><td>420</td><td><b>180</b></td><td>7</td><td>15</td><td>1.2%</td></tr><tr><td>40</td><td>3431.5</td><td>3358.3</td><td>8</td><td>9</td><td>560</td><td><b>272</b></td><td>7</td><td>16</td><td>2.1%</td></tr><tr><td>50</td><td>3425.8</td><td>3306.4</td><td>10</td><td>10</td><td>700</td><td><b>220</b></td><td>7</td><td>11</td><td>3.5%</td></tr><tr><td>60</td><td>3426.0</td><td>3543.4</td><td>12</td><td>12</td><td>840</td><td><b>264</b></td><td>7</td><td>11</td><td>3.3%</td></tr><tr><td>70</td><td>3419.7</td><td>3457.7</td><td>14</td><td>14</td><td>980</td><td><b>336</b></td><td>7</td><td>12</td><td>1.1%</td></tr><tr><td>80</td><td>3418.1</td><td>3379.3</td><td>16</td><td>16</td><td>1120</td><td><b>480</b></td><td>7</td><td>15</td><td>&lt; 1%</td></tr><tr><td>90</td><td>3416.6</td><td>3353.5</td><td>18</td><td>18</td><td>1260</td><td><b>648</b></td><td>7</td><td>18</td><td>2.3%</td></tr><tr><td>100</td><td>3413.6</td><td>3479.0</td><td>20</td><td>20</td><td>1400</td><td><b>520</b></td><td>7</td><td>13</td><td>1.9%</td></tr></tbody></table>

Figure 1 shows an example of database of brain images<fnref n="4" /> together with the optimal denoised version obtained by Algorithm 1 for single Gaussian noise estimation.

<footnote n="4">⁴ OASIS online database.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1: Sample of 5 images of a MRI brain database: original images (upper row), noisy images (middle row) and optimal denoised images (bottom row), $\lambda_S = 3280.5$.</figcaption></figure>

**Multiple noise estimation.** We consider now a more interesting application of (1.1)-(1.2) where the image is corrupted by noises with different distributions. We consider the case where a combination of Gaussian and impulse noise is present. The fidelity term for the impulse distributed component is $\phi_1(u, f_k) = |u - f_k|$, whereas, as above, for the Gaussian noise we consider the fidelity $\phi_2(u, f_k) = (u - f_k)^2$, for every $k$. Each fidelity term is weighted by a parameter $\lambda_i$, $i = 1, 2$. Thus, we aim to solve:

$$ \min_{(\lambda_1, \lambda_2), \lambda_i \ge 0} \frac{1}{2N} \sum_{k=1}^{N} \|\hat{u}_k - u_k\|_{L^2(\Omega)}^2 \quad (3.4) $$

where, for each $k$, $\hat{u}_k$ is now the solution of the regularised PDE:

$$ -\varepsilon \Delta \hat{u}_k - \operatorname{div}(h_\gamma(\nabla \hat{u}_k)) + \lambda_1 h_\gamma^1(\hat{u}_k - f_k) + \lambda_2 (\hat{u}_k - f_k) = 0, \quad k = 1, \dots, N. \quad (3.5) $$

In (3.5) the first and the second terms are as before while the third one corresponds to the Huber-type regularisation of sgn($\hat{u}_k - f_k$). The adjoint state is computed as in [5], in a similar manner as (3.3). By taking also into account equations (2.7)-(2.8), we solve (3.4)-(3.5) with $\varepsilon = 10^{-12}$, $\gamma = 100$ by means of Algorithm 1.

We take as example slices of the brain database shown in Figure 1 corrupted with both Gaussian noise distributed as $N(0, 0.005)$ and impulse noise with fraction of missing pixels $d = 5\%$, and again solve (1.1)-(1.2) by solving the PDE constraints all at once and by using Dynamic Sampling for different $N$. In Table 2 we report the results for the estimation of $\lambda_1$ and $\lambda_2$.
<table><caption>Table 2: $\hat{\lambda}_{1S}$ and $\hat{\lambda}_{2S}$ are the optimal weights for (3.4)-(3.5) estimated with Dynamic Sampling. We observe again a clear improvement in efficiency (i.e. number of PDEs solved). As above, $|S_0| = 20\%N$ and $\theta = 0.5$.</caption><thead><tr><th>N</th><th>λ̂<sub>1S</sub></th><th>λ̂<sub>2S</sub></th><th>|S<sub>0</sub>|</th><th>|S<sub>end</sub>|</th><th>eff.</th><th>eff. Dyn.S.</th><th>diff.</th></tr></thead><tbody><tr><td>10</td><td>86.31</td><td>28.43</td><td>2</td><td>7</td><td>180</td><td><b>70</b></td><td>5.2%</td></tr><tr><td>20</td><td>90.61</td><td>26.96</td><td>4</td><td>6</td><td>920</td><td><b>180</b></td><td>5.3%</td></tr><tr><td>30</td><td>94.36</td><td>29.04</td><td>6</td><td>7</td><td>2100</td><td><b>314</b></td><td>5.6%</td></tr><tr><td>40</td><td>88.88</td><td>31.56</td><td>8</td><td>8</td><td>880</td><td><b>496</b></td><td>1.2%</td></tr><tr><td>50</td><td>88.92</td><td>29.81</td><td>10</td><td>10</td><td>2200</td><td><b>560</b></td><td>&lt; 1%</td></tr><tr><td>60</td><td>89.64</td><td>28.36</td><td>12</td><td>12</td><td>1920</td><td><b>336</b></td><td>1.9%</td></tr><tr><td>70</td><td>86.09</td><td>28.09</td><td>14</td><td>14</td><td>2940</td><td><b>532</b></td><td>3.3%</td></tr><tr><td>80</td><td>87.68</td><td>29.97</td><td>16</td><td>16</td><td>3520</td><td><b>448</b></td><td>&lt; 1%</td></tr></tbody></table>

**Convergence and sensitivity.** Figure 2 shows two features of Algorithm 1 applied to solve problem (3.1)-(3.2). On the left we represent the evolution of the cost functional along the BFGS iterations. Because of the sampling strategy, in the early iterations of BFGS the problem considered varies quite a lot, thus showing oscillations. Once evolving the process, the convergence is super-linear. On the right we represent the sensitivity with respect to the accuracy parameter $\theta$ (cf. (2.9)): smaller values of $\theta$ penalise larger variances on $\nabla J_S$, thus favouring larger samples. Larger values of $\theta$ allow larger variances on $\nabla J_S$ and, consequently, smaller sample sizes. In this case, efficiency improves, but accuracy might suffer as shown in Table 3.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Left: evolution of BFGS with Dynamic Sampling along the iterations. Right: samples size changes in Algorithm 1 for different values of $\theta$. For each value of $\theta$, the result is plotted till convergence. For this example $N = 20$, $|S_0| = 2$.</figcaption></figure>

<table><caption>Table 3: As $\theta$ increases we observe improvements upon the efficiency as smaller samples are allowed. However, the relative difference with the value estimated without sampling shows that accuracy suffers.</caption><thead><tr><th>θ</th><th>efficiency</th><th>difference</th></tr></thead><tbody><tr><td>0.1</td><td>516</td><td>0.07%</td></tr><tr><td>0.3</td><td>246</td><td>4.3%</td></tr><tr><td>0.5</td><td>92</td><td>5.9%</td></tr><tr><td>0.7</td><td>68</td><td>15%</td></tr></tbody></table>
Conclusions

In this paper, we propose an efficient and competitive technique to solve numer-
ically the constrained optimisation problem (1.1)-(1.2) designed for learning the
noise model in a TV denoising framework accounting for different types of noise.
The set of nonsmooth PDE constraints resembles a large-size training database
of clean and noisy images that allows a more robust estimation of parameters.
To solve the problem, we use *Dynamic Sampling* methods, proposed in [3] for
*linear* constrained problems. The idea consists in selecting just a small sample
of the PDEs that need to be solved over the whole database and then, during
the progression of the algorithm, verify whether such a size produces approx-
imations that are accurate enough. Extended to our nonlinear framework, the
results show a remarkable improvement in efficiency, which reflects in reduced
computational times for both single noise estimations as well as for mixed ones.
Further directions for future research are an accurate analysis of convergence
properties of such a scheme as well as the design of a similar algorithm for the
case of a $L^1$-regularisation on the parameter vector.

References

1. Ambrosio, L.; Fusco, N.; Pallara, D.: Functions of Bounded Variation and Free Discontinuity problems, Oxford Mathematical Monographs. Oxford: Clarendon Press. xviii, (2000).

2. Bottou, L., Bousquet, O.: *The tradeoffs of large scale learning*, Advances in Neural Information Processing Systems, 20, 161-168 (2008).

3. Byrd R. H., Chin G. M., Nocedal J., Wu Y.: *Sample size selection in optimization methods for machine learning*, Math. Program. Ser. B, 134, 127-155 (2012).

4. Byrd R. H., Chin G. M., Neveitt W., Nocedal J.: *On the use of stochastic Hessian information in unconstrained optimisation*, SIAM J. Optim., 21(3), 977-995 (2011).

5. De Los Reyes J. C., Schönlieb C.-B.: *Image denoising: Learning noise distribution via nonsmooth PDE-constrained optimisation*, Inverse Problems and Imaging, vol. 7(4), 1183-1214 (2013).

6. Haber E., Chung M., Herrmann F., *An effective method for parameter estimation with PDE constraints with multiple right-hand sides*, SIAM J. Optim., 22(3), 739-757 (2012).

7. Mantle M.D., Sederman A.J., Gladden L.F.: *Single- and two- phase flow in fixed-bed reactors: MRI flow visualisation and lattice-Boltzmann simulations*, Chemical Engineering Science, 56, 523-529 (2001).

8. Nesterov, Y.: *Primal-dual subgradient methods for convex problems*, Math. Program., 120(1), 221-259 (2009).

9. Polyak, B., Juditsky, A.: *Acceleration of stochastic approximation by averaging*, SIAM J. Control Optim., 30, (1992).

10. Robbins, H., Monro, S.: *A stochastic approximation method*, Ann. Math. Stat., 22(3), 400-407 (1951).