Learning with Known Operators reduces Maximum Training Error Bounds

TL;DR

This paper introduces a method to incorporate known operators into machine learning models, reducing maximum training error bounds and the number of free parameters, with applications in physics and imaging.

Contribution

It presents a novel framework for embedding prior knowledge via known operators, leading to tighter error bounds and improved model efficiency.

Findings

Inclusion of known operators reduces maximum training error bounds.

Using known operators decreases the number of free parameters.

The approach is effective across diverse applications like CT reconstruction and vessel segmentation.

Abstract

We describe an approach for incorporating prior knowledge into machine learning algorithms. We aim at applications in physics and signal processing in which we know that certain operations must be embedded into the algorithm. Any operation that allows computation of a gradient or sub-gradient towards its inputs is suited for our framework. We derive a maximal error bound for deep nets that demonstrates that inclusion of prior knowledge results in its reduction. Furthermore, we also show experimentally that known operators reduce the number of free parameters. We apply this approach to various tasks ranging from CT image reconstruction over vessel segmentation to the derivation of previously unknown imaging algorithms. As such the concept is widely applicable for many researchers in physics, imaging, and signal processing. We assume that our analysis will support further investigation of known operators in other fields of physics, imaging, and signal processing.

Equations

$${\hat{u}}({\bf x})=\sum_{i=0}^{{\color[rgb]{0,0,0}N}-1}u_{i}{\varphi}({\bf w}_{i}^{\top}{\bf{x}})$$

$$|\hat{u}(\mathbf{x})-u(\mathbf{x})|\le {ϵ}_u$$

$$|U_k(\mathbf{x})-u_k(\mathbf{x})|\le {ϵ}_{u_k}.$$

$$f(\mathbf{x})=g(\mathbf{u}(\mathbf{x})):\mathcal{D}\to \mathbb{R}$$

$${\hat{f}}_g(\mathbf{x})=$$

$${\hat{f}}_u(\mathbf{x})=$$

$$\hat{f}(\mathbf{x})=$$

$${\hat{u}}_k(\mathbf{x})$$

$$|e_f|\le l_g\cdot {||{\mathbf{e}}_u||}_p$$

$$\hat{g}(\mathbf{x})$$

$${\hat{u}}_k(\mathbf{x})$$

$$|e_f|\le \sum _{j=0}^{N_{\hat{g}}-1}|g_j|\cdot l_{\phi }\cdot |e_{{\mathbf{u}}_j}|+{ϵ}_g.$$

$$\underset{\text{only dependent on }\hat{\mathbf{u}}(\mathbf{x})}{\underbrace{\sum _{j=0}^{N_{\hat{g}}-1}|g_j|\cdot l_{\phi }\cdot |e_{{\mathbf{u}}_j}|}}+\underset{\text{only dependent on }\hat{g}(\mathbf{x})}{\underbrace{{ϵ}_g}}$$

$$|e_{f,\ell ,k}|\le \sum _{{\ell }_i=1}^{\ell }||{{\mathbf{e}}_u}_{,{\ell }_i}|{|}_p\cdot l_{{\mathbf{f}}_{{\ell }_i-1}}$$

$$\mathbf{A}\mathbf{y}=\mathbf{x}$$

$$\mathbf{y}={\mathbf{A}}^{\top }(\mathbf{A}{\mathbf{A}}^{\top }{)}^{-1}\mathbf{x}$$

$$\mathbf{y}=\text{ReLu}({\mathbf{A}}^{\top }\mathbf{K}\mathbf{W}\mathbf{x})$$

$$V_{0}(\sigma)=\left\{\begin{array}[]{ll}0,&\text{ if }\lambda_{2}<0,\\ \exp(-\frac{R_{B}^{2}}{2\beta^{2}})(1-\exp(-\frac{S^{2}}{2c^{2}})),&\text{ otherwise,}\end{array}\right.$$

$${\mathbf{p}}_{\text{PB}}={\mathbf{A}}_{\text{PB}}\mathbf{v}$$

$${\mathbf{p}}_{\text{CB}}={\mathbf{A}}_{\text{CB}}\mathbf{v}$$

$$\mathbf{v}={\mathbf{A}}_{\text{PB}}^{\top }({\mathbf{A}}_{\text{PB}}{\mathbf{A}}_{\text{PB}}^{\top }{)}^{-1}{\mathbf{p}}_{\text{PB}}={\mathbf{A}}_{\text{PB}}^{\top }({\mathbf{A}}_{\text{PB}}{\mathbf{A}}_{\text{PB}}^{\top }{)}^{-1}{\mathbf{F}}^H\mathbf{x}$$

$${\mathbf{p}}_{\text{CB}}={\mathbf{A}}_{\text{CB}}{\mathbf{A}}_{\text{PB}}^{\top }({\mathbf{A}}_{\text{PB}}{\mathbf{A}}_{\text{PB}}^{\top }{)}^{-1}{\mathbf{F}}^H\mathbf{x}.$$

$${\mathbf{p}}_{\text{CB}}={\mathbf{A}}_{\text{CB}}{\mathbf{A}}_{\text{PB}}^{\top }\mathbf{W}{\mathbf{F}}^H\mathbf{C}\mathbf{F}{\mathbf{F}}^H\mathbf{x}={\mathbf{A}}_{\text{CB}}{\mathbf{A}}_{\text{PB}}^{\top }\mathbf{W}{\mathbf{F}}^H\mathbf{C}\mathbf{x}$$

Full text

1 Introduction

Pattern analysis and machine intelligence have been focussed predominantly on tasks that mimic perceptual problems. These are typically modelled as classification or regression tasks in which the actual reference stems from a human observer that defines the ground-truth. As we have only limited understanding on how these man-made classes emerge from the human mind, there is only limited knowledge available. As such, pattern recognition has relied on expert knowledge to design features that are suited towards a particular recognition task [1]. In order to alleviate the task of feature-design, researchers started also learning feature descriptors as a part of the training procedure [2]. Implementation of such on efficient hardware gave rise to first models that could outperform classical feature extraction methods significantly [3] and was one of the milestone works in the emerging field of deep learning.

With the rise of deep learning [4], researchers became aware that these methods of general function learning are applicable to a much wider range than mere perceptual tasks. Today, machine learning is applied in a much wider range of applications. Examples range from image super resolution [5], image denoising and inpainting [6], or even computed tomography [7]. In these fields, the methods from deep learning are often directly applied and often show performances that are either en par or even significantly better than results found with state-of-the-art methods. Yet, there are also reports that present surprising results in which parts of the image are hallucinated [8, 9]. In particular [9] demonstrates that mismatches in training and test data leads to dramatic changes in the produced result. Hence, blind deep learning methods have to be performed with care in order to be successful.

In this article, we explore the use of known operations within machine learning algorithms. First, we analyze the problem from a theoretical perspective and study the effect of prior knowledge in terms of maximal error bounds. This is followed by three applications in which we use prior operators to study to their effect on the respective regression or classification problem. Lastly, we discuss our observations in relation to other works in literature and give an outlook on future work. Note that some of the work presented here is based on prior conference publications [10, 11, 12, 13].

2 Known Operator Learning

The general idea of known operator learning is to embed entire operations into a learning problem. Figure 1 presents the idea graphically. We generally refer to the ${{N_{D}}}$-dimensional input of our trained algorithm as ${\bf x}^{\prime}\in\mathbb{R}^{{N_{D}}}$. In order to increase readability, we use an extended version ${\bf x}\in\mathbb{R}^{{N_{D}}+1}$ such that inner products with some weight vector $\bf w^{\prime}$ plus bias $w_{0}$ can be conveniently written, i. e. ${\bf w^{\prime}}^{\top}{\bf x^{\prime}}+w_{0}={{\bf w}}^{\top}{{\bf x}}$. Before looking into the properties of this approach and in particular maximal error bounds, we shortly summarize the Universal Approximation Theorem as it is closely related to our analysis. Note that the supplementary material to this article contains all proofs for the presented theorems in this section.

2.1 Universal Approximation

Theorem 1** (Universal Approximation Theorem).**

Let ${\varphi}(t):\mathbb{R}\rightarrow\mathbb{R}$ be a non-constant, bounded, and continuous function and ${u}({\bf x})$ be a continuous function on a compact set $\mathcal{D}\subset\mathbb{R}^{{N_{D}}+1}$. Then, there exist an integer {\color[rgb]{0,0,0}N}, weights ${\bf w}\in\mathbb{R}^{{N_{D}}+1}$, and $u_{i}\in\mathbb{R}$ that form an approximation ${\hat{u}}({\bf x})$

[TABLE]

such that the inequality

[TABLE]

holds for all ${\bf x}\in\mathcal{D}$ and $\epsilon_{u}>0$.

Theorem 1 states that for any continuous function ${u}({\bf x})$ an approximation ${\hat{u}}({\bf x})$ can be found such that the difference between true function and approximation is bounded by $\epsilon_{u}$. With increasing number of nodes {\color[rgb]{0,0,0}N}, $\epsilon_{u}$ will decrease. In literature, this result is often referred to as Universal Approximation Theorem [14, 15] and forms the fundamental result that neural networks with just a single hidden layer are general function approximators. Yet, this type of approximation may result in a very high requirement for the choice of {\color[rgb]{0,0,0}N} which is the reason why stacked layers of different type are known to be more successful [16].

We can extend Theorem 1 to vector-valued functions ${\bf u}({\bf x}):\mathcal{D}\rightarrow\mathbb{R}^{{N_{O}}}$ on $\mathcal{D}$ by postulating Theorem 1 for each of their components $k$

[TABLE]

Hence, universal approximation generally also applies to ${N_{O}}$-dimensional functions.

2.2 Known Operator Error Bounds

Knowing the limits of general function approximation, we are now interested in finding limits for mixing known and approximated operators. As previously mentioned, deep networks are never constructed out of a single layer, but rather take the form of the configuration shown in Figure 1. Hence, we need to consider layered networks to analyze the maximal error bounds. Instead of investigating entire networks, we choose to simplify our theoretical analysis to the special case

[TABLE]

with ${g}({\bf x}):S\rightarrow\mathbb{R}$, ${\bf u}({\bf x}):\mathcal{D}\rightarrow\mathcal{S}$, and compact sets $\mathcal{S}\subset\mathbb{R}^{N_{I}+1}$ and $\mathcal{D}\subset\mathbb{R}^{{N_{D}}+1}$. Note that this simplification does not limit the generality of our analysis, as we can map any knowledge on the structure of the network architecture either onto the output function ${g}({\bf x})$, the intermediate function ${\bf u}({\bf x})$, or directly as a transform of the inputs ${\bf x}$. Generalisation to ${N_{O}}$-dimensional functions is also possible following the idea shown in Eq. 3.

Previous definition of $f({\bf x})$ allows us to investigate different forms of approximation. In particular, we are able to introduce approximations ${\bf\hat{u}}({\bf x})$ and ${\hat{g}}({\bf x})$ following Theorem 1:

[TABLE]

Here $|{e}_{u}|\leq\epsilon_{u}$, $|{e}_{g}|\leq\epsilon_{g}$, and $|{e}_{f}|\leq\epsilon_{f}$ denote the errors that are introduced by respective approximation of ${\bf u}$, ${g}$, and $f$.

Next, we are interested in finding bounds on $|e_{f}|$ using above approximations. For the case of known ${\bf u}({\bf x})$, we can substitute ${\bf x}^{\ast}:={\bf u}({\bf x})$, as ${\bf u}({\bf x})$ is a fixed function. In this case Theorem 1 directly applies and a bound on $|e_{f}|$ is found as $|e_{f}|\leq\epsilon_{g}$ with $|e_{g}|\leq\epsilon_{g}$. If we would know ${g}({\bf x})$ in addition, $e_{g}$ would be 0 and the bound would shrink to the case of equality.

The case described in Eq. 5 is slightly more complicated, but we are also able to find general bounds as shown in Theorem 2.

Theorem 2** (Known Output Operator Theorem).**

Let ${\varphi}({\bf x}):\mathbb{R}\rightarrow\mathbb{R}$ be a non-constant, bounded, and continuous function and $f({\bf x})=g({\bf u}({\bf x})):\mathcal{D}\rightarrow\mathbb{R}$ be a continuous function on $\mathcal{D}\subset\mathbb{R}^{{N_{D}}+1}$. Further let $g({\bf x}):\mathcal{S}\rightarrow\mathbb{R}$ be Lipschitz-continuous function with Lipschitz constant ${l_{g}}=\textnormal{sup}\{||\nabla g({\bf x})||_{p}\}$ with $p\in\{1,2\}$ on $\mathcal{S}\subset\mathbb{R}^{N_{I}+1}$ and

[TABLE]

be a general function approximator of ${\bf u}({\bf x})$ with integer $N_{\hat{u}_{k}}$, weight ${\bf w}_{i,k}\in\mathbb{R}^{{N_{D}}+1}$, and $u_{i}\in\mathbb{R}$. Then, $e_{f}=f({\bf x})-{\hat{f}}({\bf x})$ with ${\hat{f}}({\bf x})={\hat{g}}({\bf\hat{u}}({\bf x}))={g}({\bf\hat{u}}({\bf x}))$ — as ${g}$ is known — is generally bounded for all ${\bf x}\in\mathcal{D}$ by

[TABLE]

with $e_{u}=e_{f}$ and component-wise approximation errors ${{{\bf e}_{u}}}=[{e}_{u_{0}},\ldots e_{u_{{N_{I}}}}]^{\top}$.

The bound for $|e_{f}|$ is found using a Lipschitz constant ${l_{g}}$ on ${g}({\bf x})$ which implies that the theorem will only hold, if Lipschitz-bounded functions are used for ${g}({\bf x})$. Analysis of Eq. 8 reveals that knowing ${\bf u}({\bf x})$ in this case, would imply ${{\bf e}_{u}}={\bf 0}$ which also yields equality on both sides.

We further explore this idea in Theorem 3. It describes a bound for the case that both ${g}({\bf x})$ and ${\bf u}({\bf x})$ are approximated.

Theorem 3** (Unknown Operator Theorem).**

Let ${\varphi}({\bf x}):\mathbb{R}\rightarrow\mathbb{R}$ be a non-constant, bounded, and continuous function with Lipschitz-bound ${l_{\varphi}}$ and $f({\bf x})={g}({\bf u}({\bf x})):\mathcal{D}\rightarrow\mathbb{R}$ be a continuous function on $\mathcal{D}\subset\mathbb{R}^{{N_{D}}+1}$. Further let

[TABLE]

be general function approximators of $g({\bf x}):\mathcal{S}\rightarrow\mathbb{R}$ and ${\bf u}({\bf x}):\mathcal{D}\rightarrow\mathcal{S}$ with integers $N_{\hat{g}}$ and $N_{\hat{u}_{k}}$, weights ${\bf w}_{j}\in\mathbb{R}^{N_{I}+1}$, ${\bf w}_{i,k}\in\mathbb{R}^{{N_{D}}+1}$, $g_{j},u_{i}\in\mathbb{R}$, and compact sets $\mathcal{S}\subset\mathbb{R}^{N_{I}+1}$ and $\mathcal{D}\subset\mathbb{R}^{{N_{D}}+1}$. Then, $e_{f}=f({\bf x})-{\hat{f}}({\bf x})$ with ${\hat{f}}({\bf x})={\hat{g}}({\bf\hat{u}}({\bf x}))$ is generally bounded for all ${\bf x}\in\mathcal{D}$ by

[TABLE]

where $\epsilon_{g}\geq|{e}_{g}|$, ${e_{{\bf u}_{j}}}={\bf w}_{j}^{\top}{{{\bf e}_{u}}}$, and ${{\bf e}_{u}}=[{e}_{u_{0}},\ldots e_{u_{{N_{I}}}}]^{\top}$ is the vector of errors introduced by the components of ${\bf\hat{u}}({\bf x})$.

The bound is comprised of two terms in an additive relation:

[TABLE]

where the first term vanishes, if ${\bf u}({\bf x})$ is known as $|{e_{{\bf u}_{j}}}|=0~{}\forall j$ and the second term vanishes for known ${g}({\bf x})$ as $e_{g}=0$. Hence for all of the considered cases, knowing ${g}({\bf x})$ or ${\bf u}({\bf x})$ is beneficial and allows to shrink the maximal training error bounds.

Given the previous observations, we can now also explore deeper networks that try to mimic the structure of the original function. This gives rise to Theorem 4.

Theorem 4** (Unknown Operators in Deep Networks).**

Let ${{\bf u}_{\ell}}({\bf x}_{\ell}):\mathcal{D}_{\ell}\rightarrow\mathcal{D}_{\ell-1}$ be a continuous function with Lipschitz-bound $l_{{\bf u}_{\ell}}$ on compact set $\mathcal{D}_{\ell}\subset\mathbb{R}^{N_{\ell}}$ with integer $\ell>0$. Further let ${\bf f}_{\ell}({\bf x}_{\ell}):\mathcal{D}_{\ell}\rightarrow\mathcal{D}$ be a function composed of $\ell$ layers / function blocks defined as recursion ${\bf f}_{\ell}({\bf x}_{\ell})={\bf f}_{\ell-1}({{\bf u}_{\ell}}({\bf x}_{\ell}))$ with ${\bf f}_{\ell=0}({\bf x})={\bf x}$ on compact set $\mathcal{D}\subset\mathbb{R}^{{N_{D}}+1}$ bound by Lipschitz constant $l_{{\bf f}_{\ell}}$ with $l_{{\bf f}_{\ell=0}}=1$. Recursive function ${\hat{{\bf f}}_{\ell}}({\bf x}_{\ell})={\bf\hat{f}}_{\ell-1}({\bf\hat{u}}_{\ell}({\bf x}_{\ell}))$ with ${\bf\hat{f}}_{\ell=0}({\bf x})={\bf x}$ is then an approximation of ${\bf f}_{\ell}({\bf x}_{\ell})$. Then, ${\bf e}_{f,\ell}={\bf f}_{\ell}({\bf x}_{\ell})-{\hat{{\bf f}}_{\ell}}({\bf x}_{\ell})$ is generally bounded for all ${\bf x}_{\ell}\in\mathcal{D}_{\ell}$ and for all $\ell>0$ in each component $k$ by

[TABLE]

where ${{\bf e}_{u}}_{,\ell}=[{e}_{u,\ell,{0}},\ldots e_{u,\ell,{{N_{I}}}}]^{\top}$ is the vector of errors introduced by ${\bf\hat{u}}_{\ell}({\bf x}_{\ell})$.

If we investigate Theorem 4 closely, we identify similar properties to Theorem 3. The errors of each layer / function block ${{\bf u}_{\ell}}({\bf x}_{\ell})$ are additive. If a layer is known, the respective error vector ${{\bf e}_{u}}_{,\ell}:={\bf 0}$ vanishes and the respective part of the bound cancels out. Furthermore, later layers have a multiplier effect on the error as their Lipschitz constants amplify $||{{\bf e}_{u}}_{,\ell}||_{p}$. Note that the relation $l_{{\bf f}_{\ell}}\leq\prod_{\ell_{i}=1}^{\ell}l_{{{\bf u}_{\ell}}_{i}}$ is shown in the supplementary material. A large advantage of Theorem 4 over Theorem 3 is that the Lipschitz constants $l_{{\bf f}_{\ell}}$ that appear in the error term $||{{\bf e}_{u}}_{,\ell_{i}}||_{p}\cdot l_{{{\bf f}}_{\ell_{i}-1}}$ are the ones of the true function ${\bf f}_{\ell}({\bf x}_{\ell})$. Therefore, the amplification effects are only dependent of the structure of the true function and independent of the actual choice of the universal function approximator. The approximator only influences the actual error ${{\bf e}_{u}}_{,\ell_{i}}$.

Above observations pave the way to incorporating prior operators into different architectures. In the following, we will highlight several applications in which we explore blending deep learning with prior operators.

3 Application Examples

We believe that known operators have a wide range of applications in physics and signal processing. Here, we highlight three approaches to use such operators. All three applications are from the domain of medical imaging, yet the method is applicable to many more disciplines to be discovered in the future. The results presented here are based on conference contributions [10, 11, 13]. Note that the supplementary material contains descriptions of experiments, data, and additional figures that were omitted here for brevity.

3.1 Deep Learning Computed Tomography

In computed tomography, we are interested in computing a reconstruction ${\bf y}$ from a set of projection images ${\bf x}$. Both are related by the X-ray transfrom ${\bf A}$:

[TABLE]

Solving for ${\bf y}$ requires inversion of above formula. The Moore-Penrose inverse of ${\bf A}$ yields the following solution:

[TABLE]

This type of inversion gives rise to the class of filtered back-projection methods, as it can be shown that $({\bf A}{\bf A}^{\top})^{-1}$ takes the form of a circulant matrix ${\bf K}$, i. e. ${\bf K}=({\bf A}{\bf A}^{\top})^{-1}={\bf F}^{H}{\bf C}{\bf F}$, where ${\bf F}$ denotes the Fourier transform, ${\bf F}^{H}$ its inverse, and ${\bf C}$ a diagonal matrix that corresponds to the Fourier transform of ${\bf K}$. As ${\bf K}$ typically is associated with a large receptive field, it is typically implemented in Fourier space. In order to be applicable for other geometries, such as fan-beam reconstruction additional Parker and cosine weights have to be incorporated that can elegantly be summarised in an additional diagonal matrix ${\bf W}$ to yield

[TABLE]

where $\text{ReLu}(\cdot)$ suppresses negative values as the final reconstruction algorithm.

Following the paradigm of known operator learning, Eq. 14 can also interpreted as a neural network structure as it only contains diagonal, circulant, and fully connected matrix operators displayed in Figure 2. A practical limitation of ${\bf A}$ is that it typically is a very large and sparse matrix. In practice, it is therefore never instantiated, but only evaluated on the fly using fast ray-tracing methods. For 3-D problems, the full matrix size is way beyond the memory restrictions of today’s compute systems. Furthermore, none of the parameters need to be trained as all of them are known for complete data acquisitions.

Incomplete data cannot be reconstructed with this type of algorithm and would lead to strong artifacts. We can still tackle limited data problems if we apply additional training of our network. As ${\bf A}^{\top}$ is large, we treat it as fixed during the training and only allow modification of ${\bf W}$ and ${\bf K}$. Results and experimental details are demonstrated in the supplementary material. Training of both matrices clearly improves the image reconstruction result. In particular, the trained algorithm learns to compensate for the loss of mass in areas of the reconstruction in which rays are missing.

As the trained algorithm is mathematically equivalent to the original filtered back-projection method, we are able to map the trained weights back onto their original interpretation which allows comparison to state-of-the-art weights. In Figure 3, we can see that the trained weights show similarity with the approach published by Schäfer et al. [18]. In contrast to Schäfer et al. who arrived at their weights following intuition, our approach is optimal with respect to our training data. In our present model, we have to re-train the algorithm for every new geometry. This could be avoided by modelling the weights using a continuous function which is sampled by the reconstruction network.

3.2 Learning from Heuristic Algorithms

Incorporating known operators generally allows blending of deep learning methods with traditional image processing approaches. In particular, we are able to choose heuristic methods that are well understood and augment them with neural network layers.

One example for such a heuristic method is Frangi’s vesselness [19]. The vesselness values for dark tubes are calculated using the following formula:

[TABLE]

where $|\lambda_{1}|<|\lambda_{2}|$ are the eigenvalues, $S=\sqrt{\lambda_{1}^{2}+\lambda_{2}^{2}}$ is the second order structureness, $R_{B}=\frac{\|\lambda_{1}\|}{\|\lambda_{2}\|}$ is the blobness measure, $\beta,c$ are image-dependent parameters for blobness and structureness terms, and $V_{0}$ stands for the vesselness value.

The entire multi-scale framework of Frangi filter can be mapped onto a neural network architecture [11]. In Frangi-Net, each step of the Frangi filter is replaced with a network counterpart and data normalization layers are added to enable end-to-end training. Multi-scale analysis is formed as a series of trainable filters, followed by eigenvalue computation in specialized fixed function network blocks. This is followed by another fixed function – the actual vesselness measure as described in Eq. 15.

We compare the segmentation result of the proposed Frangi-Net with the original Frangi filter, and show that the Frangi-Net outperforms Frangi filter regarding all evaluation metrics. In comparison to the state-of-the-art image segmentation model U-Net, Frangi-Net contains less than 6 % the number of trainable parameters, while achieving an AUC score around 0.960, which is only 1 % inferior to that of the U-Net. Adding a trainable guided-filter before Frangi-Net as preprocessing step yields an AUC 0.972 with only 8.5 % of the trainable parameters of U-Net which is statistically not distinguishable from U-Net’s AUC of 0.974.

Hence using our approach of known operators, we are able to augment heuristic methods by blending them with methods of deep learning saving many trainable parameters.

3.3 Deriving Networks

A third application of known operator learning that we would like to highlight in this paper, is the derivation of new network architectures from mathematical relations of the signal processing problem at hand. In the following, we are interested in hybrid imaging of magnetic resonance imaging (MRI) and X-ray imaging simultaneously. One major problem is that MRI $k$-space acquisitions typically allow parallel projection geometries, i. e. a line through the center $k$-space, while X-rays are associated with a divergent geometry such as fan- or cone-beam geometries. Both modalities allow different contrast mechanisms and simultaneous acquisition and overlay in the same image would be highly desirable for improved interventional guidance.

In the following, we assume to have sampled MRI projections ${\bf x}$ in $k$-space. By inverse Fourier Transform ${\bf F}^{H}$, they can be transformed into parallel projections ${\bf p}_{\text{PB}}={\bf F}^{H}{\bf x}$. Both parallel and cone-beam projections ${\bf p}_{\text{CB}}$ are related to the volume under consideration ${\bf v}$ by associated projection operations ${\bf A}_{\text{PB}}$ and ${\bf A}_{\text{CB}}$:

[TABLE]

As ${\bf v}$ appears in both relations, we can solve Eq. 16 for ${\bf v}$ using the Moore-Penrose Pseudo Inverse:

[TABLE]

Next, we can use ${\bf v}$ in Eq. 17 to yield

[TABLE]

Note that all operations on the path from $k$-space to ${\bf p}_{\text{CB}}$ are known. Yet, $({\bf A}_{\text{PB}}{\bf A}_{\text{PB}}^{\top})^{-1}$ is expensive to determine and may need significant amounts of memory. As we know from reconstruction theory, this matrix often takes the form of a circulant matrix, i. e. a convolution. As such, we can approximate it with the chain of operations ${\bf F}^{H}{\bf C}{\bf F}$ where ${\bf C}$ is a diagonal matrix. In order to add a few more degrees of freedom, we further add another diagonal operator in spatial domain ${\bf W}$ to yield

[TABLE]

as parallel to cone rebinning formula. In this formulation, only ${\bf C}$ and ${\bf W}$ are unknown and need to be trained. By design both matrices are diagonal and therewith only have few unknown parameters.

Even though the training was conducted merely on numerical phantoms we can apply the learned algorithm on data acquired with a real MRI system without any loss of generality. Using only 15 parallel-beam MR projections we were able to compute a stacked fan-beam projection with both approaches. In Figure 5 the results of the analytical and learned algorithms are shown. The result of the learned algorithm has much sharper visual impression compared to the analytical approach which intrinsically suffers from ray-by-ray interpolation and thus from a blurring effect. Note that additional smoothing could be incorporated into the network by regularization of the filter or additional hard-coded filter steps at request.

4 Discussion

For many applications, we do not know which operation is required in the ideal processing pipeline. Most machine learning tasks focus either on perceptual problems or man-made classes. Therefore, we only have limited knowledge on the ideal processing chain. In many cases, the human brain seems to have identified suitable solutions. Yet, our knowledge of the human brain is incomplete and search for well-suited deep architectures is a process of trial and error. Still, deep learning has shown to be able to solve tasks that were deemed as hard or close to impossible [20].

Now that deep learning also starts addressing fields of physics and classical signal processing, we are entering areas in which we have much better understanding of the underlying processes and therefore know that kind of mathematical operations need to be present in order to solve specific problems. Yet, during the derivation of our mathematical models, we often introduce simplifications that allow more compact descriptions and a more elegant solution. Still these simplifications introduce slight errors along the way and are often compensated using heuristic correction methods [21].

In this paper, we have shown that inclusion of known operators is beneficial in terms of maximal error bounds. We demonstrated that in all cases in which we are able to use partial knowledge on the function at hand, the maximal errors that may remain after training of the network are reduced even for networks of arbitrary depth. Note that in the future tighter error bounds than the ones described in this work might be identified that are independent of the use of known operators. Yet, our error analysis is still useful, as for the case of increasing number of known operations in the network, the magnitude of the bound shrinks up to the point of identity, if all operations are known. To the knowledge of the authors, this is the first paper to attempt such a theoretical analysis of the use of known operators in neural network training.

In our experiments with CT reconstruction, we could demonstrate that we are able to tackle limited angle reconstructions using a standard filtered back-projection-type of algorithm. In fact, we only adopted weights while run-time, behaviour, and computational complexity remained unchanged. As we can map the trained algorithm back onto its original interpretation, we could also investigate shape and function of the learned weights. They demonstrated similarity to a heuristic method that could previously only be explained by intuition rather than by showing optimality. For the case, of our trained weights, we can demonstrate that they are optimal with respect to the training data.

Based on Frangi’s vesselness, we could develop a trainable network for vessel detection. In our experiments, we could demonstrate that training of this net already yields improved filters for vessel detection that are close in terms of performance with a much more complex U-Net. Further inclusion of a trainable denoising step yielded an accuracy that is statistically not distinguishable from U-Net.

As last application of our approach, we investigated rebinning of MR data to a divergent beam geometry. For this kind of rebinning procedure, a fast convolution-based algorithm was previously unknown. Prior approaches relied on ray-by-ray interpolation that is typically introducing blurring. With our hypothesis that the inverse matrix operator takes the form of a circulant matrix in spatial domain in combination with an additional multiplicative weight, we could train a new algorithm successfully. The new approach is not just computationally efficient, it also features images of a degree of sharpness that was previously not reported in literature.

Although only applications from the medical domain are shown in this paper, this does not limit the generality of our theoretical analysis. Similar problems are found in many fields, e. g. computer vision [22], image super resolution [23], or audio signal processing [24].

Obviously, known operators have been embedded into neural networks already for a long time. Already, LeCun et al. [2] suggested convolution and pooling operators. Janderberg et al. introduced differentiable spatial transformations and their resampling into deep learning [25]. Lin et al. use this for image generation [26]. Kulkarni et al. developed an entire deep convolutional graphics pipeline [27]. Zhu et al. include differentiable projectors to disentangle 3D shape information from appearance [28]. Tewari et al. integrate a differentiable model-based image generator to synthesize facial images [29]. Adler et al. shows an approach to partially learn the solution for ill-posed inverse problems[30]. Ye et al. [31] introduced the Wavelet transform as multi-scale operator, Hammernik et al. [32] mapped entire energy minimization problems onto networks, and Wu et al. even included the guided filter as layer into a network [33]. As this list could be continued with many more references, we see this as an emerging trend in deep learning. In fact, any operation that allows the computation of a sub-gradient [34] is suited to be used in combination with the back-propagation algorithm. In order to integrate a new operator, only the partial derivatives / sub-gradients with respect to its inputs and its parameters have to be implemented. This allows inclusion of a wide range of operations. To the best of our knowledge, this is the first paper giving a general argument for the effectiveness of such approaches.

Next, the introduction of a known operator is also associated with a reduction of trainable parameters. We demonstrate this in this paper in all of our experiments. This allows us to work with much fewer training data and helps us to create models that can be transferred from synthetic training data to real measured data. Zarei et al. [35] drive this approach so far that they are able to train user-dependent image denoising filters using only few clicks from alternate forced-choice experiments. Thus, we believe that known operators may be a suitable approach to problems for which only limited data is available.

At present we are unaware how to predict the benefit of using known operators before the actual experiment. Our analysis only focuses on maximum error bounds. Therefore, investigation of expected errors following for example the approach of Barron seems interesting for future work [36]. Also analysis of the bias variance trade-off seems interesting. In [37, Chapter 9] Duda and Hart already hinted at the positive effect of prior knowledge on this trade-off.

Lastly, we believe that known operators may be key in order to gain better understanding of deep networks. Similar to our experiments with Frangi-Net, we can start replacing layers with known operations and observe the effect on the performance of the network. From our theoretical analysis, we expect that inclusion of a known operation will not or only insignificantly reduce the system’s performance. This may allow us to find configurations for networks that only have few unknown operations while showing large parts that are explainable and understood. Figure 6 shows a variant of this process that is inspired by [38]. Here, we offer a set of known operations in parallel and determine their optimal superposition by training of the network. In a second step, connections with low weights can be removed to iteratively determine the optimal sequence of operations. Furthermore, any known operator sequence can also be regarded as a hypothesis for a suitable algorithm for the problem at hand. By training, we are able to validate of falsify our hypothesis similar to our example of the derivation of a new network architecture.

5 Conclusion

We believe that the use of known operators in deep networks is a promising method. In this paper, we demonstrate that the use of such reduces maximal error bounds and experimentally show an reduction in the number of trainable parameters. Furthermore, we applied this to the case of learning CT reconstruction yielding networks that are interpretable and that can be analysed with classical signal processing tools. Also mixing of deep and known operator learning is beneficial, as it allows us to build smaller networks with only 6 % of the parameters of a competing U-Net while being close with respect to their performance. Lastly, the known operators can also be found using mathematical derivation of networks. While keeping large parts of the mathematical operations, we only replace inefficient or unknown operations with deep learning techniques to find entirely new imaging formulas. While all of the applications shown in this paper stem only from the medical domain, we believe that this approach is applicable to all fields of physics and signal processing which is the focus of our future work.

Acknowledgment

The research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (ERC Grant No. 810316).

Author Contributions

Andreas Maier is the main author of the paper and is responsible for the writing of the manuscript, theoretical analysis, and experimental design. Christopher Syben and Bernhard Stimpel contributed to the writing of Section 3.3 and the supporting experiments. Tobias Würfl and Mathis Hoffmann supported writing Section 3.1 and performed the experiments reported in this section. Frank Schebesch contributed to the mathematical analysis and the writing thereof. Weilin Fu conducted the experiments supporting Section 3.2 and contributed to their description. Leonid Mill, Lasse Kling, and Silke Christiansen contributed to the experimental design and the writing of the manuscript.

Data Availability Statement

All data in this publication are publicly available. The experiments in Section 3.1 use data of the low-dose CT challenge [39]. Section 3.2 uses the DRIVE database [40]. The data for Section 3.3 is available in a Code Ocean Capsule available at https://doi.org/10.24433/CO.8086142.v2 [41].

Code Availability Statement

The code and data for this article, along with an accompanying computational environment, are available and executable online as a Code Ocean Capsule. Experiments in Section 3.1 can be found at https://doi.org/10.24433/CO.2164960.v1 [42]. The code on learning vesselness in Section 3.2 are published at https://doi.org/10.24433/CO.5016803.v2 [41]. Code for Section 3.3 is available at https://doi.org/10.24433/CO.8086142.v2 [43]. The code capsules for the experiments in Section 3.1 and Section 3.3 were implemented using the open source framework PYRO-NN [44].