Rank-Deficient and Discrete Ill-Posed Problems

The front matter includes the title page, series page, copyright page, TOC, preface, and symbols and acronyms

A very large condition number of the coefficient matrix A in a linear system of equations A x =b implies that some (or all) of the equations are numerically linearly dependent. Hence, the standard advice to avoid solving such systems numerically is not bad. Indeed, it is sometimes the case that the large condition number is caused by an incorrect mathematical model which should be modified before one attempts to compute a numerical solution. Numerical “tools,” such as the singular value decomposition (SVD) (see §2.1), can identify the linear dependencies and thus help to improve the model and lead to a modified system with a better-conditioned matrix. This modified system can then be solved by standard numerical techniques [36], [154], [230].

However, there are classes of problems for which the coefficient matrix is correctly very ill conditioned, i.e., where this property is part of the formulation of the problem. Then the standard linear algebra techniques no longer apply, and the numerical treatment often becomes more difficult.

This book gives a survey of advanced numerical methods for solving such problems with ill-conditioned matrices.

1.1. Problems with Ill-Conditioned Matrices

The numerical treatment of very ill conditioned linear systems of equations is more complicated than the treatment of well-conditioned systems, for the following two reasons.

• The user should know what kind of ill conditioning to expect and how to deal with it. Is the problem rank deficient or ill posed? Is it possible to regularize, i.e., include additional information to stabilize the solution? What additional information is available and is it suited for stabilization purposes?

• The user should also know which numerical regularization method should be used to treat the problem efficiently and reliably on a computer. Are both the analysis and the solution phases important? Should one prefer a direct or an iterative method? How much stabilization should be added?

Before we turn to the numerical regularization methods in the following chapters, it is convenient to summarize the “canonical decompositions” which are important theoretical as well as computational tools in connection with rank-deficient and discrete ill-posed problems. Along with this discussion, it is natural to briefly summarize methods for transforming regularization problems in general form into standard form. We also describe how the singular value expansion (SVE) of a kernel can be computed by means of the singular value decomposition (SVD).

2.1. The SVD and Its Generalizations

The superior numerical “tools” for analysis of rank-deficient and discrete ill-posed problems are the (ordinary) SVD of A and its generalization to two matrices, the generalized SVD (GSVD) of the matrix pair (A, L); see [154, §§2.5.3–4 and 8.7.3]. The SVD reveals all the difficulties associated with the ill-conditioning of the matrix A, while the GSVD of (A, L) yields important insight into regularization problems involving both the matrix A and the regularization matrix L, such as in (1.14).

The use of the SVD and the GSVD in the analysis of discrete ill-posed problems goes back to Hanson [201] and Varah [352], [353]. The SVD has many similarities with the SVE discussed in §1.2.1, and the early history of both is described in [324].

2.1.1. The (Ordinary) SVD

Let A ∈ ℝm×n be a rectangular or square matrix, and assume for ease of presentation that m≥n .

In this chapter we discuss numerical methods that are suited for regularization of problems with a numerically rank-deficient coefficient matrix A, i.e., problems for which there is a well-determined gap between the large and small singular values of A.

Such problems come in two “flavors” : those that involve the solution of a (possibly overdetermined) system of linear equations, and those that involve the computation of a rank-k matrix approximation to a given matrix. As we shall see, these two problems are strongly connected—although the applications in which the problems arise can be very different.

We start with a discussion of the numerical rank of a matrix, as defined via the singular value decomposition (SVD). Next we describe various solutions to rank-deficient problems, defined in terms of the SVD and generalized SVD (GSVD), and we give some perturbation bounds for these solutions. Finally, we focus on algorithms that use rank-revealing decompositions as computational alternatives to the SVD, and we relate the corresponding solutions to the SVD-based solutions.

While the difference between two vectors is naturally measured as the norm of their differences, it is perhaps not intuitively clear how to compare two sub-spaces. The subspace angle is a convenient tool for measuring the difference between two subspaces. Given the two subspaces S1 and S2 of the same dimension, spanned by the orthonormal columns of the matrices V1 and V2 , the subspace angle Θ between S1 and S2 is defined via the relation sin Θ =∥ V1 V 1 T − V2 V 2 T ∥ 2 , 3.1 where Vi V i T is the orthogonal projection matrix associated with the subspace Si . We note that sin Θ is a distance function, and therefore sin Θ =0 if and only if S1 = S2 . See [154, §2.6] or [368] for more information about orthogonal projections and subspace angles.

3.1. Numerical Rank

The rank of a matrix A is defined as the number of linearly independent columns of A.

The purpose of this chapter is to summarize important results about discrete ill-posed problems, i.e., systems of equations (either square or overdetermined) derived from discretization of ill-posed problems. The main feature of these problems is that all the singular values of the coefficient matrix decay gradually to zero, with no gap anywhere in the spectrum. Whatever threshold ϵ is used in Eq. (3.3), the numerical ϵ-rank is highly ill determined, and therefore the concept of “numerical rank” is not useful for these problems.

As a consequence, the regularization of discrete ill-posed problems is more complicated than merely filtering out a cluster of small singular values. For this reason, it is convenient to have a variety of mathematical tools at hand for obtaining more insight into the problem as well as the available regularization methods. Among these tools we find the filter factors, the resolution matrix, and the L-curve, all of which are described in detail below. Numerical examples that illustrate all these tools are presented in the last section of this chapter.

4.1. Characteristics of Discrete Ill-Posed Problems

From a strictly mathematical point of view, a finite-dimensional problem always satisfies the Picard condition (1.10), the minimum-norm solution is stable, and no regularization is required. Indeed, in a purely mathematical sense the transformation of a continuous problem to a discrete problem (“regularization by discretization”) always has a regularizing effect; see, e.g., [160, Chapter 4], [226, Chapter 3], or [227, Chapter 17]. However, this point of view does not account for the disastrous effects of rounding errors when the system is solved, due to the huge condition number of the coefficient matrix; cf. [160, Eq. (4.10)].

In practical treatments of discrete ill-posed problems it is therefore necessary to incorporate some kind of regularization in the solution procedure for the discretized system A x=b or min ∥A x−b∥2 , in order to compute a useful solution. It is also convenient to introduce the concept of a discrete Picard condition, equivalent to the Picard condition described in §4.5.

From a general perspective, we can essentially choose between two classes of regularization methods: those that are based on some kind of “canonical decomposition” such as the QR factorization or the singular value decomposition (SVD), and those that avoid such decompositions. In the first class of methods we find the direct methods treated in this chapter, while the iterative methods (treated in the next chapter) belong to the second class of algorithms.

In this chapter we also briefly consider a third class of regularization problems, namely, those which lead to nonlinear optimization problems.

Since all direct regularization methods are based on standard operations and decompositions in numerical linear algebra (such as matrix multiplications, backsubstitutions, orthogonal transformations, QR factorizations, bidiagonal reductions, SVDs, etc.), the computational effort involved in a direct method can be estimated a priori.

The standard linear algebra operations and decompositions are almost always available on high-performance computers in highly efficient implementations, often based on the basic linear algebra subprograms (BLAS) and/or the LAPACK Library [4]. Efficient parallel numerical linear algebra routines [25], [83], [165] are available on most parallel computers. Hence, it is not too difficult to develop an efficient implementation of a direct regularization method.

We wish to emphasize again that there are no computational difficulties involved with the use of rectangular regularization matrices L in connection with general-form regularization, when the transformation algorithms from §2.3 are used. Moreover, as long as L is a banded matrix, the computational overhead is negligible compared to standard-form regularization.

5.1. Tikhonov Regularization

The numerical treatment of Tikhonov regularization has been described by many authors, and we include it here mainly for completeness. The method was developed independently by Phillips [280] (see also Twomey's paper [342]) and Tikhonov [336]. Some recent presentations of the method can be found in [107], [111, Chapter 5], [160], [161, §5.1], and [171].

Iterative methods for linear systems of equations and linear least squares problems are based on iteration schemes that access the coefficient matrix A only via matrix-vector multiplications with A and AT , and they produce a sequence of iteration vectors x (k) , k=1,2,… , that converge to the desired solution. Iterative methods are preferable to direct methods when the coefficient matrix is so large that it is too time-consuming or too memory-demanding to work with an explicit decomposition of A.

Obviously, we can apply iterative methods to the symmetric positive definite Tikhonov system ( AT A+ λ2 LT L) x= AT b or the equivalent least squares problem, and in this way compute a regularized solution [47], [273]. However, this approach requires a new iteration process for each new value of λ, and for general matrices it is difficult to construct efficient preconditioners.

Instead, we are interested in iterative regularization methods in which each iteration vector x (k) can be considered as a regularized solution, with the iteration number k playing the role of the regularization parameter. Hence, we need iteration schemes with the intrinsic property that they, when applied to discrete ill-posed problems, initially pick up those singular value decomposition (SVD) components ( u i T b/ σi ) vi corresponding to the largest singular values (which are the ones that are desired in a regularized solution), in such a way that the iteration number k can be considered as a regularization parameter. This phenomenon is sometimes referred to as semiconvergence [259, p. 89], because the iteration vector x (k) initially approaches a regularized solution and then, in later stages of the iterations, converges to some other undesired vector—often the least squares solution xLS = A† b .

Below, we first describe some classical iterative methods that have been studied extensively in the literature, and then we focus on the use of the conjugate gradient (CG) algorithm and related algorithms based on Lanczos bidiagonalization. We give a detailed description of the regularizing properties of these algorithms in terms of the associated filter factors, and we discuss the influence of finite-precision arithmetic. Finally, we describe how a hybrid method can be obtained by performing “inner” regularization in each step of an iterative regularization algorithm.

So far, we have discussed various algorithms for computing a regularized solution. Every algorithm has its advantages and disadvantages in terms of implementation issues, filter properties, etc. However, no regularization method is complete without a method for choosing the regularization parameter, either the continuous parameter λ or the discrete parameter k, and in this chapter we discuss several parameter-choice methods and their implementations. Surveys of parameter-choice methods are rare in the literature; we are currently aware of [171] and [333].

Whenever possible, we discuss these methods in terms of a continuous regularization parameter λ, but we emphasize that similar results hold when λ is replaced by a discrete regularization parameter k. We will also restrict our discussion to the standard-form case whenever the extensions to the general case are obvious, and we use the notation xλ = x In ,λ .

Most of the parameter-choice methods are based on residual norms and, in the case of the L-curve, also on the solution's “size” Ω, typically its seminorm. When solving general-form problems via a standard-form transformation, it is therefore important to recall the relations ∥L x ∥2 = ∥ x ¯ ∥2 and ∥A x−b ∥2 =∥ A ¯ x ¯ − b ¯ ∥2 ; cf. (2.39) and (2.45). These relations ensure that application of a norm-based parameter-choice rule to the original problem with A and b, or to the standard-form problem with A ¯ and b ¯ , yields exactly the same regularization parameter.

For certain problems it may be advantageous to use more information about the residual vector than its norm. For example, the residual's autocorrelation may be useful; see [13] for an example. Such methods do not seem to have found widespread use.

We first describe some parameter-choice methods that make explicit use of the norm of the error in the right-hand side. Then we turn to methods that do not require this information, in particular, the generalized cross-validation method and the L-curve method. We finish the chapter with two extensive numerical examples that illustrate the behavior of the parameter-choice methods used in connection with some of the direct and iterative regularization methods presented in the previous two chapters.

The software package Regularization Tools [187], Version 3.1, consists of 53 Matlab routines for analysis and solution of discrete ill-posed problems. By means of this package, the user can experiment with different regularization strategies, compare them, and draw conclusions that would otherwise require a major programming effort. In addition to the analysis and solution routines, the package also includes 12 test problems. The package, as well as the corresponding manual [186], is available from Netlib in the file numeralgo/na4, as well as from the author.

The purpose of Regularization Tools is to provide the user with easy-to-use routines, based on numerically robust and efficient algorithms, for doing experiments with regularization of discrete ill-posed problems. By means of this package, the user can experiment with different regularization strategies, compare them, and draw conclusions that would otherwise require a major programming effort. For discrete ill-posed problems, which are indeed difficult to treat numerically, such an approach is superior to a single black-box routine.

There have been other attempts to write general software for discrete ill-posed problems, all of them implementing Tikhonov regularization. They are listed here in chronological order.

• What appears to be the first package was written in the object-oriented language Simula by Eldén [99], and it concentrated on the use of bidiagonalization (§5.1.1) and standard-form transformation (§2.3.1) to solve the Tikhonov problem in its three incarnations (5.1), (5.14), and (5.15).

• The Fortran program CONTIN by Provencher [288] uses the quadrature method to discretize the integral equation (1.2). Equality and inequality constraints can be imposed on the discretized solution. The regularization parameter is chosen either by solving V (λ) = σ 0 2 (cf. (7.9)) or via a confidence region estimate. The linear systems are solved by means of orthogonal transformations and SVD, and the inequality constraints are handled by means of the LDP least distance programming routine from [231].

The back matter includes bibliography and index