Discretization of Learned NETT Regularization for Solving Inverse Problems

Deep learning based reconstruction methods deliver outstanding results for solving inverse problems and are therefore becoming increasingly important. A recently invented class of learning-based reconstruction methods is the so-called NETT (for Network Tikhonov Regularization), which contains a trained neural network as regularizer in generalized Tikhonov regularization. The existing analysis of NETT considers fixed operators and fixed regularizers and analyzes the convergence as the noise level in the data approaches zero. In this paper, we extend the frameworks and analysis considerably to reflect various practical aspects and take into account discretization of the data space, the solution space, the forward operator and the neural network defining the regularizer. We show the asymptotic convergence of the discretized NETT approach for decreasing noise levels and discretization errors. Additionally, we derive convergence rates and present numerical results for a limited data problem in photoacoustic tomography.


Introduction
In this paper, we are interested in neural network based solutions to inverse problems of the form Find x from data y δ = Ax + η .
Here A is a potentially non-linear operator between Banach spaces X and Y, y δ are the given noisy data, x is the unknown to be recovered, η is the unknown noise perturbation and δ ≥ 0 indicates the noise level. Numerous image reconstruction problems, parameter identification tasks and geophysical applications can be stated as such inverse problems [1][2][3][4]. Special challenges in solving inverse problems are the non-uniqueness of the solutions and the instability of the solutions with respect to the given data. To overcome these issues, regularization methods are needed, which select specific solutions and at the same time stabilize the inversion process.

Reconstruction with Learned Regularizers
One of the most established class of methods for solving inverse problems is variational regularization where regularized solutions are defined as minimizers of the generelaized Tikhonov functional [2,5,6] T y δ ,α : X → [0, ∞] : x → D(Ax, y δ ) + αR(x) . ( Here D is a distance like function measuring closeness of the data, R a regularization term enforcing regularity of the minimizer and α is the regularization parameter. Taking minimizers of this functional as regularized solution is also called (generalized) Tikhonov regularization. In the case that D and the regularizer are defined by the Hilbert space norms, (2) is classical Tikhonov regularization for which the theory is quite complete [1,7]. In particular, in this case, convergence rates, which name quantitative estimates for the distance between the true noise-free solution and regularized solutions from noisy data, are well known. Convergence rates for non-convex regularizers are derived in [8].
Typical regularization techniques are based on simple hand crafted regularization terms such as the total variation f TV = |∇ f | or quadratic Sobolev norms ∇ f 2 2 = |∇ f | 2 on some function space. However, these regularizers are quite simplistic and might not well reflect the actual complexity of the underlying class of functions. Therefore, recently, it has been proposed and analyzed in [9] to use machine learning to construct regularizers in a data driven manner. In particular, the strategy in [9] is to construct a data-driven regularizer via the following consecutive steps: Choose a family of desired reconstructions ( For imaging applications, the function class (Φ θ ) θ∈Θ can be chosen as convolutional neural networks which have demonstrated to give powerful classes of mappings between image spaces. The function r measures distance between a potential reconstruction x and the output of the network Φ(x), and possibly contains additional regularization [10,11]. According to the training strategy in item (T4) the value of the regularizer will be small if the reconstruction is similar to elements in (x i ) n i=1 and large for elements in (BAx i ) n i=1 . A simple example that we will use for our numerical results is the learned regularizer Convergence analysis and convergence rates for NETT (which stands for Network Tikhonov; referring to variants of (2), where the regularization term is given by a neural network) as well as training strategies have been established in [9,11,12]. A different training strategy for learning a regularizer has been proposed in [13,14]. Note that learning the regularizer first and then minimizing the Tikhonov functional is different from variational and iterative networks [15][16][17][18][19][20] where an iterative scheme is applied to enroll the functional D θ (Ax, y δ ) + αR θ (x) which is then trained in an end to end fashion. Training the regularizer first has the advantage of being more modular, sharing some similarity with plug and play techniques [21], and the network training is independent of the forward operator A. Moreover, it enables to derive a convergence analysis as the noise level tends to zero and therefore comes with theoretical recovery guarantees.

Discrete NETT
The existing analysis of NETT considers minimizers of the Tikhonov functional (2) with regularizer of the form R(x) = r(x, Φ(x)) before discretization, typically in an infinite dimensional setting. However, in practice, only finite dimensional approximations of the unknown, the operator and the neural network are given. To address these issues, in this paper, we study discrete NETT regularization which considers minimizers of T y δ ,α,n : X n → Y : x → D(A n x, y δ ) + αR n (x) .
Here (X n ) n∈N , (A n ) n∈N and (R n ) n∈N are families of subspaces of X n ⊆ X, mappings A n : X → Y and regularizers R n : X → [0, ∞], respectively, which reflect discretization of all involved operations. We present a full convergence analysis as the noise level δ converges to zero and n, α are chosen accordingly. Discretization of variational regularization has studied in [22] for the case that D is given by the norm distance and the regularizer R is taken convex and fixed. However, in the case of discrete NETT regularization it is natural to consider the case where the regularization depends on the discretization as regularization is learned in a discretized setting based on actual data. For that purpose our analysis includes non-convex regularizers that are allowed to depend on the discretization and the noise level.

Outline
The convergence analysis including convergence rates is presented in Section 2. In Section 3 we will present numerical results for a non-standard limited data problem in photoacoustic tomography that can be considered as simultaneous inpainting and artifact removal problem. We conclude the paper with a short summary and conclusion presented in Section 4.

Convergence Analysis
In this section we study the convergence of (3) and derive convergence rates.

Well-Posedness
First we state the assumptions that we will use for well-posedness (existence and stability) of minimizing NETT.
Conditions (W2)-(W5) are quite standard for Tikhonov regularization in Banach spaces to guarantee the existence and stability of minimizers of the Tikhonov functional and the given conditions are similar to [2,[8][9][10]12,23,24]. In particular, (W2) describes the properties that the distance measure D should have. Clearly, the norm distance on Y fulfills these properties. Moreover, (W2a) holds for the norm with τ = 1 since it then corresponds to the triangle inequality. Item (W2c) is the continuity of D(y, ·) while (W2d) considers the continuity of D(·, y) at y. While (W2c) is not needed for existence and convergence of NETT it is required for the stability result as shown in [10] (Example 2.7). On the other hand (W2e) implies that the Tikhonov functional is wslsc which is needed for existence. Assumption (W5) is a coercivity condition; see [9] (Remark 2.4f.) on how to achieve this for a regularizer defined by neural networks. Item (W8) poses some restrictions on the regularizers. For NETT this is not an issue as neural networks used in practice are continuous. Note that for convergence and convergence rates we will require additional conditions that concern the discretization of the reconstruction space, the forward operator and regularizer.
The references [8][9][10]23] all consider general distance measures and allow non-convex regularizers. However, existence and stability of minimizing (2) are shown under assumptions slightly different from (W1)-(W5). Below we therefore give a short proof of the existence and stability results. Theorem 1 (Existence and Stability). Let Assumption 1 hold. Then for all y δ ∈ Y, α > 0, n ∈ N the following assertions hold true: Let (y k ) k∈N ∈ Y N with y k → y and consider x k ∈ argmin T y k ,α,n .
• (x k ) k∈N has at least one weak accumulation point.
The statements in (a),(b) also hold for T y,α in place of T y,α,n , Proof. Since (W1), (W6)-(W9) for T y,α,n when n ∈ N are fixed give the same assumption as (W1), (W3)-(W5) for the non-discrete counterpart T y,α , it is sufficient to verify (a), (b) for the latter. Existence of minimizers follows from (W1), (W2e), (W3)-(W5), because these items imply that the T y,α is a wslsc coercive functional defined on a nonempty weakly sequentially closed subset of a reflexive Banach space. To show stability one notes that according to (W2a) for all x ∈ X we have According to (W2c), (W2d), (W5) there exists x ∈ X such that the right hand side is bounded, which by (W5) shows that (x k ) k has a weak accumulation point. Following the standard proof [2] (Theorem 3.23) shows the weak accumulation points (x k ) k are minimizers of T y,α . This uses the fact that the weak topology is indeed weaker than the norm topology, and that the involved functionals are wslsc.
In the following we write x δ α,n for minimizers of T y δ ,α,n . For Lemma 1 (Existence of R-minimizing solutions). Let Assumption 1 hold. For any y ∈ A(D) an R-minimizing solution of Ax = y exists. Likewise, if n ∈ N and y ∈ A n (D) an R n -minimizing solution of A n x = y exists.
Proof. Again is is sufficient the verify the claim for R-minimizing solution. Because is bounded according to (W5). By (W1) X is reflexive and therefore (x k ) k∈N has a weak accumulation point x + . From (W1), (W4), (W3) we conclude that x + is an R-minimizing solution of Ax = y. The case of R n -minimizing solutions follows analogous.

Convergence
Next we proof that discrete NETT converges as the noise level goes to zero and the discretization as well as the regularization parameter are chosen properly. We write D n,M := {x ∈ D ∩ X n | R n (x) ≤ M} and formulate the following approximation conditions for obtaining convergence.
Assumption 2 (Conditions for convergence). Element x + ∈ D satisfies the following for all M > 0: Conditions (C1) and (C3) concerns the approximation of the true unknown x with elements in the discretization space, that is compatible with the discretization of the forward operator and regularizer. Conditions (C2) and (C4) are uniform approximation properties of the operator and the regularizer on R n -bounded sets.
Theorem 2 (Convergence). Let (W1)-(W9) hold, y ∈ A(D) and let x + be an R-minimizing solution of Ax = y that satisfies (C1)-(C4). Moreover, suppose (δ k ) k∈N ∈ (0, ∞) N converges to zero and (y k ) k∈N ∈ Y N satisfies D(y, y k ) ≤ δ k . Choose (α k ) k∈N and (n k ) k∈N such that as k → ∞ we have Then for x k ∈ argmin T y k ,δ k ,n k the following hold: (a) (x k ) k∈N has a weakly convergent subsequence (x σ(k) ) k∈N (b) The weak limit of (x σ(k) ) k∈N is an R-minimizing solution of Ax = y.
If the R-minimizing solution of Ax = y is unique, then (x k ) k∈N x + .
Proof. For convenience and some abuse of notation we use the abbreviations R k := R n k , A k := A n k , a k := a n k , z k := z n k and ρ k := ρ n k . Because x k is a minimizer of the discrete NETT functional T y k ,δ k ,n k by (W2) we have According to (C1), (4), we get According to (C1), (C3), (5), (6) the right hand side in (7) converges to zero and the right hand side in (8) to R(x + ). Together with (C2) we obtain R(x k ) ≤ R k (x k ) + ρ k → R(x + ) and D(Ax k , y) ≤ D(A k x k , y) + a k ≤ τD(A k x k , y) + a k + τδ k → 0. This shows that (D(Ax k , y) + R(x k )) k∈N is bounded and by (W1), (W9) there exists a weakly convergent subsequence (x σ(k) ) k∈N . We denote the weak limit by x ∈ X. From (W2), (W4) we obtain Ax = y. The weak lower semi-continuity of R assumed in (W3) shows Consequently, x is an R-minimizing solution of Ax = y and R(x σ(k) ) → R(x ). If the R-minimizing solution is unique then x + is the only weak accumulation point of (x k ) k∈N which concludes the proof.

Convergence Rates
Next we derive quantitative error estimates (convergence rates) in terms of the absolute Bregman distance. Recall that a function R : X → [0, ∞] is Gâteaux differentiable at some x ∈ X if the directional derivative R (x )(h) := (R(x + th) − R(x ))/t exist for every h ∈ X. We denote by R (x ) the Gâteaux derivative of R at x. In [9] we intro- We write sup y δ H(y δ ) := sup{H(y δ ) | y δ ∈ X ∧ D(Ax + , y δ ) ≤ δ}. Convergence rates in terms of the Bregman distance are derived under a smoothness assumption on the true solution in the form of a certain variational inequality. More precisely we assume the following: Assumption 3 (Conditions for convergence rates). Element x + ∈ D satisfies the following for all M, δ > 0: There exist a concave, continuous, strictly increasing ϕ According to (R5) the inverse function ϕ −1 : [0, ∞) → [0, ∞) exists and is convex. We denote by ϕ − * (s) := sup{s t − ϕ −1 (t) | t ≥ 0} its Fenchel conjugate.
For the second inequality we used (C1) and (C2). We have D(A n x δ α,n , y δ ) + αR n (x δ α,n ) ≤ D(A n z n , y δ ) + αR n (z n ) and thus we get an estimate for R n (x δ α,n ) − R n (z n ) which we used for the third inequality. For the next inequality we used (R2) and (R3). Finally we used Young's inequality αϕ(τt) ≤ t + τ −1 ϕ − * (τα) for the last step.

Remark 1.
The error estimate (10) includes the approximation quality of the discrete or inexact forward operator A n and the discrete or inexact regularizer R n described by a n,δ and ρ n , respectively. What might be unexpected at first is the inclusion of two new parameters λ n and γ n,δ . These factors both arise from the approximation of X by the finite dimensional spaces X n , where γ n,δ reflects approximation accuracy in the image of the operator A and λ n approximation accuracy with respect to the true regularization functional R. Note that in the case where the forward operator, the regularizer, and the solution space X are given precisely, we have a n,δ = γ n,δ = λ n = ρ n = 0. In this particular case we recover the estimate derived for the NETT in [9].
Next we verify that a variational inequality of the form (R5) is satisfied with ϕ(t) = c √ t under a typical source like condition.
Lemma 2 (Variational inequality under source condition). Let R, A be Gâteaux differentiable at x + ∈ X, consider the distance measure D(y 1 , y 2 ) = y 1 − y 2 2 and assume there exist η ∈ X and c 1 , Using the Cauchy-Schwarz inequality and Equation (12), we can estimate Putting this together we get

Corollary 1 (Convergence rates under source condition). Let the conditions of Lemma 2 hold and suppose
Then we have the convergence rates result Proof. Follows from Theorem 3 and Lemma 2. Note that we use · in the theorem, while D(y 1 , y 1 ) = y 1 − y 2 2 uses the squared norm · 2 and thus the approximation rates for the terms concerning A n(δ) are order √ δ instead of δ as in Theorem 3.
In Corollary 1, the approximation quality of the discrete operator A n and the discrete and inexact regularization functional R n need to be of the same order.

Application to a Limited Data Problem in PAT
Photoacoustic Tomography (PAT) is an emerging non-invasive coupled-physics biomedical imaging technique with high contrast and high spatial resolution [25,26]. It works by illuminating a semi-transparent sample with short optical pulses which causes heating of the sample followed by expansion and the subsequent emission of an acoustic wave. Sensors on the outside of the sample measure the acoustic wave and these measurements are then used to reconstruct the initial pressure f : R d → R, which provides information about the interior of the object. The cases d = 2 and d = 3 are relevant for applications in PAT. Here we only consider the case d = 2 and assume a circular measurement geometry. The 2D case arises for example when using integrating line detectors in PAT [26].

Discrete Forward Operator
The pressure data p : R 2 × [0, ∞) → R satisfies the wave equation (∂ 2 t − ∆)p(r, t) = 0 for (r, t) ∈ R 2 × (0, ∞) with initial data p( · , 0) = and ∂ t p( · , 0) = 0. In the case of circular measurement geometry one assumes that f vanishes outside the unit disc D 1 := {r ∈ R 2 | r < 1} and the measurement sensors are located on the boundary ∂D 1 = S 1 . We assume that the phantom will not generate any data for some region I ⊆ D 1 , for example when the acoustic pressure generated inside I is too small to be recorded. This masked PAT problem consists in the recovery of the function f from sampled noisy measurements of g = W (1 I c f ) where W denotes the solution operator of the wave equation and 1 I c the indicator function on I c := R 2 \ I. Note that the resulting inverse problem can be seen of the combination of an inpainting problem and in inverse problems for the wave equation.
Here I m is the modified Bessel function of the first kind of order n ∈ N and the parameters γ > 0 and R denote the window taper and support radius, respectively.
. For convenience we will use a pseudo-3D approach where use the 3D solution of W ψ for which there exists an analytical representation [29]. Denote by s k uniformly spaced sensor locations on S 1 and by t j > 0 uniformly sampled measurement times in [0,2]. Define the N t N s × N 2 model matrix by W N t (k−1)+j,N(i 1 −1)+i 2 = W (ψ( · − r i ))(s k , t j ) and an N 2 × N 2 diagonal matrix by (M I ) N(i 1 −1)+i 2 ,N(i 1 −1)+i 2 = 1 if r i ∈ I c and zero otherwise. Let WM I = UΣV be the singular valued decomposition. We then consider the discrete forward matrix A = UΣ V where Σ is the diagonal matrix derived from Σ by setting singular values smaller than some σ to zero. This allows us to easily calculate A + = VΣ + U where Σ + is calculated by inverting all diagonal elements of Σ that are greater than zero. In our experiments we use N = N t = 128, N s = 150 and take I fixed as a diagonal stripe of width 0.34.

Discrete NETT
We consider the discrete NETT with discrepancy term D(Ax, y δ ) = Ax − y δ 2 2 /2 and regularizer given by where ∇x 1, := ∑ 128 version of the total variation [31] and Φ (m) is a learnable network. We take Φ (m) as the U-Net [32] with residual connection, which has first been applied to PAT image reconstruction in [33]. Here m stands for the number of down-/upsampling steps performed in the U-Net (the original one had m = 4 i.e., m ∈ N. This means that larger m yield a deeper network with more parameters. We generate training data that consist of square shaped rings with random profile and random location. See Figure 1 for an example of one such phantom (note that all plots in signal space use the same colorbar) and the corresponding data. We get a set of phantoms x 1 , . . . , x 1000 and corresponding basic reconstructions h a := A + (Ax a + η a ), where A + is the pseudo-inverse and η a is Gaussian noise with standard deviation of σ Ax a ∞ with σ = 0.01. The networks are trained by minimizing ∑ 1000 where we used the Adam optimizer with learning rate 0.01 and γ = 0.1. The considered loss is that we want the trained regularizer to give small values for x a and large values for h a . The strategy is similar to [9] but we use the final output of the network for the regularizer as proposed in [34]. To minimize (15) we use Algorithm 1 which implements a forward-backward scheme [35]. The most expensive step of this algorithm is the matrix inversion but since we use constant stepsize one also has to option to only calculate the inverse of the matrix once and reuse it. Thus one only has to perform two matrix-vector multiplications which are of the order O(N 2 N t N s ) and O(N 4 ) since N 2 is the dimension of our phantoms. On the other hand calculating the gradient has similar complexity than applying the neural network which is in the order O(F 2 LN 2 ) with F the number of convolution channels and L the number of layers.

Numerical Results
For the numerical results we train two regularizers R (1) and R (3) as described in Section 3.2. The networks are implemented using PyTorch [36]. We also use PyTorch in order to calculate the gradient ∇ x R (m) . We take N iter = 15, s = 0.25 and x 0 = Φ (m) A y in Algorithm 1 and compute the inverse (A A − s Id) −1 only once and then use it for all examples. We set α = 0.015 for the noise-free case, α = 0.016 for the low noise case and α = 0.02 for the high noise cases, respectively, and selected a fixed β = 15. We expect that the NETT functional will yield better results due to data consistency, which is mainly helpful outside the masked center diagonal.
First we use the phantom from the testdata shown in Figure 1. The results using post processing and NETT are shown in Figure 2. One sees that all results with higher noise than used during training are not very good. This indicates that one should use similar noise as in the later applications even for the NETT. Figure 3 shows the average error using 10 test phantoms similar to the on in Figure 1. Careful numerical comparison of the numerical convergence rates and the theoretical results of Theorem 1 is an interesting aspect of further research. To investigate the stability of our method with respect to phantoms that are different from the training data we create a phantom with different structures as seen in Figure 4. As expected, the post processing network Φ (3) is not really able to reconstruct the circles object, since it is quite different from the training data, but it also does not break down completely. On the other hand, the NETT approach yields good results due to data consistency.    (1) . Middle row: NETT reconstructions using R (1) . Bottom row: NETT reconstructions using R (3) . From Left to Right: Reconstructions from data without noise, low noise (σ = 0.01) and high noise (σ = 0.1).

Conclusions
We have analyzed the convergence a discretized NETT approach and derived the convergence rates under certain assumptions on the approximation quality of the involved operators. We performed numerical experiments using a limited data problem for PAT that is the combination of an inverse problem for the wave equation and an inpainting problem. To the best of our knowledge this is the first such problem studied with deep learning. The NETT approach yields better results that post processing for phantoms different from the training data. NETT still fails to recover some missing parts of the phantom in cases the data contains more noise than the training data. This highlights the relevance of using different regularizers for different noise levels. Finding ways to make the regularizers less dependent on the noise level used during training is a possible future research direction. Another interesting question is if this results can be combined with with approximation error estimates for neural networks e.g., [37,38]. It seems not obvious how these two approaches can be combined. Furthermore, studying how one can define neural network based regularizers that fulfill (12) might also be an interesting line of future research.