On Regularization via Frame Decompositions with Applications in Tomography

In this paper, we consider linear ill-posed problems in Hilbert spaces and their regularization via frame decompositions, which are generalizations of the singular-value decomposition. In particular, we prove convergence for a general class of continuous regularization methods and derive convergence rates under both a-priori and a-posteriori parameter choice rules. Furthermore, we apply our derived results to a standard tomography problem based on the Radon transform.


Introduction
In this paper, we consider linear inverse problems in the standard form Ax = y , (1.1) where A : X → Y is a bounded linear operator between real or complex Hilbert spaces X and Y . Additionally, we assume that A is compact, which implies both that solving (1.1) is an ill-posed problem, and that there exists a singular system (σ k , v k , u k ) ∞ k=1 such that A admits a singular-value decomposition (SVD) of the following form (cf. e.g., [10]): (1.2) 1. The sequence {σ 2 k } ∞ k=1 consists of the non-zero eigenvalues of A * A written down in decreasing order, taking multiplicity into account and observing σ k > 0.
3. The singular functions u k are defined by u k := (1/σ k )Av k . Hence, they satisfy AA * u k = σ 2 k u k and form a complete orthonormal system spanning R(A).
Note that due to the above definition, the singular functions u k and v k satisfy The SVD is an important tool for analysing and solving ill-posed problems. In particular, the minimum-norm least-squares solution x † of (1.1) is characterized by which is well-defined if and only if the so-called Picard condition holds: (1.6) Furthermore, given noisy data y δ , which are typically assumed to satisfy where δ denotes the noise level, one can define stable approximations x δ α of x † by where g α is a properly selected approximation of s → 1/s. If the regularization parameter α is suitably chosen, e.g., by an a-priori or an a-posteriori parameter choice rule [9], then one can prove that x δ α → x † as δ → 0. Furthermore, if, e.g., the source condition holds, then one can even prove (order-optimal) convergence rates of the form [10,20] x δ α − x † X = O δ 2µ 2µ+1 .
Most classic regularization methods can be identified with specific choices of g α , which are also called spectral filter functions. For example, for Tikhonov regularization, Landweber iteration, or the truncated singular-value decomposition (TSVD), there respectively holds [10,20] g α (s) = Tikhonov, s −1 1 − (1 − s) 1/α , Landweber, 1/s , s ≥ α , 0 , else , TSVD. (1.10) Unfortunately, for many operators no explicit representation of the SVD is known, and even if it is, its numerical computation might be infeasible (see, e.g., [22]). Additionally, generalizing an available SVD of an operator defined between Hilbert spaces over regular domains to irregular domains is often impossible, e.g. in case orthogonality is lost. The situation is somewhat better for finite-dimensional problems, where one can work with the analogously defined SVD of a matrix, but this can also become infeasible when the problem is medium-to large-scale. Hence, even though an important tool in the analysis of ill-posed problems, the SVD often only enjoys limited use in practice.
In order to remedy this situation, several researchers have studied generalizations of the SVD such as the Wavelet-Vaguelette Decomposition (WVD) [1,6,7,12,18,19] or the Frame Decomposition (FD) [8,13,15,16]. The idea is that by weakening some of the requirements of the SVD such as the eigenvalue properties and the resulting orthogonality of the functions u k and v k , one may end up with decompositions similar to (1.2) and (1.5) which are easier to derive explicitly for a given operator. For example, the WVD requires an orthonormal wavelet basis {ψ jk } j,k∈N on X and two bi-orthogonal sets (vaguelettes) {u jk } j,k∈N , {v jk } j,k∈N on Y , connected by the quasi-singular relations λ jk v jk = Aψ jk , and λ jk ψ jk = A * u jk , ∀ j, k ∈ N . (1.11) which are reminiscent of (1.4). The resulting decomposition of the operator A is similar to (1.2) and can be used for solving (1.1). For applications of the WVD to the practically relevant problems of computerized and photoacoustic tomography including different aspects of its numerical realization see e.g. [7,11,12]. Extensions of the WVD which also have applications to the two-and three-dimensional Radon transform are e.g. given by the biorthogonal curvelet and shearlet decompositions [2,4]. In contrast to the WVD, the FD does not specifically work with (orthogonal) wavelets and vaguelettes but with general frames in Hilbert spaces; cf. Section 2.1. In particular, the FD requires a frame {e k } k∈N over X and a frame {f k } k∈N over Y , connected via where λ k denotes the complex conjugate of the coefficient λ k ∈ C. The above condition is clearly connected to (1.4) and (1.11), and leads to the following decomposition: Here the functions {f k } k∈N denote the dual frame of the frame {f k } k∈N ; cf. Section 2.1. Furthermore, for y ∈ Y one can then analogously to (1.5) consider the operator where similarly to above {ẽ k } k∈N denotes the dual frame of the frame {e k } k∈N . The element Ay is well-defined if the following analog of the Picard condition (1.6) holds: The properties of the operator A and its use for obtaining (approximate) solutions of (1.1) was studied in detail in [16]. In particular, it was investigated in which cases Ay is either a minimum-coefficient or a minimum-norm (least-squares) solution of (1.1).
While we refer to [16] for details, we here only want to mention the special case that A satisfies a stability condition of the form for some constants c 1 , c 2 > 0 with Z ⊆ Y being a Hilbert space. In this case, for each y ∈ R(A) the unique solution of (1.1) is precisely given by Ay. Furthermore, in this case it is possible to give recipes for finding frames {e k } k∈N and {f k } k∈N which satisfy (1.12). For example, starting with a frame {f k } k∈N which satisfies the condition for a sequence of coefficients 0 = α k ∈ R and some constants a 1 , a 2 > 0, one can define and it then follows that {e k } k∈N forms a frame over X satisfying (1.12) with λ k = 1/α k [16]. A typical example for an operator satisfying condition (1.16) is given by the Radon transform [20,21], and (1.17) can e.g. be satisfied if Y and Z are (suitably connected) Sobolev spaces and {f k } k∈N is either an exponential or a wavelet frame/basis. Note also that condition (1.16) is satisfied for any continuously invertible operator A with Z = X, in which case also (1.17) holds for any frame {f k } k∈N over Y . Further details on all of these topics can be found in [16], which also includes a generalization of condition (1.12) that can be useful if the Hilbert spaces X and Y have a particular product structure [15,16,24]. In this paper, we focus on different aspects of regularization via FDs. In particular, analogously to (1.8) for the SVD we consider stable approximations of Ay of the form and derive convergence and convergence rate results under a-prior and a-posteriori parameter choice rules similar to those for the SVD [10]. We want to note that this work is inspired by our recent investigations of general frame decompositions in Hilbert spaces [16] and their specific application to the atmospheric tomography problem [15,24]. However, we emphasize that different aspects of regularization via WVDs and specific FDs have already been considered in the literature before [8,12,13]. For example, a WVD approach to photoacoustic tomography was developed in [12], where regularization with wavelet sparsity constraints is employed via a soft-thresholding approach. This was generalized to sparse regularization for inverse problems using FDs and nonlinear soft-thresholding in [13], where also convergence rates were proven under a-priori parameter choice rules. The analysis in [13] is based on the assumption that {f k } k∈N forms a frame over R(A) instead of over Y as in our current setting. The same is also assumed in the preprint [8], in addition to the requirement that {e k } k∈N forms a frame over N (A) ⊥ = R(A * ) and that λ k ∈ (0, ∞) for all k ∈ N. Note that under these assumptions, for all y ∈ D(A † ) there holds A † y = Ay. Working within this setting, the authors of [8] prove convergence and convergence rates under a-priori parameter choice rules for continuous regularization methods based on FDs similar to our results. As in this paper, the analysis is based on the standard approach to linear inverse problems [10].
In contrast to these papers, here we consider the general FD setup introduced above, i.e., we allow arbitrary λ k ∈ C and assume that {e k } k∈N and {f k } k∈N form frames over X and Y , respectively. This setting is beneficial, since it allows more general and potentially different frame decompositions of a given operator. Furthermore, it is often easier to find a frame over the whole spaces X and Y instead of over N (A) ⊥ and R(A), respectively. Of course, one may theoretically obtain such frames by projecting the frames {e k } k∈N and {f k } k∈N onto these subspaces. However, doing so is practically infeasible, since in most cases an explicit characterization of N (A) ⊥ and R(A), and thus of the orthogonal projectors onto these subspaces is unavailable or involves the operator A † . This situation is comparable to the SVD of a compact operator: the existence is guaranteed but explicit representations are often unavailable, which is a motivation for considering frames in the first place. Furthermore, while (1.12) implies that {f k } k∈N,λ k =0 forms a frame over R(A), the set {e k } k∈N,λ k =0 does not necessarily form a frame over N (A) ⊥ , unless further assumptions on the values λ k , the frame {f k } k∈N , and/or the mapping properties of the operator A are made. Note that by redefining the frame functions e k one could always achieve that λ k ∈ R + 0 , at the cost of potentially changing the solution properties of the element Ay; compare with the conditions in [16]. However, the fact that λ k may also be zero is a crucial difference to the setting of [8], since then (1.12) no longer implies that all frame functions e k are elements in N (A) ⊥ = R(A * ). Hence, every dual frame functionẽ k may also have some component outside of N (A) ⊥ , since by definition it depends on all frame functions e k including those which correspond to λ k = 0. This may translate to the element Ay, except in those situations characterized in [16] in which Ay = A † y ∈ N (A) ⊥ . The benefit of allowing λ k = 0 is that it can make the search for frames satisfying (1.12) easier, as is the case e.g. for the atmospheric tomography problem [15]. In contrast, under the more restrictive assumptions of [8] there always holds A † y = Ay, which is not necessarily the case in our more general setting (see above). Hence, while [8] considers the stable approximation of A † y in the presence of noisy data via its representation in terms of frames, here we consider the stable approximation of the (approximate) solution Ay. Finally, note that in contrast to [8] we also present convergence rates results for an a-posteriori parameter choice rule adapted from the discrepancy principle, as well as numerical results illustrating our derived theory on the example of a standard tomography problem based on the Radon transform.
The outline of this paper is as follows: In Section 2, after reviewing some necessary material on frames in Hilbert spaces, we consider continuous regularization methods based on frame decompositions and show convergence and convergence rates under standard assumptions, both under a-priori and a-posteriori parameter choice rules. In Section 3, we then apply our results to a standard tomography problem based on the Radon transform, providing numerical examples for specific FDs and comparing the results of different regularization methods. Section 4 then summarizes our results.

Background on Frames in Hilbert Spaces
Before deriving our results on regularization via FDs, we first recall some basic facts on frames in Hilbert spaces. This short summary, based on the seminal work [5], is adapted from our previous publications [15,16]. First, recall the definition of a frame. Definition 2.1. A sequence {e k } k∈N in a Hilbert space X is called a frame over X, if and only if there exist frame bounds 0 < B 1 , B 2 ∈ R such that for all x ∈ X there holds For a given frame {e k } k∈N one can consider the frame (analysis) operator F and its adjoint (synthesis) operator F * , which are given by Due to (2.1) there holds Furthermore, one can define the operator S := F * F , i.e., x, e k X e k , which is a bounded and continuously invertible linear operator with for all x ∈ X, and thus the set {ẽ k } k∈N also forms a frame over X with frame bounds For the corresponding operators it follows analogously to (2.3) that In particular, for any sequence of coefficients a = {a k } k∈N there holds Furthermore, it can be shown that there holds and thus any x ∈ X can be written in the form x, e k Xẽ k . (2.9) Note that for any frame {e k } k∈N the following statements are equivalent (see, e.g., [3]): no element can be deleted) . (2.10) In general there holds holds {0} ⊂ N (F * ) = N (F * ), an thus the decomposition of x given in (2.9) is not unique, which is a key differences between frames and bases. However, this decomposition can be understood as the most economical one (cf. [5]).

Continuous Regularization Methods
In this section, we consider general continuous regularization methods based on FDs similar to those based on the SVD [10]. More precisely, for α > 0 we consider the functions as approximations of Ay given noisy data y δ , as well as the functions in the noise-free case. Here, g α : C → C is a suitable approximation of s → 1/s to be specified below. We will prove convergence as well as convergence rates of the form for both a-priori and a-posterior parameter choice rules under standard assumptions. Throughout the analysis, which is based on classical arguments (see e.g. [10]), we use Assumption 2.1. The operator A : X → Y is bounded, linear, and compact between the (complex) Hilbert spaces X and Y . Furthermore, the set {e k } k∈N forms a frame over X with frame bounds B 1 , B 2 , and the set {f k } k∈N forms a frame over Y with frame bounds C 1 , C 2 . Moreover, there exist coefficients λ k ∈ C such that (1.12) holds. In addition, let g α : C → C be a parameter-dependent family of piecewise continuous, bounded functions defined for all α > 0, satisfying Additionally, assume there exists a constant C > 0 independent of α such that (2.14) Note that (2.13) and (2.14), as well as the following definitions, which we need for the upcoming analysis, are the same as those used in the classic SVD analysis [10,20]: as well as Next, concerning the well-definedness of z δ α and z α we have Lemma 2.1. Let Assumption 2.1 hold and let y, y δ ∈ Y . Then z δ α and z α as given in (2.11) and (2.12) are well-defined. Furthermore, if additionally the Picard condition (1.15) holds, then Ay given in (1.14) is well-defined.
Proof. Let y ∈ Y be arbitrary but fixed. Since the set {e k } k∈N forms a frame over X with frame bounds B 1 , B 2 , it follows from (2.8) that Hence, since we assumed that g α is bounded it follows that z α is well-defined. Analogously we can show that z δ α is well-defined. Furthermore, similarly to above we have Hence, if (1.15) is satisfied then also Ay is well-defined, which concludes the proof.
The following result establishes convergence of z α to Ay in the noise-free case. Ay − z α X = 0 .
Proof. Let y ∈ Y be arbitrary but fixed and let (1.15) hold. Then due to Lemma 2.1 both Ay and z α are well-defined. Furthermore, due to (1.14) and (2.12) there holds Since {e k } k∈N forms a frame with frame bounds B 1 , B 2 , it follows from (2.8) that Since due to (2.14) and the definition of r α there holds |r α (λ)| ≤ (C + 1), it follows that The right-hand side in the above inequality is uniformly bounded independently of α due to (1.15). Hence, we can apply the dominated convergence theorem to obtain Since due to (2.13) there holds lim α→0 r α (λ) = 0 for all λ = 0, it follows that lim α→0 Ay − z α X = 0 , which yields the assertion and thus concludes the proof.
Next, we derive an upper bound on the data-propagation error in the following Theorem 2.3. Let Assumption 2.1 hold, let y, y δ ∈ Y satisfy (1.7), and let z α , z δ α be as in (2.12), (2.11), respectively. Then with G α as defined in (2.16) there holds (1.7) be arbitrary but fixed. Then due to Lemma (2.1) both z α and z δ α are well-defined. By their definitions (2.12) and (2.11) it follows that Since {e k } k∈N forms a frame with frame bounds B 1 , B 2 , it follows from (2.8) that and thus together with the definition (2.16) of G α and (2.14) we get Hence, since {f k } k∈N forms a frame with frame bounds C 1 , C 2 , it follows that which after taking the square root on both sides yields the assertion.
Combining Theorem 2.2 and Theorem 2.3, for the total error we obtain where the first term tends to 0 for α → 0 if y satisfies (1.15), similarly as in the SVD case. Hence, if the regularization parameter α = α(δ) is chosen such that then we obtain convergence of z δ α(δ) to Ay as the noise level δ → 0.

A-priori Parameter Choice Rules
In this section we derive convergence rates results under a-priori parameter choice rules similar to those for the SVD case (cf. [10]). These typically require source conditions such as the Hölder-type condition (1.9) for some µ > 0. Alternatively, this can be rewritten as Using the SVD (1.3) of the operator A it follows and thus (1.9) is equivalent to For the upcoming analysis, we use a source condition similar to (2.21), namely which analogously to (2.22) implies the decay condition We now start our analysis by deriving a convergence rate estimate given exact data.

24)
for some constant c µ > 0. Then if the source condition (2.23) holds it follows that Proof. Let y ∈ Y be arbitrary but fixed and let (1.15) hold. Then due to Lemma 2.1 both Ay and z α are well-defined. Then due to (2.18) there holds which together with the source condition (2.23) implies Hence, since {e k } k∈N forms a frame with bounds B 1 , B 2 , it follows with (2.8) that Together with (2.24) this implies which yields the assertion.
Next, we derive convergence rate estimates also in the case of inexact data.
Remark 2.1. Note that (2.24) could be replaced by the more general condition for some function ω µ : (0, α 0 ) → R. With this, one can derive similar results as in Theorem 2.4 and Theorem 2.5 also under more general source conditions than (2.23).

A-posteriori Parameter Choice Rule
In this section we derive convergence rates results under an a-posteriori parameter choice rule similar to the SVD case. More precisely, we consider a variant of the well-known discrepancy principle, which defines the regularization parameter α DP (δ, y δ ) via where the constant τ DP is chosen such that Since from the properties of the SVD it follows (after some computation) that where Q is the orthogonal projector onto R(A), it follows that (2.27) is equivalent to This motivates the following a-posteriori parameter choice rule for the FD case: 3. Let γ be as in (2.17) and let the parameter τ > 0 be such that where as before C 2 denotes the upper frame bound of {f k } k∈N . Then we define In case that α(δ, y δ ) = +∞, z δ α(δ,y δ ) is understood in the sense of a limit, i.e., Remark 2.2. While the SVD is generally not known explicitly and thus mostly used as a theoretical tool, we here assume that the FD is known explicitly. Hence, the sum in (2.29) can be computed and thus our a-posteriori rule can also be used in practice.
Proof. Note first that due to (2.1) and (2.17), for all α > 0 there holds Hence, we can apply the dominated convergence theorem to obtain Hence, for each ε > 0 there exists an α(ε) such that and thus D defined in (2.31) is non-empty and consequently α(δ, y δ ) is a well-defined element in (0, ∞]. Next, note that since we assumed that for all λ ∈ C the functional α → g α (λ) is continuous from the left, it follows that the same is also true for the functional Hence, the supremum in (2.29) is attained and thus (2.32) holds whenever α(δ, y δ ) < ∞. For the limit case α(δ, y δ ) = ∞ observe that due to (2.8) and (2.16) there holds Hence, together with (2.1) and (2.33) we obtain and thus taking the limit we obtain z δ ∞ = lim α→∞ z δ α = 0, which concludes the proof.
We are now able to derive convergence rate estimates in the presence of noisy data.

34)
for some constant c µ+1/2 > 0. Furthermore, let (2.33) and the source condition (2.23) hold, and let the function α → g α (λ) be continuous from the left for all λ ∈ C. Then for α = α(δ, y δ ) chosen via the a-posteriori stopping rule (2.28), (2.29), it follows that Proof. First, assume that there are sequences δ n → 0 and y δn ∈ Y satisfying such that for all n there holds α n := α(δ n , y δn ) = ∞. Then due to (2.29) there holds for all α > 0. Hence, together with (2.1) and the definition (2.17) of γ we obtain Letting n → ∞ in the above inequality we thus obtain that for all α > 0 there holds Hence, using the dominated convergence theorem we obtain Now since from the definition (2.16) of G α and (2.33) there follows Since this implies y, f k Y = 0 for all λ k = 0, it follows from the definition of Ay that Ay = 0 = z δn αn . Hence, for the remainder of this proof we can assume that α(δ, y δ ) < ∞ for all y δ satisfying (1.7) with δ sufficiently small. Next, note that in (2.19) we have shown that Since by the Hölder inequality with p = 2µ + 1 and q = (2µ + 1)/(2µ) there follows we thus obtain that .

(2.36)
We now consider each of these sums separately, where we have already shown in (2.35) that the second sum in (2.36) can be estimated by Concerning the first sum in (2.36), note that due to the source condition (2.23), Together with the definition of γ and since {e k } k∈N forms a frame, we obtain Hence, inserting (2.35) and (2.37) into (2.36) we obtain which simplifies to and thus provides an (optimal) convergence rate given exact data. Next, note that due to the reverse triangle inequality there holds By the definition (2.29) of our a-posteriori parameter choice rule we obtain Thus, rearranging (2.39), together with the definition of γ, (1.7), and (2.3) we obtain Now since due to (2.28) there holds Since due to (2.34) and the definition of γ there holds , it follows together with the source condition (2.23) that Hence, together with (2.1) we obtain (cδ) 2 ≤ c 2 µ+1/2 (2α) 2µ+1 B 2 w 2 X , and thus we find that Now, note that due to (2.20) and (2.33) there holds (2.41) Inserting (2.38) and (2.40) into this inequality yields , which now yields the assertion.

Application to the Radon Transform
In this section, we illustrate our theoretical results on continuous regularization via FDs by applying them to a standard tomography problem based on the Radon transform [20,21], which in 2D is given by where ω(ϕ) = (cos(ϕ), sin(ϕ)) T for ϕ ∈ [0, 2π) and s ∈ R. After recalling some recently derived FDs of the Radon transform on which we base the numerical illustration of our theoretical results derived above, we consider the details of the implementation of these decompositions and provide a number of numerical examples. For comparison, we also present numerical results based on the SVD of the Radon transform; cf. [20,21].

Frame Decompositions of the Radon Transform
A number of different decompositions of the Radon transform fitting into the category of frame decompositions have been studied in the past. These include e.g. the WVD [7] or the the biorthogonal curvelet/shearlet decompositions [2,4], for which also efficient implementations are available. For the numerical examples presented in this paper, which mainly serve to illustrate the theoretical results derived above, we focus on a class of FDs recently introduced in [16]. These are based on a general strategy for deriving FDs for arbitrary bounded linear operators in Hilbert spaces satisfying a stability condition of the form (1.16). In contrast to the classic SVD, available for the setting where Ω D := {x ∈ R 2 | |x| ≤ 1} and Ω S := R × [0, 2π), these FDs of the Radon transform are also available for the very general settings where Ω S := [−1, 1] × [0, 2π). While we refer to [16] for the most general version of these FDs, here we focus only on two special cases based on wavelets and exponentials.  1). Furthermore, let {ψ j,k } j,k∈Z be an orthonormal wavelet basis corresponding to an r-regular multiresolution analysis of L 2 (R) with r > β, let {w l } l∈N be an orthonormal basis of L 2 (0, 2π), and define f j,k,l (s, ϕ) := ψ j,k (s)w l (ϕ) , and e j,k,l := 1 + 2 −2j(β+1/2) 1/2 A * f j,k,l . x, A * f j,k,l f j,k,l .
Furthermore, for any y ∈ R(A) the unique solution of Ax = y is given by A possible choice for the orthonormal basis {w l } l∈N in the above theorem is e.g., which was also used throughout the numerical experiments presented below.
(3.5) Then the set {v j,k } j,k∈Z forms a frame over H β 0 (Ω D ) and Furthermore, for any y ∈ R(A) the unique solution of Ax = y is given by Ay = j,k∈Z 1 + |j| 2 (β+1/2)/2 y, w j,k L 2 (Ω S )ṽ j,k . (3.6) Note that the FDs and the corresponding formulas for computing Ay given in Theorem 3.1 and 3.2 can be efficiently implemented by approximating the involved inner products via fast Fourier and wavelet transforms. While this can entail a loss of accuracy compared to higher order integration methods, the increased efficiency is beneficial for applications with real-time requirements such as atmospheric tomography [16].

Implementation and Computational Aspects
In this section, we consider the implementation of the FDs of the Radon transform presented in Section 3.1, focusing on the relevant special case β = 0 in (3.2), i.e., The key difference in implementation between the FDs given in Theorem 3.1 and Theorem 3.2 lies in the different sets of frame functions and their corresponding properties. Note that if the dual frame functionsẽ j,k,l orṽ j,k are pre-computed and stored, then for each right-hand side y the computation of Ay amounts only to the computation of either the inner products y, f j,k,l L 2 (Ω S ) or y, w j,k L 2 (Ω S ) , and a corresponding summation according to (3.4) or (3.6), respectively. As mentioned above, these inner products can be efficiently implemented using fast Fourier and wavelet transforms.

Problem Discretization and Computational Environment
For the discretization of the problem we have used the AIR Tools II toolbox [14], which is based on a piecewise constant discretization of both the definition and the image space of A. More precisely, a density function x ∈ L 2 (Ω D ) is approximated by a piecewise constant function with values given on a uniform N × N pixel grid. Similarly, sinogram data y are also considered as piecewise constant functions on a uniform p×N θ pixel grid, where p denotes the number of equidistant, parallel lines on [−1, 1], and N θ denotes the number of different angles θ n . Hence, the discretized problem can be written as For all tests presented below, we have used the choice N = p = 60, and N θ = 180, with uniformly spaced angles θ n = nπ/180 for n = 0 , . . . , 179, which due to Ax(s, ϕ) = Ax(−s, ϕ + π) amounts to a full-angle tomography problem. The infinite sums in formulas (3.4) and (3.6) for computing Ay have been replaced by finite sums as detailed below, and the involved integrals have been approximated using the trapezoidal rule. All computations have been performed using Matlab 2019b on a desktop computer running on Windows 10 with a 4 core processor (Intel Core i5-6500 CPU@3.20GHz) and 16GB RAM, except the computation of the dual frames, which is discussed in detail below, and which was performed on the high performance computing cluster Radon1 [17].

Computation of the Dual Frame Functions
Next, we consider the computation of the dual frame functionsẽ k . While it is generally not possible to give explicit expressions, one can use the recursive approximation [5] e k ≈ẽ M k = 2 where B 1 , B 2 denote the frame bound of {e k } k∈N and S is as in (2.4). The approximation error can be estimated by [5]: Alternatively, one can compute the dual frame functionsẽ k directly via their definitionẽ k = S −1 e k . Discretizing as above, one has to solve the system of equations whereẽ k and e k are piecewise constant discretizations ofẽ k and e k , respectively, and where {φ i } i∈N 2 denotes the piecewise constant pixel basis on Ω D used for discretization.  (2.9) in terms of the frame {v j,k } j,k∈Z , with the functions v j,k as defined in (3.5) using the explicitly (middle) and the recursively (right) computed dual frames. Note that these computations are based on a coarser discretization with N θ = 90 instead of N θ = 180 as explained in the text.
Using an LU-decomposition for the matrix S, the equation Sẽ k = e k can be solved efficiently for all k. However, our matrices S are ill-conditioned, with condition-numbers κ wav = λ max λ min ≈ 7.7 · 10 5 , and where λ max and λ min denote the maximal and minimal singular value of S, respectively. Hence, to obtain stable approximations ofẽ k we used Tikhonov regularization, i.e., We found optimal results for α = 0.01 for the wavelet-based FD and α = 2 for the exponential-based FD. For these particular choices of α, this explicit approach outperformed the recursive approximation. The computed dual frame functions were verified by implementing the reconstruction formula (2.9), an example of which is shown in Figure 3.1. Note that for obtaining these results, a coarser discretization of N θ = 90 had to be used in the numerical computation of the frame functions v j,k defined via (3.5), since the computation of the recursively approximated dual framesṽ j,k already takes about 8 hours in this setup, while the full setup with N θ = 180 angles is computationally infeasible. The reconstructions shown in Figure 3.1 have a relative error of 6.98% in the case of explicitly computed dual frames, and an error of 9.07% in the case of recursively approximated dual frames. Thus, not only is the explicit computation approach faster than the recursive approach, but it also outperforms it in terms of reconstruction quality. Hence, all numerical results presented below are based on the explicitly computed dual frame functions using the full setup with N θ = 180 angels. We now consider the setup from Theorem 3.1. In particular, for {ψ j,k } j,k∈Z we consider inhomogeneous orthonormal wavelet bases of the form {ψ 1 jmax,k } k∈Z ∩{ψ 2 j,k } j≤jmax,k∈Z . These are defined via ψ 1 j,k (s) = 2 −j/2 ψ 1 (2 −j s − k) and ψ 2 j,k (s) = 2 −j/2 ψ 2 (2 −j s − k) based on suitable scaling-and wavelet function ψ 1 and ψ 2 , respectively (cf. [5] for details). A typical example of such wavelet bases are the different Daubechies wavelets [5], for which the corresponding scaling-and wavelet functions are available in Matlab. Using such an inhomogeneous wavelet basis, any x ∈ L 2 (R) can be written in the form [5]:
Remark 3.1. Note that since due to (1.12) for all x ∈ L 2 (Ω D ) there holds α j x, e j,k,l L 2 (Ω D ) = Ax, f j,k,l L 2 (Ω S ) , in the the computation of the dual frame functionsẽ j,k,l the computation of the inner products x, e j,k,l L 2 (Ω D ) can be replaced by the computation of the inner products Ax, f j,k,l L 2 (Ω S ) . This is advantageous, since the support of f j,k,l is known explicitly (cf. Figure 3.2) and thus the integration domain can be restricted.

Implementation of the Exponential-based Frame Decomposition
For the setup as in Theorem 3.2 we replace the infinite sum in (3.6) by using j min = − p/2 = −30, j max = p/2 − 1 = 29, and α j = 1 + |j| 2 −1/4 , which follows from our choice of β = 0. Analogously, for the approximate solutions z δ α we use with the filter functions of Tikhonov regularization and Landweber iteration as above.

Numerical Results
In the following, we present some numerical results for the FDs of the Radon transform discussed above. The numerical experiments involve two regularization methods combined with an a-priori and an a-posteriori parameter choice rule for different noise levels. We also present equivalent results based on the SVD of the Radon transform, in order to compare the performance and behaviour of these methods. Note that since the SVD is a special case of the FD, the parameter choice rules discussed in Section 2 can be used for both. The truncation value for the indices in the SVD-implementation was chosen such that the number of singular functions coincides with the number of frame functions in the exponential-based setup.
For our numerical experiments we use the Shepp-Logan phantom x SL depicted in Figure 3.1 (left) as the exact solution, i.e., x † := x SL . For measuring the quality of the obtained reconstructions, we use both the relative L 2 -error z δ α − x SL / x SL and the structural similarity index measure (SSIM), which is defined in [23]. The SSIM is a value in [0, 1], with higher values indicating a stronger structural similarity.   Table 3.1, which also includes results for our a-posteriori parameter choice rule (2.29). However, since these results are visually very similar to those obtained with the a-priori parameter choice rule, we decided not to include them  in this paper. Note that since in our experiments the frame functions f k are orthonormal, there holds C 2 = 1, and thus (2.28) provides a computable lower bound for τ in our a-priori parameter choice rule (2.29). However, we empirically found the choice τ = 0.1 for the SVD case and τ = 20 for the FD cases to lead to much better results,  and thus used them in all of the presented numerical experiments.
Comparing the obtained results as summarized in Table 3.1, we see that for a relative noise level of 1%, additional regularization beyond the truncation inherent in the discretization is only beneficial for the SVD case. However, for a noise level of 15% noise, additional regularization via the Tikhonov and Landweber filter functions is often beneficial, which can be seen from the error measures, in particular from the SSIM. It appears that the additional regularization has the most impact on the SVD, while having less impact on the exponential-and the wavelet-based FDs. This can be explained by the fact that for the chosen range of the indices j we have α j ∈ [0.71, 0.94] for the wavelet-based FD and α j ∈ [0.15, 1] for the exponential-based FD. In contrast, the singular values σ m of the SVD lie within the interval [0.27, 3.54] in our chosen range of indices m. However, note that since these index ranges were chosen such that in both the SVD and the exponential-based FD case the same number of singular/frame functions are used, we find that the FDs lead to more accurate reconstructions in the case of low noise levels than those obtained via the SVD, while the SVD shows better stability and regularization properties in the case of high noise levels, at a comparable computational cost.
Finally, Figure 3.6 illustrates the dependence of the error-measures on the regularization parameter. The marks indicate the parameters selected by the a-priori and a-posteriori parameter choice rules. In particular, note that the optimal (maximal) value of the SSIM is reached at a larger value of α than the optimal (minimal) L 2 -error.

Conclusion
In this paper, we considered general continuous regularization methods based on FDs for linear ill-posed problems in Hilbert spaces. In particular, we proved convergence and convergence rates results under a-priori and a-posteriori parameter choice rules analogous to those for SVD-based regularization methods. Furthermore, we applied our results to a standard tomography problem based on the Radon transform, using specific FDs based on wavelets and exponential functions. The obtained results demonstrate that FDs are a viable approach for efficiently solving linear ill-posed problems.

Support
SH and RR were funded by the Austrian Science Fund (FWF): F6805-N36. LW was supported by the strategic program "Innovatives OÖ 2010 plus" by the Upper Austrian Government and by the Austrian Science Fund (FWF): W1214-N15, project DK8