Learning Filter Functions in Regularisers by Minimising Quotients

TL;DR

This paper introduces a method for learning regularisation functions in inverse problems by quotient minimisation, extending previous models to higher dimensions and multiple data types, and proposing novel non-derivative regularisers.

Contribution

It extends quotient minimisation to include higher-dimensional filter functions and multiple training data types, and introduces new non-derivative regularisers.

Findings

Behaves like total variation in 1D.

Learns scales and geometric properties.

Proposes novel non-derivative regularisers.

Abstract

Learning approaches have recently become very popular in the field of inverse problems. A large variety of methods has been established in recent years, ranging from bi-level learning to high-dimensional machine learning techniques. Most learning approaches, however, only aim at fitting parametrised models to favourable training data whilst ignoring misfit training data completely. In this paper, we follow up on the idea of learning parametrised regularisation functions by quotient minimisation as established in [3]. We extend the model therein to include higher-dimensional filter functions to be learned and allow for fit- and misfit-training data consisting of multiple functions. We first present results resembling behaviour of well-established derivative-based sparse regularisers like total variation or higher-order total variation in one-dimension. Our second and main contribution is the introduction of novel families of non-derivative-based regularisers. This is accomplished by learning favourable scales and geometric properties while at the same time avoiding unfavourable ones.

Figures

Equations

$$\hat{h}\in \underset{\begin{array}{c}\Vert h{\Vert }_2=1\\ \text{mean}(h)=0\end{array}}{\mathrm{argmin}}\frac{J(u^{+};h)}{J(u^{-};h)},J(u;h)=\Vert u\ast h{\Vert }_1,$$

$$\hat{h}\in \underset{\begin{array}{c}\Vert h{\Vert }_2=1\\ \text{mean}(h)=0\end{array}}{\mathrm{argmin}}\frac{\frac{1}{M}\sum _{i=1}^M\sum _{k=1}^{K}J(u_i^{+};h_k)}{\frac{1}{N}\sum _{j=1}^N\sum _{k=1}^{K}J(u_j^{-};h_k)},J(u;h)=\Vert u\ast h{\Vert }_1.$$

$$\hat{h}\in \underset{h}{\mathrm{argmin}}\left\{\frac{F(h)}{G(h)}\text{subject to}\Vert h{\Vert }_2=1\text{and}\text{mean}(h)=0\right\}.$$

$$\{\begin{array}{ll}{h}^{k+\frac{1}{2}}& =\underset{\text{mean}(h)=0}{\mathrm{argmin}}\left\{F(h)-{\mu }^k⟨h-h^k,s^k⟩+{\Vert h-h^k\Vert }_2^2\right\}\\ {\mu }^{k+1}& =\frac{F(h^{k+\frac{1}{2}})}{G(h^{k+\frac{1}{2}})}\\ s^{k+1}& \in \partial G(h^{k+\frac{1}{2}})\\ h^{k+1}& =\frac{h^{k+\frac{1}{2}}}{{\Vert h^{k+\frac{1}{2}}\Vert }_2}\end{array}.$$

$$\hat{u}=\underset{u\in {\mathbb{R}}^m}{\mathrm{argmin}}J(u;\hat{h})\text{subject to}\Vert u-f{\Vert }_2\le \eta \sigma \sqrt{m},$$

$${\mu }^{k+\frac{1}{2}}+\frac{1}{G(h^{k+\frac{1}{2}})}\Vert h^{k+\frac{1}{2}}-h^k{\Vert }^2\le {\mu }^k,$$

$$F(h^{k+\frac{1}{2}})+\Vert h^{k+\frac{1}{2}}-h^k{\Vert }_2^2$$

$$\le F(h^k)+{\mu }^k\left(G(h^{k+\frac{1}{2}})-G(h^k)\right)$$

$$={\mu }^{k}G(h^{k+\frac{1}{2}}),$$

$${\mu }^{k+\frac{1}{2}}+\frac{1}{G(h^{k+\frac{1}{2}})}\Vert h^{k+\frac{1}{2}}-h^k{\Vert }^2\le {\mu }^k,$$

$${\varphi }_{\gamma }(x)=\{\begin{array}{ll}\frac{x^2}{2},& |x|\le \gamma \\ \gamma \left(|x|-\frac{\gamma }{2}\right),& |x|>\gamma \end{array},$$

$$\Vert r^{k+\frac{1}{2}}-{\mu }^{k+\frac{1}{2}}\nabla G(x^{k+\frac{1}{2}}){\Vert }_2\le (2+C^{k+\frac{1}{2}}L)\Vert h^{k+\frac{1}{2}}-h^k{\Vert }_2,$$

$$r^{k+\frac{1}{2}}-{\mu }^k\nabla G(h^k)=2(h^k-h^{k+\frac{1}{2}}),$$

$$\Vert r^{k+\frac{1}{2}}-{\mu }^{k+\frac{1}{2}}\nabla G(x^{k+\frac{1}{2}}){\Vert }_2$$

$$=\Vert 2(h^k-h^{k+\frac{1}{2}})+{\mu }^k\nabla G(x^k)-{\mu }^{k+\frac{1}{2}}\nabla G(x^{k+\frac{1}{2}}){\Vert }_2$$

$$\le 2\Vert h^{k+\frac{1}{2}}-h^k{\Vert }_2+C^{k+\frac{1}{2}}\Vert \nabla G(x^{k+\frac{1}{2}})-\nabla G(x^k){\Vert }_2$$

$$\Vert r^{k+\frac{1}{2}}-{\mu }^{k+\frac{1}{2}}\nabla G(x^{k+\frac{1}{2}}){\Vert }_2$$

$$=(2+C^{k+\frac{1}{2}}L)\Vert h^{k+\frac{1}{2}}-h^k{\Vert }_2.$$

$$\frac{\Vert h_1\ast v{\Vert }_1+\Vert h_2\ast w{\Vert }_1}{\Vert h_1\ast u^{-}{\Vert }_1+\Vert h_2\ast u^{-}{\Vert }_1},$$

Full text

Keywords: Regularisation Learning, Non-linear Eigenproblem, Sparse Regularisation, Generalised Inverse Power Method

1 Introduction

Learning approaches for variational regularisation models constitute an active area of current research. In so-called bi-level learning approaches [8, 14], for instance, one seeks to minimise a cost functional subject to a variational minimisation problem usually consisting of a data fidelity term and a regularisation term. Application of such models range from learning of suitable regularisation parameters to learning the correct operator or entire model, strongly dependent on the type of the underlying problem. In [9], the authors compared performance of Total Variation (TV), Infimal Convolution TV (ICTV) and second-order Total Generalised Variation (TGV2) regularisers combined with both $L_{1}$ and $L_{2}$ cost functions for denoising of 200 images of the BSDS300 dataset measured by SSIM, PSNR and an objective value. There was no unique regulariser that always performed best. The images in the above-mentioned dataset differ from each other significantly enough such that advantages of the different regularisers become apparent for images with different prominent features such as sharp edges or piecewise linear regions.

Another approach to variational regularisation learning is dictionary learning [7]. In this approach the basic paradigm is that local image regions (patches) can be composed based on a linear combination of very few atoms from some dictionary. The dictionary could be global, for example wavelets or DCT-based, and in those cases a basis. However, it was revealed that tailored dictionaries to the specific image (or class of images), which are overcomplete, outperform the global dictionaries. These dictionaries are typically learned from the noisy image itself using algorithms such as K-SVD [1] based on orthogonal-matching-pursuit (OMP). Recent studies have shown relations between convolutional neural nets and convolutional sparse coding [17]. One can conceptually perceive our proposed convolution filter set $\{h_{i}\}$, which defines the target-specific regulariser, as a small dictionary which is learned based on a few positive and negative image examples.

Learning approaches for variational regularisation models are aiming to design appropriate regularisation and customise it to particular structures present in the image. A somewhat separate route is the mathematical analysis of model-based regularisation, aiming at understanding the main building-blocks of existing regularisers to pave the path for designing new ones. In [2], for instance, the concept of ground states, singular values and singular vectors of regularisation functionals has been introduced, enabling the computation of solutions of variational regularisation schemes that can be reconstructed perfectly (up to a systematic bias).

In [3] a new model motivated by generalised, non-linear Eigenproblems has been proposed to learn parametrised regularisation functions. As a novelty, both wanted and unwanted outcomes are incorporated in the model by integrating the former in the numerator and the latter in the denominator:

[TABLE]

where $h$ is a parametrisation of a regularisation functional $J$ and $u^{+}$ and $u^{-}$ are desired and undesired input signals, respectively.

The basic idea underlying the model in [3] is optimisation of a quotient with respect to a convolution kernel parametrising certain regularisation functionals. In the paper, the authors investigate the same regularisation function, which is the one-norm of a signal convolved with a kernel $h$ both in the numerator and denominator. In the former, the input is a desirable signal, i.e. a function, which is preferred to be sparse once convolved with the kernel, whereas the latter is a signal to be avoided, yielding a large one-norm once convolved. In [3], the undesirable signal has only been chosen to be pure noise or a noisy version of the desired input signal. In this work, however, we are also going to use clean signals as undesirable signals, with specific geometric properties or scales one wants to avoid. In Section 4 we are going to see that this will enable us to derive tailored filters superior to those derived merely from desirable fitting data.

2 The Proposed Learning Model

In order to be able to incorporate multiple input functions, different regularisation functionals and multi-dimensional filter functions, we generalise the model in [3] as follows:

[TABLE]

Now, $\hat{h}=(\hat{h}_{1},\dots,\hat{h}_{K})$, where $\hat{h}_{k}\in\mathbb{R}^{n}$ for all $k\in\{1,\dots,K\}$, is a combination of multiple filter functions. The signals $u^{+},u^{-}\in\mathbb{R}^{m}$ are one-dimensional or two-dimensional images written as a column vector. In the following section we want to describe how we want to solve (2) numerically.

2.1 Numerical Implementation

Viewing the quotients in (1) and (2) as generalised Rayleigh quotients, we observe that we deal with the (numerical) solution of generalised Eigenvalue problems. In order to solve (1) and (2) with the same algorithm, we write down an abstract algorithm for the solution of

[TABLE]

The optimality condition of (3) is given via $0\in\partial F(\hat{h})-\hat{\mu}\,\partial G(\hat{h})$, where $\partial F(\hat{h})$ and $\partial G(\hat{h})$ denote the subdifferential of $F$ and $G$ at $\hat{h}$, respectively, and $\hat{\mu}=F(\hat{h})/G(\hat{h})$. Note that the Lagrange multipliers for the constraints are zero by the same argumentation as in [18, Section 2], and can therefore be omitted.

In [13] the authors have proposed a generalised inverse power method to tackle problems of the form (3). We, however, follow [5] and use a modification with added penalisation of the squared two-norm between $h$ and the previous iterate, to guarantee coercivity (and therefore existence and uniqueness of the solution) of the main update. The proposed algorithm for solving (3) therefore reads as

[TABLE]

Similar to [3] we are using the CVX MATLAB® software for disciplined convex programming [12]. Due to the non-convexity of the overall problem and the resulting dependence on random initialisations of the filter, we re-initialise $h$ and iterate (4) 100 times. As also explained in [3], reconstruction of a noisy signal in order to test the behaviour of the optimal filter is obtained by solving the following constrained optimisation problem:

[TABLE]

where $f$ is the sum of $u^{+}$ and Gaussian noise with zero mean and variance $\sigma^{2}$, $\eta$ is a weighting factor and $m$ is the number of elements of $u^{+}$.

Remark.

In the setting of [3], there is indeed an even more efficient way of finding suitable filter functions $h$. Simplifying the model to a variant without the need of having a negative input function $u^{-}$ yields the same results, which is a clear indicator that in the above-mentioned framework the numerator plays a dominant role. In fact, varying model (1) by replacing the one-norm in the denominator by $\|h\|_{2}$ returns exactly the same solutions. However, we would like to stress that the denominator is going to play a more important role in our extended model, since we are able to incorporate more than one input function $u^{-}$, especially ones which are different from pure noise. In fact, one can think of a large variety of undesired input signals such as specific textures and shapes.

Despite existing convergence results previously stated in [13] and [5] we want to briefly state a simplified convergence result for global convergence of Algorithm 4 in the following.

2.2 A brief convergence analysis

Following [4, Section 3.2], we show two results that are essential for proving global convergence of Algorithm (4): a descent lemma and a bound of the subgradient by the iterates gap. We start with the sufficient decrease property of the objective.

Lemma 1**.**

Let $F$ and $G$ be proper, lower semi-continuous and convex functions. Then the iterates of Algorithm (4) satisfy

[TABLE]

if we further assume $G(h^{k+\frac{1}{2}})\neq 0$ for all $k\in\mathbb{N}$.

Proof.

From the first equation of Algorithm (4) we observe

[TABLE]

due to the convexity of $G$. If we divide by $G(h^{k+\frac{1}{2}})$ on both sides of the equation, we obtain

[TABLE]

which concludes the proof. ∎

In order to further prove a bound of the subgradient by the iterates gap, we assume that $G$ is smooth and further has a Lipschitz-continuous gradient $\nabla G$. We want to point out that this excludes choices for $G$ such as in (1) and (2), as the one-norm is neither smooth nor are its subgradients Lipschitz-continuous. A remedy here is the smoothing of the one-norms in (1) and (2). If we replace the one-norm(s) in the denominator with Huber one-norms, i.e. we replace the modulus in the one-norm with the Huber function

[TABLE]

we can achieve smoothness and Lipschitz-continuity of the gradient, where the Lipschitz parameter depends on the smoothing parameter $\gamma$. We want to note that for $\gamma$ small enough we have not seen any significant difference in numerical performance between using the one-norm or its Huber counterpart.

Lemma 2**.**

Let $F$ and $G$ be proper, lower semi-continuous and convex functions, and let $G$ be differentiable with $L$-Lipschitz-continuous gradient, i.e. $\|\nabla G(h_{1})-\nabla G(h_{2})\|_{2}\leq L\|h_{1}-h_{2}\|_{2}$ for all $h_{1}$ and $h_{2}$ and a fixed constant $L$. Then the iterates of Algorithm (4) satisfy

[TABLE]

for some constant $C^{k+\frac{1}{2}}$, $r^{k+\frac{1}{2}}\in\partial F(h^{k+\frac{1}{2}})$ and $\mu^{k+\frac{1}{2}}:=\mu^{k+1}=F(h^{k+\frac{1}{2}})/G(h^{k+\frac{1}{2}})$.

Proof.

This follows almost instantly from the optimality condition and the Lipschitz-continuity of $\nabla G$. We obtain

[TABLE]

for $r^{k+\frac{1}{2}}\in\partial F(h^{k+\frac{1}{2}})$, as the optimality condition of the first sub-problem of (4) - note that we can omit the zero-mean constraint with a similar argumentation as earlier. Hence, we obtain

[TABLE]

This concludes the proof. ∎

Under the additional assumption that the function $F/G$ satisfies the Kurdyka-Łojasiewicz property (cf. [16, 15]) we can now use Lemma 1 and Lemma 2 to show finite length of the iterates (4) similar to [4, Theorem 1], following the general recipe of [4, Section 3.2]. Note that we further have to substitute $h^{k+\frac{1}{2}}=\|h^{k+\frac{1}{2}}\|_{2}h^{k+1}$ in order to also show global convergence of the normalised iterates.

3 Reproducing Standard Sparse Penalties

In this section we want to demonstrate that we are able to reproduce standard first- and second-order total variation regularisation penalties in 1D.

Figure 1 shows results for different sizes of the kernel $h$. In all experiments the filter function is indeed resembling a two-point stencil functioning as a finite differences discretisation of TV. This is expected as the desired input function is a TV Eigenfunction.

In Figure 2 (a) we can reproduce a filter resembling a second-order derivative. This is indeed expected as we choose three different piecewise-linear functions as desired input signals. The reconstruction in (b) is performed according to (5).

In a more sophisticated example, we mimic a TV-TV2 infimal convolution model, where we are given a known decomposition $u^{+}=v+w$, i.e. $u^{+}$ consists of a smooth part $v$ and a piecewise constant part $w$. When minimising

[TABLE]

with respect to $h_{1}$ and $h_{2}$, we indeed obtain two filters resembling a second- and first-order derivative, respectively (cf. Figure 3).

4 Novel Sparse Filters

In this section we derive a new family of regularisers not necessarily related to derivatives in contrast to the total variation. They have the interesting property of reconstructing piecewise-constant both vertical and horizontal lines in the corresponding null-spaces. Consequently, we are able to almost perfectly reconstruct those types of images and obtain better denoising results compared to standard TV denoising. In [10], a definition of desirable features of a regulariser, which is adapted for a specific type of images, is given. It is suggested that in the ideal case, all instances belonging to the desired clean class should be in the null-space of the regulariser (see [10, Section 2] and also compare [6]). This is exactly what we obtain in the following.

In Figure 4, a new family of diagonal regularisers is established for piecewise-constant images with stripes in both vertical and horizontal direction. For denoising purposes, those filters yield superior results over TV denoising as they additionally avoid loss of contrast, which would occur when performing TV denoising for these examples. The reason for that is simply that if we consider a $2\times 2$ diagonal-shaped filter $h=[1,-1;-1,1]$ in the variational problem (5), we can expect both horizontal and vertical stripes to be in its null-space. Therefore, we obtain perfect shape preservation for any regularisation parameter $\eta$. Note that 1-pixel-thick stripes are in the null-space as well (cf. Figure 5 (left), where $f-u=0$).

In Figure 5 (centre) we can observe that rectangles are also well preserved with this filter. For better performance, however, one could additionally use a contrast preserving mechanism such as Bregman iteration or as in our case low-pass filtering. On the right, Figure 5 illustrates how the diagonal filter is capable of removing diagonal structures in an image.

However, for denoising tasks this filter is not optimal. For images that do not consist exclusively of horizontal or vertical stripes it produces undesired artefacts in the reconstruction such as additional thin stripes in horizontal and vertical direction. Also, we obtain the same filter despite having used quite different training data (rotated by 90°) in both experiments.

As a consequence, we are going to focus in the following on distinguishing between different shapes or scales.

Figure 6 shows a variety of experiments aiming at finding an optimal filter function $\hat{h}$, which favours the specific texture, orientation or scale in input image $u^{+}$ and disfavours the one present in $u^{-}$. Again, we perform reconstruction according to (5) choosing $\eta=1$ and $\sigma^{2}=0.005$, where the given image $f$ is a noise-free combination of $u^{+}$ and $u^{-}$. We state the functional evaluations of $J(u^{+};\hat{h})$ and $J(u^{-};\hat{h})$ below the respective figures. It can be observed that the former are in most cases significantly smaller than the latter, confirming the usefulness of model (1) and Algorithm (4). In (a), we can clearly see that the left-hand-side can be almost perfectly reconstructed in $\hat{u}$ while undesired artefacts are occurring for the unfavourable horizontal stripes. Similar results can be observed in (b), where the optimal filter is sparser than the one in (a). In (c), we have diagonal stripes in different angles as input images. Again, the left-hand side is almost perfectly reconstructed whereas the stripes on the right-hand side appear blurred. In (d), the diagonal stripes are only one pixel thick and hence the appearance of the filter changes significantly. The regulariser is able to reconstruct the left-hand side of $f$ very well. We exchange $u^{+}$ and $u^{-}$ in (e) and (f). First, the circle is the desired and the square is the undesired input signal. The opposite case holds true for example (f). Here, again the diagonal-shaped filter performs best and blurs the circular structure enforcing edges in vertical and horizontal direction.

Sometimes it is necessary to find suitable filters by increasing their size, as can be seen in Figure 7 (a). On the right-hand side we can see the optimal filter calculated assuming a size of $5\times 5$. By increasing the size slightly to $7\times 7$, we obtain a much more reasonable filter being able to recover the circles. The filter itself almost resembles a circle with a radius of three pixels. In (b), given the ground truth image on the top left and four rotated versions of it in angles between 0 and 45 degrees as favourable input as well as five noise signals as negative input, we obtain the filter on the bottom right, which performs surprisingly well at denoising the image on the top right.

In Figure 8, convincing denoising results can be achieved for examples (a) and (c) in Figure 6, adding noise to the ground truth image $f$ and using the calculated optimal filter. One can clearly see that the structures on the left-hand sides are denoised better.

We would like to remark that the setup with one filter is indeed rather simple and cannot mimic 2D differential-based filters like TV or alike, but therefore it is even more surprising that the filters presented above perform really well.

5 Conclusions and Outlook

Starting from the model in [3], we derived a more generalised formulation suitable for minimisation with respect to multi-dimensional filter functions. In addition, our flexible framework allows for multiple desired and undesired input signals.

We were able to reproduce different common first- and second-order regularisers such as TV and TV2 in the 1D case. Furthermore, we created a new family of non-derivative-based regularisers suitable for specific types of images. Also, we showed that specific shapes such as diagonal stripes can be eliminated while applying such parametrised regularisers. In addition, we believe that our learning approach is suitable for distinguishing between different shapes, scales and textures. A great advantage and novelty is that we are able to include both favourable and unfavourable input signals in our framework.

Regarding numerical implementation, we would like to stress that for computational simplicity in combination with the CVX framework we have only considered Dirichlet boundary conditions so far, but will use different boundary conditions (like the more suitable Neumann boundary conditions) in the future.

We further assume that the ansatz of filter functions respectively convolutions as parametrisations for the regularisers is too generic especially for denoising tasks. It will be interesting to look into dictionary-based sparsity approaches, and to then learn basis functions with the presented quotient model.

Moreover, future work might include applications in biomedical imaging such as reconstruction in CT or MRI as well as denoising or object detection in light microscopy images. In [19], the authors present a multiscale segmentation method for circulating tumour cells, where they are able to detect cells of different sizes. Using our model, we believe that shape or texture priors incorporated in sparsity-based regularisers could be well improved. One possible application could be mitotic cell detection (cf. [11]).

6 Acknowledgements

MB acknowledges support from the Leverhulme Trust early career fellowship ”Learning from mistakes: a supervised feedback-loop for imaging applications” and the Newton Trust. GG acknowledges support by the Israel Science Foundation (grant 718/15). JSG acknowledges support by the NIHR Cambridge Biomedical Research Centre. CBS acknowledges support from Leverhulme Trust project ’Breaking the non-convexity barrier’, EPSRC grant ’EP/M00483X/1’, EPSRC centre ’EP/N014588/1’, the Cantab Capital Institute for the Mathematics of Information, and from CHiPS (Horizon 2020 RISE project grant).

Data Statement.

The corresponding MATLAB® code is publicly available on Apollo - University of Cambridge Repository (https://doi.org/10.17863/CAM.8419).