Our world is full of physics-driven data where effective mappings between data manifolds are desired. There is an increasing demand for understanding combined model-based and data-driven methods. We propose a nonlinear, learned singular value decomposition (L-SVD), which combines autoencoders that simultaneously learn and connect latent codes for desired signals and given measurements. We provide a convergence analysis for a specifically structured L-SVD that acts as a regularisation method. In a more general setting, we investigate the topic of model reduction via data dimensionality reduction to obtain a regularised inversion. We present a promising direction for solving inverse problems in cases where the underlying physics are not fully understood or have very complex behaviour. We show that the building blocks of learned inversion maps can be obtained automatically, with improved performance upon classical methods and better interpretability than black-box methods.

1. Introduction. We are living in a world full of physics-driven data with an increasing demand for combining model-based and data-driven approaches in areas of science, industry and society. In many cases it is essential to reliably recover hidden multi-dimensional model parameters (signals) from noisy indirect observations (measurements) in imaging or sensing technology in medicine, engineering, astronomy or geophysics. These , are often ill-posed, suffering from non-uniqueness and instability in direct inversion. Classical model-based research on inverse problems has focused on variational regularisation methods to guarantee existence and stable approximation of solutions under uncertainty like noise in the measurements [19, 7]. For linear inverse problems the singular value decomposition (SVD) [23] is a classical tool to directly construct a regularised inverse, e.g. in the sense of Tikhonov regularisation [19]. Recent research in inverse problems has focused on combining deep learning with model-based approaches based on knowledge of the underlying physics [4]. Precise knowledge is often not available; for now we rely mainly on empirical evidence that such approaches can still be applied when one makes use of inexact operators that approximate the exact physical process [25]. The main limitation of such methods are that they require an iterative application of expensive, possibly nonlinear mappings. Moreover, they are hard to interpret due to the lack of connection between data structure and structure of the mapping.

Figure 1.1: L-SVD learns the inversion mapping via a hybrid nonlinear data manifold learning.

We propose the ‘learned singular value decomposition’ (L-SVD): a direct method that provides the inversion procedure with an explainable connection between measurements and signals. The method does not rely on an iterative application of expensive mappings; it does not need any information on the mapping at all. L-SVD makes use of two connected autoencoders: the first one encodes measurement to latent code , while the second one encodes signal x to latent code in a nonlinear way; both latent codes are connected with a linear ‘scaling’ layer. The training of all parameters is done simultaneously, which enforces the latent codes to preserve as much information on the measurements and signals as possible, while making sure that the codes have very similar structure. After training, a reconstruction is obtained by consecutive application of encoding, scaling and decoding (see Figure 1.1).

In case a forward mapping is available, it can be used to extract an effective encoding and decoding, e.g. via the SVD. In such cases, the L-SVD can be shaped to a data-driven Tikhonov regularisation method, for which we prove convergence with respect to noise. In case a forward mapping is not available, L-SVD allows to learn the nonlinear inversion dynamics via two autoencoders and a linear scaling layer. This has the advantage that advances in interpretable autoencoding have a direct effect on the L-SVD method, making it easier to analyse than other fully learned inversion methods. Assuming that the autoencoders can be trained with high accuracy, finding the connection between both codes is a much lower dimensional problem than finding a nonlinear map between the measurements and signals directly. Moreover, in a semi-supervised scenario where not all training data is available in supervised pairs, the autoencoders still allow to learn from all samples, which is not possible in most supervised neural networks.

Current research in inverse problems is focused at developing theory for combining data-driven deep learning with model-based approaches [4]. Recently developed methods with theoretical guarantees such as convergence and stability include those in which a regularisation term is explicitly learned [37, 35] and those where an initial imperfect reconstruction is post-processed by a neural network [46, 10]. These methods are two-step procedures and require knowledge of the forward mapping, while L-SVD provides a one-step procedure without this knowledge. A series of works about the optimal regularised inverse matrix (ORIM) investigated the problem of finding a linear inverse matrix for a noisy inverse problem, both for the case where the linear forward matrix is known [16, 17] and for the case where it is not known [16], i.e. a fully learned scenario. Similar to the L-SVD, the idea of data-driven approximation of nonlinear mappings via model reduction was proposed in [8]. In that paper, the interest is in approximating a forward mapping, while we are interested in solving an inverse problem. The idea of connecting two autoencoders was exploited in [53] and [24], where the authors solved a superresolution and a deconvolution problem with a patch-based method. In our work, we consider a more general method that does not assume identical domains for measurement and signal. An extensive literature review of similar methods and the embedding of L-SVD in its research context is given in Section 5.

1.1. Contribution. This paper proposes the learned singular value decomposition (L- SVD), a general data-driven method that nonlinearly encodes (compresses) data in two vector spaces and connects them in an easy-to-understand way. The contributions of our method can be seen as the extension of existing methods in the following way:

1. Data-driven solution of inverse problems: L-SVD is a nonlinear generalisation of Tikhonov regularisation in Bayesian inverse problems [20, 28, 50] and piecewise linear estimates [52], which are linear data-driven variants of classical SVD approaches. We show that L-SVD can be shaped to a data-driven Tikhonov regularisation for which we provide a convergence analysis. In general, L-SVD requires no a-priori information on the forward mapping to achieve good reconstruction quality.

2. Improved generalisation via nonlinear hybrid encoding: Autoencoders show that nonlinear encoding provides better encodings than linear encoding [26]. L-SVD shows that this is also the case when encodings are used to solve inverse problems. An autoencoder can act as a regulariser when attached to a supervised neural network that is trained for a supervised task, e.g. classification [54, 31]. L-SVD makes use of two autoencoders for the task of solving an inverse problem, which gives improved generalisation performance. Moreover, it enables high-quality reconstruction in a semi-supervised setup.

1.2. Overview of the paper. Throughout the paper, we make use of the above two perspectives to show the advantage of using the L-SVD method for addressing the previously mentioned limitations. In Section 2, a brief overview of the classical SVD and inversion methods is given, which serves as motivation for the L-SVD method. Next, a precise definition of L-SVD is provided in Section 3, accompanied with various architecture choices that exploit its potential. One of these choices results in a data-driven Tikhonov regularisation, while a second one results in a fully learned L-SVD method. In Section 4, we analyse the L-SVD method by showing its connection with Bayesian inverse problems, by showing a convergence analysis of the data-driven Tikhonov regularisation and providing a stability and error estimate for the fully learned L-SVD method. After that, in Section 5, the connection of our work with several fields of research are discussed. In Section 6, we explain simulation experiments that show the transition from non-learned to learned, linear to nonlinear, and single to hybrid encoding. Results are provided in Section 7, where we visualise these transitions by looking at the latent space, the decoded representations of the latent space and the dependency of L-SVD on noise. The section is completed with a comparison between L-SVD and state-of-the art reconstruction methods applied on a biomedical CT data set. In Section 8, we conclude with some remarks and outlook for future work.

2. Motivation: SVD and inversion methods. The motivation of our L-SVD method can be found in the application of classical SVD and its variants in inversion methods. In our work we consider the finite dimensional version of the equation introduced in Section 1. That is, we make use of the ‘first discretise, then optimise’ approach. We define the inverse problem as

(2.1)

where we wish to reconstruct the signals from measurements corrupted by additive noise are Banach spaces. The mapping is in general a nonlinear one. For this section however, we assume a linear operator that we call A. Any A can be written in its singular value decomposition: are unitary matrices and is a diagonal matrix with nonnegative real numbers (singular values) on the diagonal. We now summarise well-known inversion methods that can be written as the application of an SVD [23].

2.1. Maximum likelihood estimator (MLE). The MLE is defined via the maximisation over x given measurements ]. Its solution is obtained by applying the MoorePenrose inverse (to the measurements

2.2. Tikhonov regularisation. Tikhonov regularisation is a method which puts a uniform variance prior on the desired solution x. It solves an -weighted minimisation problem that can be solved directly via its regularised Moore-Penrose inverse:

The diagonal elements of are defined as ), which means that for smaller scales , the new inverse scale goes to zero as gets larger. The optimal depends on the type of noise; usually increases with noise level

For a specific model choice (see Section 3.1.1), L-SVD becomes a data-driven Tikhonov regularisation method, for which we provide a convergence analysis in Section 4.2.

2.3. Truncated SVD. The best Frobenius approximation of A with rank r is given [18] by the truncated SVD (T-SVD):

Here we made use of the ‘thin’ representation, where the top r rows of U and V respectively. is a diagonal square matrix that consists of the largest r singular values of A. With the thin representation, we lose one desirable property, namely that of unitary matrices: while still holds, generally

T-SVD can be applied in an inversion method instead of the standard SVD for noisy measurements becomes small for i large, noise is amplified by 1problem is mitigated by solving instead. It has the additional benefit that the thin decomposition is smaller than the full SVD, requiring less memory and computation time.

3. The learned singular value decomposition. In this section, we provide the general L-SVD method for solving inverse problems. It aims to solve the inverse problem as defined in (2.1), where the forward mapping A may be nonlinear. L-SVD can be seen as a nonlinear learned variant of the inversion methods in Section 2, where is replaced by a nonlinear encoder and V by a nonlinear decoder.

3.1. Model statement. The L-SVD model (Figure 3.1) is a trained neural network that consists of a measurement autoencoder (green), a signal autoencoder (blue) and a reconstruction component (red). Reconstruction ˆx from measurement is obtained via the latent representations , which are part of the autoencoders. The latent space is a low-dimensional space, i.e. . A more formal definition is given as:

Figure 3.1: Schematic of the L-SVD method. Green: autoencoder for measurement autoencoder for signal x. Red: reconstruction procedure. The standard network does not use the gray connections for training. Note that is used multiple times, but has shared weights.

The L-SVD model is obtained by minimising a neural network loss function that consists of three parts:

where we minimise over all trainable parameters parand over all samples (i) in the training set. The distance functions ) can be any metric; often used in neural networks are (Wasserstein) metrics. The L-SVD model encodes measurements representation that contains sufficient information to approximately reconstruct the clean data y, while being able to map to an encoded representation that can approximately be decoded to the desired signal x. Since the output of the data autoencoder is the clean data y, instead of the corrupted measurements , it can be seen as a denoising autoencoder (DAE) [51]. This means that noise will not necessarily be represented in the latent variable

3.1.1. Opportunities by various model choices. Below some specific model choices and variations on the standard model are discussed, which establish certain capabilities of the L-SVD model:

Data-driven Tikhonov regularisation: With the SVD, a linear encoder and decoder can be derived from the operator A (see Section 2). If this is combined with the a nonlinear function, a data-driven Tikhonov functional is obtained. A convergence analysis is given in Section 4.2.

Linearly connected nonlinear representations: Results from autoencoding [26] motivate the search for a nonlinear encoding and decoding such that a nonlinear representation of signals and measurements is obtained. The expressiveness of the nonlinear encoder and decoder allows us to restrict the scaling layer to be a square matrix Σ , either full or diagonal. Generic stability and reconstruction estimates are given in Section 4.3. They depend on the reconstruction quality of the autoencoders, which can be freely chosen depending on the application at hand. This includes regularised autoencoders, such as sparse [40] or contractive autoencoders [41] of a fully-connected or convolutional type.

The autoencoder on the measurement side can be chosen to be a regular autoencoder instead of a denoising autoencoder, with the effect that noise is represented in the encoded version . This means that the latent dimension should be large enough, since unstructured noise can not be compressed. Moreover, this means that Σ should be able to remove (part of) the noise, since the latent variable noise-free.

Structured latent space: No specific structures of the latent spaces imposed. If control on these spaces is desired, one could sample from a desired set in the latent space and add one of the following losses to (3.1):

This means the sampled latent code is decoded to a signal x (gray in Figure 3.1), after which it takes either the blue autoencoder path or the red reconstruction path, without the final decoder step . Although not guaranteed, it is likely that due to this additional loss, the encoder will map all samples in the training set to this subset E. If this is the case, it means that we have control over the latent space . Moreover it turns out that having a bound on enables us to compute a uniform error bound for the reconstruction procedure (see Section 4.3).

4. Analysis. In this section, we first consider the case where the encoder and decoder are derived from the SVD of A. In Section 4.1, we show that training this L-SVD model with a linear scaling coincides with learning the covariance matrix of a prior in Bayesian inverse problems. In Section 4.2, we provide a convergence analysis with respect to noise in case a nonlinear scaling is used. Finally, in Section 4.3, we provide a stability and error estimate for the L-SVD model with nonlinear encoding and decoding.

4.1. Connection with Bayesian inverse problems. Here we provide an explicit connection between a linear L-SVD model and the solution of a Bayesian inverse problem with Gaussian noise, Gaussian prior and known forward operator A. For an introduction to statistical and Bayesian inverse problems, we refer to [20, 28, 50].

Banach spaces. Consider the inverse problem

where has full row-rank, i.e. rank(, with thin SVD decomposition . Moreover, let be a Gaussian prior measure on x. We define and restrict the covariance matrix to be of rank n that can be written as is positive definite.

Then the maximum a posteriori (MAP) estimate written as an SVD inversion method in the following way:

For the proof we refer to Appendix A. The connection between (4.1) and L-SVD is clear if we define the linear measurement encoder to be , the linear signal decoder to be and we assume the noise covariance matrix B to be known. Then it can be seen in (4.1) that learning a linear Σ is equivalent to optimising over the prior covariance matrix , defined via

4.1.2. Scale dependency on Gaussian noise level. For many inverse problems one assumes an additive noise term that originates from the Gaussian distribution where the noise level is either known or estimated. If the data covariance matrix B is replaced with ) is simplified to

This implies a stronger regularisation (by an increased noise level . Thus, by learning the scales Σ, we learn the prior distribution on x which regularises our inverse problem. For a prior distribution Id, it is easily shown that we get the formulation for classical Tikhonov-regularisation (2.3) back:

4.2. Data-driven Tikhonov regularisation with L-SVD. In case the forward operator A is known, it is possible to use the linear SVD encoder and decoder and train a nonlinear scaling Σ. By doing this in a structured manner, a data-driven Tikhonov regularisation can be learned, where the scales Σ nonlinearly depend on the measurements , potentially via

are Hilbert spaces. Let be a nonlinear function s.t. for all be the singular value decomposition as defined in Section 2. We define

L-SVD provides a direct reconstruction through the encoder, scaling and decoder. The function N is a trained neural network that can be bounded by construction. In case it has the form of (4.3), the acquired solution is the unique minimiser of a convex functional:

be defined as in (4.3). Then unique minimiser of the functional

(4.4)

is strictly convex. Moreover, limhas a unique minimiser, which can be found by checking the first-order optimality condition. This yields the expression of (4.3).

L-SVD in the form of (4.3) provides a regularisation method: in the noisy case mines which scales should be regularised more and which less. In the limit of noise decreasing to zero, the solution converges to the unregularised solution, as shown in Theorem 4.4.

defined as in (4.3). If

(4.6) lim

Proof. Define the sequence With we define the functional (4.4) with variables as defined above. By Theorem 4.3, is the unique minimiser of . Hence with

From this, we obtain

Because is bounded, it has a convergent subsequence

Since the bounded linear operator A is sequentially continuous,

Again by Theorem 4.3, we obtain that

Together with (4.8) this implies

(4.10) Av = y.

Since any minimiser of is in ker(and therefore ker(. Together with (4.10), by [19, Theorem 2.5], this implies that . By applying the same argument to all subsequences, we obtain

(4.11)

Since the sequence is arbitrarily chosen s.t. 0, the desired expression (4.6) follows.

Finally, the convergence rate of (4.3) is derived. The proof of Theorem 4.5 makes use of several theorems in [19].

. Then, the optimal parameter choice is , which provides the convergence rate

(4.12)

This function meets the assumptions of [19, Theorem 4.1]:

and

for all ]. Next, for 0 we define

Furthermore, we define

(4.14)

Finally, we define for 0

(4.15)

For this 1, the requirement in [19, Theorem 4.3] holds:

(4.16)

An expanded derivation of (4.16) is provided in Appendix B. For Corollary 4.4], the parameter choice is of optimal order in

(. The best possible convergence rate is obtained with . This provides the convergence rate (4.12).

4.3. Fully learned L-SVD. In Section 4.2, it was shown that a data-driven Tikhonov regularisation is obtained by giving L-SVD the structure as shown in (4.3). The linear encoder and decoder are obtained from the SVD of A, which allows for the convergence analysis that was provided. In the current section we analyse the more general case where the L-SVD is fully learned, which means that the SVD of A can not be used. Furthermore, we consider the L-SVD with nonlinear encoder and decoder and diagonal matrix Σ, i.e. the second model choice of Section 3.1.1.

Definition 4.6 (Nonlinear L-SVD with diagonal scaling). Define the nonlinear functions and variables as in Definition 3.1. Let Σ : be a diagonal matrix Σ with on the diagonal.

Theorem 4.7 (Stability estimate). Let the variables and functions in L-SVD be defined as in Definition 4.6. Define two different measurements Then there is an M > 0 that depends on the weights and nonlinearities in the L-SVD network such that

(4.17)

Proof. From Definition 4.6 and on its turn Definition 3.1, we compute the bound

(4.18)

where represents the operator norm. After training, the weights of the L-SVD model are fixed, which provides the desired bound (4.17).

consisting of L layers with weight matrices and pointwise nonlinearities , the following bound is achieved: . In case the nonlinearity is a ReLU, leaky ReLU or hyperbolic tangent, C = 1 and thus the norm only depends on the norms of the weight matrices. In this case the error estimate (4.18) can be written as

The error estimate can thus be controlled by assuring small -norms of the weight matrices, which can be achieved by adding weight regularisation in the neural network objective function.

Next, we prove that an error estimate on the difference between any reconstructed signal and the true solution exists, provided that its associated measurement maps to a ball in the latent space. Before we provide this error estimate, we first prove a supporting Lemma.

be a continuous function. Assume for some given 0 is the closed ball centered at 0 with radius r.

Proof. Let us first define a scaled function this scaled function,

The closed unit ball is a contractible space. By definition of a contractible space [22]: continuous such that for all where is the contraction point. Since ), it holds that are continuous, its composition is continuous. From this it follows that ) is a contractible space, which implies that ) is a contractible space. In other words, ) does not have any ‘holes’.

Left to show is that the boundary of ) lies outside , which implies ) is contractible. For this we make use of [43, Theorem 4.22]: since F is a continuous mapping of a metric space () into a metric space (), and the boundary of the unit ball (i.e. ) is a connected subset of implies ) is connected. Moreover, . This implies that the boundary of the unit ball lies completely outside , which completes the proof.

Theorem 4.10 (Reconstruction error estimate). Let the variables and functions in L-SVD be defined as in Definition 4.6. Assume that for some 0 . Then there is an M > 0 that depends on the weights and nonlinearities in the L-SVD network such that for all . Moreover, for all , we have the error estimate

Proof. The first part of the proof can be obtained by combining the operator norm of the decoding function in x and the given bound:

For the second part of the proof, we make use of Lemma 4.9. We define Σ)), which is a continuous function from . Since for all , we know that all elements in the ball are in the range of ). Therefore, for all , there exists a , such that the same bound (4.19) holds.

Note that the reconstruction error estimate depends on . The former can be kept small by regularising the weights of the decoder in the training phase. The latter can be kept small by including in the cost function, as described in the last point of Section 3.1.1. After training, points in the unit ball can be sampled uniformly and passed through the network to examine an actual value for

5. Research context. Our paper presents the L-SVD method for solving inverse problems via hybrid autoencoding. Within the method, a low-dimensional (i.e. sparse) representation or manifold is explicitly learned. This method has connections with many research fields, which are pointed out in this section.

5.1. Combining data and models for solving inverse problems. Recent research in inverse problems seeks to develop a mathematically coherent foundation for combining data-driven deep learning with model-based approaches based on physical-analytical domain knowledge [4]. A first class of methods are partially learned variational and iterative methods [1, 30, 25, 12]. These methods can be seen as a learned variant of gradient, proximal or primal-dual methods. They require less iterations than their non-learned counterparts, but the demand on training time is substantial, while the mathematical analysis of these methods is limited. A second approach is to learn an explicit regularisation term [15, 2, 37, 35]. Signals affected by artefacts are penalised, while the desired signals are not. Reconstructions are of higher quality compared to classical regularisation choices, but their computation time is of the same order. A third approach is to perform learned post-processing of initial reconstructions obtained by classical methods, which may be affected by artefacts [27, 46]. Data-consistent reconstructions can be obtained without an iterative procedure [46]. However, the quality of the reconstructions heavily depends on the initial reconstruction, which is often obtained by applying a pseudo-inverse to the data.

The above methods depend on precise knowledge of the physical process, modelled by a forward mapping, which is not always available. However, emperically it was shown that learned iterative methods can still be used in case of inexact forward operators [25]. A recent work [36] aims to improve an inexact forward operator (linear mapping) explicitly with a neural network, after which it is applied in a variational framework. It was proven that small training losses ensure that the optimisation procedure finds a solution close to the one that would have been found with the true forward operator.

The optimal regularised inverse matrix (ORIM) method [16, 17] is a data-driven method that aims to find a linear inverse matrix for a noisy inverse problem. A global minimiser of the Bayes’ risk is found when there is knowledge about the forward operator and about the probability distribution of the signals. A fully learned variant where the forward operator is not known is also available [16], albeit at a higher computational cost.

5.2. Fully learned image reconstruction. Our work is closely related to the work of Zeng et al. [53]. In this paper, the task of superresolution is solved by autoencoding patches of both a low- and high-resolution image and finding a nonlinear mapping between them. Gupta et al. [24] used this method for the task of removing motion blur, which is a specific case of a deconvolution problem. In both cases, a one-layered autoencoder was applied to patches of the distorted image (measurement) and desired image (signal). This is only possible if measurement and signal lie in the same domain and if the forward mapping has very little effect outside the local patch. We consider a more general method that does not work patch-based and therefore does not assume identical domains for measurement and signal. As a result, the forward mapping may have a global behaviour.

A different fully learned reconstruction method without this restriction is proposed by Zhu et al. [55], where the problem of finding a reconstruction from undersampled MRI data is considered. Their neural network consists of three fully-connected layers, followed by three convolutional layers, which maps Fourier measurements directly to the desired image. A joint low-dimensional manifold is learned implicitly, since there is no explicit low-dimensional representation within the network architecture. In our work, an explicit representation of a joint manifold is learned in the form of two linearly connected latent codes. Results in [55] show a clear improvement over non-learned state-of-the-art methods for in-vivo data.

These works display the potential for fully learned methods in image reconstruction and inverse problems in general: high quality reconstructions are obtained, while no exact knowledge of underlying physics or specifics of the measurement system is required.

5.3. Manifold learning. Many relevant inverse problems in medicine, engineering, astronomy or geophysics are large-scale in both signal and measurement space. However, seen from a statistical point of view, probability mass concentrates near manifolds of a much lower dimension than the original data spaces [5]. To detect linear manifolds, principal component analysis is a suitable and simple method. However, since manifolds for real data are expected to be strongly nonlinear [5], one needs to make use of nonlinear techniques. One of the best known methods that achieves this is kernel PCA [45]. In this method, data is mapped to a reproducing kernel Hilbert space by applying a nonlinear kernel, after which the linear PCA is applied. Other methods are principal geodesic analysis [21], which can be applied for Riemannian manifolds and geodesic PCA [9], which acts in a Wasserstein space, which is nonlinear.

Another approach to learn nonlinear manifolds is to use autoencoders, which can also be seen as a generalisation of linear PCA [26]. Autoencoders have shown to learn explicit representations of nonlinear manifolds [41] and provide better low-dimensional latent code in terms of clustering and reconstruction performance [26] than their linear counterpart. For the inverse problems (2.1) considered in our paper, there is an explicit relationship between signals and measurement via the forward mapping that models the physics. This means that signals and measurements share one data manifold that is learned by connecting two autoencoders.

The observation of a shared manifold, or “an unknown underlying relationship between two domains” [56], is also the idea driving the cycle GAN. Unlike in our paper, the goal of the cycle GAN is not to find a unique and supervised one-to-one mapping from one domain to the other, but to identify the shared parameters and add such elements so that the output is realistic in its respective domain. This cycle GAN was applied for inverse problems in different forms [47, 48]. In these works, the manifold is implicitly learned, unlike the explicitly learned representation that we study in our paper.

5.4. Transfer learning with autoencoders. Transfer learning is used to exploit similarities between different tasks to share information necessary for both tasks. Representation learning, such as manifold learning, has a strong influence in transfer learning scenarios [5] since the learned representation can guide the supervised reconstruction task. Recent research has shown that autoencoders can be used as a regulariser for a supervised training task, such as classification [54, 31]. Such networks, coined supervised autoencoders (SAE), help to generalise the supervised problem. They are specifically useful in a semi-supervised scenario, where a lot of unsupervised training data is available, but supervised training pairs are scarce.

L-SVD profits from the same regularisation and generalisation effect by attaching two autoencoders to the supervised reconstruction problem. For inverse problems a semi-supervised scenario is also often encountered: imagine a training set of undersampled medical MRI or CT data sets and a set of high-quality reconstructions; not all pairs are available because not all patients have had a fully sampled scan that is needed for a high-quality reconstruction.

5.5. Model reduction and learning. Model reduction is a mathematical and computational field of study that derives low-dimensional models of complex systems [14, 6]. Via projections and decompositions it is possible to represent approximations of large-scale high-fidelity computational models resulting from discretization of partial differential equations. Recent developments focus on data-driven learning of governing equations [13, 44, 39] and learned model reduction [38].

Our work focuses on learning the inverse map for problems that are often physics-driven. The nonlinear equations or parameters of this map are implicitly learned through the latent representations via autoencoders.

5.6. Bayesian inversion and sparsity. The goal of the Bayesian approach to solving inverse problems is to find the posterior measure, given sampled data and a prior measure [50]. The posterior measure will contain information about the relative probability of different outputs, given the data. Often the posterior is too complex to recover and the goal is shifted to finding a maximum a posteriori (MAP) estimate. In Section 4.1 we showed that a linear variant of L-SVD coincides with learning the covariance matrix of the prior measure.

Yu et al. [52] developed piecewise linear estimates (PLE), a method based on Gaussian mixture models (GMM), which are estimated via the MAP expectation-maximization algorithm (MAP-EM). The method makes use of GMM as prior measures on local patches, which results in a linear reconstruction model for each patch. One could think of this procedure as finding a patchwork of locally linear tangent spaces which approximate a nonlinear manifold. If measurement and signal domain are the same, it can be shown that PLE is equivalent to learning a linear L-SVD method for a group of similar patches.

6. Experiments and implementation. In this section, we explain three simulation ex- periments: two that demonstrate the contributions as stated in Section 1.1 in a relatively low-resolution scenario and one that shows the application of L-SVD to a biomedical data-set in a higher resolution. The forward operator in all experiments is chosen to be the Radon transform [29], which is a nonlocal linear operator that produces an ill-posed inverse problem. We will refer to the measurements as ‘sinograms’ and to the signals as ‘images’ in their respective spaces Y and X. In Appendix C, the neural network architectures and their parameter choices are provided, together with details of the training set.

6.1. Experiment 1: from model-based to data-driven. The goal of this simulation ex- periment is to demonstrate the first perspective of Section 1.1: data-driven solution of inverse problems. This is done by comparing classical non-learned methods to learned methods. For fair comparison and clarity, only linear inversion methods are considered.

This experiment makes use of the MNIST data set [32], after rescaling it to 64We apply the Radon transform with 64 uniformly samples angles, i.e. a ‘full-angle’ Radon operator A, with a discretisation of 64 for each angular view. Moreover, there is no bottleneck latent space. This results in equally large spaces . For training, the ‘clean’ simulated full-angle measurements are normalised, after which Gaussian noise with a noise level 05 is added. The following linear reconstruction methods are compared:

(a) Tikhonov-regularised reconstruction with optimally chosen (see Section 2.2). (b) Truncated SVD reconstruction with optimal truncation number r (see Section 2.3). (c) Optimal regularised inverse matrix (ORIM) [16, Theorem 2.1]. (d) Data-driven Tikhonov regularisation, i.e. U and V are from the SVD, while Σ is a learned diagonal matrix with the structure of (4.3) (see Section 4.2).

(e) Reconstruction from a learned covariance matrix of the prior, i.e. U and V are from the SVD, while Σ is a full matrix that is learned (see Section 4.1.1).

(f) Fully learned L-SVD: nonlinear L-SVD as in Definition 4.6 (see Section 4.3). A regular autoencoder is used on the sinogram side, which means that noise should be reconstructed after the sinogram is encoded and decoded. By training on sinograms with noise levels between 0 and 0.2 instead of the aforementioned 0.05, the effect of noise on the matrix Σ is investigated. Instead of training for several noise levels individually, only one training set is created, where each sinogram has a randomly chosen noise level 2]. To process this data set adequately, the static scaling matrix Σ is exchanged for a noise-aware component that depends on the noisy input data (c.f. Section 3.1.1). This component is a nonlinear fully-connected network which takes as input and provides the diagonal scales () as output, which are multiplied with

6.2. Experiment 2: from linear autoencoding to hybrid nonlinear autoencoding. The goal of this simulation experiment is to demonstrate the second perspective of Section 1.1: nonlinear encoding is more effective than linear encoding; moreover, combining two nonlinear autoencoders has a regularising effect on the reconstruction and gives a more insightful latent representation than one autoencoder. For this experiment, we make use of the fully learned L-SVD.

This experiment again makes use of the MNIST data set [32], after rescaling it to 64 pixels. A limited-angle Radon transform of 8 uniformly sampled angles is applied, with a discretisation of 64 for each angular view. The latent space is chosen to have 64 dimensions, which means that it acts as a bottleneck. This results in the spaces , providing a dimensionality reduction of 12.5% and 1.56% compared to sinogram and image space respectively. Gaussian noise with a noise level 05 is added to the limited-angle measurements. We analyse the following methods:

(a) linear autoencoder;

(b) nonlinear autoencoder;

(c) linear L-SVD;

(d) nonlinear L-SVD (

(e) nonlinear L-SVD.

Here the first two methods are only applied on the image side and not on the sinogram side, meaning that the autoencoder only connects . As can be seen in Figure 3.1, L-SVD connects the sinogram side with the image side, where a denoising autoencoder (DAE) is used on the sinogram side. This means that noise should not be reconstructed after decoding to Y , which allows for a dimensionality reduction. The fourth method has the same network structure of nonlinear L-SVD, but without the autoencoders on either side (i.e. in (3.1)), which allows us to investigate the regularising effect of the autoencoders. It is obvious that this experiment investigates the second perspective from Section 1.1, since the transition from linear to nonlinear is made, as well as the transition from a single to hybrid autoencoder.

Finally, we also compare nonlinear L-SVD for the following three cases:

(a) All supervised training pairs () are available.

(b) Only 10% of the training samples is available, all in pairs (

(c) All training samples are available, but only 10% is paired. In the third case, if the training data is unpaired, then the encoders and decoders are only trained by the autoencoders. That is, the reconstruction loss ) is set to zero. We investigate how L-SVD can exploit the semi-supervised case in which additional unpaired training data is available.

6.3. Experiment 3: L-SVD on human chest CT images. The third simulation experiment demonstrates the capacity of L-SVD to reconstruct biomedical images without any knowledge on the forward operator. It exploits the Low-Dose Parallel Beam CT (LoDoPaB) data set [33], which is based on the LIDC/IDRI data set [3]. In our paper, we only make use of the high quality CT reconstructions in the LoDoPaB-CT data set that we use as ‘ground truth’ for our setup. The images are scaled to 128 128 pixels, after which a Radon transform of 36 uniformly sampled angles is applied with a discretisation of 192 for each angular view. The latent space of L-SVD is chosen to have 2048 dimensions. This results in the spaces , providing a dimensionality reduction of around 29.6% and 12.5% compared to sinogram and image space respectively. Gaussian noise with an signal-to-noise ratio (SNR) of 40dB is added to the measurements. The following reconstruction methods are compared:

(a) T-SVD with optimal truncation number

(b) optimal regularised inverse matrix (ORIM) [16];

(c) total variation regularisation with optimal regularisation parameter

(d) fully learned nonlinear L-SVD.

Truncated SVD can be seen as the linear counterpart of L-SVD and makes use of the forward operator A. To obtain at least the same dimensionality reduction as L-SVD, the truncation number r of T-SVD is chosen smaller than 2048, but such that it gives the best reconstruction quality. ORIM is a linear data-driven method, for which we use the same training set as for L-SVD. For computational reasons, we make use of [16, equation (12)]: the ORIM method that requires knowledge about the mean and covariance of signals Moreover, the forward operator A is needed for this method. Total variation (TV) [42] is an adequate nonlinear regularised reconstruction method for this problem, since the ground truth CT-images consist of almost entirely piecewise constant structures. This method also needs the forward operator A: we minimise the functional min), implemented as described in [11].

7. Numerical results. In this section, the results of the simulation experiments explained in Section 6 are shown and discussed.

7.1. Experiment 1: from model-based to data-driven. The results of the first simulation experiments for a randomly chosen sample in the test set are shown in Figure 7.1. For this sample, Gaussian noise with a noise level 05 was added. Visually, the reconstruction improves gradually as we move from model-based methods with only one tunable parameter (b,c) to combinations of model- and data-driven methods with tunable scaling matrix (d,e) to fully data-driven methods (f,g,h). This is most noticeable in the background of the reconstruction, which should be constant. The visual improvement is confirmed by the peak signal-to-noise ratio (PSNR), displayed above each reconstruction. These values represent the mean and standard deviation of the PSNR values for the first 1000 images in the test set.

Figure 7.1: Reconstructions ˆ(see Figure 3.1) for different models ranging from fully model-based (top left) to more data-driven (bottom right) with increasing amount of learning.

Next, we analyse the effect of noise on the diagonal scaling Σ with the network component as explained at the end of Section 6.1. We compare the networks where encoder and decoder are fixed as with the fully learned networks. All networks are compared with Tikhonov regularisation, in which only one tunable parameter is chosen such that the smallest MSE is obtained. Gaussian noise is added with 6 different noise levels to all sinograms in the test set. The average scales per noise level are shown in Figure 7.2, where (a-c) have the same ordering as the classical SVD and (d) is ordered from large to small scales. For visualisation purposes, the graphs were smoothed by a Gaussian filter with scale 10.

Figure 7.2: Comparison of the scales depending on noise level for methods with increasing amount of learning. Although similar, methods (b-d), where the scales are learned individually, show a greater noise dependency than (a) Tikhonov regularisation. Note that methods (a-c) use the SVD ordering, while (d) is ordered from large to small scales at the highest noise level.

All methods show a decay of scales as the noise level grows. Tikhonov regularisation shows a very similar decay over all latent dimensions, while the decay is much more dimension dependent in (b-d). Here, the scales that coincide with large in the SVD case, i.e. the first dimensions, are relatively noise independent. The scales that coincide with small SVD case, i.e. dimensions 500-3000, show a relatively large decay, meaning that they are more affected by noise. Finally, the last part of all graphs show more than a thousand dimensions that have a small scale for all noise levels. To sum up, the behaviour of the scales to noise is similar for all methods in which the scales are learned individually, regardless of the encoding and decoding used. Moreover, most structural information is encoded in a limited number of latent dimensions, while the other dimensions encode of a substantial amount of noise: a compression is learned.

7.2. Experiment 2: from linear autoencoding to hybrid nonlinear autoencoding. The results of the second simulation experiments for a randomly chosen sample in the test set are shown in Figure 7.3. For the conciseness of this section, only ˆis shown for all variants of L-SVD: not much difference between ˆwas observed. For the autoencoders ˆshown, since ˆis not available there. Note that for ˆ, no noise was added to their inputs x, while for the reconstruction outputs ˆ, noise with a noise-level of 0.05 was added to their inputs

Figure 7.3: Comparison of the outputs ˆ(see Figure 3.1) for linear and nonlinear variants of the AE and fully learned L-SVD network.

It can be seen that the linear methods produce side-lobes to the main intensities, which are often observed in frequency-based compressions. The nonlinear methods provide a more homogeneous background in the reconstruction, especially L-SVD.

Table 7.1: Generalisation performance analysis by comparing train and test losses (MSE). Nonlinear L-SVD with hybrid autoencoder shows a smaller difference between train and test loss than other nonlinear methods, indicating better generalisation performance.

In Table 7.1, the train and test losses of all compared methods are shown. The first thing that can be seen is that for the linear networks, errors are larger then for the nonlinear networks. The second thing is that there is no significant difference between train and test loss for the linear networks, which indicates that they generalise well. This difference is present in the nonlinear networks, but for L-SVD not to the same extent as for AE en LSVD. From the autoencoder point of view, it seems that L-SVD benefits from the sinogram branch of the network in terms of generalisation. From a reconstruction point of view, L-SVD benefits from the incorporation of the two autoencoders, which contribute to its generalisation capacity in the reconstruction output ˆ. Finally, since the test loss for nonlinear L-SVD is smaller than for L-SVD, we conclude that the two autoencoders act as regularisers for the reconstruction.

Next, canonical basis vectors in the latent space are decoded to the image space, to compose a ‘dictionary’ of elements. While in a linear case all outputs could be reconstructed from these elements in a linear way, this is not true for the nonlinear case. This means that this dictionary only gives a partial view on the decoder.

Figure 7.4 shows four selected elements that are exemplary for the complete dictionaries, which are provided in Appendix E. Because the elements of linear AE are very similar to linear L-SVD, they are not shown here. In the top left, linear L-SVD shows an element that is very similar to the Euclidean mean of all training samples, while the top right and

Figure 7.4: Selected elements in the latent space , decoded to the image space X. Nonlinear L-SVDand nonlinear L-SVD learn a more interpretable representation than other methods by combining features from sinogram and image space in their ‘dictionary’. Due to its similarity to linear L-SVD, linear AE is not shown here. All 64 elements are shown in Appendix E.

bottom images show low and high frequency components that are only active in the middle. Nonlinear AE shows much smaller structures in its elements, which do not seem to have a visual coherent structure. The nonlinear networks L-SVDand L-SVD do provide this visually coherent structure, where L-SVD seems to provide somewhat ‘smoother’ and better connected structures than L-SVD. Their dictionaries consist of various elements, of which one is similar to the Euclidean mean (top left), some are similar to digits (top right), some that show a combination of line segments (bottom left) and some with high-frequency components (bottom right) which were also visible in the linear dictionary. With this diversity, fully learned L-SVD combines information from images, sinograms and operator.

Table 7.2: Comparison of test losses (MSE) for supervised L-SVD and semi-supervised L-SVD, where 90% of the data is only trained by the autoencoders incorporated in L-SVD.

Finally, the capacity of L-SVD in a semi-supervised setup is shown in Table 7.2. Only 10% of the training samples are available in pairs (i.e. supervised), and the other 90% are available, but not in pairs (i.e. unsupervised). Table 7.2 shows that the semi-supervised setup provides a highly increased performance over the supervised case where only the pairs (i.e. 10%) are used. This demonstrates the advantage of adding the autoencoders in L-SVD: they help to efficiently encode and decode, even if pairs between data and signal are not available.

7.3. Experiment 3: L-SVD on human chest CT images. Here we demonstrate the capacity of L-SVD to reconstruct biomedical images without any knowledge on the forward operator. We compare it to T-SVD, ORIM and TV, which are all methods that make use of the forward operator A. For a randomly chosen sample in the test set, the outputs of all methods are shown in Figure 7.5. It can be seen that none of the methods can reconstruct all details that are apparent in the ground truth. This is probably due to the limited number of angles that were used and the noise added to the sinogram. The T-SVD and ORIM reconstructions give typical ringing artefacts, while the TV reconstruction has the typical staircase behaviour with clusters of piecewise constant structures. The L-SVD reconstruction is smoother than the other methods, which can be seen in the large piecewise constant structure in the bottom left of the image, as well as close to the boundary of the imaged body, where a low contrast between tissues should be reconstructed. To get a better idea of the reconstruction biases of all methods, three additional samples from the test set are provided in Appendix F.

Figure 7.5: Comparison between fully learned L-SVD and reconstruction methods that make use of the forward operator A.

In Table 7.3, the peak signal-to-noise ratio (PSNR) and structural similarity (SSIM) of all compared methods are given.

truncated SVD (T-SVD) 33optimal regularised inverse matrix (ORIM) 34total variation (TV) 35learned SVD (L-SVD) 36

Table 7.3: Comparison in PSNR and SSIM between L-SVD and non-learned reconstruction methods that make use of the forward operator A.

Firstly, it can be seen that TV and L-SVD, both nonlinear methods, provide a higher quality than T-SVD and ORIM, linear reconstruction methods. Secondly, from the methods that contain a bottleneck in their architecture, L-SVD is superior to T-SVD. This can be due to its nonlinearity or the fact that it is a learned method. Thirdly, from the data-driven methods, L-SVD shows better results than ORIM, which is probably due to the nonlinearity of L-SVD. Finally, TV and L-SVD are comparable in terms of PSNR and SSIM. However, to use TV as a reconstruction method, it is necessary to know the forward operator A exactly, which is not always possible. Moreover, this quality of TV comes at the expense of computation power: L-SVD is a direct reconstruction method, while TV needs close to a 1000 iterations to converge to its solution. On the other hand, L-SVD requires a large training set and training time, while there is only one parameter to tune for TV.

8. Conclusion and outlook. We proposed the learned SVD for inverse problems: a recon- struction method that connects low-dimensional representations of corrupted measurements and desired signals. When the forward mapping is known, the L-SVD can be shaped into a data-driven Tikhonov regularisation method, for which we provided a convergence analysis. When the forward mapping is unknown, nonlinear representations of measurements and signals can be fully learned from data, with a connecting layer that is fully connected or sparse, linear or nonlinear, noise dependent or independent. One specific choice is to incorporate the necessary nonlinearity of the learned inverse function in the autoencoders, while the sparse diagonal scaling layer is chosen to be linear, making the connection between measurement and signal manifold easy to understand. In simulation experiments, it was shown that this nonlinear reconstruction gives superior performance to other methods, while providing interpretable autoencoding. Moreover, since the reconstruction error estimate depends on the autoencoding quality, L-SVD can benefit from general advances in nonlinear autoencoding.

Results show that L-SVD makes use of information from both measurements and signals; by doing so, it learns elements of the physics operator, although not explicitly provided. Therefore the method is especially promising in applications where the forward physics are not completely understood or computationally expensive to simulate. Learning a joint manifold by two connected autoencoders also enables the possibility of a semi-supervised setup: the autoencoders provide regularisation for reconstruction of the signal.

Due to its generic formulation, L-SVD is very flexible for other architecture choices in autoencoding. Therefore, future efforts will lie in investigating other architectures, such as convolutional autoencoders and ladder variational autoencoders [49], for their inclusion in L-SVD for large scale inverse problems. Furthermore, other loss functions such as the Wasserstein loss or a learned discriminator could be investigated. Finally, we will focus on other ways of incorporating (partial) information of the forward mapping in the architecture of the L-SVD method, to create a regularisation method with a convergence analysis that benefits from the advantages of nonlinear autoencoding.

Acknowledgements. We thank Srirang Manohar, Johannes Schwab, Kathrin Smetana and Leonie Zeune for valuable discussions on the topic. We thank Johannes Schwab for providing his numerical implementation of the Radon transform. The collaboration project is co-funded by the PPP allowance made available by HealthHolland, Top Sector Life Sciences & Health, to stimulate public private partnerships. CB acknowledges support by the 4TU programme Precision Medicine and the European Union Horizon 2020 research and innovation programme under the Marie Skodowska-Curie grant agreement No. 777826 NoMADS.