Networks for Nonlinear Diffusion Problems in Imaging

A multitude of imaging and vision tasks have seen recently a major transformation by deep learning methods and in particular by the application of convolutional neural networks. These methods achieve impressive results, even for applications where it is not apparent that convolutions are suited to capture the underlying physics. In this work, we develop a network architecture based on nonlinear diffusion processes, named DiffNet. By design, we obtain a nonlinear network architecture that is well suited for diffusion-related problems in imaging. Furthermore, the performed updates are explicit, by which we obtain better interpretability and generalisability compared to classical convolutional neural network architectures. The performance of DiffNet is tested on the inverse problem of nonlinear diffusion with the Perona–Malik filter on the STL-10 image dataset. We obtain competitive results to the established U-Net architecture, with a fraction of parameters and necessary training data.


Introduction
We are currently undergoing a paradigm shift in imaging and vision tasks from classical analytic to learning and data based methods.In particular by deep learning and the application of convolutional neural networks (CNN).Whereas highly superior results are obtained, interpretability and analysis of the involved processes is a challenging and ongoing task [18,22].
In a general setting, Deep Learning proposes to develop a non-linear mapping A Θ : X → Y between elements of two spaces X, Y (which may be the same) parametrised by a finite set of parameters Θ, which need to be learned: This learning based approach is in marked contrast to classical methods with a physical interpretation of the process (1.1) for both computational modelling of data given a physical model (which we call a forward problem), and the estimation of parameters of a physical model from measured, usually noisy, data (which we call an inverse problem).Several recent papers have discussed the application of Deep Learning in both forward [21,27,32,33,36,37] and inverse [19,20,25,38,39] problems.Several questions become essential when applying learned models to both forward and inverse problems including: • how and to what extent can learned models replace physical models?
• how do learned models depend on training protocols and how well do they generalise?• what are appropriate architectures for the learned models, what is the size of the parameter set Θ that needs to be learned, and how can these be interpreted?
Date: November 30, 2018.S. Arridge and A. Hauptmann are with the Department of Computer Science; University College London, London, United Kingdom.
In this paper we aim to answer some of these questions, by taking a few steps back and looking at an analytic motivation for network architectures.Here we consider in particular mappings between images u in dimension d, i.e.X = Y = L p (Ω ⊂ R d ), for which there exist several widely used linear and nonlinear mappings A defined by differential and/or integral operators.For example a general linear integral transform such as (1.2) u obs (x) = (Au true )(x) = Ω K(x, y)u true (y)dy includes stationary convolution as a special case if the kernel is translation invariant, i.e.K(x, y) ≡ K(x − y).If the kernel depends on the image, i.e.K(x, y) ≡ K(x, y; u) then (1.2) becomes nonlinear.Alternatively, the forward model may be modelled as the end-point of an evolution process which becomes nonlinear if the kernel depends on the image state K(x, y, t) ≡ K(x, y; u(t)); see eq.(2.7) below for a specific example.Furthermore, a widely studied class of image mappings is characterised by defining the evolution through a partial differential equation (PDE) where flow is generated by the local structure of u [23,30,34], such that (1.3) u t = F u, ∇u, ∂ 2 u ∂x i ∂x j , . . . .These problems in general do not admit an explicit integral transform representation and are solved instead by numerical techniques, but since they depend on a physical model, described by the underlying PDE, they can be analysed and understood thoroughly [35].This also includes a statistical interpretation [17] as hyperpriors in a Bayesian setting [4].Furthermore, models like (1.3) can be used to develop priors for nonlinear inverse problems, such as optical tomography [7,13] and electrical impedance tomography [11].Motivated by the success of these analytic methods to imaging problems in the past, we propose to combine physical models with data driven methods to formulate network architectures for solving both forward and inverse problems that take the underlying physics into account.We limit ourselves to the case where the physical model is of diffusion type, although more general models could be considered in the future.The leading incentive is given by the observation that the underlying processes in a neural network do not need to be limited to convolutions.
Similar ideas of combining partial differential equations with deep learning have been considered earlier.For instance learning coefficients of a PDE via optimal control [24], as well as deriving CNN architectures motivated by diffusion processes [5], deriving stable architectures by drawing connections to ordinary differential equations [9] and constraining CNNs [29] by the interpretation as a partial differential equation.Another interpretation of our approach can be seen as introducing the imaging model into the network architecture; such approaches have led to a major improvement in reconstruction quality for tomographic problems [2,12,15].This paper is structured as following.In section 2 we review some theoretical aspects of diffusion processes for imaging and the inversion based on theory of partial differential equations and differential operators.We formulate the underlying conjecture for our network architecture, that the diffusion process can be inverted by a set of local non-stationary filters.In the following we introduce the notion of continuum networks in section 3 and formally define the underlying layer operator needed to formulate network architectures in a continuum setting.We draw connections to the established convolutional neural networks in our continuum setting.We then proceed to define the proposed layer operators for diffusion networks in section 4 and derive an implementable architecture by discretising the involved differential operator.In particular we derive a network architecture that is capable of reproducing inverse filtering with regularisation for the inversion of nonlinear diffusion processes.We examine the reconstruction quality of the proposed DiffNet in the following section 5 for an illustrative example of deconvolution and the challenging inverse problem of inverting nonlinear diffusion with the Perona-Malik filter.We achieve results that are competitive to popular CNN architectures with a fraction of the amount of parameters and training data.Furthermore, all computed components that are involved in the update process are interpretable and can be analysed empirically.In section 6 we examine the generalisablity of the proposed network with respect to necessary training data.Additionally, we empirically analyse the obtained filters and test our underlying conjecture.Section 7 presents some conclusions and further ideas.

Diffusion and flow processes for imaging
In the following we want to explore the possibility to include a nonlinear process as an underlying model for network architectures.Specifically, the motivation for this study is given by diffusion processes that have been widely used in imaging and vision.Let us consider here the general diffusion process in R d ; then on a fixed time interval with some diffusivity γ, to be defined later, we have Remark 2.1.When considering bounded domains Ω ⊂ R d we will augment (2.1) with boundary conditions on ∂Ω.We return to this point in section 4.
In the following we denote the spatial derivative by Let us first consider the isotropic diffusion case, then the differential operator becomes the spatial Laplacian L(γ = 1) = ∆.In this case, the solution of (2.1) at time T is given by convolution with a Green's function where 4T in dimension d and we recall that the convolution of two functions f, g Definition 2.2.The Green's operator is defined with the Green's function as its kernel In the general case for an anisotropic diffusion flow (ADF) we are interested in a scalar diffusivity γ ∈ [0, 1] that depends on u itself, i.e. (2.6) This is now an example of a non-linear evolution where K ADF (x, y, u(x, t)) is now a non-stationary, non-linear and time-dependent kernel.In general there is no explicit expression for K ADF and numerical methods are required for the solution of (2.6).
Remark 2.3.Not considered here, but a possible extension to (2.6) is where γ is a tensor, which, for d = 2 takes the form Furthermore, extensions exist for the case where u is vector or tensor-valued.We do not consider these cases here, see [34] for an overview.
2.1.Forward Solvers.First of all, let us establish a process between two states of the function u.Integrating over time from t = t 0 to t = t 1 = t 0 + δt yields Note, that the left hand side can be expressed as t1 t0 ∂ t u(x, t) dt = u(x, t 1 )−u(x, t 0 ) and we denote the right hand side by an integral operator A δt (γ), such that (2.8) In the following we denote the solution of (2.1) at time instances t n as u (n) = u(x, t n ).Then we can establish a relation between two time instances of u by (2.9) where Id denotes the identity and t n+1 = t n + δt.
Since we cannot compute u (n+1) by (2.9) without the explicit knowledge of the (possibly time-dependent) diffusivity γ, it is helpful to rather consider a fixed diffusivity at each time instance γ (n) = γ(x, t = t n ), or γ (n) = γ(u(x, t = t n )) in the nonlinear case; then by using the differential operator (2.2) we have an approximation of (2.8) by We can now solve (2.6) approximately by iterating for time steps δt using either an explicit scheme (2.10) and coincide with the integral formulation of (2.9).It's also useful to look at the Green's functions solutions.
Lemma 2.4.Consider the isotropic case with γ ≡ 1.Then we may write, with time steps δt = T /N , (2.13) and use the convolution theorem to give Let us also note that in Fourier domain, by Taylor series expansion we have and therefore in the spatial domain the finite difference step and the Gaussian convolution step are the same 2.2.Inverse Filtering.Let us now consider the inverse problem of reversing the diffusion process.That is we have u T and aim to recover the initial condition u 0 .This is a typical ill-posed problem as we discuss in the following.
2.2.1.Isotropic case γ ≡ 1.As the forward problem is represented as convolution in the spatial domain, the inverse mapping u T → u 0 is a (stationary) deconvolution.
We remind that ûT = e −k 2 T û0 (k), then the inversion is formally given by division in the Fourier domain as However we note: i.) The factor e k 2 T is unbounded and hence the equivalent convolution kernel in the spatial domain does not exist.ii.) (2.15) is unstable in the presence of even a small amount of additive noise and hence it has to be regularised in practice.
Nevertheless, let's consider formally with e k 2 T = e k 2 δt N that by Taylor series we get ) .Motivated by this we define an operator for the inversion process δt coincides with the inverse of the implicit update in (2.11), and (2.17) is an estimate for the deconvolution problem which (in the absence of noise) is correct in the limit 1 Note the definition of Fourier Transform is chosen to give the correct normalisation so that 2.2.2.Anisotropic case.In this case the diffused function is given by (2.7).Following Lemma 2.4 we may put )u 0 and we also have Remark 2.6 (Unsharp Masking).We recall that a simple method for "deconvolution" is called Unsharp Masking which is usually considered as for some blur value σ and sufficiently small .By similar methods as above, we find We may choose to interpret the operators E δt (ζ) as a kind of "non-stationary unsharp masking".

Discretisation. We introduce the definition of a sparse matrix operator representing local non-stationary convolution
Definition 2.7.W is called a Sparse Sub-Diagonal (SSD) matrix if its non-zero entries are all on sub-diagonals corresponding to the local neighbourhood of pixels on its diagonal.In the following we restrict ourselves to the typical 4-connected neighbourhood of pixels in dimension d = 2.For the numerical implementation we have the Laplacian stencil from which we have that L ∆ is zero-mean.Similarly we will have for the numerical approximation of E iso δt the matrix operator Further we conjecture that in the numerical setting E δt (ζ) is approximated by the sum of identity plus a SSD matrix operator as (2.23) In the presence of noise, the inverse operator E δt (ζ) can be explicitly regularised by addition of a smoothing operation Whereas in classical approaches to inverse filtering the regularisation operator would be defined a priori, the approach in this paper is to learn the operator W and interpret it as the sum of a differentiating operator L and a (learned) regulariser S.This is discussed further in section 4 below

Continuum Networks
Motivated by the previous section, we aim to build network architectures based on diffusion processes.We first discuss the notion of (neural) networks in a continuum setting for which we introduce the concept of a continuum network as a mapping between function spaces.That is, given a function on a bounded domain Ω ⊂ R d with f ∈ L p (Ω), we are interested in finding a non-linear parametrised operator H Θ : L p (Ω) → L p (Ω) acting on the function f .We will consider in the following the case p ∈ {1, 2}; extensions to other spaces depend on the involved operations and will be the subject of future studies.
We will proceed by defining the essential building blocks of a continuum network and thence to discuss specific choices to obtain a continuum version of the most common convolutional neural networks.Based on this we will then introduce our proposed architecture as a diffusion network in the next chapter.

3.1.
Formulating a continuum network.The essential building blocks of a deep neural network are obviously the several layers of neurons, but since these have a specific notion in classical neural networks, see for instance [31], we will not use the term of neurons to avoid confusion.We rather introduce the concept of layers and channels as the building blocks of a continuum network.In this construction, each layer consists of a set of functions on a product space and each function represents a channel.Definition 3.1 (Layer and channels).
Then we call: F k the layer k with I channels and corresponding index set I.
The continuum network is then built by defining a relation or operation between layers.In the most general sense we define the concept of a layer operator for this task.Definition 3.2 (Layer operator).Given two layers F k and F t , k = t, with channel index set I, J, respectively, we call the mapping H : I L p (Ω) → J L p (Ω) with HF k = F t a layer operator.If the layer operator depends on a set of parameters Θ, then we write H Θ .
We note that for simplicity, we will not index the set of parameters, i.e Θ generally stands for both all involved parameters of each layer separately, or the whole network.The classical structure for layer operators follows the principle of affine linear transformations followed by a nonlinear operation.Ideally the affine linear transformation should be parameterisable by a few parameters, whereas the nonlinear operation is often fixed and acts pointwise.A popular choice is the maximum operator also called the "Rectified Linear Unit": The continuum network is then given by the composition of all involved layer functions.For example in monochromatic imaging applications we typically have an input image f 0 and a desired output f K with several layer functions inbetween, that perform a specific task such as denoising or sharpening.In this case the input and output consists of one channel, i.e.

Continuum convolutional networks.
Let us now proceed to discuss a specific choice for the layer operator, namely convolutions.Such that we will obtain a continuum version of the widely used convolutional neural networks, which we will call here a continuum convolutional network, to avoid confusion with the established convolutional neural networks (CNN).We note that similar ideas have been addressed as well in [2].
Let us further consider linearly ordered network architectures, that means each layer operator maps between consecutive layers.The essential layer operator for a continuum convolutional network is then given by the following definition.Definition 3.3 (Convolutional layer operator).Given two layers F k−1 and F k with channel index set I, J, respectively.Let b j ∈ R and ω i,j ∈ L p (Ω), with compact support in Ω, be the layer operator's parameters for all i ∈ I, j ∈ J.We call C (k) Θ,ϕ the convolutional layer operator for layer k, if for each output channel with a point-wise nonlinear operator ϕ : L p (Ω) → L p (Ω).
If the layer operator does not include a nonlinearity, we write C Θ,Id .Now we can introduce the simplest convolutional network architecture by applying K ≥ 1 convolutional layer operators consecutively.Definition 3.4 (K-layer Continuum Convolutional Network).Let K ≥ 1, then we call the composition of K convolutional layer operator, denoted by C K Θ , a K-layer Continuum Convolutional Network, such that In the following we will also refer to a K-layer CNN as the practical implementation of a K-layer continuum convolutional network.A popular network architecture that extends this simple idea is given by a resdiual network (ResNet) [16], that is based on the repetition of a 2-layer CNN with a residual connection, that consists of addition.That is, the network learns a series of additive updates to the input.The underlying structure in ResNet is the repeated application of the following residual block, given by Θ,ϕ + Id, R Θ,ϕ F 0 = F 2 .Note that the second convolutional layer does not include a nonlinearity.Furthermore, it is necessary that |F 0 | = |F 2 |, but typically it is often chosen such that The full continuum ResNet architecture can then be summarized as follows.Let K ≥ 1, then the composition of K residual blocks, denoted by Θ,ϕ .Note that the two additional convolutional layers in the beginning and end are necessary to raise the cardinality of the input/output layer to the cardinality needed in the residual blocks.A complete K-block ResNet then consists of 2(K + 1) layers.Note that in the original work [16], the network was primarily designed for an image classification task rather than image-to-image mapping that we consider here.

DiffNet: discretisation and implementation
Next we want to establish a layer operator based on the diffusion processes discussed in chapter 2. This means that we now interpret the layers F k of the continuum network as time-states of the function u : Ω × R + → R, where u is a solution of the diffusion equation (2.1).In the following we assume single channel networks, i.e. |F k | = 1 for all layers.Then we can associate each layer with the solution u such that F k = u (k) = u(x, t = t k ).To build a network architecture based on the continuum setting, we introduce the layer operator versions of (2.10), and (2.21): Definition 4.1 (Diffusion and filtering layer operator).Given two layers F k and , then a diffusion layer operator D Θ , with parameters Θ = {γ, δt}, is given by Similarly, an inverse filtering layer operator, with parameters Θ = {ζ, δt} is given by Note that this formulation includes a learnable time-step and hence the time instances that each layer represents changes.That also means that a stable step size is implicitly learned, if there are enough layers.In the following we discuss a few options on the implementation of the above layer operator, depending on the type of diffusivity.4.1.Discretisation of a continuum network.Let us briefly discuss some aspects on the discretisation of a continuum network; let us first start with affine linear networks, such as the convolutional networks discussed in section 3.2.Rather than discussing the computational implementation of a CNN, (see e.g. the comprehensive description in [8]), we concentrate instead on an algebraic matrix-vector formulation that serves our purposes.
For simplicity we concentrate on the two-dimensional d = 2 case here.Let us then assume that the functions f i in each layer are represented as a square n-by-n image and we denote the vectorised form as f ∈ R n 2 .Then any linear operation on f can be represented by some matrix A; in particular we can represent convolutions as a matrix.
We now proceed by rewriting the convolutional operation (3.1) in matrix-vector notation.Given two layers F k and F k−1 with channel index set I, J, respectively, then we can represent the input layer as vector F k−1 ∈ R Jn 2 and similarly the output layer F k ∈ R In 2 .Let A i ∈ R n 2 ×Jn 2 represent the sum of convolutions in (3.1), then we can write the layer operator in the discrete setting as matrix-vector operation by , where Id denotes the identity matrix of suitable size.Furthermore, following the above notation, we can introduce a stacked matrix A ∈ R In 2 ×Jn 2 consisting of all A i and a matrix B ∈ R In 2 ×Jn 2 by stacking the diagonal matrices b i Id.Then we can represent the whole convolutional layer in the discrete setting as Now the parameters of each layer are contained in the two matrices A and B.

4.2.
Learned Forward and Inverse Operators.Let us now discuss a similar construction for the diffusion layers.For the implementation of the diffusion network, we consider the explicit formulation in (2.10) with the differential operator L(γ (k) ) = ∇ • γ (k) ∇ approximated by the stencil (including the time-step δt) The basis of learning a diffusion network is now given as estimating the diagonals of L δt (γ) and the time-step δt.This can be done either explicitly as for the CNN or indirectly by an estimator network, as we will discuss next.

Formulating DiffNet.
Let us first note, that if we construct our network by strictly following the update in (4.6), we restrict ourselves to the forward diffusion.
To generalise to the inverse problem we consider Additionally, there are two fundamentally different cases for the diffusivity γ we need to consider before formulating a network architecture to capture the underlying behaviour.These two cases are i.) Linear diffusion; spatially varying and possible time dependence, γ = γ(x, t).
ii.) Nonlinear diffusion; diffusivity depending on the solution u, γ = γ(u(x, t)).In the first case, we could simply to try learn the diffusivity explicitly, to reproduce the diffusion process.In the second case, this is not possible and hence an estimation of the diffusivity needs to be performed separately in each time step from the image itself, before the diffusion step can be performed.This leads to two conceptually different network architectures.
The linear case i.) corresponds to the diffusion layer operator (4.1) and is aimed to reproduce a linear diffusion process with fixed diffusivity.Thus, learning the mean-free filter suffices to capture the physics.The resulting network architecture is outlined in Figure 1.Here, the learned filters can be directly interpreted as the diffusivity of layer k and are then applied to F k−1 to produce F k .
Nonlinear Diffusion Network (DiffNet) ζ (1)   W δt (ζ (2) ) W δt (ζ (3) ) Here the filters ζ are implicitly estimated by a small k-layer CNN and then applied to the image in the filtering layer.
In the nonlinear case ii.) we follow the same update structure, but now the filters are not learned explicitly, they are rather estimated from the input layer itself, as illustrated in Figure 2. Furthermore, since this architecture is designed for inversion of the nonlinear diffusion process, we employ the generalised stencil W δt (ζ).Then, given layer F k , the filters ζ are estimated by a small CNN from F k−1 , which are then applied following an explicit update as in (4.6) to produce F k .Note that the diagonals in the update matrix are produced by the estimation CNN.We will refer to this nonlinear filtering architecture as the DiffNet under consideration in the rest of this study.4.3.1.Implementation.The essential part for the implementation of a diffusion network is to perform the update (4.6) with either L δt (γ) or W δt (ζ).For computational reasons it is not practical to build the sparse diagonal matrix and evaluate (4.6), we rather represent the filters γ and ζ as an n × n-image and apply the filters as pointwise matrix-matrix multiplication to a shifted and cropped image, according to the position in the stencil.This way, the zero Neumann boundary condition (4.5) is also automatically incorporated.
For the linear diffusion network, we would need to learn the parameter set Θ, consisting of filters and time steps, explicitly.This has the advantage of learning a global operation on the image where all parameters are interpretable, but it comes with a few disadvantages.First of all, in this form we are limited to linear diffusion processes and a fixed image size.Furthermore, the parameters grow with the image size, i.e. for an image of size n × n we need 5n 2 parameters per layer.Thus, applications may be limited.
For the nonlinear architecture of DiffNet, where the filters depend on the image at each time step, we introduced an additional estimator consisting of a K-layer CNN.This CNN gets an image, given as layer F k , as input and estimates the filters ζ.The architecture for this K-layer CNN as estimator is chosen to be rather simplistic, as illustrated in Figure 3.The input F k consists of one channel, which is processed by K − 1 convolutional layers with 32 channels and a ReLU nonlinearity, followed by the last layer without nonlinearity and 5 channels for each filter, which are represented as matrices of the same size as the input F k .In particular for the chosen filter size of 3 × 3 we have exactly 9 • (32 + 32 • 5 + 32 2 • (K − 2)) convolutional parameters and 32 • (K − 1) + 5 biases per diffusion layer.That is for a 5-layer CNN 29.509 parameters independent of image size.K-layer CNN for filter estimation Figure 3. Architecture of the K-layer CNN used as diffusivity estimator in the nonlinear diffusion network (DiffNet).

Computational experiments
In the following we will examine the reconstruction capabilities of the proposed DiffNet.The experiments are divided into a simple case of deconvolution, where we can examine the learned features, and a more challenging problem of recovering an image from its nonlinear diffused and noise corrupted version.The forward problem is given by (2.1) with zero Neumann boundary condition (4.5) and constant diffusivity γ ≡ 1.For the experiment we choose T = 1, which results in a small uniform blurring, as shown in Figure 4. We remind that for the isotropic diffusion, the forward model is equivalent to convolution in space with the kernel G √ 2T , see (2.3).As it is also illustrated in Figure 4, convolution in the spatial domain is equivalent to multiplication in Fourier domain.In particular, high frequencies get damped and the convolved image is dominated by low frequencies.Hence, for the reconstruction task without noise, we essentially need to recover the high frequencies.
The training and test data for DiffNet consists of simple disks of varying radius and contrast.The training set consists of 1024 samples and the test set of an additional 128, each of size 64 × 64.The network architecture is chosen following the schematic in Figure 2, with 3 diffusion layers and a final projection to the positive numbers by a ReLU layer.The filter estimator is given by a 4-layer CNN, as described in section 4.3.1.All networks were implemented in Python with TensorFlow [1].
The input to the network is given by the convolved image without noise and we have minimised the 2 -loss of the output to the ground-truth image.The optimisation is performed for about 1000 epochs in batches of 16 with the Adam algorithm and initial learning rate of 4 • 10 −4 and a gradual decrease to 10 −6 .Training on a single Nvidia Titan Xp GPU takes about 24 minutes.The final training and test error is both at a PSNR of 86.24, which corresponds to a relative 2 -error of 2.5 • 10 −4 .We remind that this experiment was performed without noise.
The result of the network and intermediate updates for one example from the test data are illustrated in Figure 5.We also show the filters ζ (k) computed as the output of the trained CNN in each layer, k = 1 . . .3. The output of the last diffusion layer F 3 is additionally processed by a ReLU layer to enforce positivity in the final result.It can be seen that the network gradually reintroduces the high frequencies in the Fourier domain; especially the last layer mainly reintroduces the high frequencies to the reconstruction.It is interesting to see, that the learned filters follow indeed the convention that the central filter is of different sign than the directional filters.This enforces the assumption that the filter consists of a mean free part and a regularising part, which should be small in this case, since we do not have any noise in the data.Lastly, we note that the final layer, before projection to the positive numbers, has a clearly negative part around the target, which will be cut off resulting in a sharp reconstruction of the ball.

Nonlinear diffusion.
Let us now consider the nonlinear diffusion process with the Perona-Malik filter function [26] for (2.1) with zero Neumann boundary condition (4.5).In this model the diffusivity is given as a function of the gradient with contrast parameter λ > 0. We mainly concentrate here on the inverse problem of restoring an image that has been diffused with the Perona-Malik filter and contaminated by noise.
For the experiments we have used the test data from the STL-10 database [6], which consists of 100,000 RGB images with resolution 96 × 96.These images have been converted to grayscale and divided to 90,000 for training and 10,000 for testing.The obtained images were then diffused for 4 time steps with δt = 0.1 and λ = 0.2.A few sample images from the test data with the result of the diffusion are displayed in Figure 6.The task is then to revert the diffusion process with additional regularisation to deal with noise in the data.
For all experiments we have used the same network architecture of DiffNet using the architecture illustrated in Figure 2. By performing initial tests on the inversion without noise, we have found that 5 diffusion layers with a 4-layer CNN, following performance.Additionally, we have used a ReLU layer at the end to enforce nonnegativity of the output, similarly to the last experiment.We emphasise that this architecture was used for all experiments and hence some improvements for the high noise cases might be expected with more layers.All networks were trained for 18 epochs, with a batch size of 16, and 2 -loss.For the optimisation we have used the Adam algorithm with initial learning rate of 2 • 10 −3 and a gradual decrease to 4 • 10 −6 .Training on a single Nvidia Titan Xp GPU takes about 75 minutes.
As benchmark, we have performed the same experiments with a widely used network architecture known as U-Net [28].This architecture has been widely applied in inverse problems [3,10,19,20] and is mainly used to post-process initial directly reconstructed images from undersampled or noisy data.For instance by filtered back-projection in X-ray tomography or the inverse Fourier transform in magnetic resonance imaging [14].The network architecture we are using follows the publication [19] and differs from the original mainly by a residual connection at the end.That means the network is trained to remove noise and undersampling artefacts from the initial reconstruction.In our context, the network needs to learn how to remove noise and reintroduce edges.For training we have followed a similar protocol as for DiffNet.The only difference is that we started with an initial learning rate of 5 The reconstruction results, for some samples of the test data, with DiffNet can be seen in Figure 6 for the case without noise and in Figure 7 for 1% noise on the diffused images.A comparison of the reconstruction results with U-Net and DiffNet are shown in Figure 8 for the test without noise and in Figure 8 for 1% noise.Qualitatively, the reconstructed images are very similar, as can be seen in the residual images in the last column.The leftover noise pattern for both networks is concentrated on the fine structures of the ship.Quantitatively, for the noise free experiment DiffNet has an increase of 4dB in PSNR compared to the result of U-Net, 65.34 (DiffNet) compared to 61.08 (U-Net).For the case with 1% noise, the quantitative measures are very similar.Here U-Net has a slightly higher PSNR with 35.27 compared to DiffNet with 34.96.A thorough study of reconstruction quality of both networks follows in the next section as well as some interpretation of the learned features in DiffNet.

Discussion
First of all we note that the updates in DiffNet are performed explicitly and that the CNN in the architecture is only used to produce the filters ζ.This means that DiffNet needs to learn a problem specific processing, in contrast to a purely data driven processing in a CNN.Consequently, the amount of necessary learnable parameters is much lower.For instance the 5-layer DiffNet used for inversion of the nonlinear diffusion in section 5. ∼ 0.3% of parameters compared to U-Net and hence the learned features can be seen to be much more explicit.In the following we discuss some aspects that arise from this observation, such as generalisability and interpretability.Results for the four scenarios are shown in Figure 10.Most notably DiffNet outperforms U-Net clearly for the forward problem and the noise free inversion, by 4dB and 3dB, respectively.For the noisy cases both networks perform very similar for the full training data size of 90,000.The biggest difference overall is that DiffNet achieves its maximum test error already with 500-1,000 samples independent of the noise case, whereas the U-Net test error saturates earlier with higher noise.In conclusion we can say, that for the noisy cases both networks are very comparable in reconstruction quality, but for small amounts of data the explicit nature of DiffNet is clearly superior.6.2.Interpretation of learned filters.Since all updates are performed explicitly with the output from the filter estimation CNN, we can interpret some of the learned features.For this purpose we show the filters for the ship image from section 5.2 for the three inversion scenarios under consideration.In Figure 11 we show the sum of all learned filters, i.e.Additionally, we also show the diagonal filters ζ 5 in Figure 12.
We would expect that with increasing noise level, the filters will incorporate more smoothing to deal with the noise; this implies that the edges get wider with increasing noise level.This can be nicely observed for the diagonal filters in Figure 12.For the smoothing in Figure 11, we see that the first layer consists of broader details and edges that are refined in the noise free case for increasing layers.In the noisy case the later layers include some smooth features that might depict the regularisation necessary in the inversion procedure.It is generally interesting to observe, that the final layer shows very fine local details, necessary to restore fine details for the final output.
Finally, we have computed training data of a wider noise range to examine the regularisation properties of the learned network.For this we have taken the full 90,000 training samples and contaminated the diffused image with noise ranging from 0.01% to 10% noise.As we conjectured in section 2.3, the learned update the deviation from the mean-free part.Furthermore, we interpret the smoothing level α as the mean of the absolute value of S. The resulting graph of smoothing versus noise level is shown in Figure 13.As we have conjectured, the estimate of α increases clearly with the noise level and hence we believe our interpretation of the learned filters as the composition of a mean-free part and a smoothing necessary for ill-posed inverse problems is valid.
Increasing smoothing with noise level

Conclusions
In this paper we have explored the possibility to establish novel network architectures based on physical models other than convolutions; in particular we concentrated here on diffusion processes.As main contributions, we have introduced some non-linear forward mappings, modelled through learning rather than just through PDEs or integral transforms.We have reviewed (regularised) inverse diffusion processes for inverting such maps.In particular, we have conjectured that these inverse diffusion processes can be represented by local non-stationary filters, which can be learned in a network architecture.More specific, these local filters can be represented by a sparse sub-diagonal (SSD) matrix and hence efficiently used in the discrete setting of a neural network.We emphasise that even though we have concentrated this study on a specific structure for these SSD matrices based on diffusion, other (higher order) models can be considered.
We obtain higher interpretability of the network architecture, since the image processing is explicitly performed by the application of the SSD matrices.Consequently, this means that only a fraction of parameters is needed in comparison to classical CNN architectures to obtain similar reconstruction results.We believe that the presented framework and the proposed network architectures can be useful for learning physical models in the context of imaging and inverse problems, especially, where a physical interpretation of the learned features is crucial to establish confidence in the imaging task.

Furthermore
we are going to consider that a class of SSD matrices W(ζ) with learned parameters ζ can be decomposed as W(ζ) = S(ζ)+L(ζ) where S is smoothing and L(ζ) is zero-mean; i.e.L(ζ) has one zero eigenvalue such that its application to a constant image gives a zero valued image L(ζ)1 = 0

Remark 4 . 2 .
The assumption of a single channel network, i.e. |F k | = 1 for all k, can be relaxed easily.Either by assuming |F k | = m for some m ∈ N and all layers, or by introducing a channel mixing as in the convolutional operator (3.1).

Figure 2 .
Figure 2. Illustration for a nonlinear 3-layer Diffusion network.Here the filters ζ are implicitly estimated by a small k-layer CNN and then applied to the image in the filtering layer.

5. 1 .Figure 4 .
Figure 4. Illustration of the deconvolution problem for a simple ball.Top row shows the image space and the bottom row shows the corresponding absolute value of the Fourier coefficients.All images are plotted on their own scale.

Figure 5 .Figure 6 .Figure 7 .
Figure 5. Illustration of the deconvolution process with three layers of DiffNet.The left column shows the processed image and intermediate steps.The right column shows the corresponding absolute value of Fourier coefficients.All images are plotted on their own scale.

Figure 8 .
Figure 8.Comparison of reconstruction quality for reconstruction from nonlinear diffused image without noise.Both networks are trained on the full set of 90,000 images.PNSR: DiffNet 65.34, U-Net 61.08

Figure 9 .
Figure 9.Comparison of reconstruction quality for reconstruction from 1% noise contaminated nonlinear diffused image.Both networks are trained on the full set of 90,000 images.PNSR: DiffNet 34.96, U-Net 35.27

Figure 10 .
Figure 10.Generalisation plot for the forward and inverse problem of nonlinear diffusion and varying noise levels.Test error depending on the amount of training data, for both DiffNet and U-Net.

4 i=1 ζ i − ζ 5 .
If the network would only learn the mean-free differentiating part, then these images should be zero.That implies that the illustrated filters in Figure11can be related to the learned regularisation S(ζ).

Figure 11 .Figure 12 .
Figure 11.Filter updates 4 i=1 ζ i − ζ 5 for different noise levels.Each image is displayed on its own scale.

Figure 13 .
Figure 13.Estimate of the smoothing level α for increasing noise in the inverse problem.Computed over a sample of 32 images from the test data.
To test the generalisation properties of the proposed DiffNet we have performed similar experiments as in section 5.2 for nonlinear diffusion, but with increasing amounts of training data.Under the assumption that DiffNet learns a more explicit update than a classic CNN, we would expect also to require less training data to achieve a good test error.To certify this assumption, we have examined 4 settings of nonlinear diffusion with the Perona-Malik filter: learning the forward model, learning to reconstruct from the diffused image without noise, as well as with 0.1% and 1% noise.We then created training data sets of increasing size from just 10 samples up the full size of 90,000.For all scenarios we have trained DiffNet and U-Net following the training protocol described in 5.2.Additionally, we have aborted the procedure when the networks started to clearly overfit the training data. 6.1.Generalisability.