GMRES, L-Curves, and Discrete Ill-Posed Problems

The GMRES method is a popular iterative method for the solution of large linear systems of equations with a nonsymmetric nonsingular matrix. This paper discusses application of the GMRES method to the solution of large linear systems of equations that arise from the discretization of linear ill-posed problems. These linear systems are severely ill-conditioned and are referred to as discrete ill-posed problems. We are concerned with the situation when the right-hand side vector is contaminated by measurement errors, and we discuss how a meaningful approximate solution of the discrete ill-posed problem can be determined by early termination of the iterations with the GMRES method. We propose a termination criterion based on the condition number of the projected matrices defined by the GMRES method. Under certain conditions on the linear system, the termination index corresponds to the vertex of an L-shaped curve.     

An extrapolated TSVD method for linear discrete ill-posed problems with Kronecker structure

This paper describes a new numerical method for the solution of large linear discrete ill-posed problems, whose matrix is a Kronecker product. Problems of this kind arise, for instance, from the discretization of Fredholm integral equations of the first kind in two space-dimensions with a separable kernel. The available data (right-hand side) of many linear discrete ill-posed problems that arise in applications is contaminated by measurement errors. Straightforward solution of these problems generally is not meaningful because of severe error propagation. We discuss how to combine the truncated singular value decomposition (TSVD) with reduced rank vector extrapolation to determine computed approximate solutions that are fairly insensitive to the error in the data. Exploitation of the structure of the problem keeps the computational effort quite small.     

约束Fractional Tikhonov正则化的模迭代方法

反问题是现在数学物理研究中的一个热点问题,而反问题求解面临的一个本质性困难是不适定性。求解不适定问题的普遍方法是:用与原不适定问题相“邻近”的适定问题的解去逼近原问题的解,这种方法称为正则化方法.如何建立有效的正则化方法是反问题领域中不适定问题研究的重要内容.当前,最为流行的正则化方法有基于变分原理的Tikhonov正则化及其改进方法,此类方法是求解不适定问题的较为有效的方法,在各类反问题的研究中被广泛采用,并得到深入研究.     

On the regularizing properties of the GMRES method

Summary. The GMRES method is a popular iterative method for the solution of large linear systems of equations with a nonsymmetric nonsingular matrix. However, little is known about the behavior of this method when it is applied to the solution of nonsymmetric linear ill-posed problems with a right-hand side that is contaminated by errors. We show that when the associated error-free right-hand side lies in a finite-dimensional Krylov subspace, the GMRES method is a regularization method. The iterations are terminated by a stopping rule based on the discrepancy principle. Received November 10, 2000 / Revised version received April 11, 2001 / Published online October 17, 2001     

Estimation of the L-Curve via Lanczos Bidiagonalization

The L-curve criterion is often applied to determine a suitable value of the regularization parameter when solving ill-conditioned linear systems of equations with a right-hand side contaminated by errors of unknown norm. However, the computation of the L-curve is quite costly for large problems; the determination of a point on the L-curve requires that both the norm of the regularized approximate solution and the norm of the corresponding residual vector be available. Therefore, usually only a few points on the L-curve are computed and these values, rather than the L-curve, are used to determine a value of the regularization parameter. We propose a new approach to determine a value of the regularization parameter based on computing an L-ribbon that contains the L-curve in its interior. An L-ribbon can be computed fairly inexpensively by partial Lanczos bidiagonalization of the matrix of the given linear system of equations. A suitable value of the regularization parameter is then determined from the L-ribbon, and we show that an associated approximate solution of the linear system can be computed with little additional work.     

A new L-curve for ill-posed problems

The truncated singular value decomposition is a popular method for the solution of linear ill-posed problems. The method requires the choice of a truncation index, which affects the quality of the computed approximate solution. This paper proposes that an L-curve, which is determined by how well the given data (right-hand side) can be approximated by a linear combination of the first (few) left singular vectors (or functions), be used as an aid for determining the truncation index.     

A generalized global Arnoldi method for ill-posed matrix equations

This paper discusses the solution of large-scale linear discrete ill-posed problems with a noise-contaminated right-hand side. Tikhonov regularization is used to reduce the influence of the noise on the computed approximate solution. We consider problems in which the coefficient matrix is the sum of Kronecker products of matrices and present a generalized global Arnoldi method, that respects the structure of the equation, for the solution of the regularized problem. Theoretical properties of the method are shown and applications to image deblurring are described.     

A truncated projected SVD method for linear discrete ill-posed problems

Truncated singular value decomposition is a popular solution method for linear discrete ill-posed problems. However, since the singular value decomposition of the matrix is independent of the right-hand side, there are linear discrete ill-posed problems for which this method fails to yield an accurate approximate solution. This paper describes a new approach to incorporating knowledge about properties of the desired solution into the solution process through an initial projection of the linear discrete ill-posed problem. The projected problem is solved by truncated singular value decomposition. Computed examples illustrate that suitably chosen projections can enhance the accuracy of the computed solution.     

An L-ribbon for large underdetermined linear discrete ill-posed problems

The L-curve is a popular aid for determining a suitable value of the regularization parameter when solving ill-conditioned linear systems of equations with a right-hand side vector, which is contaminated by errors of unknown size. However, for large problems, the computation of the L-curve can be quite expensive, because the determination of a point on the L-curve requires that both the norm of the regularized approximate solution and the norm of the corresponding residual vector be available. Recently, an approximation of the L-curve, referred to as the L-ribbon, was introduced to address this difficulty. The present paper discusses how to organize the computation of the L-ribbon when the matrix of the linear system of equations has many more columns than rows. Numerical examples include an application to computerized tomography.     

Deflated GMRES for systems with multiple shifts and multiple right-hand sides

We consider solution of multiply shifted systems of nonsymmetric linear equations, possibly also with multiple right-hand sides. First, for a single right-hand side, the matrix is shifted by several multiples of the identity. Such problems arise in a number of applications, including lattice quantum chromodynamics where the matrices are complex and non-Hermitian. Some Krylov iterative methods such as GMRES and BiCGStab have been used to solve multiply shifted systems for about the cost of solving just one system. Restarted GMRES can be improved by deflating eigenvalues for matrices that have a few small eigenvalues. We show that a particular deflated method, GMRES-DR, can be applied to multiply shifted systems.In quantum chromodynamics, it is common to have multiple right-hand sides with multiple shifts for each right-hand side. We develop a method that efficiently solves the multiple right-hand sides by using a deflated version of GMRES and yet keeps costs for all of the multiply shifted systems close to those for one shift. An example is given showing this can be extremely effective with a quantum chromodynamics matrix.     

An iterative method for linear discrete ill-posed problems with box constraints

Many questions in science and engineering give rise to linear discrete ill-posed problems. Often it is desirable that the computed approximate solution satisfies certain constraints, e.g., that some or all elements of the computed solution be nonnegative. This paper describes an iterative method of active set-type for the solution of large-scale problems of this kind. The method employs conjugate gradient iteration with a stopping criterion based on the discrepancy principle and allows updates of the active set by more than one index at a time.     

Algorithms for range restricted iterative methods for linear discrete ill-posed problems

Range restricted iterative methods based on the Arnoldi process are attractive for the solution of large nonsymmetric linear discrete ill-posed problems with error-contaminated data (right-hand side). Several derivations of this type of iterative methods are compared in Neuman et al. (Linear Algebra Appl. in press). We describe MATLAB codes for the best of these implementations. MATLAB codes for range restricted iterative methods for symmetric linear discrete ill-posed problems are also presented.     

Implementations of range restricted iterative methods for linear discrete ill-posed problems

This paper is concerned with iterative solution methods for large linear systems of equations with a matrix of ill-determined rank and an error-contaminated right-hand side. The numerical solution is delicate, because the matrix is very ill-conditioned and may be singular. It is natural to require that the computed iterates live in the range of the matrix when the latter is symmetric, because then the iterates are orthogonal to the null space. Computational experience indicates that it can be beneficial to require that the iterates live in the range of the matrix also when the latter is nonsymmetric. We discuss the design and implementation of iterative methods that determine iterates with this property. New implementations that are particularly well suited for use with the discrepancy principle are described.     

Analysis and solution of the Cauchy problem for linear Volterra integro-differential equations with functional delay

We study how to transform Cauchy problems for Volterra integro-differential equations with functional delays to resolving Volterra integral equations with conventional argument by using a modification of a function of flexible structure. We show that such a transformation is possible for all linear Volterra integro-differential equations of retarded type. There exists a unique solution of the resolving equation provided that the kernels and the right-hand side are bounded in the closed square. The presence of parameters in the expression for the function of flexible structure permits one to choose these parameters in an optimal way in the course of the solution of the problem so as to represent the solution in closed form or, if this is difficult, optimize an approximate solution method. The accuracy of the approximate solutions is estimated.     

Optimal control of systems with discontinuous differential equations

In this paper we discuss the problem of verifying and computing optimal controls of systems whose dynamics is governed by differential systems with a discontinuous right-hand side. In our work, we are motivated by optimal control of mechanical systems with Coulomb friction, which exhibit such a right-hand side. Notwithstanding the impressive development of nonsmooth and set-valued analysis, these systems have not been closely studied either computationally or analytically. We show that even when the solution crosses and does not stay on the discontinuity, differentiating the results of a simulation gives gradients that have errors of a size independent of the stepsize. This means that the strategy of “optimize the discretization” will usually fail for problems of this kind. We approximate the discontinuous right-hand side for the differential equations or inclusions by a smooth right-hand side. For these smoothed approximations, we show that the resulting gradients approach the true gradients provided that the start and end points of the trajectory do not lie on the discontinuity and that Euler’s method is used where the step size is “sufficiently small” in comparison with the smoothing parameter. Numerical results are presented for a crude model of car racing that involves Coulomb friction and slip showing that this approach is practical and can handle problems of moderate complexity.     

Cascadic multilevel methods for ill-posed problems

Multilevel methods are popular for the solution of well-posed problems, such as certain boundary value problems for partial differential equations and Fredholm integral equations of the second kind. However, little is known about the behavior of multilevel methods when applied to the solution of linear ill-posed problems, such as Fredholm integral equations of the first kind, with a right-hand side that is contaminated by error. This paper shows that cascadic multilevel methods with a conjugate gradient-type method as basic iterative scheme are regularization methods. The iterations are terminated by a stopping rule based on the discrepancy principle.     

A stopping criterion for the iterative solution of an overdetermined system of linear algebraic equations

For an overdetermined system of linear algebraic equations, systems obtained by introducing independent random errors into the original right-hand side are examined. Under certain assumptions on how these random variables are distributed, a practical stopping criterion is proposed for an iterative process that minimizes the sum of the squares of the residuals for the above systems. Numerical results demonstrating the efficiency of this criterion for some ill-conditioned problems are presented.     

Sparse approximate inverse preconditioning for dense linear systems arising in computational electromagnetics

We investigate the use of sparse approximate inverse preconditioners for the iterative solution of linear systems with dense complex coefficient matrices arising in industrial electromagnetic problems. An approximate inverse is computed via a Frobenius norm approach with a prescribed nonzero pattern. Some strategies for determining the nonzero pattern of an approximate inverse are described. The results of numerical experiments suggest that sparse approximate inverse preconditioning is a viable approach for the solution of large-scale dense linear systems on parallel computers. This revised version was published online in August 2006 with corrections to the Cover Date.     

Approximation of solutions to second order nonlinear Picard problems with Carathéodory right-hand side

We present an approximation method for Picard second order boundary value problems with Carathéodory righthand side. The method is based on the idea of replacing a measurable function in the right-hand side of the problem with its Kantorovich polynomial. We will show that this approximation scheme recovers essential solutions to the original BVP. We also consider the corresponding finite dimensional problem. We suggest a suitable mapping of solutions to finite dimensional problems to piecewise constant functions so that the later approximate a solution to the original BVP. That is why the presented idea may be used in numerical computations.     

Fast solution methods for fredholm integral equations of the second kind

Summary The main purpose of this paper is to describe a fast solution method for one-dimensional Fredholm integral equations of the second kind with a smooth kernel and a non-smooth right-hand side function. Let the integral equation be defined on the interval [–1, 1]. We discretize by a Nyström method with nodes {cos(j/N)} _j ^=0/N . This yields a linear system of algebraic equations with an (N+1)×(N+1) matrixA. GenerallyN has to be chosen fairly large in order to obtain an accurate approximate solution of the integral equation. We show by Fourier analysis thatA can be approximated well by , a low-rank modification of the identity matrix. ReplacingA by in the linear system of algebraic equations yields a new linear system of equations, whose elements, and whose solution , can be computed inO (N logN) arithmetic operations. If the kernel has two more derivatives than the right-hand side function, then is shown to converge optimally to the solution of the integral equation asN increases.We also consider iterative solution of the linear system of algebraic equations. The iterative schemes use bothA andÃ. They yield the solution inO (N ²) arithmetic operations under mild restrictions on the kernel and the right-hand side function.Finally, we discuss discretization by the Chebyshev-Galerkin method. The techniques developed for the Nyström method carry over to this discretization method, and we develop solution schemes that are faster than those previously presented in the literature. The schemes presented carry over in a straightforward manner to Fredholm integral equations of the second kind defined on a hypercube.