Error Estimates for Arnoldo-Tikhonov Regularization for Ill-Posed Operator Equations

Most of the literature on the solution of linear ill-posed operator equations, or their discretization, focuses only on the infinite-dimensional setting or only on the solution of the algebraic linear system of equations obtained by discretization. This paper discusses the influence of the discretization error on the computed solution. We consider the situation when the discretization used yields an algebraic linear system of equations with a large matrix. An approximate solution of this system is computed by first determining a reduced system of fairly small size by carrying out a few steps of the Arnoldi process. Tikhonov regularization is applied to the reduced problem and the regularization parameter is determined by the discrepancy principle. Errors incurred in each step of the solution process are discussed. Computed examples illustrate the error bounds derived.


Introduction. Let
A : X → Y be an injective linear operator between the Hilbert spaces X and Y, and assume that A is not continuously invertible. We are concerned with the solution of operator equations of the form Ax = y, x ∈ X , y ∈ Y. (1.1) Let · X and · Y denote the norms of the spaces X and Y, respectively. We will assume that equation (1.1) is consistent and are interested in determining the solution of minimal norm. We denote this solution by x. The solution x might not depend continuously on y. Therefore its computation is an ill-posed problem. The right-hand side y of (1.1) is assumed not to be available; only an errorcontaminated approximation y δ ∈ Y of y is known. We assume that y δ satisfies with a known bound δ > 0. The solution of the equation obtained by replacing y by y δ in (1.1), if it exists, generally, is not a meaningful approximation of the desired solution x since A is not continuously invertible. In fact, equation (1.3) might not have a solution even when equation (1.1) does. A regularization method, which replaces the operator A by a nearby operator, such that the solution of the modified equation so obtained exists and is less sensitive to the error in y δ , has to be used to obtain a meaningful approximation of x. The numerical solution of (1.3) requires discretization at a certain stage of the process. In general, this can be done in two ways: (i) Regularize then discretize: In this approach, the infinite-dimensional ill-posed problem is transformed into a well-posed problem, e.g., by means of Tikhonov regularization. Then the well-known error estimates for the regularized solution, see, e.g., [7] for regularization in Hilbert spaces, can be applied. In a second step the now well-posed regularized equation is discretized, and available error estimators for well posed-problems, e.g., from the theory of finite elements, can be used. This approach has been followed in, e.g., [20,16,10,4].
(ii) Discretize then regularize: The discretization of the ill-posed operator equation (1.3) yields a linear system of algebraic equations A n x n = y δ n (1.4) with an ill-conditioned, possibly singular, matrix A n ∈ R n,n , and vectors x n , y δ n ∈ R n . Well-known methods from linear algebra are used for its solution; see, e.g., [3,11,13]. The difficulty with this approach is to obtain convergence and convergence rate results for the distance between the solution of the infinite-dimensional problem (1.1) and its discretized counterpart x n fulfilling (1.4), see, e.g., [18,5]. As mentioned above, the approach (i) works particularly well for different variants of Tikhonov regularization. Iterative methods, however, require frequent application of the operator, and maybe of its adjoint, which is only possible in a discretized form. Iterative methods therefore belong to category (ii). An analysis of approach (ii) has been carried out for an adaptive version of Landweber iteration [24], but to the best of our knowledge this approach has not been investigated for Krylov subspace methods. Additionally, methods that work exceptionally well in finite dimensions but have no infinite-dimensional counterpart, or for which an error analysis is missing in infinite dimensions, belong to category (ii). The latter holds for the method discussed in the present paper.
In this paper, we start with two continuous linear operator equations, (1.1) and (1.3), and discretize the latter to obtain the linear system of algebraic equations (1.4). We are concerned with the situation when the matrix A n is large and, in particular, when A n is too large to make the computation of its singular value decomposition attractive. Then we apply the Arnoldi process to compute an approximation of fairly low rank of the matrix A n in (1.4). We replace A n in (1.4) by this low-rank approximation, and compute an approximate solution of the linear system of equations with the aid of Tikhonov regularization. The replacement of A n by a low-rank approximation reduces the computational effort required for Tikhonov regularization when the matrix A n is large, which is the situation of interest to us. Our approach allows us to solve problems with a matrix A n that is too large to make the use of direct solution methods, which require factorization of a large matrix, e.g., of A n or a related matrix, too expensive to be attractive or feasible. We will discuss the effect on the computed solution of discretization errors that stem from replacing the operator A by the matrix A n , as well as the effect of the error in the right-hand side y δ . Moreover, we are concerned with the influence on the computed solution of the replacement of the matrix A n in the linear system (1.4) by a low-rank matrix determined by the Arnoldi process. We remark that Tikhonov regularization based on partial Arnoldi decomposition, and some variations thereof, have been described in [3,6,8,14,15] and in references therein. The contribution of this paper is to provide an error analysis. This paper is organized as follows. Section 2 discusses results by Natterer [18] on the discretization of integral operators. Discretization yields the linear system of algebraic equations (1.4). We assume that the matrix A n determined by discretization is so large that factorization is unattractive or unfeasible. Section 3 reviews the Arnoldi process for computing an approximation of fairly low rank of the matrix A n in (1.4). We use this low-rank approximation in Tikhonov regularization and obtain a quite efficient solution method. To analyze the performance of this solution approach, we have to take into account the discretization error as well as the error incurred by approximating the matrix A n by the low-rank approximation furnished by the Arnoldi process. Section 4 applies bounds due to Neubauer [20] to the computed solution obtained by the Tikhonov regularized problem that uses the approximation of the matrix A n in (1.4) computed with the Arnoldi process. We remark that while the bounds provided by Natterer [18] shed light on the influence of the discretization error on the computed solution, they are not useful for assessing the effect of approximating the matrix A n by a low-rank approximation determined by the Arnoldi process. We will comment further on this issue in Remark 4.10 of Section 4. A few computed examples that illustrate the theory are presented in Section 5, and concluding remarks can be found in Section 6.

Discretization of the operator equation.
To be able to numerically compute an approximate solution of equation (1.3) in the infinite-dimensional Hilbert space X , the equation first has to be discretized. This results in the finite-dimensional equation (1.4). We introduce a discretization and define a finite-dimensional leastsquares problem similarly as Natterer [18], who investigated regularization properties of projection methods.
Introduce the finite-dimensional subspaces and define projectors P n : X → X n and Q n : Y → Y n . The space X n is chosen for its convenience to use in applications and for the approximation properties of its elements. For instance, X n may be a space of piece-wise polynomials or finite elements. Consider the linear system of equations Q n AP n x = Q n y δ . (2.1) We identify the matrix A n and vector y δ n in (1.4) with the finite-dimensional operator Q n AP n and the right-hand side Q n y δ in (2.1), respectively. The unique least-squares solution of minimal norm of equation (2.1) is given by x n := A † n y δ n , where A † n denotes the Moore-Penrose pseudoinverse of the matrix A n . We identify this solution of (2.1) with the solution in R n of (1.4).
Let {e j } n j=1 form a convenient basis for X n , such as a basis of piece-wise polynomials or finite elements with local support. Consider the representation of an element x n ∈ X n . We identify the function x n with the vector To shed light on how x n X relates to x n 2 , we introduce an orthonormal basis { e j } n j=1 for X n . There is a nonsingular matrix M n = [m ij ] ∈ R n×n such that (e 1 , e 2 , . . . , e n ) = ( e 1 , e 2 , . . . , e n )M n , i.e., e j = n i=1 m i,j e i for j = 1, 2, . . . , n. For instance, when the basis { e j } n j=1 is determined from {e j } n j=1 by the Gram-Schmidt process, the matrix M n is upper triangular.
It follows that where σ max (M n ) denotes the largest singular value of the matrix M n . Let σ min (M n ) stand for the smallest singular value of M n . Then we obtain analogously to (2.3) that We will assume that there are constants c min and c max (independent of n) such that Then Thus, the norms · X and · 2 are equivalent. We will therefore simply write x n as x n . The equivalence will be explicitly used in Section 4.
The solution x n ∈ X n of (2.1) might not be a useful approximation of the desired solution x of (1.1) due to a large propagated error stemming from the error in the available data y δ n . We therefore would like to determine a bound for x − x n X . This is generally not possible without some additional assumptions on the solution x of (1.1); in particular, it is not sufficient that A and A n be close.
Let Ω ⊆ R N , X = L 2 (Ω), and define the Sobolev spaces H l = H l (Ω) for l ∈ R. Assume that holds for all x ∈ H −l and some 0 < l < ∞, i.e., the operator A : H −l → Y is continuously invertible. The theory developed by Natterer [18] requires that (2.5) holds for a finite value of l.
Let ω ⊥ be a unit vector perpendicular to ω, and define the Radon transform Then (2.5) holds with l = 1/2; see [19]. ✷ Example 2.2. Let for some kernel function k ∈ L 2 (R d ). If the Fourier transformk of k satisfies |k(ξ)| ∼ (1 + |ξ| 2 ) −β/2 , then (2.5) holds with l = β, see, e.g., [9]. ✷ Assume that the minimal norm solution x of (1.1) lives in H k . Natterer [18] shows that if the operator A is injective and the subspaces X n , n = 1, 2, . . . , are chosen so that an inverse estimate is fulfilled (see [18, eq. (4.1)-(4.5)] for details on the latter), then one obtains the bound for some constant C that can be chosen independently of h(n), x, and δ. Here, h = h(n) > 0 is a discretization parameter that depends on the approximation property of the subspaces X n , n = 1, 2, . . . , i.e., on how well x can be approximated by an element in X n ; in particular, h ց 0 as n → ∞. The parameter δ > 0 in (2.6) is the bound (1.2); see [18]. An optimal dimension of the discretized problem is given by and yields the bound for some constant C ′ > 0; see Natterer [18] for details. For instance, spline and finite element approximation spaces X n allow for bounds of the type (2.6) and (2.8).
Natterer [18] proposes that the dimension n of the solution subspace of the discretized problem (1.4) be chosen according to (2.7). This choice provides regularization of the operator equation (1.3) and no further regularization is necessary. We note that the use of wavelet-based projection methods also has been investigated for the solution of ill-posed problems. Regularization properties of wavelet methods have been shown by Dicken and Maaß [5].
Convergence rates analogous to (2.8), when h is chosen according to (2.7), also can be established in a different setting; see Mathé and Pereverzev [17]. They assume that the operator A is continuously invertible in Hilbert scales (which resembles the condition (2.5)), and show convergence rates in a stochastic noise setting with respect to norms of the relevant Hilbert scales; see [17,Theorem 6.3].
We conclude this section with a comment on condition (1.2). Let y n = Q n y and y δ n = Q n y δ . We will assume that y − y δ Y ≈ y n − y δ n Y . Then (1.2) translates to It is convenient to replace the norm in (2.9) by the Euclidean norm. This can be achieved analogously as in the beginning of this section: Let {f j } n j=1 form a convenient basis for Y n , such as f j = Ae j , where {e j } n j=1 is a basis for X n . Represent the element y n ∈ Y n as and define the vector y n = (y We would like to bound y n Y in terms of y n 2 . Introduce an orthonormal basis { f j } n j=1 for Y n . There is a nonsingular matrix N n ∈ R n×n such that It follows similarly as (2.3) that where σ max (N n ) denotes the largest singular value of N n . We will assume that there is an upper bound d max , independent of n, such that In computations, we will apply the discrepancy principle based on the inequality 3. Arnoldi decomposition of a matrix. Let A n and y δ n be as in (1.4), and assume that n is large. The Arnoldi process is a popular approach to reduce a large matrix to a small one by evaluating matrix-vector products with the large matrix and applying Gram-Schmidt orthogonalization. The small matrix, denoted by H ℓ+1,ℓ below, is an orthogonal projection of A n . Application of 1 ≤ ℓ ≪ n steps of the Arnoldi process to the matrix A n with initial vector y δ n gives the decomposition where the columns of the matrix V n,ℓ+1 ∈ R n,ℓ+1 form an orthonormal basis for the Krylov subspace K ℓ+1 (A n , y δ n ) := span{y δ n , A n y δ n , . . . , A ℓ n y δ n } with respect to the inner product and associated norm see, e.g., [26]. We also will denote the spectral norm of a matrix by · 2 . The matrix V n,ℓ in (3.1) is made up of the first ℓ columns of V n,ℓ+1 , and H ℓ+1,ℓ ∈ R ℓ+1,ℓ is an upper Hessenberg matrix, i.e., all entries below the subdiagonal vanish. We assume the generic situation that the subspace K ℓ+1 (A n , y n ) is of dimension ℓ+1 for all ℓ ≥ 0, otherwise the computations simplify; see below.
Algorithm 1 The Arnoldi process 1: Input: An ∈ R n,n , y δ n ∈ R n \{0}, number of steps ℓ. hj+1,j := w 2; vj+1 := w/hj+1,j Algorithm 1 describes the Arnoldi process for computing the decomposition (3.1). The algorithm is said to break down in iteration j if h k+1,k > 0 for 1 ≤ k < j, and h j+1,j = 0 in line 9. Then the decomposition (3.1) simplifies to A n V n,j = V n,j H j,j and the solution of (1.4) lives in the Krylov subspace K j (A n , y δ n ) if the matrix H j,j is nonsingular. This is secured, e.g., if the matrix A n is nonsingular. A discussion on how to continue the Arnoldi process in case of breakdown when H j,j is singular is provided in [25].
We remark that the Arnoldi process simplifies to the Lanczos process when the matrix A n is symmetric; see [26,Chapter 6].
4. The Arnoldi-Tikhonov method. The results of Section 2 suggest that the discretized system can be solved without further regularization if the discretization is carried out on a suitably (but not too) fine grid. However, numerical realization of regularization by discretization only often leads to difficulties, because an appropriate value of the dimension n of the solution subspace generally is not known before the computations are begun. For instance, when X n is a finite element space, it may be necessary to determine several discretizations and associated solutions (for different values of n) to find a suitable n-value. We therefore prefer to first discretize the spaces X and Y to obtain n-dimensional subspaces X n and Y n , respectively, that allow approximation of elements of X and Y with sufficient accuracy, and then regularize (1.4) by Tikhonov's method. In the remainder of this section, we identify the spaces X n and Y n with R n and the finite-dimensional operator Q n AP n in (2.1) with the matrix A n ∈ R n,n in (1.4).
The solution method considered consists of three steps: 1. Discretization of the (infinite-dimensional) operator equation. This requires an estimate of the distance between the solution of the infinitedimensional system and the solution of its finite-dimensional approximation.
2. Definition of a regularized finite-dimensional system. Estimate the distance between the solution of the finite-dimensional system and its regularized version.

3.
Compute an approximate solution of the regularized solution. Estimate the distance between the solution of the regularized finite-dimensional system and its computed approximation.
The error of the computed solution is bounded by the sum of the norms of these three errors. We will discuss each one of these errors separately.
Let the Arnoldi decomposition (3.1) be available, and introduce the approximation In what follows, we need to compute an estimate for the distance between A n and A n . It has recently been shown that a few of the largest singular values of a large matrix that stems from the discretization of a linear ill-posed problem can be computed quite inexpensively; see [22] for discussions and illustrations. Alternatively, we may use the easily computable Frobenius norm, where A n = [a ij ] n i,j=1 , and apply the bound

Assume that
for some scalar h ℓ > 0 and define the Tikhonov functional where α > 0 is a regularization parameter. We will solve the minimization problem x δ α,n,ℓ := arg min xn∈R n {J α,n,ℓ (x n )} .  The proper choice of the parameter pair {ℓ, α} has been studied by Neubauer [20], who considers the computation of an approximate solution of an operator equation where X and Y are Hilbert spaces, by first discretizing and then solving the discretized equation using Tikhonov regularization, Here T h,ℓ denotes a discretization and modification of T (see below), and T * h,ℓ is the adjoint of T h,ℓ . Neubauer [20] requires the operator T h,ℓ to satisfy where R ℓ is an orthogonal projector onto an ℓ-dimensional subspace W ℓ ⊂ Y to be specified below. The dimension ℓ is finite and typically quite small. Moreover, T h is a discretization of T and T h,ℓ is a modification of T h determined by R ℓ .
In our application of the results of Neubauer [20], we let T := A n and X := Y := R n . Thus, we set denotes the orthogonal projector onto the (closure of the) range of n . The operator T h is not important to us; we only use T h,ℓ . We are in a position to show the following results.
Proposition 4.1. Assume that the Arnoldi process does not break down. With the operators defined as above, we have where the bound (4.9) is inspired by (2.10).
Proof. First note that the ranges of the operators (matrices) A n and A (ℓ) n are closed as they are maps between finite-dimensional spaces. It follows from (3.1) and (4.1) that and, therefore, i.e., property (4.7) holds. Furthermore, This establishes (4.8). Finally, we have It remains to show (4.10). According to (4.7), we have R(A (ℓ) n ) ⊂ R(A n ). We will show that for every y n ∈ R(A n ), there exists an ℓ ≥ 1 such that y n ∈ R(A (ℓ) n ). Let x n ∈ R n and define y n = A n x n . Note that V n,ℓ V * n,ℓ is an orthogonal projector onto the space R(V n,ℓ ). Assuming that the Arnoldi process does not break down, there is an ℓ ≥ 1 (in the worst case, ℓ = n) such that x n ∈ R(V n,ℓ ) and, therefore, V n,ℓ V * n,ℓ x n = x n . It follows from (4.11) that i.e., y n ∈ R(A (ℓ) n ) and, consequently, R ℓ y n = y n . This shows the point-wise convergence of the projector R ℓ to I as ℓ increases.
Thus, the requirements of Neubauer [20, Assumption 2.3] are fulfilled, and we get the following result from [20, Proposition 2.6 and Theorem 3.1]: Proposition 4.2. Let x n be an approximate solution of (2.1) such that for some constant ρ ≥ 0 independent of n, and assume that y n − Q n y δ 2 ≤ δ. Let the regularization parameter α > 0 satisfy 14) where the inner product ·, · is defined by (3.2) and the constants C > 1 and E > 3 x n 2 are chosen such that Then the associated solution of (4.4) satisfies Additionally, O has to be replaced by O for ν = 1.

Remark 4.3.
The smoothness condition (4.12) is for infinite-dimensional problems a fairly strong restriction. In finite dimension, we observe that (4.12) implies that x n ∈ N (A n ) ⊥ . Therefore, there is a unique v n ∈ N (A n ) ⊥ such that x n = (A * n A n ) ν v n . However, the uniform boundedness of v n , cf. inequality (4.13), generally remains an open problem; see Proposition 4.5 below.
Remark 4.4. The quantities in (4.16) may depend on n. Generally, h ℓ does not vary much as ℓ is kept fixed and n is increased; see Section 5 for illustrations. When n is fixed and ℓ increases, h ℓ decreases. We are interested in choosing ℓ large enough so that both terms in the right-hand side of (4.16) are sufficiently small; see Corollary 4.6 below. Also, the condition (4.15) requires ℓ to be large enough.
Let us now give an example where the uniform boundedness of the source elements v n , required in Proposition 4.2, can be guaranteed: Proposition 4.5. Let the conditions of Proposition 4.2 except for condition (4.13) hold. Assume that A is self-adjoint, fulfilling (2.5), A n = P n AP n , and the solution x ∈ H k of the equation Ax = y fulfills a source condition with ν = 1/2 and source element v ∈ Hk. If A n is injective, then also the solutions of the equations A n x n = y n fulfill a source condition with ν = 1/2, and the associated source elements v n are uniformly bounded.
Proof. For ν = 1/2 and self-adjoint operator A, the source condition transfers to As A n is also self-adjoint, finite-dimensional, and injective, x n also fulfills a source condition, see Remark 4.3, x n = (A * n A n ) 1/2 v n = A n v n with a unique v n . As A fulfills (2.5), the distance between x and x n can be bounded by n , see (2.8), and δ n = y − y n ց 0 as n → ∞. Using again (2.8) for solving (4.17) with i.e., v n → v and consequently v n is uniformly bounded. The best convergence rates can be achieved for ν = 1: Corollary 4.6. Assume that the conditions of Proposition 4.2 hold, let ν = 1, and letα solve (4.14). Then for ℓ such that we have x δ α,n,ℓ − x n 2 = O δ 2/3 as δ ց 0.
Proof. The first term on the right-hand side of (4.16) behaves like O(δ 3/2 ) if h ℓ ∼ δ. For the second term, we have Since R ℓ → I as ℓ → n, we can choose ℓ large enough such that (I − R ℓ )A (ℓ) n 2 ≤ δ and obtain as δ ց 0.
With the same argument we achieve optimal convergence rates for each ν ∈ (0, 1) if p(l, ν) = O(δ 2ν/(2ν+1) ), which holds for l small enough. Now let us further specify the orthogonal projector R ℓ .

Moreover,
The singular value decomposition of H ℓ+1,ℓ yields Now using (4.19) and (4.20), we obtain Finally, when H ℓ+1,ℓ is of rank q ≤ ℓ, we have with I q being the q × q identity matrix.
The use of the discrepancy principle requires the solution of equation (4.14). The following result is concerned with the evaluation of the left-hand side of this equation.
Since V n,ℓ+1 V * n,ℓ+1 and I − V n,ℓ+1 V * n,ℓ+1 are complementary orthogonal projectors, it follows that Introduce the vector Then n * + αI 3 z δ n = R ℓ y δ n and, therefore, This shows (4.22). Equation (4.23) now follows by substituting (4.18) into (4.22). In actual computations, the matrix H ℓ+1,ℓ typically is small; see Section 4 for illustrations. The singular value decomposition of H ℓ+1,ℓ therefore is quite inexpensive to compute and the left-hand side of (4.22) easily can be evaluated. Corollary 4.9. Let the conditions in Section 2 hold and choose n according to (2.7). Assume that 1 ≤ ℓ ≤ n is large enough so that (4.14) has a solution, which we denote byα. Consider the regularized solution x δ α,n,ℓ , defined by (4.4) with α =α, an element in X n . Assume that the conditions of Corollary 4.6 hold. Then for a suitable constant C ′ > 0 with the parameter l the same as in (2.5).
Proof. Let x n ∈ X n be the minimal-norm solution (2.2) of (2.1) with n chosen according to (2.7). Then we obtain by the triangle inequality and (2.8) that Now considering x n and x δ α,n,ℓ elements in R n , we obtain from (2.4) that The inequality (4.24) now follows from Corollary 4.6. Remark 4.10. We conclude this section with a comment on why our analysis requires results by both Natterer [18] and Neubauer [20], because it may appear more natural to choose T h,ℓ = A (ℓ) n and apply Neubauer's result, only, without invoking those of Natterer. Our reason for using the bounds provided by Natterer is that in order to be able to use the results of Neubauer, without applying those of Natterer, we need a bound for where we consider A (ℓ) n an operator from X to Y and · denotes the appropriate operator norm. For many standard discretizations with suitable basis functions such a bound can be determined. However, this is not the case for the Arnoldi approximation A (ℓ) n , as the Arnoldi process depends on the starting vector. Therefore, we need a discrete approximation A n of A so that we are able to evaluate n is considered a matrix. The application of the Arnoldi process to A n gives an approximation of the solution of the discretized equation. Natterer's bounds are required to bound the distance to the solution of the infinite-dimensional problem.

Computed examples.
We apply the Arnoldi-Tikhonov method to a few ill-posed operator equations and illustrate the influence of different discretizations. All computations were carried out using MATLAB with about 15 significant decimal digits. n ℓ h ℓ α x δ α,n,ℓ − x n 2 / x n 2 x δ α,n,ℓ − x δ α,n 2 / x n 2 1000 20 1.14 · 10 −1 4.90 2.28 · 10 −1 3.64 · 10 −4 1000 30 1.13 · 10 −1 4.96 2.28 · 10 −1 3.60 · 10 −4 1000 40 1.    where the solution x(t), kernel κ(s, t), and right-hand side y(s) are given by We discretize this integral equation by a Nyström method based on the composite trapezoidal rule with n nodes. This yields a nonsymmetric matrix A n ∈ R n,n . The vector x n ∈ R n is a discretization of the exact solution (5.2). We define the associated right-hand side y n = A n x n , which is assumed to be unknown. An associated contaminated right-hand side, y δ n ∈ R n , which is assumed to be known, is obtained by adding a vector e n ∈ R n with normally distributed random entries with mean zero, that models "noise," to y n . The noise vector e n is scaled to correspond to a prescribed noise level ν = e n 2 y n 2 .
We will use δ = ν y n 2 when determining the regularization parameter α by solving (4.14). Application of ℓ steps of the Arnoldi process to the matrix A n with initial vector v 1 = y δ / y δ 2 yields the decomposition (3.1), as well as the low-rank approximation A (ℓ) n of A n defined by (4.1).
cf. (4.2). We determine the regularization parameter α by solving (4.14) with E = 3 x n 2 and C = 1, as suggested by Proposition 4.2, and then solve the regularized problem (4.4) with the low-rank matrix A (ℓ) n for x δ α,n,ℓ . The inequality (4.15) holds for all examples in this section. Table 5.1 shows the relative error x δ α,n,ℓ − x n 2 / x n 2 . This error depends both on the error in y δ n and on the approximation error (5.4). For fixed n, the approximation error (5.4) is seen to decrease as ℓ increases in Table 5.1.
Let x δ α,n denote the solution of the regularized problem (4.5) with the matrix A n . We are interested in how much the replacement of A n by the low-rank approximation A (ℓ) n affects the quality of the computed solution. Therefore, we tabulate the normalized difference x δ α,n,ℓ − x δ α,n 2 / x n 2 . Table 5.1 shows this difference to be much smaller than x δ α,n,ℓ − x n 2 / x n 2 . Hence, the use of A (ℓ) n instead of A n , with a fixed value of α, does not affect the quality of the computed solution significantly. Table 5.1 shows results for different problem sizes, n ∈ {1000, 2000, 4000}, and noise level 1%. The quality of the computed solution x δ α,n,ℓ is seen not to be very sensitive to the problem size n or to the number of steps ℓ carried out with the Arnoldi process. For n fixed, Table 5.1 shows h ℓ to decrease as ℓ increases. Also the relative error x δ α,n,ℓ − x n 2 / x n 2 can be seen to decrease slowly as ℓ increases. Moreover, the error decreases when n increases and ℓ is kept fixed.
The quality of the computed solution is, of course, sensitive to the noise level. This is illustrated by Table 5.2, which shows results for noise level 0.1%. The α-values of Table 5.2 are smaller than of Table 5.1, as can be expected. Moreover, the relative errors x δ α,n,ℓ − x n 2 / x n reported in Table 5.2 are smaller than the corresponding errors of  Table  5.1. The size of the largest singular value is seen to be independent of ℓ, while the smallest singular value decreases slowly as ℓ and n increase. ✷ Example 5.2. This example also considers the integral equation (5.1), but uses a different discretization. The discretization is computed with the MATLAB function phillips from Regularization Tools by Hansen [12]. This function uses a Galerkin method with n orthonormal box functions as test and trial functions and yields a symmetric indefinite matrix A n ∈ R n×n . The vector x n ∈ R n is a scaled discretization n ℓ σ   of the exact solution (5.2). Since the matrix A n is symmetric, the Arnoldi process (Algorithm 1) simplifies to the Lanczos process. Table 5.4 is analogous to Table  5.1 and shows results for the noise level (5.3) 1%. Results for noise level 0.1% are displayed in Table 5.5, which is analogous to Table 5.2. Due to the different scaling of matrices and right-hand sides in this and the previous examples, the quantities h ℓ and α will differ. However, the relative errors tabulated in the last two columns are comparable, and it is clear that the Galerkin method of the present example furnishes more accurate approximations of x n than the Nyström discretization of Example 5.1. The exact solution x n and the computed approximation x δ α,n,ℓ for n = 2000, ℓ = 30, and ν = 1 · 10 −2 , are shown in Figure 5.1. ✷ Example 5.3. We turn to the Fredholm integral equation of the first kind discussed by Baart [1], π 0 κ(s, t)x(t)dt = g(s), where κ(s, t) = exp(s cos(t)) and g(s) = 2 sinh(s)/s. The solution is given by x(t) = sin(t). We discretize this integral equation by a Galerkin method using n orthonormal box functions as test and trial functions. The discretization is computed with the MATLAB function baart from [12] and gives a nonsymmetric matrix A n ∈ R n,n and n ℓ h ℓ α x δ α,n,ℓ −xn 2 xn 2 x δ α,n,ℓ −x δ α,n 2 xn 2  a vector x n ∈ R n that is a scaled discretization of the exact solution. Similarly as in Example 5.1, we define the "unknown" exact right-hand side by y n = A n x n , and obtain the associated contaminated right-hand side y δ n ∈ R n , which is assumed to be known, by adding a vector e n ∈ R n with normally distributed entries with zero mean to y n . The vector e n is scaled to correspond to a prescribed noise level. A few computed results are displayed in Table 5.6. The table shows the relative error x δ α,n,ℓ − x n 2 / x n 2 to be independent of n for n large, and to decrease as the noise level (5.3) decreases.
We remark that since the singular values of A decrease exponentially with their index number, the condition (2.5) is not valid for any finite l. Nevertheless, this example illustrates that the approximation method described in this paper also can be applied in this situation. ✷ 6. Conclusion and extensions. The paper presents an analysis of the influence of discretization and truncation errors on the computed approximate solution. These errors are caused by replacing an operator A first by a large matrix A n , which in turn is approximated by a matrix A (ℓ) n of rank at most ℓ ≪ n. The choice of the regularization parameter in Tikhonov regularization is discussed. Computed example illustrate the theory.
The matrix A (ℓ) n is determined by the application of ℓ steps of the Arnoldi process to the large matrix A n . Other approaches to determine low-rank approximations are available, such as methods based on Golub-Kahan bidiagonalization or block Golub-Kahan bidiagonalization; see, e.g., Bentbib et al. [2] and Gazzola et al. [8]. The analyses for these methods differ from the one of this paper and are presently being pursued.