Networks for Nonlinear Diffusion Problems in Imaging
Abstract.
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 tested on the inverse problem of nonlinear diffusion with the PeronaMalik filter on the STL10 image dataset. We obtain competitive results to the established UNet architecture, with a fraction of parameters and necessary training data.
1. 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 nonlinear mapping between elements of two spaces , (which may be the same) parametrised by a finite set of parameters , which need to be learned:
(1.1) 
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?
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 in dimension , i.e. , for which there exist several widely used linear and nonlinear mappings defined by differential and/or integral operators. For example a general linear integral transform such as
(1.2) 
includes stationary convolution as a special case if the kernel is translation invariant, i.e. . If the kernel depends on the image, i.e. then (1.2) becomes nonlinear. Alternatively, the forward model may be modelled as the endpoint of an evolution process which becomes nonlinear if the kernel depends on the image state ; 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 [23, 30, 34], such that
(1.3) 
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 nonstationary 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 PeronaMalik 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.
2. 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 ; then on a fixed time interval with some diffusivity , to be defined later, we have
(2.1) 
Remark 2.1.
In the following we denote the spatial derivative by
(2.2) 
Let us first consider the isotropic diffusion case, then the differential operator becomes the spatial Laplacian . In this case, the solution of (2.1) at time is given by convolution with a Green’s function
(2.3) 
where in dimension and we recall that the convolution of two functions is defined by
(2.4) 
Definition 2.2.
The Green’s operator is defined with the Green’s function as its kernel
(2.5) 
by which we have
In the general case for an anisotropic diffusion flow (ADF) we are interested in a scalar diffusivity that depends on itself, i.e.
(2.6) 
This is now an example of a nonlinear evolution
(2.7) 
where is now a nonstationary, nonlinear and timedependent kernel. In general there is no explicit expression for and numerical methods are required for the solution of (2.6).
Remark 2.3.
2.1. Forward Solvers
First of all, let us establish a process between two states of the function . Integrating over time from to yields
Note, that the left hand side can be expressed as and we denote the right hand side by an integral operator , such that
(2.8) 
In the following we denote the solution of (2.1) at time instances as . Then we can establish a relation between two time instances of by
(2.9) 
where denotes the identity and .
Since we cannot compute by (2.9) without the explicit knowledge of the (possibly timedependent) diffusivity , it is helpful to rather consider a fixed diffusivity at each time instance , or 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 using either an explicit scheme
(2.10) 
or an implicit scheme
(2.11) 
Whereas (2.10) is stable only if CFL conditions are satisfied and (2.11) is unconditionally stable, they are both only accurate for sufficiently small steps . In fact, by the Neumann series, the schemes are equivalent in the limit
(2.12) 
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 . Then we may write, with time steps ,
(2.13) 
Proof.
Take the Fourier Transform^{1}^{1}1Note the definition of Fourier Transform is chosen to give the correct normalisation so that 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.14) 
2.2. Inverse Filtering
Let us now consider the inverse problem of reversing the diffusion process. That is we have and aim to recover the initial condition . This is a typical illposed problem as we discuss in the following.
2.2.1. Isotropic case
As the forward problem is represented as convolution in the spatial domain, the inverse mapping is a (stationary) deconvolution. We remind that , then the inversion is formally given by division in the Fourier domain as
(2.15) 
However we note:

The factor is unbounded and hence the equivalent convolution kernel in the spatial domain does not exist.

(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 that by Taylor series we get
Motivated by this we define an operator for the inversion process
(2.16) 
Clearly 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
(2.18) 
2.2.2. Anisotropic case
In this case the diffused function is given by (2.7). Following Lemma 2.4 we may put
(2.19) 
and we also have
(2.20) 
Conjecture 2.5.
There exists a set of local (nonstationary) filters where
(2.21) 
and where has only local support and such that
(2.22) 
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 as a kind of "nonstationary unsharp masking".
2.3. Discretisation
We introduce the definition of a sparse matrix operator representing local nonstationary convolution
Definition 2.7.
is called a Sparse SubDiagonal (SSD) matrix if its nonzero entries are all on subdiagonals corresponding to the local neighbourhood of pixels on its diagonal.
Furthermore we are going to consider that a class of SSD matrices with learned parameters can be decomposed as where is smoothing and is zeromean; i.e. has one zero eigenvalue such that its application to a constant image gives a zero valued image
In the following we restrict ourselves to the typical 4connected neighbourhood of pixels in dimension . For the numerical implementation we have the Laplacian stencil
from which we have that is zeromean. Similarly we will have for the numerical approximation of the matrix operator
Further we conjecture that in the numerical setting is approximated by the sum of identity plus a SSD matrix operator as
(2.23) 
In the presence of noise, the inverse operator can be explicitly regularised by addition of a smoothing operation
(2.24) 
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 and interpret it as the sum of a differentiating operator and a (learned) regulariser . This is discussed further in section 4 below
3. 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 with , we are interested in finding a nonlinear parametrised operator acting on the function . We will consider in the following the case ; 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).
For , let be a set of functions for , . Then we call: the layer with channels and corresponding index set .
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 and , , with channel index set , respectively, we call the mapping with
a layer operator. If the layer operator depends on a set of parameters , then we write .
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 and a desired output 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. ; consequently for colour images (in RGB) we have .
3.2. 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 and with channel index set , respectively. Let and , with compact support in , be the layer operator’s parameters for all . We call the convolutional layer operator for layer , if for each output channel
(3.1) 
with a pointwise nonlinear operator .
If the layer operator does not include a nonlinearity, we write . Now we can introduce the simplest convolutional network architecture by applying convolutional layer operators consecutively.
Definition 3.4 (layer Continuum Convolutional Network).
Let , then we call the composition of convolutional layer operator, denoted by , a Klayer Continuum Convolutional Network, such that
(3.2) 
In the following we will also refer to a layer CNN as the practical implementation of a 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 2layer 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
(3.3) 
Note that the second convolutional layer does not include a nonlinearity. Furthermore, it is necessary that , but typically it is often chosen such that . The full continuum ResNet architecture can then be summarized as follows. Let , then the composition of residual blocks, denoted by , defines a block continuum ResNet
(3.4) 
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 Kblock ResNet then consists of layers. Note that in the original work [16], the network was primarily designed for an image classification task rather than imagetoimage mapping that we consider here.
4. 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 of the continuum network as timestates of the function , where is a solution of the diffusion equation (2.1). In the following we assume single channel networks, i.e. for all layers. Then we can associate each layer with the solution such that . 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 and , such that and , then a diffusion layer operator , with parameters , is given by
(4.1) 
Similarly, an inverse filtering layer operator, with parameters is given by
(4.2) 
Note that this formulation includes a learnable timestep 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.
Remark 4.2.
The assumption of a single channel network, i.e. for all , can be relaxed easily. Either by assuming for some and all layers, or by introducing a channel mixing as in the convolutional operator (3.1).
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 matrixvector formulation that serves our purposes.
For simplicity we concentrate on the twodimensional case here. Let us then assume that the functions in each layer are represented as a square by image and we denote the vectorised form as . Then any linear operation on can be represented by some matrix ; in particular we can represent convolutions as a matrix.
We now proceed by rewriting the convolutional operation (3.1) in matrixvector notation. Given two layers and with channel index set , respectively, then we can represent the input layer as vector and similarly the output layer . Let represent the sum of convolutions in (3.1), then we can write the layer operator in the discrete setting as matrixvector operation by
where denotes the identity matrix of suitable size. Furthermore, following the above notation, we can introduce a stacked matrix consisting of all and a matrix by stacking the diagonal matrices . Then we can represent the whole convolutional layer in the discrete setting as
(4.3) 
Now the parameters of each layer are contained in the two matrices and .
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 approximated by the stencil (including the timestep )
(4.4) 
We use zero Neumann boundary conditions on the domain boundary
(4.5) 
Then we can represent (4.1) by
(4.6) 
The basis of learning a diffusion network is now given as estimating the diagonals of and the timestep . This can be done either explicitly as for the CNN or indirectly by an estimator network, as we will discuss next.
4.3. 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
(4.7) 
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

Linear diffusion; spatially varying and possible time dependence, .

Nonlinear diffusion; diffusivity depending on the solution , .
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 meanfree 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 and are then applied to to produce .
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 . Then, given layer , the filters are estimated by a small CNN from , which are then applied following an explicit update as in (4.6) to produce . 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 or . 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 image and apply the filters as pointwise matrixmatrix 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 we need 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 layer CNN. This CNN gets an image, given as layer , as input and estimates the filters . The architecture for this layer CNN as estimator is chosen to be rather simplistic, as illustrated in Figure 3. The input consists of one channel, which is processed by 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 . In particular for the chosen filter size of we have exactly convolutional parameters and biases per diffusion layer. That is for a 5layer CNN 29.509 parameters independent of image size.
5. 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.
5.1. Deconvolution with DiffNet
We first examine a simple deconvolution experiment to determine what features the DiffNet learns in an inverse problem. For this task we will only consider deconvolution without noise.
The forward problem is given by (2.1) with zero Neumann boundary condition (4.5) and constant diffusivity . For the experiment we choose , 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 , 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 . 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 4layer 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 loss of the output to the groundtruth image. The optimisation is performed for about 1000 epochs in batches of 16 with the Adam algorithm and initial learning rate of and a gradual decrease to . 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 error of . 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 computed as the output of the trained CNN in each layer, . The output of the last diffusion layer 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.
5.2. Nonlinear diffusion
Let us now consider the nonlinear diffusion process with the PeronaMalik 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
(5.1) 
with contrast parameter . We mainly concentrate here on the inverse problem of restoring an image that has been diffused with the PeronaMalik filter and contaminated by noise.
For the experiments we have used the test data from the STL10 database [6], which consists of 100,000 RGB images with resolution . 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 and . 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 4layer CNN, following the architecture in 3, gave the best tradeoff between reconstruction quality and network size. Increasing the amount of either layers, led to minimal increase in 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 loss. For the optimisation we have used the Adam algorithm with initial learning rate of and a gradual decrease to . 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 UNet [28]. This architecture has been widely applied in inverse problems [3, 10, 19, 20] and is mainly used to postprocess initial directly reconstructed images from undersampled or noisy data. For instance by filtered backprojection in Xray 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 with a gradual decrease to . Training of UNet takes about 3 hours.
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 UNet 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 UNet, 65.34 (DiffNet) compared to 61.08 (UNet). For the case with 1% noise, the quantitative measures are very similar. Here UNet 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.
6. 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 5layer DiffNet used for inversion of the nonlinear diffusion in section 5.2 has 101,310 parameters, whereas the used UNet with a filter size of has a total of 34,512,705 parameters; i.e DiffNet uses only of parameters compared to UNet 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.
6.1. Generalisability
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 PeronaMalik 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 UNet following the training protocol described in 5.2. Additionally, we have aborted the procedure when the networks started to clearly overfit the training data.
Results for the four scenarios are shown in Figure 10. Most notably DiffNet outperforms UNet 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 5001,000 samples independent of the noise case, whereas the UNet 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. . If the network would only learn the meanfree differentiating part, then these images should be zero. That implies that the illustrated filters in Figure 11 can be related to the learned regularisation . Additionally, we also show the diagonal filters 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 to 10% noise. As we conjectured in section 2.3, the learned update filters can be decomposed to a meanfree operation and a smoothing part . This implies that the magnitude of has to increase with higher noise. To examine this conjecture, we have taken (fixed) 32 samples from the reconstructed test data for each noise level and computed an estimate of as the sum , i.e. the deviation from the meanfree part. Furthermore, we interpret the smoothing level as the mean of the absolute value of . 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 meanfree part and a smoothing necessary for illposed inverse problems is valid.
7. 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 nonlinear 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 nonstationary filters, which can be learned in a network architecture. More specific, these local filters can be represented by a sparse subdiagonal (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.
Acknowledgements
We thank Jonas Adler and Sebastian Lunz for valubale discussions and comments. We acknowledge the support of NVIDIA Corporation with one Titan Xp GPU.
References
 [1] M. Abadi et al., TensorFlow: Largescale machine learning on heterogeneous systems, 2015. Software available from https://www.tensorflow.org/.
 [2] J. Adler and O. Öktem, Solving illposed inverse problems using iterative deep neural networks, Inverse Problems, 33 (2017), p. 124007.
 [3] S. Antholzer, M. Haltmeier, and J. Schwab, Deep learning for photoacoustic tomography from sparse data, Inverse Problems in Science and Engineering, (2018), pp. 1–19.
 [4] D. Calvetti and E. Somersalo, Hypermodels in the Bayesian imaging framework, Inverse Problems, 24 (2008), p. 034013.
 [5] Y. Chen, W. Yu, and T. Pock, On learning optimized reaction diffusion processes for effective image restoration, in Proceedings of the IEEE conference on computer vision and pattern recognition, 2015, pp. 5261–5269.
 [6] A. Coates, A. Ng, and H. Lee, An analysis of singlelayer networks in unsupervised feature learning, in Proceedings of the fourteenth international conference on artificial intelligence and statistics, 2011, pp. 215–223.
 [7] A. Douiri, M. Schweiger, J. Riley, and S. Arridge, Local diffusion regularization method for optical tomography reconstruction by using robust statistics, Opt. Lett., 30 (2005), pp. 2439–2441.
 [8] I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning, MIT Press, 2016. http://www.deeplearningbook.org.
 [9] E. Haber and L. Ruthotto, Stable architectures for deep neural networks, Inverse Problems, 34 (2017), p. 014004.
 [10] S. J. Hamilton and A. Hauptmann, Deep dbar: Real time electrical impedance tomography imaging with deep neural networks, IEEE Transactions on Medical Imaging, (2018).
 [11] S. J. Hamilton, A. Hauptmann, and S. Siltanen, A datadriven edgepreserving Dbar method for electrical impedance tomography, Inverse Problems and Imaging, 8 (2014), pp. 1053–1072.
 [12] K. Hammernik, T. Klatzer, E. Kobler, M. P. Recht, D. K. Sodickson, T. Pock, and F. Knoll, Learning a variational network for reconstruction of accelerated MRI data, Magnetic resonance in medicine, 79 (2018), pp. 3055–3071.
 [13] A. Hannukainen, L. Harhanen, N. Hyvönen, and H. Majander, Edgepromoting reconstruction of absorption and diffusivity in optical tomography, Inverse Problems, 32 (2015), p. 015008.
 [14] A. Hauptmann, S. Arridge, F. Lucka, V. Muthurangu, and J. Steeden, Realtime cardiovascular MR with spatiotemporal artifact suppression using deep learningproof of concept in congenital heart disease., Magnetic resonance in medicine, (2018).
 [15] A. Hauptmann, F. Lucka, M. Betcke, N. Huynh, J. Adler, B. Cox, P. Beard, S. Ourselin, and S. Arridge, Modelbased learning for accelerated, limitedview 3d photoacoustic tomography, IEEE transactions on medical imaging, 37 (2018), pp. 1382–1393.
 [16] K. He, X. Zhang, S. Ren, and J. Sun, Deep residual learning for image recognition, in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
 [17] T. Helin and M. Lassas, Hierarchical models in statistical inverse problems and the Mumford–Shah functional, Inverse problems, 27 (2010), p. 015008.
 [18] C. Hofer, R. Kwitt, M. Niethammer, and A. Uhl, Deep learning with topological signatures, in Advances in Neural Information Processing Systems, 2017, pp. 1634–1644.
 [19] K. H. Jin, M. T. McCann, E. Froustey, and M. Unser, Deep convolutional neural network for inverse problems in imaging, IEEE Transactions on Image Processing, 26 (2017), pp. 4509–4522.
 [20] E. Kang, J. Min, and J. C. Ye, A deep convolutional neural network using directional wavelets for lowdose Xray CT reconstruction, Medical Physics, 44 (2017).
 [21] Y. Khoo, J. Lu, and L. Ying, Solving parametric PDE problems with artificial neural networks, ArXiv, 1707.03351v2 (2017).
 [22] B. Kim, M. Wattenberg, J. Gilmer, C. Cai, J. Wexler, F. Viegas, et al., Interpretability beyond feature attribution: Quantitative testing with concept activation vectors (tcav), in International Conference on Machine Learning, 2018, pp. 2673–2682.
 [23] R. Kimmel, Numerical geometry of images: Theory, algorithms, and applications, Springer Science & Business Media, 2003.
 [24] R. Liu, Z. Lin, W. Zhang, and Z. Su, Learning pdes for image restoration via optimal control, in European Conference on Computer Vision, Springer, 2010, pp. 115–128.
 [25] T. Meinhardt, M. Moeller, C. Hazirbad, and D. Cremers, Learning proximal operators: Using denoising networks for regularizing inverse imaging problems, in International Conference on Computer Vision, 2017, pp. 1781–1790.
 [26] P. Perona and J. Malik, Scalespace and edge detection using anisotropic diffusion, IEEE Transactions on Pattern Analysis and Machine Intelligence, 12 (1990), pp. 629–639.
 [27] M. Raissi and G. E. Karniadakis, Hidden physics models: Machine learning of nonlinear partial differential equations, ArXiv, 1708.00588v2 (2017).
 [28] O. Ronneberger, P. Fischer, and T. Brox, Unet: Convolutional networks for biomedical image segmentation, in International Conference on Medical Image Computing and ComputerAssisted Intervention, Springer, 2015, pp. 234–241.
 [29] L. Ruthotto and E. Haber, Deep neural networks motivated by partial differential equations, arXiv preprint arXiv:1804.04272, (2018).
 [30] G. Sapiro, Geometric partial differential equations and image analysis, Cambridge university press, 2006.
 [31] S. ShalevShwartz and S. BenDavid, Understanding machine learning: From theory to algorithms, Cambridge university press, 2014.
 [32] J. Sirignano and K. Spiliopoulos, DGM: A deep learning algorithm for solving partial differential equations, ArXiv, 1708.07469v1 (2017).
 [33] J. Tompson, K. Schlachter, P. Sprechmann, and K. Perlin, Accelerating Eulerian fluid simulation with convolutional networks, ArXiv, 1607.03597v6 (2017).
 [34] J. Weickert, Anisotropic Diffusion in Image Processing, Teubner Stuttgart, 1998.
 [35] J. Weickert, B. T. H. Romeny, and M. A. Viergever, Efficient and reliable schemes for nonlinear diffusion filtering, IEEE transactions on image processing, 7 (1998), pp. 398–410.
 [36] E. Weinan, H. Jiequn, and J. Arnulf, Deep learningbased numerical methods for highdimensional parabolic partial differential equations and backward stochastic differential equations, ArXiv, 1706.04702v1 (2017).
 [37] Y. Wu, P. Zhang, H. Shen, , and H. Zhai, Visualizing neural network developing perturbation theory, ArXiv, 1802.03930v2 (2018).
 [38] J. C. Ye, Y. Han, and E. Cha, Deep convolutional framelets: A general deep learning framework for inverse problems, SIAM Journal on Imaging Sciences, 11 (2018), pp. 991–1048.
 [39] B. Zhu, J. Z. Liu, S. F. Cauley, B. R. Rosen, and M. S. Rosen, Image reconstruction by domaintransform manifold learning, Nature, 555 (2018), pp. 487–489.