A regularizing L-curve Lanczos method for underdetermined linear systems

https://doi.org/10.1016/S0096-3003(99)00262-3Get rights and content

Abstract

Many real applications give rise to the solution of underdetermined linear systems of equations with a very ill conditioned matrix A, whose dimensions are so large as to make solution by direct methods impractical or infeasible. Image reconstruction from projections is a well-known example of such systems. In order to facilitate the computation of a meaningful approximate solution, we regularize the linear system, i.e., we replace it by a nearby system that is better conditioned. The amount of regularization is determined by a regularization parameter. Its optimal value is, in most applications, not known a priori. A well-known method to determine it is given by the L-curve approach. We present an iterative method based on the Lanczos algorithm for inexpensively evaluating an approximation of the points on the L-curve and then determine the value of the optimal regularization parameter which lets us compute an approximate solution of the regularized system of equations.

Introduction

In this article, we describe a new regularizing iterative method for the solution of ill conditioned underdetermined linear systems of equationsAx=b,A∈Rm×n,xRn,bRm,m<n.Linear systems of this kind arise, for example, in the reconstruction of images from projections. In this application, each column corresponds to one pixel of the image to be reconstructed, and each row to a “ray” from an emitter to a detector. The matrices that arise are typically large and have many more columns than rows. Further, they are generally very ill conditioned, because the linear system is obtained by discretizing an ill posed problem, the inversion of a Radon transform. The ill conditioning makes a solution of (1) very sensitive to perturbations in the right-hand side vector. Such perturbations arise, e.g., from measurement errors [1], [2].

In order to be able to compute a meaningful approximate solution of (1), the system has to be regularized, i.e., the linear system (1) has to be replaced by a linear system of equations with a not very ill conditioned matrix, such that the unique solution of this system yields an acceptable approximate solution of (1). We introduce the matrixB:=AAT.The solution of the minimal norm of (1) can be written as x:=ATy, where yRm solves the linear systemBy=b.We regularize this system by adding a multiple of the identity matrix to B and obtain(B+λI)y=b,where λ>0 is a regularization parameter. We denote the solution of (3) by yλ. The smaller the value of λ, the closer the solution yλ of (3) is to the minimal norm solution of (2), and the more ill conditioned the matrix B+λI. We wish to find a value of λ in (3), such that the matrix B+λI is not very ill conditioned, and the solution of (3) is an acceptable approximate solution of (2). Generally, it is not known a priori for which values of λ these requirements are satisfied.

The approximate solution xλ=ATyλ of (1), obtained by first solving (3) for yλ, can also be computed by standard Tikhonov regularization of system (1). This approach, however, as has been shown in [3], requires more computer storage and arithmetic work when the linear system has many more columns than rows.

Let us consider the minimization problemminyRm{∥Byb2+λ∥ATy2},where ∥·∥ denotes the Euclidean norm here and in the sequel. The solution yλ of (4) satisfies the equation(BTB+λAAT)y=BTb.

The comparison of , , shows the equivalence of problem (3) and (5) when B is a nonsingular matrix and this justifies what follows.

Let components of the right-hand side vector b be contaminated by errors, e.g., due to inaccurate measurements, and bexact and e denote the unperturbed right-hand side and the error, respectively. If the norm of the error e is known, then the value of λ can be chosen so that the associated solution satisfies the Morozov discrepancy principle [3]. However, in many applications the norm of e is not explicitly known, and other approaches to determine a value of the regularization parameter λ have to be employed.

We will show that the plot of the curve (∥Byλb∥,∥ATyλ) is shaped roughly like an “L”, and the optimal value of the regularization parameter λ corresponds to the “vertex” point of the L-curve. This idea of L-curve was observed by Lawson and Hanson [4] to determine the regularizing parameter λ in the Tikhonov approach, while the choice of the optimal regularization parameter as the vertex point was suggested recently in [5], [6], and successively considered in [7].

The computation of the L-curve is quite costly for large problems; in fact, 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, which we will refer to as the discrete L-curve, are used to determine a value of the regularization parameter.

In [8], a new method for approximately determining the points on the L-curve is proposed. They showed how rectangular regions that contain points on the L-curve (L-ribbon) can be determined fairly inexpensively for several values of the regularization parameter without having to solve the corresponding regularized minimization problem but using Gauss and Gauss–Radau quadrature rules. When the L-ribbon is narrow, it is possible to infer the approximate location of the vertex of the L-curve from the shape of the L-ribbon.

In order to solve problem (4), or equivalently (3), we cannot follow the Calvetti et al. [8] approach because their assumptions are not satisfied due to the term AT in the regularizing functional.

We present an iterative procedure that determines successive encapsulated intervals containing the optimal regularization parameter corresponding to the vertex of the L-curve.

At each step, some new points on the L-curve are computed by a few steps of the Lanczos bidiagonalization algorithm, while a new interval for λ is determined by a minimum distance from the local origin criterion.

The final value obtained for the regularization parameter is used to compute an approximation of the solution of the regularized system (3).

This article is organized as follows. In Section 2, we give a brief description of properties of the L-curve and we show how to use it to solve problem (4). Section 3 explains the relation between the bidiagonalization Lanczos algorithm and the approximated L-curve and how to update the information extracted from the Lanczos process to compute a new point. Section 4 describes the regularizing algorithm for the solution of the linear system (3). Section 5 contains several computed examples which illustrate the performance of the algorithm.

Section snippets

L-curve for the choice of regularization parameter value

A recent method to choose the regularization parameter is the so-called L-curve. The L-curve is a parametric plot of (ρ(λ),η(λ)), where η(λ) and ρ(λ) measure the size of the regularized solution and the corresponding residual. The underlying idea is that a good method for choosing the regularization parameter for discrete ill posed problems must incorporate information about the solution size in addition to using information about the residual size. This is indeed quite natural, because we are

The Lanczos algorithm and regularization

In this section, we briefly review the bidiagonalization Lanczos process [9], and, we explore its relation with the functions ρk(λ) and ηk(λ) introduced in Section 2.

Given a symmetric matrix B=AAT and a nonvanishing vector b, if we execute k iterations of the Lanczos bidiagonalization Algorithm 1 in exact arithmetic, we get the orthogonal m×k matrixVk=[v0,v1,…,vk−1]and the orthogonal n×k matrixWk=[w0,w1,…,wk−1].

The matrices Vk and Wk are connected by the two equationsAWk=VkLkkvkekT,ATVk=WkLkT,

A regularizing algorithm based on the Lanczos process

In this section, our method for the solution of underdetermined ill conditioned linear systems (1) is described.

The evaluation of points on the L-curve requires that the exact solution of the linear system (3) is available. However, as already mentioned, computation of the exact or a highly accurate solution of (3) for each value of λ is expensive, and therefore, we replace the equations for ρ(λ) and η(λ) by , , respectively. It follows from Proposition 2 and (28) that the norm of the residual

Numerical results

The purpose of this section is to illustrate the performance of Algorithm 2 in some numerical examples obtained on a Sun Ultra workstation with a unit roundoff of ϵ=2−52≈2.2204×10−16. We implemented our method to compute an approximation of the L-curve in FORTRAN77 language, and we used MATLAB commercial software package [11] to obtain the exact computation of points on the L-curve.

In order to obtain small and large scale test problems, we generated an m-by-n matrix A defined by its singular

Conclusions

We have demonstrated that the approximated L-curve is a versatile and inexpensive tool for the solution of large-scale underdetermined linear ill posed problems. Our numerical examples clearly illustrate how Algorithm 2 reconstructs with a good accuracy the discrete L-curve providing a suitable value for the regularization parameter also in more general problems. At the same time, the proposed algorithm allows us to compare results related to different values of the regularization parameter

References (11)

  • D. Calvetti et al.

    A regularizing Lanczos iteration method for undetermined linear systems

    J. Comp. Appl. Math.

    (2000)
  • M. Hanke et al.

    Regularization methods for large-scale problems

    Surv. Math. Ind.

    (1993)
  • G.T. Herman

    Image Reconstruction from Projections

    (1980)
  • C.L. Lawson et al.

    Solving Least Square Problems

    (1995)
  • P.C. Hansen et al.

    The use of the L-curve in the regularization of discrete ill posed problems

    SIAM J. Sci. Comput.

    (1993)
There are more references available in the full text version of this article.

Cited by (22)

  • Inverse Unit Load Method for Full-Field Reconstruction of Bending Stiffness in Girder Bridges

    2023, ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part A: Civil Engineering
  • A Novel Recurrent Neural Network for Robot Control

    2023, SpringerBriefs in Computer Science
View all citing articles on Scopus
1

This research was supported by the University of Bologna, funds for research topics.

View full text