A regularizing L-curve Lanczos method for underdetermined linear systems
Introduction
In this article, we describe a new regularizing iterative method for the solution of ill conditioned underdetermined linear systems of equationsLinear 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 matrixThe solution of the minimal norm of (1) can be written as , where solves the linear systemWe regularize this system by adding a multiple of the identity matrix to B and obtainwhere λ>0 is a regularization parameter. We denote the solution of (3) by . The smaller the value of λ, the closer the solution 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 of (1), obtained by first solving (3) for , 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 problemwhere ∥·∥ denotes the Euclidean norm here and in the sequel. The solution of (4) satisfies the equation
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 be contaminated by errors, e.g., due to inaccurate measurements, and and denote the unperturbed right-hand side and the error, respectively. If the norm of the error 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 () 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 matrixand the orthogonal n×k matrix
The matrices Vk and Wk are connected by the two equations
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 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)
- et al.
A regularizing Lanczos iteration method for undetermined linear systems
J. Comp. Appl. Math.
(2000) - et al.
Regularization methods for large-scale problems
Surv. Math. Ind.
(1993) Image Reconstruction from Projections
(1980)- et al.
Solving Least Square Problems
(1995) - et al.
The use of the L-curve in the regularization of discrete ill posed problems
SIAM J. Sci. Comput.
(1993)
Cited by (22)
Distributed bending stiffness estimation of bridges using adaptive inverse unit load method
2023, Engineering Structures3D reconstruction of nonuniform and unsteady flow velocity distribution using nonlinear acoustic ray tracing and tomography
2023, International Journal of Heat and Mass TransferA study on the influence of measurement location and regularization on the evaluation of boundary tractions by inverse finite element method
2009, Finite Elements in Analysis and DesignInverse 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 EngineeringA Novel Recurrent Neural Network for Robot Control
2023, SpringerBriefs in Computer Science
- 1
This research was supported by the University of Bologna, funds for research topics.