RTV 4 Manual - Regu Tools

Regularization Tools
A Matlab Package for

Analysis and Solution of Discrete Ill-Posed Problems
Version 4.1 for Matlab 7.3
Per Christian Hansen

Informatics and Mathematical Modelling
Building 321, Technical University of Denmark
DK-2800 Lyngby, Denmark
pch@imm.dtu.dk
http://www.imm.dtu.dk/~pch
March 2008
The software described in this report was originally published in

Numerical Algorithms 6 (1994), pp. 1–35.
The current version is published in Numer. Algo. 46 (2007), pp. 189–194,

and it is available from www.netlib.org/numeralgo
and www.mathworks.com/matlabcentral/fileexchange
Contents
Changes from Earlier Versions 3
1 Introduction 5
2 Discrete Ill-Posed Problems and their Regularization 9

2.1 Discrete Ill-Posed Problems . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Regularization Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3 SVD and Generalized SVD . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3.1 The Singular Value Decomposition . . . . . . . . . . . . . . . . 13
2.3.2 The Generalized Singular Value Decomposition . . . . . . . . . 14
2.4 The Discrete Picard Condition and Filter Factors . . . . . . . . . . . . 16
2.5 The L-Curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.6 Transformation to Standard Form . . . . . . . . . . . . . . . . . . . . 21
2.6.1 Transformation for Direct Methods . . . . . . . . . . . . . . . . 21
2.6.2 Transformation for Iterative Methods . . . . . . . . . . . . . . 22
2.6.3 Norm Relations etc. . . . . . . . . . . . . . . . . . . . . . . . . 24
2.7 Direct Regularization Methods . . . . . . . . . . . . . . . . . . . . . . 25
2.7.1 Tikhonov Regularization . . . . . . . . . . . . . . . . . . . . . . 25
2.7.2 Least Squares with a Quadratic Constraint . . . . . . . . . . . 25
2.7.3 TSVD, MTSVD, and TGSVD . . . . . . . . . . . . . . . . . . 26
2.7.4 Damped SVD/GSVD . . . . . . . . . . . . . . . . . . . . . . . 27
2.7.5 Maximum Entropy Regularization . . . . . . . . . . . . . . . . 28
2.7.6 Truncated Total Least Squares . . . . . . . . . . . . . . . . . . 29
2.8 Iterative Regularization Methods . . . . . . . . . . . . . . . . . . . . . 29
2.8.1 Conjugate Gradients and LSQR . . . . . . . . . . . . . . . . . 29
2.8.2 Bidiagonalization with Regularization . . . . . . . . . . . . . . 32
2.8.3 The ν-Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.8.4 Extension to General-Form Problems . . . . . . . . . . . . . . . 33
2.9 Methods for Choosing the Regularization Parameter . . . . . . . . . . 34
2.10 New Functions in Version 4.1 . . . . . . . . . . . . . . . . . . . . . . . 36
2 CONTENTS
3 Regularization Tools Tutorial 39

3.1 The Discrete Picard Condition . . . . . . . . . . . . . . . . . . . . . . 39
3.2 Filter Factors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.3 The L-Curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.4 Regularization Parameters . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.5 Standard Form Versus General Form . . . . . . . . . . . . . . . . . . . 43
3.6 No Square Integrable Solution . . . . . . . . . . . . . . . . . . . . . . . 46
4 Regularization Tools Reference 47

Routines by Subject Area . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
The Test Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
Alphabetical List of Routines . . . . . . . . . . . . . . . . . . . . . . . . . . 51
Bibliography 121
Changes from Earlier Versions
The following is a list of the major changes since Version 2.0 of the package.
• Replaced gsvd by cgsvd which has a different sequence of output arguments.
• Removed the obsolete function csdecomp (which replaced the function csd)
• Deleted the function mgs.
• Changed the storage format of bidiagonal matrices to sparse, instead of a dense

matrix with two columns.
• Removed the obsolete function bsvd.
• Added the function regutm that generates random test matrices for regulariza-
tion methods.
• Added the blur test problem.
• Functions tsvd and tgsvd now allow k = 0, and functions tgsvd and tikhonov now
allow a square L.
• Added output arguments rho and eta to functions dsvd, mtsvd, tgsvd, tikhonov,
and tsvd.
• Added a priori guess x 0 as input to tikhonov.
• Corrected get l such that the sign of L*x is correct.
• Added MGS reorthogonalization of the normal equation residual vectors in the

two functions cgls and pcgls.
• Added the method ’ttls’ to the function fil fac.
• More precise computation of the regularization parameter in gcv, lcurve, and

quasiopt.
• Changed heb new and newton to work with λ instead of λ2 .
• Added legend to lagrange and picard.

4 CONTENTS
The following major changes were made since Version 3.0 of the package.
• Renamed lsqr and plsqr to lsqr b and plsqr b, respectively, and removed the
option reorth = 2.
• Corrected the routines to work for complex problems.

• Changed eta to seminorm in tgsvd, and in dsvd and tikhonov for the general-form
case.
• Changed cgsvd, discrep, dsvd, lsqi, tgsvd, and tikhonov to allow for an underde-
termined A matrix.
• Added new test problems gravity and tomo.
• Added new parameter-choice methods corner and ncp.
• Added new iterative regularization methods art, mr2, pmr2, prrgmres, rrgmres,
and splsqr.
• Changed l curve and l corner to use the new function corner if the Spline Toolbox
is not available.
• Renamed ilaplace to i laplace (to avoid name overlap with the Symbolic Math
Toolbox).
1. Introduction
Ill-posed problems—and regularization methods for computing stabilized solutions to

the ill-posed problems—occur frequently enough in science and engineering to make
it worthwhile to present a general framework for their numerical treatment. The
purpose of this package of Matlab routines is to provide the user with easy-to-use
routines, based on numerically robust and efficient algorithms, for doing experiments
with analysis and solution of discrete ill-posed problems by means of regularization
methods.
The theory for ill-posed problems is well developed in the literature. We can
easily illustrate the main difficulties associated with such problems by means of a
small numerical example. Consider the following least squares problem
min kA x − bk2
x
with coefficient matrix A and right-hand side b given by

   
0.16 0.10 0.27
A =  0.17 0.11  , b =  0.25  .
2.02 1.29 3.33
Here, the right-hand side b is generated by adding a small perturbation to an exact

right-hand side corresponding to the exact solution x̄T = (1 1):
   
0.16 0.10 µ ¶ 0.01
1.00
b =  0.17 0.11  +  −0.03  .
1.00
2.02 1.29 0.02
The difficulty with this least squares problem is that the matrix A is ill-conditioned;
its condition number is 1.1·103 . This implies that the computed solution is potentially
very sensitive to perturbations of the data. Indeed, if we compute the ordinary least-
squares solution xLSQ by means of a QR factorization of A, then we obtain
µ ¶
7.01
xLSQ = .
−8.40
This solution is obviously worthless, and something must be done in order to compute
a better approximation to the exact solution x̄T = (1 1).
6 Chapter 1. Introduction
The large condition number implies that the columns of A are nearly linearly
dependent. One could therefore think of replacing the ill-conditioned matrix A =
(a1 a2 ) with either (a1 0) or (0 a2 ), both of which are well conditioned. The two
corresponding so-called basic solutions are
µ ¶ µ ¶
(1) 1.65 (2) 0.00
xB = , xB = .
0.00 2.58
Although these solutions are much less sensitive to perturbations of the data, and
although the corresponding residual norms are both small,
(1) (2)
kA xB − bk2 = 0.031 , kA xB − bk2 = 0.036 ,
the basic solutions nevertheless have nothing in common with x̄T = (1 1).
A major difficulty with the ordinary least squares solution xLSQ is that its norm
is significantly greater than the norm of the exact solution. One may therefore try
another approach to solving the least squares problem by adding the side constraint
that the solution norm must not exceed a certain value α,
min kA x − bk2 subject to kxk2 ≤ α .

x
The such computed solution xα depends in a non-linear way on α, and for α equal to
0.1, 1, 1.385, and 10 we obtain
µ ¶ µ ¶ µ ¶ µ ¶
0.08 0.84 1.17 6.51
x0.1 = , x1 = , x1.385 = , x10 = .
0.05 0.54 0.74 −7.60
We see that by a proper choice of α we can indeed compute a solution x1.385 which
is fairly close to the desired exact solution x̄T = (1 1). However, care must be taken
when choosing α, and the proper choice of α is not obvious.
Although the above example is a small one, it highlights the three main difficulties
associates with discrete ill-posed problems:
1. the condition number of the matrix A is large
2. replacing A by a well-conditioned matrix derived from A does not necessarily

lead to a useful solution
3. care must taken when imposing additional constraints.
The purpose of numerical regularization theory is to provide efficient and numerically

stable methods for including proper side constraints that lead to useful stabilized
solutions, and to provide robust methods for choosing the optimal weight given to
these side constraints such that the regularized solution is a good approximation to
the desired unknown solution.
Introduction 7
The routines provided in this package are examples of such procedures. In ad-
dition, we provide a number of utility routines for analyzing the discrete ill-posed
problems in details, for displaying these properties, and for easy generation of simple
test problems. By means of the routines in Regularization Tools, the user can—
at least for small to medium-size problems—experiment with different regularization
strategies, compare them, and draw conclusions from these experiments that would
otherwise require a major programming effort. For discrete ill-posed problems, which
are indeed difficult to treat numerically, such an approach is certainly superior to a
single black-box routine.
The package was mainly developed in the period 1990–1992 at UNI•C and to
some extent it reflects the author’s own work. Prof. Dianne P. O’Leary and Dr.
Martin Hanke helped with the iterative methods. Prof. Lars Eldén—and his 1979
Simula package [23] with the same purpose as Regularization Tools—provided
a great source of inspiration. The package was also inspired by a paper by Natterer
[68] where a “numerical analyst’s toolkit for ill-posed problems” is suggested. The
Fortran programs by Drake [20], te Riele [76], and Wahba and her co-workers [6] also
deserve to be mentioned here.
The package has been officially updated twice since its release in 1994. The first
update [50] from 1999 incorporated several revisions that were made necessary by
changes in Matlab itself, as well as changes suggested by users. The second update is
the present version 4.0 from 2007. The modification that was most often requested
by the users was to allow for under-determined problems, and this is allowed in the
current version, which also incorporates two new test problems, two new parameter-
choice methods, and several new iterative regularization methods.
Acknowledgement. I am indebted to the many users of Regularization Tools

who have provided criticism, feedback and suggestions to the package since its release!
8 Chapter 1. Introduction
2. Discrete Ill-Posed Problems and
their Regularization
In this chapter we give a brief introduction to discrete ill-posed problems, we discuss

some numerical regularization methods, and we introduce several numerical “tools”
such as the singular value decomposition, the discrete Picard condition, and the L-
curve, which are suited for analysis of the discrete ill-posed problems. A more com-
plete treatment of all these aspects is given in [49].
2.1. Discrete Ill-Posed Problems

The concept of ill-posed problems goes back to Hadamard in the beginning of this
century, cf. e.g. [35]. Hadamard essentially defined a problem to be ill-posed if the
solution is not unique or if it is not a continuous function of the data—i.e., if an
arbitrarily small perturbation of the data can cause an arbitrarily large perturbation
of the solution. Hadamard believed that ill-posed problems were “artificial” in that
they would not describe physical systems. He was wrong, though, and today there is
a vast amount of literature on ill-posed problems arising in many areas of science and
engineering, cf. e.g. [15, 16, 17, 33, 66, 71, 73, 79, 89].
The classical example of an ill-posed problem is a Fredholm integral equation of
the first kind with a square integrable kernel [32],
Z b
K(s, t) f (t) dt = g(s) , c≤s≤d, (2.1)
a
where the right-hand side g and the kernel K are given, and where f is the unknown
solution. If the solution f is perturbed by
∆f (t) = ² sin(2πp t) , p = 1, 2, . . . , ² = constant
then the corresponding perturbation of the right-hand side g is given by
Z b
∆g(s) = ² K(s, t) sin(2πp t) dt , p = 1, 2, . . .
a
and due to the Riemann-Lebesgue lemma it follows that ∆g → 0 as p → ∞ [32, p. 2].

Hence, the ratio k∆f k/k∆gk can become arbitrary large by choosing the integer p
10 DISCRETE ILL-POSED PROBLEMS
large enough, thus showing that (2.1) is an ill-posed problem. In particular, this
example illustrates that Fredholm integral equations of the first kind with square
integrable kernels are extremely sensitive to high-frequency perturbations.
Strictly speaking, ill-posed problems must be infinite dimensional—otherwise the
ratio k∆f k/k∆gk stays bounded, although it may become very large. However, cer-
tain finite-dimensional discrete problems have properties very similar to those of ill-
posed problems, such as being highly sensitive to high-frequency perturbations, and
it is natural to associate the term discrete ill-posed problems with these problems. We
can be more precise with this characterization for linear systems of equations
Ax = b , A ∈ Rn×n (2.2)
and linear least-squares problems
min kA x − bk2 , A ∈ Rm×n , m>n. (2.3)

x
We say that these are discrete ill-posed problems if both of the following criteria are
satisfied:
1. the singular values of A decay gradually to zero
2. the ratio between the largest and the smallest nonzero singular values is large.
Singular values are discussed in detail in Section 2.3. Criterion 2 implies that the ma-
trix A is ill-conditioned, i.e., that the solution is potentially very sensitive to pertur-
bations; criterion 1 implies that there is no “nearby” problem with a well-conditioned
coefficient matrix and with well-determined numerical rank.
The typical manifestations of discrete ill-posed problems are systems of linear
equations and linear least-squares problems arising from discretization of ill-posed
problems. E.g., if a Galerkin-type method [3] is used to discretize the Fredholm
integral equation (2.1), then a problem of the form (2.2) or (2.3) arises—depending
on the type of collocation method used—with the elements aij and bi of the matrix
A and the right-hand side b given by
Z bZ d Z d
aij = K(s, t) φi (s) ψj (t) ds dt , bi = φi (s) g(s) ds , (2.4)
a c c
where φi and ψj are the particular basis functions used in the Galerkin method. For
such problems, the close relationship between the ill-posedness of the integral equation
and the ill-conditioning of the matrix A are well understood [1, 41, 88]. In particular,
it can be shown that the singular values of A decay in such a way that both criteria
1 and 2 above are satisfied.
An interesting and important aspect of discrete ill-posed problems is that the ill-
conditioning of the problem does not mean that a meaningful approximate solution
cannot be computed. Rather, the ill-conditioning implies that standard methods in
2.2. Regularization Methods 11
numerical linear algebra [9, 30] for solving (2.2) and (2.3), such as LU, Cholesky,
or QR factorization, cannot be used in a straightforward manner to compute such
a solution. Instead, more sophisticated methods must be applied in order to ensure
the computation of a meaningful solution. This is the essential goal of regularization
methods.
The package Regularization Tools provides a collection of easy-to-use Matlab
routines for the numerical treatment of discrete ill-posed problems. The philosophy
behind Regularization Tools is modularity and regularity between the routines.
Many routines require the SVD of the coefficient matrix A—this is not necessarily
the best approach in a given application, but it is certainly well suited for Matlab [63]
and for this package.
The numerical treatment of integral equations in general is treated in standard
references such as [4, 5, 14, 18, 19], and surveys of regularization theory can be found
in, e.g., [7, 10, 32, 33, 48, 49, 57, 60, 75, 83, 89].
2.2. Regularization Methods

The primary difficulty with the discrete ill-posed problems (2.2) and (2.3) is that
they are essentially underdetermined due to the cluster of small singular values of A.
Hence, it is necessary to incorporate further information about the desired solution
in order to stabilize the problem and to single out a useful and stable solution. This
is the purpose of regularization.
Although many types of additional information about the solution x is possible in
principle, the dominating approach to regularization of discrete ill-posed problems is
to require that the 2-norm—or an appropriate seminorm—of the solution be small.
An initial estimate x∗ of the solution may also be included in the side constraint.
Hence, the side constraint involves minimization of the quantity
Ω(x) = kL (x − x∗ )k2 . (2.5)
Here, the matrix L is typically either the identity matrix In or a p × n discrete
approximation of the (n−p)-th derivative operator, in which case L is a banded matrix
with full row rank. In some cases it is more appropriate that the side constraint be a
Sobolev norm of the form
q
X
Ω(x)2 = α02 kx − x∗ k22 + αi2 kLi (x − x∗ )k22 ,
i=1
where Li approximates the ith derivative operator. Notice that this Ω can always
be written
Pq in the form (2.5) by setting L equal to the Cholesky factor of the matrix
α02 In + i=1 αi2 LTi Li . By means of the side constraint Ω one can therefore control
the smoothness of the regularized solution.
When the side constraint Ω(x) is introduced, one must give up the requirement
that A x equals b in the linear system (2.2) and instead seek a solution that provides
a fair balance between minimizing Ω(x) and minimizing the residual norm kA x −
bk2 . The underlying idea is that a regularized solution with small (semi)norm and
a suitably small residual norm is not too far from the desired, unknown solution to
the unperturbed problem underlying the given problem. The same idea of course also
applies to the least squares problem (2.3).
Undoubtedly, the most common and well-known form of regularization is the one
known as Tikhonov regularization [72, 77, 78]. Here, the idea is to define the reg-
ularized solution xλ as the minimizer of the following weighted combination of the
residual norm and the side constraint
© ª
xλ = argmin kA x − bk22 + λ2 kL (x − x∗ )k22 , (2.6)
where the regularization parameter λ controls the weight given to minimization of

the side constraint relative to minimization of the residual norm. Clearly, a large λ
(equivalent to a large amount of regularization) favors a small solution seminorm at
the cost of a large residual norm, while a small λ (i.e., a small amount of regular-
ization) has the opposite effect. As we shall see in Eq. (2.14), λ also controls the
sensitivity of the regularized solution xλ to perturbations in A and b, and the per-
turbation bound is proportional to λ−1 . Thus, the regularization parameter λ is an
important quantity which controls the properties of the regularized solution, and λ
should therefore be chosen with care. In Section 2.9 we return to numerical methods
for actually computing λ.
We remark that an underlying assumption for the use of Tikhonov regularization
in the form of Eq. (2.6) is that the errors in the right-hand side are unbiased and that
their covariance matrix is proportional to the identity matrix. If the latter condition
is not satisfied one should incorporate the additional information and rescale the
problem or use a regularized version of the general Gauss-Markov linear model:
© ª
min kuk22 + λ2 kL xk22 subject to Ax + C u = b , (2.7)
where C is the Cholesky factor of the covariance matrix. The latter approach using
(2.7) must be used if the covariance matrix is rank deficient, i.e., if C is not a square
matrix. For a discussion of this approach and a numerical algorithm for solving (2.7),
cf. [90].
Besides Tikhonov regularization, there are many other regularization methods
with properties that make them better suited to certain problems or certain comput-
ers. We return to these methods in Sections 2.7 and 2.8, but first it is convenient
to introduce some important numerical “tools” for analysis of discrete ill-posed prob-
lems in Sections 2.3–2.5. As we shall demonstrate, getting insight into the discrete
ill-posed problem is often at least as important as computing a solution, because the
regularized solution should be computed with such care. Finally, in Section 2.9 we
shall describe some methods for choosing the regularization parameter.
2.3. SVD and Generalized SVD 13
2.3. SVD and Generalized SVD

The superior numerical “tools” for analysis of discrete ill-posed problems are the
singular value decomposition (SVD) of A and its generalization to two matrices, the
generalized singular value decomposition (GSVD) of the matrix pair (A, L) [30, §2.5.3
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 the regularization
problem involving both the coefficient matrix A and the regularization matrix L, such
as in (2.6).
2.3.1. The Singular Value Decomposition

Let A ∈ Rm×n be a rectangular matrix with m ≥ n. Then the SVD of A is a
decomposition of the form
n
X
T
A = U ΣV = ui σi viT , (2.8)
i=1
where U = (u1 , . . . , un ) and V = (v1 , . . . , vn ) are matrices with orthonormal columns,

U T U = V T V = In , and where Σ = diag(σ1 , . . . , σn ) has non-negative diagonal ele-
ments appearing in non-increasing order such that
σ1 ≥ . . . ≥ σn ≥ 0 . (2.9)
The numbers σi are the singular values of A while the vectors ui and vi are the left
and right singular vectors of A, respectively. The condition number of A is equal to
the ratio σ1 /σn .
From the relations AT A = V Σ2 V T and A AT = U Σ2 U T we see that the SVD
of A is strongly linked to the eigenvalue decompositions of the symmetric positive
semi-definite matrices AT A and A AT . This shows that the SVD is unique for a given
matrix A—except for singular vectors associated with multiple singular values. In
connection with discrete ill-posed problems, two characteristic features of the SVD of
A are very often found.
• The singular values σi decay gradually to zero with no particular gap in the
spectrum. An increase of the dimensions of A will increase the number of small
singular values.
• The left and right singular vectors ui and vi tend to have more sign changes in
their elements as the index i increases, i.e., as σi decreases.
Although these features are found in many discrete ill-posed problems arising in prac-
tical applications, they are unfortunately very difficult—or perhaps impossible—to
prove in general.
To see how the SVD gives insight into the ill-conditioning of A, consider the
following relations which follow directly from Eq. (2.8):
¾
A vi = σi ui
i = 1, . . . , n . (2.10)
kA vi k2 = σi
We see that a small singular value σi , compared to kAk2 = σ1 , means that there
exists a certain linear combination of the columns of A, characterized by the elements
of the right singular vector vi , such that kA vi k2 = σi is small. In other words, one
or more small σi implies that A is nearly rank deficient, and the vectors vi associated
with the small σi are numerical null-vectors of A. From this and the characteristic
features of A we conclude that the matrix in a discrete ill-posed problem is always
highly ill-conditioned, and its numerical null-space is spanned by vectors with many
sign changes.
The SVD also gives important insight into another aspect of discrete ill-posed
problems, namely the smoothing effect typically associated with a square integrable
kernel. Notice that as σi decreases, the singular vectors ui and vi become more and
more oscillatory. Consider
Pn now the mapping A x of an arbitrary vector x. Using the
SVD, we get x = i=1 (viT x) vi and
n
X
Ax = σi (viT x) ui .
i=1
This clearly shows that the due to the multiplication with the σi the high-frequency
components of x are more damped in A x than then low-frequency components. More-
over, the inverse problem, namely that of computing x from A x = b or min kA x−bk2 ,
must have the opposite effect: it amplifies the high-frequency oscillations in the right-
hand side b.
2.3.2. The Generalized Singular Value Decomposition

The GSVD of the matrix pair (A, L) is a generalization of the SVD of A in the sense
that the generalized singular values of (A, L) are the square roots of the generalized
eigenvalues of the matrix pair (AT A, LT L). In order to keep our exposition simple,
we assume that the dimensions of A ∈ Rm×n and L ∈ Rp×n satisfy m ≥ n ≥ p, which
is always the case in connection with discrete ill-posed problems. Then the GSVD is
a decomposition of A and L in the form
µ ¶
Σ 0
A=U X −1 , L = V (M , 0) X −1 , (2.11)
0 In−p
where the columns of U ∈ Rm×n and V ∈ Rp×p are orthonormal, X ∈ Rn×n is

nonsingular, and Σ and M are p × p diagonal matrices:
Σ = diag(σ1 , . . . , σp ), M = diag(µ1 , . . . , µp ).
2.3. SVD and Generalized SVD 15
Moreover, the diagonal entries of Σ and M are non-negative and ordered such that
0 ≤ σ1 ≤ . . . ≤ σp ≤ 1 , 1 ≥ µ1 ≥ . . . ≥ µp > 0 ,
and they are normalized such that
σi2 + µ2i = 1 , i = 1, . . . , p .
Then the generalized singular values γi of (A, L) are defined as the ratios
γi = σi /µi , i = 1, . . . , p , (2.12)
and they obviously appear in non-decreasing order. For historical reasons, this order-
ing is the opposite of the ordering of the ordinary singular values of A.
For p < n the matrix L ∈ Rp×n always has a nontrivial null-space N (L). E.g.,
if L is an approximation to the second derivative operator on a regular mesh, L =
tridiag(1, −2, 1), then N (L) is spanned by the two vectors (1, 1, . . . , 1)T and (1, 2, . . . , n)T .
In the GSVD, the last n − p columns xi of the nonsingular matrix X satisfy
L xi = 0 , i = p + 1, . . . , n (2.13)
and they are therefore basis vectors for the null-space N (L).
There is a slight notational problem here because the matrices U , Σ, and V in the
GSVD of (A, L) are different from the matrices with the same symbols in the SVD
of A. However, in this presentation it will always be clear from the context which
decomposition is used. When L is the identity matrix In , then the U and V of the
GSVD are identical to the U and V of the SVD, and the generalized singular values
of (A, In ) are identical to the singular values of A—except for the ordering of the
singular values and vectors.
In general, there is no connection between the generalized singular values/vectors
and the ordinary singular values/vectors. For discrete ill-posed problems, though, we
can actually say something about the SVD-GSVD connection because L is typically
a reasonably well-conditioned matrix. When this is the case, then it can be shown
that the matrix X in (2.11) is also well-conditioned. Hence, the diagonal matrix Σ
must display the ill-conditioning of A, and since γi = σi (1 − σi2 )1/2 ≈ σi for small σi
the generalized singular values must decay gradually to zero as the ordinary singular
values do. Moreover, the oscillation properties (i.e., the increase in sign changes) of
the right singular vectors carries over to the columns of X in the GSVD: the smaller
the γi the more sign changes in xi . For more specific results, cf. [43].
As an immediate example of the use of GSVD in the analysis of discrete regular-
ization problems, we mention the following perturbation bound for Tikhonov regular-
ization derived in [42]. Let E and e denote the perturbations of A and b, respectively,
and let x̄λ denote the exact solution to the unperturbed problem; then the relative
error in the perturbed solution xλ satisfies

kxλ − x̄λ k2 kAk2 kXk2 λ−1
≤ ×
kx̄λ k2 1 − kEk2 kXk2 λ−1
µ ¶
kEk2 kek2 kEk2 kXk2 krλ k2
(1 + cond(X)) + + (2.14)
kAk2 kbλ k2 λ kbλ k2
where we have defined bλ = A xλ and rλ = b − bλ . The important conclusion we
can make from this relation is that for all reasonable λ the perturbation bound for
the regularized solution xλ is proportional to λ−1 and to the norm of the matrix X.
The latter quantity is analyzed in [43] where it is shown that kXk2 is approximately
bounded by kL† k2 , i.e., by the inverse of the smallest singular value of L. Hence,
in addition to controlling the smoothness of the regularized solution, λ and L also
control its sensitivity to perturbations of A and b.
The SVD and the GSVD are computed by means of routines csvd and cgsvd in
this package.
2.4. The Discrete Picard Condition and Filter Factors

As we have seen in Section 2.3, the integration in Eq. (2.1) with a square integrable
kernel K (2.1) has a smoothing effect on f . The opposite operation, namely, that
of solving the first kind Fredholm integral equation for f , therefore tends to amplify
oscillations in the right-hand side g. Hence, if we require that the solution f be a
square integrable solution with finite L2 -norm, then not all functions are valid as
right-hand side g. Indeed, g must be sufficiently smooth to “survive” the inversion
back to f . The mathematical formulation of this smoothness criterion on g—once the
kernel K is given—is called the Picard condition [32, §1.2].
For discrete ill-posed problems there is, strictly speaking, no Picard condition
because the norm of the solution is always bounded. Nevertheless, it makes sense to
introduce a discrete Picard condition as follows. In a real-world application, the right-
hand side b is always contaminated by various types or errors, such as measurement
errors, approximation errors, and rounding errors. Hence, b can be written as
b = b̄ + e , (2.15)
where e are the errors, and b̄ is the unperturbed right-hand side. Both b̄ and the
corresponding unperturbed solution x̄ represent the underlying unperturbed and un-
known problem. Now, if we want to be able to compute a regularized solution xreg
from the given right-hand side b such that xreg approximates the exact solution x̄,
then it is shown in [46] that the corresponding exact right-hand side b̄ must satisfy a
criterion very similar to the Picard condition:
The discrete Picard condition. The unperturbed right-hand side b̄ in a discrete

ill-posed problem with regularization matrix L satisfies the discrete Picard condition if
2.4. The Discrete Picard Condition and Filter Factors 17
the Fourier coefficients |uTi b̄| on the average decay to zero faster than the generalized
singular values γi .
The discrete Picard condition is not as “artificial” as it first may seem: it can be
shown that if the underlying integral equation (2.1) satisfies the Picard condition,
then the discrete ill-posed problem obtained by discretization of the integral equation
satisfies the discrete Picard condition [41]. See also [83, 84].
The main difficulty with discrete ill-posed problems is caused by the errors e in the
given right-hand side b (2.15), because such errors typically tend to have components
along all the left singular vectors ui . For example, if kek2 = ² and if the elements of
e are unbiased and uncorrelated, then the expected value of the Fourier coefficients
of e satisfy ¡ ¢ 1
E |uTi e| = m− 2 ² , i = 1, . . . , n . (2.16)
As a consequence, the Fourier coefficients |uTi b| of the perturbed right-hand side level
off at approximately m−1/2 ² even if the unperturbed right-hand side b̄ satisfies the
discrete Picard condition, because these Fourier coefficients are dominated by |uTi e|
for large i.
Consider now the linear system (2.2) and the least squares problem (2.3), and
assume for simplicity that A has no exact zero singular values. Using the SVD, it is
easy to show that the solutions to both systems are given by the same equation:
Xn
uTi b
xLSQ = vi . (2.17)
i=1
σi
This relation clearly illustrates the difficulties with the standard solution to (2.2)
and (2.3). Since the Fourier coefficients |uTi b| corresponding to the smaller singular
values σi do not decay as fast as the singular values—but rather tend to level off—the
solution xLSQ is dominated by the terms in the sum corresponding to the smallest
σi . As a consequence, the solution xLSQ has many sign changes and thus appears
completely random.
With this analysis in mind, we can see that the purpose of a regularization method
is to dampen or filter out the contributions to the solution corresponding to the
small generalized singular values. Hence, we will require that a regularization method
produces a regularized solution xreg which, for x∗ = 0, can be written as follows
n
X uTi b
xreg = fi vi if L = In (2.18)
i=1
σi
Xp Xn
uTi b
xreg = fi xi + (uTi b) xi if L 6= In . (2.19)
i=1
σi i=p+1
Here, the numbers fi are filter factors for the particular regularization method. The
filter factors must have the important property that as σi decreases, the correspond-
ing fi tends to zero in such a way that the contributions (uTi b/σi ) xi to the solution
from the smaller σi are effectively filtered out. The difference between the various reg-
ularization methods lies essentially in the way that these filter factors fi are defined.
Hence, the filter factors play an important role in connection with regularization the-
ory, and it is worthwhile to characterize the filter factors for the various regularization
methods that we shall present below.
For Tikhonov regularization, which plays a central role in regularization theory,
the filter factors are either fi = σi2 /(σi2 + λ2 ) (for L = In ) or fi = γi2 /(γi2 + λ2 ) (for
L 6= In ), and the filtering effectively sets in for σi < λ and γi < λ, respectively. In
particular, this shows that discrete ill-posed problems are essentially unregularized by
Tikhonov’s method for λ > σ1 and λ > γp , respectively.
Filter factors for various regularization methods can be computed by means of
routine fil fac in this package, while routine picard plots the important quantities σi ,
|uTi b|, and |uTi b|/σi if L = In , or γi , |uTi b|, and |uTi b|/γi if L 6= In .
2.5. The L-Curve
Perhaps the most convenient graphical tool for analysis of discrete ill-posed problems
is the so-called L-curve which is a plot—for all valid regularization parameters—of
the (semi)norm kL xreg k2 of the regularized solution versus the corresponding residual
norm kA xreg − bk2 . In this way, the L-curve clearly displays the compromise between
minimization of these two quantities, which is the heart of any regularization method.
The use of such plots in connection with ill-conditioned least squares problems goes
back to Miller [64] and Lawson & Hanson [61].
For discrete ill-posed problems it turns out that the L-curve, when plotted in log-
log scale, almost always has a characteristic L-shaped appearance (hence its name)
with a distinct corner separating the vertical and the horizontal parts of the curve.
To see why this is so, we notice that if x̄ denotes the exact, unregularized solution
corresponding to the exact right-hand side b̄ in Eq. (2.15), then the error xreg − x̄
in the regularized solution consists of two components, namely, a perturbation error
from the error e in the given right-hand side b, and a regularization error due to the
regularization of the error-free component b̄ in the right-hand side. The vertical part
of the L-curve corresponds to solutions where kL xreg k2 is very sensitive to changes
in the regularization parameter because the perturbation error e from dominates xreg
and because e does not satisfy the discrete Picard condition. The horizontal part of
the L-curve corresponds to solutions where it is the residual norm kA xreg − bk2 that
is most sensitive to the regularization parameter because xreg is dominated by the
regularization error—as long as b̄ satisfies the discrete Picard condition.
We can substantiate this by means of the relations for the regularized solution
xreg in terms of the filter factors. For general-form regularization (L 6= In ) Eq. (2.19)
2.5. The L-Curve 19
less filtering
log || L x || 2
more filtering 1
-------------------------------------------_/
-- --- - -,"'------:.:::---,
log || A x – b || 2
Figure 2.1: The generic form of the L-curve.
yields the following expression for the error in xreg :

 
X p T Xn Xp
u e uT b̄
xreg − x̄ =  fi i xi + (uTi e) xi  + (fi − 1) i xi . (2.20)
i=1
σi i=p+1 i=1
σi
Here, the term in the parenthesis is the perturbation error due to the perturbation
e, and the second term is the regularization error caused by regularization of the
unperturbed component b̄ of the right-hand side. When only little regularization is
introduced, most of the filter factors fi are approximately one and the error xreg − x̄ is
dominated by the perturbation error. On the other hand, with plenty of regularization
most filter factors are small, fi ¿ 1, and xreg − x̄ is dominated by the regularization
error.
In [47, 55] Eq. (2.20) was used to analyze the relationship between the error in xreg
and the behavior of the L-curve. The result is that if b̄ satisfies the discrete Picard
condition, then the horizontal part of the curve corresponds to solutions where the
regularization error dominates—i.e., where so much filtering is introduced that the
solution stays very smooth and kL xreg k2 therefore only changes a little with the
regularization parameter. In contrast, the vertical part of the L-curve corresponds to
solutions that are dominated by the perturbation error, and due to the division by the
small σi it is clear that kL xreg k2 varies dramatically with the regularization parameter
while, simultaneously, the residual norm does not change much. Moreover, it is shown
in [55] that a log-log scale emphasizes the different appearances of the vertical and
the horizontal parts. In this way, the L-curve clearly displays the trade-off between
minimizing the residual norm and the side constraint.
We note in passing that the L-curve is a continuous curve when the regularization
parameter is continuous as in Tikhonov regularization. For regularization methods
with a discrete regularization parameter, such as truncated SVD, we plot the L-curve
as a finite set of points. How to make such a “discrete L-curve” continuous is discussed
in [55]. In this reference, alternative norms for plotting the L-curve, depending on
the regularization method, are also discussed.
The L-curve for Tikhonov regularization plays a central role in connection with
regularization methods for discrete ill-posed problems because it divides the first
quadrant into two regions. It is impossible to construct any solution that corresponds
to a point below the Tikhonov L-curve; any regularized solution must lie on or above
this curve. The solution computed by Tikhonov regularization is therefore optimal in
the sense that for a given residual norm there does not exist a solution with smaller
seminorm than the Tikhonov solution—and the same is true with the roles of the
norms interchanged. A consequence of this is that one can compare other regulariza-
tion methods with Tikhonov regularization by inspecting how close the L-curve for
the alternative method is to the Tikhonov L-curve. If b̄ satisfies the discrete Picard
condition, then the two L-curves are close to each other and the solutions computed
by the two regularization methods are similar [47].
For a given fixed right-hand side b = b̄ + e, there is obviously an optimal reg-

ularization parameter that balances the perturbation error and the regularization
error in xreg . An essential feature of the L-curve is that this optimal regularization
parameter—defined in the above sense—is not far from the regularization parameter
that corresponds to the L-curve’s corner [47]. In other words, by locating the cor-
ner of the L-curve one can compute an approximation to the optimal regularization
parameter and thus, in turn, compute a regularized solution with a good balance
between the two error types. We return to this aspect in Section 2.9; suffice it here
to say that for continuous L-curves, a computationally convenient definition of the
L-curve’s corner is the point with maximum curvature in log-log scale.
In the Regularization Tools package, routine l curve produces a log-log plot of

the L-curve and—if required—also locates the corner and identifies the corresponding
regularization parameter. Given a discrete set of values of kA xreg −bk2 and kL xreg k2 ,
routine plot lc plots the corresponding L-curve, while routine l corner locates the L-
curve’s corner.
2.6. Transformation to Standard Form 21
2.6. Transformation to Standard Form

A regularization problem with side constraint Ω(x) = kL (x − x∗ )k2 (2.5) is said to
be in standard form if the matrix L is the identity matrix In . In many applications,
regularization in standard form is not the best choice, i.e., one should use some L 6= In
in the side constraint Ω(x). The proper choice of matrix L depends on the particular
application, but often an approximation to the first or second derivative operator
gives good results.
However, from a numerical point of view it is much simpler to treat problems
in standard form, basically because only one matrix, A, is involved instead of the
two matrices A and L. Hence, it is convenient to be able to transform a given
regularization problem in general form into an equivalent one in standard form by
means of numerically stable methods. For example, for Tikhonov regularization we
want a numerically stable method for transforming the general-form problem (2.6)
into the following standard-form problem
© ª
min kĀ x̄ − b̄k22 + λ2 kx̄ − x̄∗ k22 , (2.21)
where the new matrix Ā, the new right-hand side b̄, and the vector x̄∗ are derived from
the original quantities A, L, b, and x∗ . Moreover, we also want a numerically stable
scheme for transforming the solution x̄λ to (2.21) back to the general-form setting.
Finally, we prefer a transformation that leads to a simple relationship between the
SVD of Ā and the GSVD of (A, L), for then we have a perfect understanding of the
relationship between the two regularization problems.
For the simple case where L is square and invertible, the transformation is obvious:
Ā = A L−1 , b̄ = b, x̄∗ = L x∗ , and the back-transformation simply becomes xλ =
L−1 x̄λ .
In most applications, however, the matrix L is not square, and the transformation
becomes somewhat more involved than just a matrix inversion. It turns out that it is a
good idea to distinguish between direct and iterative regularization methods—cf. the
next two sections. For the direct methods we need to be able to compute the matrix
Ā explicitly by standard methods such as the QR factorization. For the iterative
methods, on the other hand, we merely need to be able to compute the matrix-
vector product Ā x̄ efficiently. Below, we describe two methods for transformation
to standard form which are suited for direct and iterative methods, respectively. We
assume that the matrix L ∈ Rp×n has full row rank, i.e., the rank of L is p.
2.6.1. Transformation for Direct Methods

The standard-form transformation for direct methods described here was developed
by Eldén [22], and it is based on two QR factorizations. The description of this trans-
formation is quite algorithmic, and it is summarized below (where, for convenience,
the subscripts p, o, and q denote matrices with p, n − p, and m − (n − p) columns,
respectively). First, compute a QR factorization of LT ,

µ ¶
T Rp
L = K R = (Kp Ko ) . (2.22)
0
We remark that since L has full rank, its pseudoinverse is simply L† = Kp Rp−T .
Moreover, the columns of Ko are an orthonormal basis for the null space of L. Next,
form the “skinny” matrix A Ko ∈ Rm×(n−p) and compute its QR factorization,
µ ¶
To
A Ko = H T = (Ho Hq ) . (2.23)
0
Then the transformed quantities Ā and b̄ are given by the following identities
Ā = HqT A L† = HqT A Kp Rp−T , b̄ = HqT b , (2.24)
and we stress that the most efficient way to compute Ā and b̄ is to apply the or-
thogonal transformations that make up K and H “on the fly” to A and b as the QR
factorizations in (2.22) and (2.23) are computed. When (2.21) has been solved for x̄,
the transformation back to the general-form setting then takes the form
x = L† x̄ + Ko To−1 HoT (b − A L† x̄) . (2.25)
The SVD of Ā is related to the GSVD of (A, L) as follows: let Ā = Ū Σ̄ V̄ T denote

the SVD of Ā, and let Ep denote the p × p exchange matrix Ep = antidiag(1, . . . , 1).
Also, let U , V , Σ, M , and X denote the GSVD matrices from Eq. (2.11). Then
U = (Hq Ū Ep , Ho ) , Σ M −1 = Ep Σ̄ Ep (2.26)
µ −1 T ¶−1
M V L
V = V̄ Ep , X = . (2.27)
HoT A
Moreover, the last n−p columns of X are given by Ko To−1 . For proofs of (2.26)–(2.27)
and an investigation of the accuracy of the GSVD matrices computed this way, cf.
[43, 45].
2.6.2. Transformation for Iterative Methods

For the iterative methods the matrix Ā is never computed explicitly. Instead, one
merely needs to be able to pre-multiply a vector with Ā and ĀT efficiently. If Ko is
an orthonormal basis for the null space of L, e.g. computed by (2.22), then y = L† x̄
and ŷ = (L† )T x, both of which are used in the iterative procedures, can easily by
computed by the following algorithms:
µ ¶−1 µ ¶
In−p 0 0 x ← x − Ko KoT x
y ← µ ¶ µ ¶−T
L x̄ ŷ L (2.28)
← x.
y ← y − Ko KoT y z 0 In−p
2.6. Transformation to Standard Form 23
By means of these simple algorithms, which are described in [8], the above standard-
transformation method can also be used for iterative methods. Notice that a basis for
the null space of L is required—often, the basis vectors can be computed explicitly,
or they can be computed from L by, say, a rank revealing factorization [13].
The algorithm from §2.6.1 can be reformulated in such a way that the pseudoin-
verse L† is replaced by a weaker generalized inverse, using an idea from [24] (and
later advocated in [36, 37, 39]). This reformulation has certain advantages for itera-
tive methods, as we shall see in §2.8.4. Define the A-weighted generalized inverse of
L as follows µ −1 ¶
† M
LA = X VT , (2.29)
0
where we emphasize that L†A is generally different from the pseudoinverse L† when
p < n. Also, define the vector
n
X
x0 = (uTi b) xi , (2.30)
i=p+1
which is the unregularized component of the regularized solution xreg , cf. Eq. (2.19),
i.e., x0 is the component of xreg that lies in the null space of L. Then the standard-
form quantities Ā and b̄ in the alternative version of the algorithm are defined as
follows
Ā = A L†A , b̄ = b − A x0 . (2.31)
Moreover, the transformation back to the general-form setting takes the simple form
x = L†A x̄ + x0 . (2.32)
This backtransformation is mathematically equivalent to the one in (2.25) since we

have
L†A = (In − Ko To−1 HoT A) L† and x0 = Ko To−1 HoT b . (2.33)
To use this alternative formulation of the standard-form transformation, we need
to compute x0 as well as the matrix-vector products L†A x̄ and (L†A )T x efficiently.
Given a basis W for N (L), the vector x0 is computed by the following algorithm:
S ← (A W )†
(2.34)
x0 ← W Sb .
This algorithm involves O(mn(n − p)) operations. To compute L†A x̄ and (L†A )T x
efficiently (we emphasize that L†A is never computed explicitly), we partition L =
(L11 , L12 ), T = (T11 , T12 ) and xT = (xT1 , xT2 ) such that L11 and T11 have p columns
and x1 has p elements. For efficiency, we also need to compute the (n − p) × n matrix
T ←SA . (2.35)
Then y = L†A x̄ and ŷ = (L†A )T x are computed by:
y L−1
← µ x̄
11 ¶ T
x ← x − T11 WTx
y −T (2.36)
y ← − W T11 y ŷ ← L11 x .
0
In the above formulas, W need not have orthonormal columns, although this is the
best choice from a numerical point of view. For more details about these algorithms,
cf. [49, Section 2.3.2].
For the latter standard-form transformation there is an even simpler relation
between the SVD of Ā and part of the GSVD of (A, L) than in §2.6.1 because
Pp
A L†A = T
i=i ui γi vi . I.e., except for the ordering the GSVD quantities ui , γi ,
and vi are identical to the similar SVD quantities, and with the same notation as in
Eq. (2.26) we have
(u1 . . . up ) = Ū Ep , V = V̄ Ep , Σ M −1 = Ep Σ̄ Ep . (2.37)
This relation is very important in connection with the iterative regularization meth-
ods.
2.6.3. Norm Relations etc.

From the above relations (2.26) and (2.37) between the SVD of Ā and the GSVD of
(A, L) we obtain the following very important relations between the norms related to
the two regularization problems
kL (x − x∗ )k2 = kx̄ − x̄∗ k2 , kA x − bk2 = kĀ x̄ − b̄k2 , (2.38)
where x denotes the solution obtained by transforming x̄ back to the general-form

setting. These relations are very important in connection with methods for choosing
the regularization parameter because they show that any parameter-choice strategy
based on these norms will yield the same regularization parameter when applied to
the general-form problem and the transformed standard-form problem.
Several routines are included in Regularization Tools for computations related
to the standard-form transformations. First of all, the routines std form and gen form
perform both transformations to standard from and back to general form. These
routines are mainly included for pedagogical reasons. Routine pinit computes the
vector x0 in Eq. (2.30) as well as the matrix T A by means of Algorithm (2.34).
Finally, the two routines lsolve and ltsolve compute L†A x̄ and (L†A )T x by the algorithms
in (2.36).
Regarding the matrix L, discrete approximations to derivative operators on a
regular mesh can be computed by routine get l which also provides a matrix W with
orthonormal basis vectors for the null space of L.
2.7. Direct Regularization Methods 25
2.7. Direct Regularization Methods

In this and the next section we shall briefly review the regularization methods for
numerical treatment of discrete ill-posed problems included in the Regularization
Tools package. Our aim is not to compare these and other methods, because that
is outside the scope of this paper. In fact, very little has been done in this area,
cf. [2, 39, 44, 62]. This section focuses on methods which are essentially direct,
i.e., methods where the solution is defined by a direct computation (which may still
involve an iterative root-finding procedure, say), while regularization methods which
are intrinsically iterative are treated in the next section.
2.7.1. Tikhonov Regularization

Tikhonov’s method is of course a direct method because the regularized solution xλ ,
defined by Eq. (2.6), is the solution to the following least squares problem
°µ ¶ µ ¶°
° A b °
min °° λL x − °
∗ ° , (2.39)
λL x 2
and it is easy to see that xλ is unique if the null spaces of A and L intersect trivially
(as they usually do in practice). The most efficient algorithm for numerical treatment
of Tikhonov’s method for a general regularization matrix L consists of three steps [22].
First, the problem is transformed into standard form by means of Eqs. (2.22)–(2.24)
from §2.6.1 (Ā = A if L = In ), then the matrix Ā is transformed into a p × p upper
bidiagonal matrix B̄ by means of left and right orthogonal transformations,
Ā = Ū B̄ V̄ T ,
and finally the resulting sparse problem with a banded B̄ is solved for V̄ T x̄λ and the
solution is transformed back to the original setting by means of (2.25).
In this package we take another approach to solving (2.39), namely, by using the
filter factors and the GSVD explicitly (or the SVD, if L = In ), cf. Eqs. (2.18) and
(2.19). This approach, which is implemented in routine tikhonov, is more suited to
Matlab’s coarse granularity. For pedagogical reasons, we also include a routine bidiag
for bidiagonalization of a matrix.
2.7.2. Least Squares with a Quadratic Constraint

There are two other regularization methods which are almost equivalent to Tikhonov’s
method, and which can be treated numerically by essentially the same technique as
mentioned above involving a transformation to standard form followed by bidiagonal-
ization of the coefficient matrix. These two methods are the following least squares
problems with a quadratic inequality constraint
min kA x − bk2 subject to kL (x − x∗ )k2 ≤ α (2.40)
∗
min kL (x − x )k2 subject to kA x − bk2 ≤ δ , (2.41)
where α and δ are nonzero parameters each playing the role as regularization pa-
rameter in (2.40) and (2.41), respectively. The solution to both these problems is
identical to xλ from Tikhonov’s method for suitably chosen values of λ that depend
in a nonlinear way on α and δ. The solution to the first problem (2.40) is computed
as follows: if kL (xLSQ − x0 )k2 ≤ α then λ ← 0 and xλ ← xLSQ , else use an iterative
scheme to solve
min kA xλ − bk2 subject to kL (xλ − x∗ )k2 = α
for λ and xλ . Similarly, the solution to the second problem (2.41) is computed as
follows (where x0 is given by Eq. (2.30)): if kA x0 − bk2 ≤ δ then λ ← ∞ and
xλ ← x0 , else use an iterative scheme to solve
min kL (xλ − x∗ )k2 subject to kA xλ − bk2 = δ
for λ and xλ . In Regularization Tools, routines lsqi and discrep solve (2.40)
and (2.41), respectively. The name “discrep” is related to the discrepancy principle
for choosing the regularization parameter, cf. Section 2.9. An efficient algorithm for
solving (2.40) when A is large and sparse, based on Gauss quadrature and Lanczos
bidiagonalization, is described in [31].
2.7.3. TSVD, MTSVD, and TGSVD

A fundamental observation regarding the abovementioned methods is that they cir-
cumvent the ill conditioning afµA by¶introducing a new problem (2.39) with a well-
A
conditioned coefficient matrix with full rank. A different way to treat the
λL
ill-conditioning of A is to derive a new problem with a well-conditioned rank defi-
cient coefficient matrix. A fundamental result about rank deficient matrices, which
can be derived from the SVD of A, is that the closest rank-k approximation Ak to
A—measured in the 2-norm—is obtained by truncating the SVD expansion in (2.8)
at k, i.e., Ak is given by
k
X
Ak = ui σi viT , k≤n. (2.42)
i=1
The truncated SVD (TSVD) [82, 40, 44] and the modified TSVD (MTSVD) [56]
regularization methods are based on this observation in that one solves the problems
min kxk2 subject to min kAk x − bk2 (2.43)
min kL xk2 subject to min kAk x − bk2 , (2.44)
where Ak is the rank-k matrix in Eq. (2.42). The solutions to these two problems are
given by
Xk
uTi b
xk = vi (2.45)
i=1
σi
2.7. Direct Regularization Methods 27
xL,k = xk − Vk (L Vk )† L xk , (2.46)
where (L Vk )† is the pseudoinverse of L Vk , and Vk ≡ (vk+1 , . . . , vn ). In other words,
the correction to xk in (2.46) is the solution to the following least squares problem
min k(L Vk ) z − L xk k2 .
We note in passing that the TSVD solution xk is the only regularized solution which
has no component in the numerical null-space of A, spanned by the columns of Vk .
All other regularized solutions, exemplified by the MTSVD solution xL,k , has some
component in A’s numerical null space in order to achieve the desired properties of
the solution, as controlled by the matrix L.
As an alternative to the abovementioned MTSVD method for general-form prob-
lems one can generalize the TSVD method to the GSVD setting [43, 46]. The result-
ing method, truncated GSVD (TGSVD), is easiest to introduce via the standard-form
transformation from §2.6.2 with Ā = A L†A , b̄ = b − A x0 , and x = L†A x̄ + x0 ⇒ L x =
x̄. In analogy with the TSVD method we now introduce a rank-k approximation Āk to
Ā via its SVD. Due to the SVD-GSVD relations between Ā and (A, PpL), computation
of the matrix Āk is essentially a “truncated GSVD” because Āk = i=p−k+1 ui γi viT .
Then we define the truncated GSVD (TGSVD) solution as x̂L,k = L†A x̄k + x0 , where
x̄k solves the problem
min kx̄k2 subject to min kĀk x̄ − b̄k2 . (2.47)
Definition (2.47) together with the GSVD of (A, L) then immediately lead to the
following simple expression of the TGSVD solution
p
X Xn
uTi b
x̂k,L = xi + (uTi b) xi , (2.48)
σi i=p+1
i=p−k+1
where the last term is the component x0 (2.30) in the null space of L. Defined this
way, the TGSVD solution is a natural generalization of the TSVD solution xk . The
TGSVD method is also a generalization of TSVD because both xk and x̂k,L can be
derived from the corresponding Tikhonov solutions (2.18) and (2.19) by substituting
0’s and 1’s for the Tikhonov filter factors fi .
The TSVD, MTSVD, and TGSVD solutions are computed by the routines with
the obvious names tsvd, mtsvd, and tgsvd.
2.7.4. Damped SVD/GSVD

A less know regularization method which is based on the SVD or the GSVD is the
damped SVD/GSVD (damped SVD was introduced in [21], and our generalization to
damped GSVD is obvious). Here, instead of using filter factors 0 and 1 as in TSVD
and TGSVD, one introduces a smoother cut-off by means of filter factors fi defined
as
σi σi
fi = (for L = In ) and fi = (for L 6= In ) . (2.49)
σi + λ σi + λ µi
These filter factors decay slower than the Tikhonov filter factors and thus, in a sense,
introduce less filtering. The damped SVD/GSVD solutions are computed by means
of routine dsvd.
2.7.5. Maximum Entropy Regularization

This regularization method is frequently used in image reconstruction and related
applications where a solution with positive elements is sought. In maximum entropy
regularization, the following nonlinear function is used as side constraint:
n
X
Ω(x) = xi log(wi xi ) , (2.50)
i=1
where xi are the positive elements of the vector x, and w1 , . . . , wn are n weights.
Notice that −Ω(x) measures the entropy of x, hence the name of this regularization
method. The mathematical justification for this particular choice of Ω(x) is that it
yields a solution x which is most objective, or maximally uncommitted, with respect
to missing information in the right-hand side, cf. e.g. [74].
Maximum entropy regularization is implemented in Regularization Tools in
the routine maxent which uses a nonlinear conjugate gradient algorithm [29, §4.1]
with inexact line search to compute the regularized solution. The typical step in this
method has the form
x(k+1) ← x(k) + αk p(k)
(2.51)
p(k+1) ← −∇F (x(k+1) ) + βk p(k)
in which F is the function to be minimized,
n
X
F (x) = kA x − bk22 + λ2 xi log(wi xi ) ,
i=1
and F ’s gradient is given by

 
1 + log(w1 x1 )
 .. 
∇F (x) = 2 AT (A x − b) + λ2  .  .
1 + log(wn xn )
In Algorithm (2.51), the step-length parameter αk miminizes F (x(k) + αk p(k) ) with
the constraint that all elements of x(k) + αk p(k) be positive, and it is computed by
means of an inexact line search. Then βk is computed by
βk = (∇F (x(k+1) ) − ∇F (x(k) ))T ∇F (x(k+1) ) / k∇F (x(k+1) )k22 .
2.8. Iterative Regularization Methods 29
This choice of βk has the potential advantage that it gives “automatic” restart to the
steepest descent direction in case of slow convergence.
2.7.6. Truncated Total Least Squares

The last direct regularization method included in Regularization Tools is trun-
cated total least squares (TTLS). For rank deficient matrices, total least squares [27]
takes its basis in an SVD of the compound matrix (A b) = Ũ Σ̃ Ṽ T with the matrix
Ṽ ∈ R(n+1)×(n+1) partitioned such that
µ ¶
Ṽ11 Ṽ12
Ṽ = , Ṽ11 ∈ Rn×k , Ṽ22 ∈ R1×(n+1−k) , (2.52)
Ṽ21 Ṽ22
where k is the numerical rank of A. Then the TLS solution is given by

† T
x̃k = −Ṽ12 Ṽ22 = −Ṽ12 Ṽ22 / kṼ22 k22 . (2.53)
The TLS solution is robust to perturbations of A because inaccuracies in A are explic-

itly taken into account in the TLS method. Therefore, for discrete ill-posed problems
with no gap in the singular value spectrum of A, it makes sense to define a truncated
TLS solution by means of (2.53) where k then plays the role of the regularization
parameter; see [28] for more details. The truncated TLS solution is computed by
means of routine ttls.
2.8. Iterative Regularization Methods

This section describes the iterative regularization methods included in Regulariza-
tion Tools. We stress that our Matlab routines should be considered as model
implementations; real implementations should incorporate any sparsity and/or struc-
ture of the matrix A. We shall first describe standard-form versions of the methods
and then describe the extension necessary for treating general-form problems. For
more details about these and other iterative methods, cf. [39, Chapter 6–7].
2.8.1. Conjugate Gradients and LSQR

The conjugate gradient (CG) algorithm is a well-known method for solving sparse sys-
tems of equations with a symmetric positive definite coefficient matrix. In connection
with discrete ill-posed problems, it is an interesting fact that when the CG algorithm
is applied to the unregularized normal equations AT A x = AT b (implemented such
that AT A is not formed) then the low-frequency components of the solution tend
to converge faster than the high-frequency components. Hence, the CG process has
some inherent regularization effect where the number of iterations plays the role of
the regularization parameter. The kth step of the CG process essentially has the form
βk ← kq(k−1) k22 /kq(k−2) k22

p(k) ← q(k−1) + βk p(k−1)
αk ← kq(k−1) k22 /kAT A p(k) k22 (2.54)
x(k) ← x(k−1) + αk p(k)
q(k) ← q(k−1) − αk AT A p(k)
where x(k) is the approximation to x after k iterations, while p(k) and q(k) are two
auxiliary iteration vectors of length n.
To explain this regularizing effect of the CG method, we introduce the Krylov
subspace
Kk (AT A, AT b) = span{AT b, AT A AT b, . . . , (AT A)k−1 AT b}
associated with the kth step of the CG algorithm applied to AT A x = AT b. It is also

convenient to introduce the Ritz polynomial Pk associated with step k:
k (k)
Y (θj )2 − σ 2
Pk (σ) = (k)
. (2.55)
j=1 (θj )2
(k)
Here, (θj )2 are the Ritz values, i.e., the k eigenvalues of AT A restricted to the
Krylov subspace Kk (AT A, AT b). The large Ritz values are approximations to some
of the large eigenvalues σi2 of the cross-product matrix AT A. Then the filter factors
associated with the solution after k steps of the CG algorithm are given by
(k)
fi = 1 − Pk (σi ) , i = 1, . . . , k . (2.56)
As k increases, and the Ritz values converge to some of the eigenvalues of AT A, then
(k)
for selected i and j we have θj ≈ σi . Moreover, as k increases these approximations
improve while, simultaneously, more eigenvalues of AT A are being approximated by
the additional Ritz values.
Equations (2.55) and (2.56) for the CG filter factors shed light on the regularizing
(k)
property of the CG method. After k iterations, if all the largest Ritz values (θj )2
have converged to all the largest eigenvalues σi2 of AT A, then the corresponding
Pk (σi ) ≈ 0 and the filter factors associated with these σi will therefore be close to
one. On the other hand, for all those eigenvalues smaller than the smallest Ritz value,
the corresponding filter factors satisfy
k
X ³ ´
(k) (k) (k) (k)
fi = σi2 (θj )−2 + O σi4 (θk )−2 (θk−1 )−2 ,
j=1
showing that these filter factors decay like σi2 for σi < θk (k).
From this analysis of the CG filter factors we see that the CG process indeed
has a regularizing effect if the Ritz values converge to the eigenvalues of AT A in
their natural order, starting with the largest. When this is the case, we are sure
that the CG algorithm is a regularizing process with the number of iterations k as
the regularization parameter. Unfortunately, proving that the Ritz values actually
converge in this order is a difficult task. For problems with a gap in the singular
value spectrum of A it is proved in [49, §6.4] that all the large eigenvalues of AT A
will be approximated by Ritz values before any of the small eigenvalues of AT A
get approximated. The case where the singular values of A decay gradually to zero
with no gap in the spectrum is more difficult to analyze—but numerical examples
and model problems [80, 81] indicate that the desired convergence of the Ritz values
actually holds as long as the discrete Picard condition is satisfied for the unperturbed
component of the right-hand side and there is a good separation among the large
singular values of A.
To put the CG method into the common framework from the previous section, we
notice that the solution x(k) after k CG steps can be defined as
min kA x − bk2 subject to x ∈ Kk (AT A, AT b) , (2.57)
where Kk (AT A, AT b) is the Krylov subspace associated with the normal equations.
Thus, we see that CG replaces the side constraint Ω(x) = kxk2 with the side constraint
x ∈ Kk (AT A, AT b). Obviously, if the Ritz values converge as desired, then the Krylov
subspace satisfies Kk (AT A, AT b) ≈ span{v1 , . . . , vk } indicating that the CG solution
x(k) is similar to, say, the TSVD solution xk .
Even the best implementation of the normal-equation CG algorithm suffers from
some loss of accuracy due to the implicit use of the cross-product matrix AT A. An
alternative iterative algorithm that avoids AT A completely is the algorithm LSQR
[70]. This algorithm uses the Lanczos bidiagonalization algorithm [30, §9.3.4] to build
up a lower bidiagonal matrix and, simultaneously, updates a QR factorization of this
bidiagonal matrix. The QR factorization, in turn, is used to update the LSQR solution
in each step. If Bk denotes the (k + 1) × k bidiagonalization matrix generated in the
(k)
kth step of LSQR, then the quantities θj in Eq. (2.55) are the singular values of this
Bk . Hence, the LSQR algorithm is mathematically equivalent to the normal-equation
CG algorithm in that the kth iteration vectors x(k) in CG and LSQR are identical in
exact arithmetic.
In real computations, the convergence of CG and LSQR is delayed due to the
influence of the finite precision arithmetic, and the dimension of the subspace in which
x(k) lies does not increase in each step. As a consequence, x(k) typically stays almost
unchanged for a few steps, then changes to a new vector and stays unchanged again for
some steps, etc. (The underlying phenomenon is related to CG and LSQR computing
“ghost” eigenvalues and singular values, respectively, cf. [30, §9.2.5]). In LSQR, it
is possible to store the so-called Lanczos vectors generated during the process and
apply some reorthogonalization scheme to them, which prevents the abovementioned
delay—but in practice it is usually less computationally demanding just to use LSQR

without any reorthogonalization. Orthogonalization can also be applied to the normal
equation residual vectors AT (A x(k) − b) in the CG algorithm.
There are several ways to implement the CG algorithm for the normal equations
in a numerically stable fashion. The one used in the routine cgls in Regularization
Tools is from [9, p. 289]. The implementation lsqr is identical to the original LSQR
algorithm from [70].
Regarding the filter factors for CG and LSQR, we have found that the expression
(2.56) using the Ritz polynomial is extremely sensitive to rounding errors. Instead, we
(k)
compute the filter factors fi by means of numerically more robust recursions derived
in [86] (see also [49]). Notice that the exact singular values σi of A are required to
compute the filter factors; hence this option is mainly of pedagogical interest.
2.8.2. Bidiagonalization with Regularization
It is possible to modify the LSQR algorithm and derive a hybrid between a direct
and an iterative regularization algorithm. The idea is to use the abovementioned
Lanczos algorithm to build up the bidiagonal matrix Bk sequentially, and in each
step to replace LSQR’s QR factorization of Bk with a direct regularization scheme
such as Tikhonov regularization or TSVD. These ideas are outlined in [8, 69]; see also
the discussion [39, Chapter 7]. The work involved in the direct regularization process
is small compared to the work in the iterative process because of the bidiagonal form
of Bk . Again, reorthogonalization of the Lanczos vectors is possible but rarely used
in practice.
One rationale behind this “hybrid” algorithm is that if the number k of Lanc-
zos bidiagonalization steps is so large that Bk becomes ill-conditioned and needs
(k)
regularization—because the singular values θk of Bk start to approximate some of
the smaller singular values of A—then hopefully all the large singular values of A are
approximated by singular values of Bk . When this is the case, then we are ensured
that the “hybrid” method computes a proper regularized solution, provided of course
that the explicit regularization in each step properly filters out the influence of the
small singular values.
The second, and perhaps most important, rational behind the “hybrid” algorithm
is that it requires a different stopping criterion which is not as dependent on choosing
the correct k as the previously mentioned methods. Provided again that the explicit
regularization scheme in each step is successful in filtering out the influence of the
small singular values, then after a certain stage k the iteration vector x(k) of the “hy-
brid” algorithm will hardly change. This is so because all the components associated
with the large singular values have been captured while the components associated
with the small singular values are filtered out. Thus, the stopping criteria should
now be based on the relative change in x(k) . With this stopping criteria, we see that
taking too many steps in the “hybrid” algorithm will not deteriorate the iteration
vector, but merely increase the computational effort.

For convenience, Regularization Tools provides the routine lanc b for com-
puting the lower bidiagonal matrix Bk as well as the corresponding left and right
Lanczos vectors. The routine csvd computes the SVD of Bk . The user can then
combine the necessary routines to form a specific “hybrid” algorithm.
2.8.3. The ν-Method

Both CG and LSQR converge rather fast to a regularized solution with damped high-
frequency components, and if the iterations are continued then the high-frequency
components very soon start to dominate the iteration vector. For some methods for
choosing the regularization parameter, i.e., the number of iterations k (such as the L-
curve criterion described in Section 2.9), this is a favorable property. However, there
are other circumstances in which is is more desirable to have an iterative scheme that
converges slower and thus is less sensitive to the choice of k.
This is exactly the philosophy behind the ν-method [11] which is similar to the CG
method except that the coefficients αk and βk used to update the iteration vectors in
Algorithm (2.54) are problem independent and depend only on the iteration number
k and a prespecified constant ν satisfying 0 < ν < 1:
(k + ν)(k + ν + 12 ) (k + ν)(k + 1)(k + 21 )

αk = 4 , βk = . (2.58)
(k + 2ν)(k + 2ν + 12 ) (k + 2ν)(k + 2ν + 21 )(k + ν + 1)
(It can be shown that the filter factors can be expressed in terms of Jacobi polyno-
mials).
A slight inconvenience with the ν-method is that it requires the problem to be
scaled such that kAk2 is slightly less than one, otherwise the method either diverges
or converges too slow. A practical way to treat this difficulty [39, §6.3] is use the
(k)
Lanczos bidiagonalization algorithm to compute a good approximation θ1 = kBk k2
(k)
to kAk2 and then rescale A and b by 0.99/θ1 . Usually a few Lanczos steps are
sufficient. This initialization process can also be used to provide the ν-method with
a good initial guess, namely, the iteration vector x(k) after a few LSQR steps.
The routine nu in Regularization Tools implements the ν-method as described
in [11] with the abovementioned rescaling. The LSQR start-vector is not used, but
can be supplied by the user if desired.
2.8.4. Extension to General-Form Problems

So far we have described several iterative methods for treating regularization prob-
lems in standard form; we shall now briefly describe the necessary extension to these
methods for treating problems in general form. The idea presented here is originally
from [36, 37]. From the discussion of the standard-form transformation for iterative
methods in §2.6.2, we see that essentially we must apply the above standard-form
iterative methods to a transformed problem with Ā and b̄. I.e., according to (2.57)
we must compute the solution x̄(k) to the problem
min kĀ x̄ − b̄k2 subject to x̄ ∈ Kk (ĀT Ā, ĀT b̄) . (2.59)
If we use the alternative formulation from §2.6.2 with Ā = A L†A and b̄ = b−A x0 ,
then the standard-form transformation can be “built into” the iterative scheme. In
this way, we work directly with x(k) and avoid the back-transformation from x̄(k) to
x(k) . To derive this technique, consider the side constraint in (2.59) which implies
that there exist constants ξ0 , . . . , ξk−1 such that
k−1
X
x̄(k) = ξi (ĀT Ā)i ĀT b̄ .
i=0
If we insert Ā = A L†A and b̄ = b − A x0 into this relation, we obtain

k−1
X ³ í
x̄(k) = ξi (L†A )T AT A L†A (L†A )T AT (b − A x0 ) .
i=0
Using Eqs. (2.31) and (2.30) for L†A and x0 together with the GSVD it is straight-
forward to show that (L†A )T AT A x0 = 0. Thus, by inserting the above expression for
x̄(k) into the back-transformation x(k) = L†A x̄(k) + x0 , we obtain
k−1
X ³ í
x(k) = ξi L†A (L†A )T AT A L†A (L†A )T AT b + x0 . (2.60)
i=0
From this relation we see that we can consider the matrix L†A (L†A )T a “preconditioner”
for the iterative methods, and we stress that the purpose of the “preconditioner” is
not to improve the condition number of the iteration matrix but rather to ensure
that the “preconditioned” iteration vector x(k) lies in the correct subspace and thus
minimizes kL x(k) k2 . Minimization of the correct residual norm is ensured by Eq.
(2.38).
“Preconditioning” is easy to build into CG, LSQR, and the ν-method by means
of the algorithms in (2.36) from §2.6.2. The special “preconditioned” versions are
implemented as routines pcgls, plsqr, and pnu in Regularization Tools.
2.9. Methods for Choosing the Regularization Parameter

No regularization package is complete without routines for computation of the reg-
ularization parameter. As we have already discussed in Section 2.5 a good regular-
ization parameter should yield a fair balance between the perturbation error and the
2.9. Methods for Choosing the Regularization Parameter 35
regularization error in the regularized solution. Throughout the years a variety of

parameter-choice strategies have been proposed. These methods can roughly be di-
vided into two classes depending on their assumption about kek2 , the norm of the
perturbation of the right-hand side. The two classes can be characterized as follows:
1. Methods based on knowledge, or a good estimate, of kek2 .
2. Methods that do not require kek2 , but instead seek to extract the necessary
information from the given right-hand side.
For many of these methods, the convergence rate for the solution as kek2 → 0 has
been analyzed [26, 34, 85]. Four parameter-choice routines are included in Regular-
ization Tools, one from class 1 and three from class 2.
The only method belonging to class 1 is the discrepancy principle [65, §27] which,
in all simplicity, amounts to choosing the regularization parameter such that the
residual norm for the regularized solution satisfies
kA xreg − bk2 = kek2 . (2.61)
When a good estimate for kek2 is known, this method yields a good regularization
parameter corresponding to a regularized solution immediately to the right of the L-
curve’s corner. Due to the steep part of the L-curve we see that an underestimate of
kek2 is likely to produce an underregularized solution with a very large (semi)norm.
On the other hand, an overestimate of kek2 produces an overregularized solution with
too large regularization error.
The three methods from class 2 that we have included in Regularization Tools
are the L-curve criterion, generalized cross-validation, and the quasi-optimality crite-
rion. The L-curve criterion has already been discussed in connection with the intro-
duction of the L-curve in Section 2.5. Our implementation follows the description in
[55] closely. For a continuous regularization parameter λ we compute the curvature
of the curve
(log kA xλ − bk2 , log kL xλ k2 )
(with λ as its parameter) and seek the point with maximum curvature, which we
then define as the L-curve’s corner. When the regularization parameter is discrete
we approximate the discrete L-curve in log-log scale by a 2D spline curve, compute
the point on the spline curve with maximum curvature, and define the corner of the
discrete L-curve as that point which is closest to the corner of the spline curve.
Generalized cross-validation (GCV) is based on the philosophy that if an arbitrary
element bi of the right-hand side b is left out, then the corresponding regularized so-
lution should predict this observation well, and the choice of regularization parameter
should be independent of an orthogonal transformation of b; cf. [87, Chapter 4] for
more details. This leads to choosing the regularization parameter which minimizes
the GCV function
kA xreg − bk22
G≡ , (2.62)
(trace(Im − A AI ))2
where AI is a matrix which produces the regularized solution xreg when multiplied
with b, i.e., xreg = AI b. Note that G is defined for both continuous and discrete
regularization parameters. The denominator in (2.62) can be computed in O(n) op-
erations if the bidiagonalization algorithm from Section 2.7 is used [25]. Alternatively,
the filter factors can be used to evaluate the denominator by means of the simple ex-
pression
Xp
trace(Im − A AI ) = m − (n − p) − fi . (2.63)
i=1
This is the approach used in routine gcv. In [47] it is illustrated that the GCV method
indeed seeks to balance the perturbation and regularization errors and thus, in turn,
is related to the corner of the L-curve.
The final method included in Regularization Tools is the quasi-optimality
criterion [65, §27]. This method is, strictly speaking, only defined for a continuous
regularization parameter λ and amounts to minimizing the function
° ° Ã p µ ¶2 !1/2
° dxλ ° X u T
b
Q ≡ λ° °
° dλ ° = fi (1 − fi ) i . (2.64)
2 i=1
γi
As demonstrated in [47], under certain assumptions the approach also corresponds to

finding a good balance between perturbation and regularization errors in xλ . For a
discrete regularization parameter k, we use λ = γk and the approximations
° °
° dxλ ° uTk b
° ° ≈ k∆xk k2 , k∆x k = , ∆λ = γk+1 − γk ≈ γk
° dλ ° |∆λ|
k 2
γk
2
to obtain the expressions

uTk b uTk b
Q≈ if L = In , Q≈ if L 6= In . (2.65)
σk γk
The discrepancy principle is implemented in routine discrep described in §2.7.2 in
connection with direct regularization methods. The L-curve criterion is implemented
in the two routines l curve and l corner, while GCV is is provided by routine gcv.
Finally, the quasi-optimality criterion is implemented in routine quasiopt.
2.10. New Functions in Version 4.1

There are two major changes in Version 4.1 of Regularization Tools, compared to
earlier versions: the package now also allows under-determined problems, and several
new functions were added. The previous sections of this chapter are kept unaltered,
but it is emphasized that functions cgsvd, discrep, dsvd, lsqi, tgsvd, and tikhonov now
also accept under-determined problems, i.e., coefficient matrices with fewer rows than
columns. This was perhaps the change that was most often requested by the users.
Note that for general-form problems, the dimensions must satisfy m + p ≥ n ≥ p.
2.10. New Functions in Version 4.1 37
New Test Problems

The test problem gravity is from the paper [51] on deconvolution and regularization
with Toeplitz matrices, and models a 1-D gravity surveying problem in which the
kernel is the vertical component of the gravity field from a point source located at
depth d:
¡ ¢−3/2
K(s, t) = d d2 + (s − t)2 .
The constant d controls the decay of the singular values (the larger the d, the faster
the decay). Three right-hand sides are available.
The test problem tomo represents a 2-D tomography problem in which the elements
of the right-hand b side are line integrals along straight rays penetrating a rectangular
domain, which is discretized into N2 cells, each cell with its own intensity stored as the
elements of the solution vector x. We use the standard ordering where the elements
in x corresponding to a stacking of the columns of the image. The elements of the
coefficient matrix A are given by
(
ìj , pixelj ∈ rayi
aij =
0 else,
where ìj is the length of the ith ray in pixel j. The exact solution is shown as
reshape (x,N,N) and it is identical to the exact image in the test problem blur. The
rays are placed randomly inside the domain, the number of rays can be specified, and
the coefficient matrix A is stored in sparse format.
New Iterative Regularization Methods

The function art implements the method by Kaczmarz, also known as the Algebraic
Reconstruction Technique (ART) [67]. Each iteration involves a sweep over all the
rows of the coefficient matrix A:
bi − aTi x
x←x+ ai , i = 1, . . . , m.
kai k22
Here bi is the ith component of the right-hand side b and aTi is the ith row of A.
We emphasize that since this method works on the individual rows of the coefficient
matrix, our implementation is much slower than cgls and other functions that use
matrix multiplications.
The function splsqr implements the subspace preconditioned LSQR (SP-LSQR)
algorithm [58] that uses LSQR to compute (an approximation to) the standard-form
Tikhonov solution xλ for a fixed value of the regularization parameter λ. The pre-
conditioner is based on the subspace spanned by the columns of an n × ksp matrix
Vsp whose columns should be “smooth” and preferably chosen such that a significant
component of the exact solution lies in the range of Vsp . For example, one can use
the first ksp columns of the DCT matrix dct (eye (n))0 .
It is emphasized that splsqr is a model implementation of the SP-LSQR algorithm.

In a real implementation the Householder transformations should use LAPACK rou-
tines, and if Vsp represents a fast transformation (such as the DCT) then explicit
storage of Vsp should be avoided. See [58] for implementation details.
Recently there has been interest in the use of minimum-residual methods that
avoid the use of AT , for regularization of problems with square matrices. It was
demonstrated in [59] that the iterative methods MR-II [38] and RRGMRES [12] –
which are variants of MINRES and GMRES, respectively, with starting vector Ab
instead of b – can have regularizing properties. In both methods, the kth iterate
solves the problem minx kA x − bk2 s.t. x ∈ span{Ab, A2 b, . . . , Ak b}. The absence of
the vector b from the Krylov subspace is crucial for the filtering of the noise in the
data [59].
The functions mr2 and rrgmres provide implementations of the MR-II and RRGM-
RES algorithms, respectively. The former has a coupled two-term recursion, while
the latter (similar to GMRES) needs to store the Arnoldi vectors and update a QR
factorization of the Hessenberg matrix. The storage requirements are thus essentially
the same as in MINRES and GMRES.
The standard-form transformations used in Regularization Tools for general-form
problems, with the smoothing (semi)norm kLxk2 , involve working with either AL† or
AL†A . Since these matrices are typically not square, the transformations cannot be
used in combination with MR-II and RRGMRES. Instead we use the smoothing-norm
preconditioning developed in [52], and the functions pmr2 and prrgmres implement this
approach for MR-II and RRGMRES.
New Parameter-Choice Methods

The function corner from [53] uses a so-called adaptive pruning algorithm to find the
corner of a discrete L-curve, e.g., produce by tsvd, tgsvd, or an iterative method. It
is now used in l curve and l corner. The advantages of the new method, compared to
the algorithm used previously in Regularization Tools, are the robustness of the
method (it avoids the “fudge factors” used in l corner) and the fact that it does not
rely on the Spline Toolbox. For backward compatibility, the spline-based algorithm
is still used if the Spline Toolbox is available.
The function ncp uses a statistical tool called the normalized cumulative peri-
odogram (NCP) to find the regularization parameter for which the corresponding
residual vector most resembles white noise [54]. The NCP is basically a cumulated
power spectrum, and its computation requires one FFT per residual vector r:
p = abs(fft(r)).^2; NCP = cumsum(p(2:1))/sum(p(2:q));
where q = floor(length(r)/2). The expected NCP for white noise is a straight line,
and we use the minimum distance from the NCP to a straight line as the criterion to
choose the regularization parameter.
3. Regularization Tools Tutorial
The purpose of this section is to give a brief introduction to the use of the routines
in Regularization Tools by means of some fairly simple examples. In particular,
we show how to compute regularized solutions and how to evaluate these solutions
by various graphical tools. Although the examples given below do not touch upon
all the features of Regularization Tools, they illustrate the fundamental ideas
underlying the package, namely, modularity and regularity between the routines. For
convenience, the examples are also available in the annotated script regudemo, with
appropriate pause statements added.
3.1. The Discrete Picard Condition

We shall first illustrate the use of the routine picard for visually checking the discrete
Picard condition, cf. Section 2.4. First, we generate a discrete ill-posed problem using
one of the many built-in test problems; the one used here is shaw which is a one-
dimensional model of an image restoration problem. Then we add white noise to the
right-hand side, thus producing a more “realistic” problem.
Before performing the analysis of the problem, we compute the SVD of the coeffi-
cient matrix; this is the typical situation in Regularization Tools, since most of
the routines make use of either the SVD (for standard-form problems) or the GSVD
(for general-form problems). We then use picard to plot the singular values and the
Fourier coefficients for both the unperturbed and the perturbed problem, see Fig. 3.1.
[A,b bar,x] = shaw (32);

randn (’seed’,41997);
e = 1e-3∗randn (size (b bar)); b = b bar + e;
[U,s,V] = csvd (A);
subplot (2,1,1); picard (U,s,b bar);
subplot (2,2,2); picard (U,s,b);
Clearly, most of the Fourier coefficients for the unperturbed problem satisfy the
discrete Picard condition—although eventually, for large i, both the singular values
and the Fourier coefficients become dominated by rounding errors. For the “noisy”
test problem, we see that the Fourier coefficients for the right-hand side become
dominated by the perturbation for i much smaller than before. We also see that
40 TUTORIAL
Picard plot
10
10
σ
i
0
10 |uTi b|
|uTi b|/σi
−10
10
−20
10
0 5 10 15 20 25 30 35
i
Picard plot
20
10
10
σi
10 T
|ui b|
0 |uTb|/σ
10 i i
−10
10
−20
10
0 5 10 15 20 25 30 35
i
Figure 3.1: Output from picard for the “pure” and the “noisy” test problems.
in the left part of curve the Fourier coefficients still decay faster than the singular
values, indicating that the unperturbed right-hand side satisfies the discrete Picard
condition. To regularize this problem, we should preferably dampen the components
for which the perturbation dominates, and leave the rest of the components (for small
i) intact.
3.2. Filter Factors

In this part of the tutorial we consider only the “noisy” test problem from Example
3.1, we compute regularized solutions by means of Tikhonov’s method and LSQR,
and we compute the filter factors for two regularization methods.
lambda = [1,3e-1,1e-1,3e-2,1e-2,3e-3,1e-3,3e-4,1e-4,3e-5];
X tikh = tikhonov (U,s,V,b,lambda);
F tikh = fil fac (s,lambda);
iter = 30; reorth = 0;
[X lsqr,rho,eta,F lsqr] = lsqr b (A,b,iter,reorth,s);
subplot (2,2,1); surf (X tikh), axis (’ij’), title (’Tikhonov solutions’)
subplot (2,2,2); surf (log10 (F tikh)), axis (’ij’), title (’Tikh filter factors, log scale’)
subplot (2,2,3); surf (X lsqr (:,1:17)), axis (’ij’), title (’LSQR solutions’)
subplot (2,2,4); surf (log10 (F lsqr)), axis (’ij’), title (’LSQR filter factors, log scale’)
3.3. The L-Curve 41
Tikhonov solutions Tikh filter factors, log scale
4 0
2
−20
0
−2 −40
0 0
10 10
20 20
5 5
40 0 40 0
LSQR solutions LSQR filter factors, log scale
5 20
0
0
−20
−5 −40
0 0
20 20
20 20
10 10
40 0 40 0
Figure 3.2: Regularized solutions and filter factors for the “noisy” test problem.
We use the command mesh to plot all the regularized solutions, cf. Fig. 3.2. This
is a very convenient to show the dependence of the solution on the regularization
parameter. The same technique is used to display the the corresponding filter factors.
We see the typical situation for regularization methods: first, when we apply much
regularization, the solution is overregularized, i.e., it is too smooth; then it becomes
a better approximation to the underlying, unperturbed solution; and eventually the
solution becomes underregularized and starts to be dominated by the perturbation
errors, and its (semi)norm “explodes”.
3.3. The L-Curve
The L-curve analysis plays an important role in the analysis phase of regularization
problems. The L-curve displays the tradeoff between minimizing the two quantities in
the regularization problem, namely, the residual norm and the solution (semi)norm,
and it shows how these quantities depend on the regularization parameter. In ad-
dition, the L-curve can be used to compute the “optimal” regularization parameter
as explained in Section 2.9 (we return to this aspect in the next example). The L-
curve can also be used to investigate the similarity between different regularization
methods—if their L-curves are close, then the regularized solutions are similar, cf. [47].
42 TUTORIAL
L−curve L−curve
3 3
10 10
2 2
10 10
5.4947e−006
solution norm || x ||2
10
1
10
1
0.00015565
0.0044091 0.1249
0 0
10 10
−3 −2 −1 0 −3 −2 −1 0
10 10 10 10 10 10 10 10
residual norm || A x − b || residual norm || A x − b ||2
2
Figure 3.3: The L-curves for Tikhonov regularization and for LSQR.
subplot (1,2,1); l curve (U,s,b); axis ([1e-3,1,1,1e3])

subplot (1,2,2); plot lc (rho,eta,’o’), axis ([1e-3,1,1,1e3]
For the particular problem considered here, the “noisy” test problem from Ex-
ample 3.1, the L-curves for both Tikhonov regularization and for LSQR, shown in
Fig. 3.3, have a particularly sharp corner. This corner corresponds to a regularized
solution where the perturbation error and the regularization error are balanced. The
similarity of the L-curves for the two methods indicate the similarity between the
methods.
3.4. Regularization Parameters

We shall now use the L-curve criterion and the GCV method to compute “optimal”
regularization parameters for two methods: Tikhonov regularization and truncated
SVD. This is very easy to do by means of the routines l curve and gcv. We then
use our knowledge about the exact solution to compute the relative errors in the
four regularized solutions. The best result may depend on the particular computer’s
floating point arithmetic—for the Toshiba PC used here, the combination of truncated
SVD and the L-curve criterion gives the best result.
3.5. Standard Form Versus General Form 43
L−curve, Tikh. corner at 0.00087288 L−curve, TSVD corner at 9

3 3
10 10
12
2 2
10 10
5.4947e−006
2
solution norm || x ||
1 1
10 10
0.00015565
0.0044091 0.1249 ---------,9 6 3
0 0
10 10
−3 −2 −1 0 −3 −2 −1 0
10 10 10 10 10 10 10 10
residual norm || A x − b ||2 residual norm || A x − b ||
2
GCV function, minimum at λ = 0.004616 −1

GCV function, minimum at k = 7
−1
10 10
−2 −2
10 10
10
10
−3
−4
/ 10
10
−3
−4
G(λ)
G(k)
−5 −5
10 10
−6 −6
10 10
−7 −7
10 10
−8 −8
10 10
−9
10 −9
−6 −5 −4 −3 −2 −1 0 10
10 10 10 10 10 10 10 0 2 4 6 8 10 12 14 16 18 20
λ k
Figure 3.4: Comparison of the L-curve criterion and the GCV method for computing
the “optimal” regularization parameter for Tikhonov’s method and for TSVD.
lambda l = l curve (U,s,b); axis ([1e-3,1,1,1e3])

k l = l curve (U,s,b,’tsvd’); axis ([1e-3,1,1,1e3])
lambda gcv = gcv (U,s,b); axis ([-6,0,-9,-1])
k gcv = gcv (U,s,b,’tsvd’); axis ([0,20,1e-9,1e-1])
x tikh l = tikhonov (U,s,V,b,lambda l);
x tikh gcv = tikhonov (U,s,V,b,lambda gcv);
x tsvd l = tsvd (U,s,V,b,k l);
x tsvd gcv = tsvd (U,s,V,b,k gcv);
[norm (x−x tikh l),norm (x−x tikh gcv),. . .
norm (x−x tsvd l),norm (x−x tsvd gcv)]/norm (x)
3.5. Standard Form Versus General Form

In this example we illustrate the difference between regularization in standard and
general form. In particular, we show that regularization in general form is sometimes
44 TUTORIAL
necessary to ensure the computation of a satisfactory solution. Unfortunately, this

judgement cannot be automated, but requires insight from the user about the desired
solution.
The test problem used here is the inverse Laplace transformation, and we generate
a problem whose exact solution is f (t) = 1 − exp (−t/2). This solution obviously
satisfies f → 1 for t → ∞, and the horizontal asymptotic part in the discretized
solution for n = 16 is visible for indices i > 8.
n = 16; [A,b,x] = i laplace (n,2);

b = b + 1e-4∗randn (b);
L = get l (n,1);
[U,s,V] = csvd (A); [UU,sm,XX] = cgsvd (A,L);
I = 1; for i=[3,6,9,12]
subplot (2,2,I); plot (1:n,V(:,i)); axis ([1,n,-1,1])
xlabel ([’i = ’,num2str(i)]), I = I + 1;
end
subplot (2,2,1), text (12,1.2,’Right singular vectors V(:,i)’)
I = 1; for i=[n-2,n-5,n-8,n-11]
subplot (2,2,I); plot (1:n,XX(:,i)); axis ([1,n,-1,1])
xlabel ([’i = ’,num2str(i)]), I = I + 1;
end
subplot (2,2,1), text (10,1.2,’Right generalized singular vectors XX(:,i)’)
k tsvd = 7; k tgsvd = 6;
X I = tsvd (U,s,V,b,1:k tsvd);
X L = tgsvd (UU,sm,XX,b,1:k tgsvd);
subplot (2,1,1); plot (1:n,X I,1:n,x,’x’), axis ([1,n,0,1.2]), xlabel (’L = I’)
subplot (2,2,2), plot (1:n,X L,1:n,x,’x’), axis ([1,n,0,1.2]), xlabel (’L \neq I’)
In this example we choose minimization of the first derivative as side constraint in

the general-form regularization, i.e., we use L = tridiag(−1, 1). It is very instructive
to first inspect the right singular vectors vi and xi from the SVD and the GSVD,
respectively, cf. Fig. 3.5. Notice that none of the SVD vectors vi include the above-
mentioned asymptotic part; hence, these vectors are not suited as basis vectors for
a regularized solution. The GSVD vectors xi , on the other hand, do posess the the
necessary asymptotic part, and they are therefore much better suited as basis vectors.
For a thorough discussion of these aspects, cf. [84].
Let us now use the truncated SVD and the truncated GSVD to compute regular-
ized solutions to this problem. The first 7 TSVD solutions and the first 6 TGSVD
solutions are shown in Fig. 3.6. From our investigation of the right singular vectors, it
is no surprise that TSVD cannot reproduce the desired solution, while TGSVD indeed
succeeds in doing so. The optimal regularization parameter for TGSVD is k = 5.
We notice that without our a priori knowledge part of the solution, we could not
immediately discard the TSVD solutions!
3.5. Standard Form Versus General Form 45
Right singular vectors V(:,i) Right generalized singular vectors XX(:,i)

1 1 1 1
0.5 0.5 0.5 0.5
0 0 0 0
−0.5 −0.5 −0.5 −0.5
−1 −1 −1 −1
5 10 15 5 10 15 5 10 15 5 10 15
i=3 i=6 i = 14 i = 11
1 1 1 1
0.5 0.5 0.5 0.5
0 0 0 0
−0.5 −0.5 −0.5 −0.5
−1 −1 −1 −1
5 10 15 5 10 15 5 10 15 5 10 15
i=9 i = 12 i=8 i=5
Figure 3.5: Comparison of the right singular vectors for SVD and GSVD for the
inverse Laplace transformation test-problem.
0.8
0.6
0.4
0.2
0
2 4 6 8 10 12 14 16
L=I
1 I L xxxxxxxxx
0.8
K
,/
0.6
0.4
0.2
0
2 4 6 8 10 12 14 16
L≠I
Figure 3.6: The first 7 TSVD solutions and the first 6 TGSVD solutions for the same
test problem as in Fig. 3.5. The exact solution is shown with ×-markers.
46 TUTORIAL
Picard plot
10
10
σi
|uTi b|
T
5
10 |ui b|/σi
0
10
−5
10
−10
10
−15
10
0 5 10 15 20 25 30 35
i
Figure 3.7: The output from picard for the test problem that does not satisfy the
discrete Picard condition.
3.6. No Square Integrable Solution

In this tutorial’s final example we use the routine ursell to generate a test problem
arising from discretization of a Fredholm integral equation of the first kind (2.1) with
no square integrable solution [19, p. 6], i.e., the integral equation does not satisfy the
Picard condition. The integral equation is
Z 1
1
f (t) dt = 1 , 0≤s≤1.
0 s + t+1
When we use picard to plot the singular values and the Fourier coefficients for the
discrete problem, see Fig. 3.7, we immediately see that discrete ill-posed problem
does not satisfy the discrete Picard condition, which indicates trouble! We stress,
however, that such an analysis cannot show whether the trouble comes from the
integral equation itself, from the discretization, or possibly from other sources.
[A,b] = ursell (32); [U,s,V] = csvd (A); picard (U,s,b);

4. Regularization Tools Reference
SVD- and GSVD-BASED REGULARIZATION ROUTINES

discrep Minimizes the solution (semi)norm subject to an upper
bound on the residual norm (discrepancy principle)
dsvd Computes a damped SVD/GSVD solution
lsqi Minimizes the residual norm subject to an upper bound
on the (semi)norm of the solution
mtsvd Computes the modified TSVD solution
tgsvd Computes the truncated GSVD solution
tikhonov Computes the Tikhonov regularized solution
tsvd Computes the truncated SVD solution
ttls Computes the truncated TLS solution
ITERATIVE REGULARIZATION ROUTINES

art Algebraic reconstruction technique (Kaczmarz’s method)
cgls Computes the least squares solution based on k steps
of the conjugate gradient algorithm
lsqr b Computes the least squares solution based on k steps
of the LSQR algorithm
maxent Computes the maximum entropy regularized solution
mr2 Solution of symmetric indefinite problems by MR-II
nu Computes the solution based on k steps of Brakhage’s
iterative ν-method
pcgls Same as cgls, but for general-form regularization
plsqr b Same as lsqr b, but for general-form regularization
pmr2 Same as mr2, but for general-form regularization
pnu Same as nu, but for general-form regularization
prrgmres Same as rrgmres, but for general-form regularization
rrgmres Range-restricted GMRES algorithm for square systems
splsqr Computes an approximate standard-form Tikhonov solution
via the subspace preconditioned LSQR algorithm (SP-LSQR)
48 Chapter 4. Regularization Tools Reference
UTILITY ROUTINES
bidiag Bidiagonalization of a matrix by Householder transformations
cgsvd Computes the compact generalized SVD of a matrix pair
csvd Computes the compact SVD of an m × n matrix
get l Produces a (n − d) × n matrix which is the discrete
approximation to the dth order derivative operator
lanc b Lanczos bidiagonalization process with/without reorthogonalization
regutm Generates random test matrices for regularization methods
ANALYSIS ROUTINES
corner Locates the corner of a discrete L-curve
fil fac Computes filter factors for some regularization methods
gcv Plots the GCV function and computes its minimum
lagrange Plots the Lagrange function kA x − bk22 + λ2 kL xk22 and its
derivative
l corner Locates the L-shaped corner of the L-curve
l curve Computes the L-curve, plots it, and computes its corner
ncp Plots the NCPs and locates the one closest to a straight line
picard Plots the (generalized) singular values, the Fourier coefficients
for the right-hand side, and the solution’s Fourier-coefficients
plot lc Plots an L-curve
quasiopt Plots the quasi-optimality function and computes its minimum
TEST PROBLEMS
baart First kind Fredholm integral equation
blur Image deblurring test problem
deriv2 Computation of second derivative
foxgood Severely ill-posed test problem
gravity One-dimensional gravity surveying problem
heat Inverse heat equation
i laplace Inverse Laplace transformation
parallax Stellar parallax problem with real observations
phillips Phillips’ “famous” test problem
shaw One-dimensional image restoration model
spikes Test problem with a “spiky” solution
tomo Two-dimensional tomography problem with sparse matrix
ursell Integral equation with no square integrable solution
wing Test problem with a discontinuous solution
Regularization Tools Reference 49
STANDARD-FORM TRANSFORMATION
gen form Transforms a standard-form solution back into the
general-form setting
std form Transforms a general-form problem into one in
standard form
AUXILIARY ROUTINES
app hh l Applies a Householder transformation from the left
gen hh Generates a Householder transformation
lsolve Inversion with A-weighted generalized inverse of L
ltsolve Inversion with transposed A-weighted generalized inverse of L
pinit Initialization for treating general-form problems
regudemo Tutorial introduction to Regularization Tools
spleval Computes points on a spline or spline curve
50 Chapter 4. Regularization Tools Reference
The Test Problems

There are 14 built-in test problems in Regularization Tools. Thirteen of them
are taken from the literature (cf. the following manual pages for references) while
the remaining, spikes, is “cooked up” for this package. All of them have in common
that they are easy to generate, and they share the characteristic features of discrete
ill-posed problems mentioned in Section 2.3.
All the test problems are derived from discretizations of a Fredholm integral equa-
tion of the first kind. Two different discretization techniques are used: the quadrature
method and the Galerkin method with orthonormal basis functions. In the quadrature
method [19, Chapter 6], the integral is approximated by a weighted sum, i.e.,
Z b n
X
K(s, t) f (t) dt ≈ In (s) = wj K(s, tj ) f (Tj ) .
a i=1
In particular, for the midpoint rule we use wj = (b − a)/n and tj = (j − 12 )(b − a)/n,
j = 1, . . . , n. Collocation in the n points s1 , . . . , sn then leads to the requirements
In (si ) = g(si ), i = 1, . . . , n. Hence, we obtain a system of linear algebraic equations
A x = b with elements given by aij = wj K(si , tj ) and bi = g(si ) for i, j = 1, . . . , n.
If the solution f is known then we represent it by x with elements xj = f (tj ),
j = 1, . . . , n. In the Galerkin method, we choose the following orthonormal box
functions as basis functions:
( (
− 12 −1
ψi (s) = h s , s ∈ [si−1 , s i ] , φi (t) = ht 2 , t ∈ [ti−1 , ti ]
0 , elsewhere 0 , elsewhere
in which hs = (d − c)/n, ht = (b − a)/n, and si = ihs , ti = iht , i = 0, . . . , n. Then the

Galerkin method [3] leads to a linear system of equations A x = b with elements given
by Eq. (2.4). Similarly, we represent the solution f by the vector x with elements
Rb
xj = a φj (t) f (t) dt, j = 1, . . . , n. We stress that for both methods the product A x
is, in general, different from b. The table below gives an overview of the 12 test
problems, while graphs of x for n = 100 are given in the individual manual pages.
problem discretization Ax = b problem discretization Ax = b

baart Galerkin no parallax Galerkin no x
blur – yes phillips Galerkin no
deriv2 Galerkin yes shaw quadrature yes
foxgood quadrature no spikes “cooked up” yes
gravity quadrature yes tomo quadrature yes
heat quadrature yes ursell Galerkin no x
ilaplace quadrature yes wing Galerkin no
app hh 51
app hh
Purpose:
Apply a Householder transformation.
Synopsis:
A = app hh (A,beta,v)
Description
app hh applies the Householder transformation, defined by the vector v and the
scaler beta, to the matrix A; i.e.,
A ← (In − beta v vT ) A .
See also:
gen hh
52 art
art
Purpose:
Algebraic reconstruction technique (Kaczmarz’s method).
Synopsis:
[X,rho,eta] = art (A,b,k)
Description:
Classical Kaczmarz iteration, or ART (algebraic reconstruction technique), applied
to the system A x = b. The number of iterations is k and the iterates are stored as
the columns of X. The vectors rho and eta hold the norms of the iterates and the
corresponding residual vectors.
In this so-called row action method, the jth iteration consists of a “sweep” through
the rows aTi = A(i, : ) of A, in which the current solution vector x[j] (for j = 0, 1, 2, . . .)
is updated as follows:
(0)
x[j ] = x[j]
for i = 1, . . . , m
(i−1)
(i) (i−1) bi − aTi x[k ]
x[j ]
= x[j ]
+ ai
kai k22
end
(m)
x[j+1] = x[j ]
Here bi is the ith component of the right-hand side b and ai is the ithe row turned
into a column vector.
Examples:
Generate a “noisy” 16 × 16 tomography problem of size n = 256, compute 50 ART
iterates, and show the last iterate:
[A,b,x] = tomo (16); b = b + 1e-3∗randn (size(b));
X = art (A,b,50); imagesc (reshape (X(:,end),16,16)); axis image
See also:
cgls, lsqr b, mr2, nu, rrgmres
References:
1. F. Natterer and F. Wübbeling, Mathematical Methods in Image Reconstruction,
SIAM, Philadelphia, 2001.
baart 53
baart
Purpose:
Test problem: Fredholm integral equation of the first kind.
Synopsis:
[A,b,x] = baart (n)
Description:
Discretization of an artificial Fredholm integral equation of the first kind (2.1)
with kernel K and right-hand side g given by
sin s
K(s, t) = exp(s cos t) , g(s) = 2 ,
s
and with integration intervals s ∈ [0, π2 ] and t ∈ [0, π]. The solution is given by
f (t) = sin t .
The size of the matrix A is n × n.
Examples:
Generate a “noisy” problem of size n = 32:
[A,b,x] = baart (32); b = b + 1e-3∗randn (size(b));
Limitations:
The order n must be even.
References:
1. M. L. Baart, The use of auto-correlation for pseudo-rank determination in noisy
ill-conditioned linear least-squares problems, IMA J. Numer. Anal. 2 (1982),
241–247.
0.18
0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0
0 10 20 30 40 50 60 70 80 90 100
54 bidiag
bidiag
Purpose:
Bidiagonalization of an m × n matrix with m ≥ n.
Synopsis:
B = bidiag (A)
[U,B,V] = bidiag (A)
Description:
If A is an m × n matriz with m ≥ n then bidiag uses Householder transformations
to compute a bidiagonalization of A:
A = U B VT ,
where B ∈ Rn×n is upper bidiagonal,
 
b11 b12
 b22 b23 
B=  .. ..  ,
. .
bnn
and the matrices U ∈ Rm×n and V ∈ Rn×n have orthonormal columns. The bidiagonal
matrix B is stored as a sparse matrix.
Examples:
Compute the bidiagonalization of A and compare the singular value of A with
those of B (they should be identical to about machine precision):
B = bidiag (A); [svd (A),csvd (B)]
Algorithm:
Alternating left and right Householder transformations are used to bidiagonalized
A. If U and V are also required, then the Householder transformations are accumu-
lated.
Limitations:
The case m < n is not allowed.
See also:
bsvd, lanc b
References:
1. L. Eldén, Algorithms for regularization of ill-conditioned least-squares problems,
BIT 17 (1977), 134-145.
blur 55
blur
Purpose:
Test problem: deblurring of images degraded by atmospheric turbulence blur.
Synopsis:
[A,b,x] = blur (N,band,sigma)
Description:
This image deblurring problem arises in connection with the degradation of digital
images by atmospheric turbulence blur, modelled by a Gaussian point-spread function:
µ ¶
1 x2 + y 2
h(x, y) = exp − .
2 π sigma2 2 sigma2
The matrix A is a symmetric N2 × N2 doubly Toeplitz matrix, stored in sparse format,

and given by A = (2 π sigma2 )−1 T ⊗T , where T is an N×N symmetric banded Toeplitz
matrix whose first row is given by
z = [exp(−([0 : band − 1].^2)/(2 ∗ sigma^2)); zeros(1, N − band)].
Only elements within a distance band − 1 from the diagonal are stored; i.e., band is
the half-bandwidth of the matrix T . If band is not specified, then band = 3 is used.
The parameter sigma controls the shape of the Gaussian point spread function
and thus the amount of smoothing (the larger the sigma, the wider the function, and
the less ill posed the problem). If sigma is not specified, then sigma = 0.7 is used.
The vector x is a columnwise stacked version of a simple test image, while b holds
a columnwise stacked version of the blurrred image; i.e, b = A*x.
Limitations:
The integer N should not be too small; we recommend N ≥ 16.
Reference:
1. M. Hanke & P. C. Hansen, Regularization methods for large-scale problems,
Surv. Math. Ind. 3 (1993), 253–315.
56 cgls
cgls
Purpose:
Conjugate gradient algorithm applied implicitly to the normal equations.
Synopsis:
[X,rho,eta,F] = cgls (A,b,k,reorth,s)
Description:
Performs k steps of the conjugate gradient algorithm applied implicitly to the
normal equations AT A x = AT b.
The routine returns all k solutions, stored as columns of the matrix X. The corre-
sponding solution norms and residual norms are returned in eta and rho, respectively.
If the singular values s are also provided, cgls computes the filter factors associated
with each step and stores them columnwise in the matrix F.
Reorthogonalization of the normal equation residual vectors AT (A X(: , i) − b),
i = 1, . . . , k is controlled by means of reorth as follows:
reorth = 0 : no reorthogonalization
reorth = 1 : reorthogonalization by means of MGS.
No reorthogonalization is assumed if reorth is not specified.
A “preconditioned” version of cgls for the general-form problem, where one mini-
mizes kL xk2 instead of kxk2 in each step, is implemented in routine pcgls.
Examples:
Perform 25 iterations and plot the corresponding L-curve:
[X,rho,eta] = cgls (A,b,25); plot lc (rho,eta,’o’);
Algorithm:
The algorithm, which avoids explicit formation of the normal-equation matrix
AT A, is described in [1]. The computation of the filter factors, using the recurrence
relations for the solution and the residual vector, is described in [2].
See also:
art, lsqr b, mr2, nu, pcgls, plsqr b, rrgmres, splsqr
References:
1. Å. Björck, Numerical Methods for Least Squares cProblems, SIAM, Philadel-
phia, 1996.
2. C. R. Vogel, Solving ill-conditioned linear systems using the conjugate gradient
method, Report, Dept. of Mathematical Sciences, Montana State University,
1987.
cgsvd 57
cgsvd
Purpose:
Compute the compact generalized SVD of a matrix pair A ∈ Rm×n and L ∈ Rp×n .
Synopsis:
sm = cgsvd (A,L)
[U,sm,X,V] = cgsvd (A,L) , where sm = [sigma,mu]
Description:
Computes the generalized singular value decomposition of the matrix pair (A, L).
If m ≥ n then the GSVD has the form
µ ¶
diag(sigma) 0
A=U X−1 , L = V (diag(mu) 0) X−1 ,
0 In−p
where the matrices U ∈ Rm×n and V ∈ Rp×p have orthonormal columns, and the
matrix X ∈ Rn×n is nonsingular. Moreover, the vectors sigma and mu have length p,
and sigma./mu are the generalized singular values of (A, L).
If m < n then it is required that m + p ≥ n, and the GSVD takes the form
µ ¶ µ ¶
0 diag(sigma) 0 −1 In−m 0 0
A=U X , L=V X−1 ,
0 0 In−p 0 diag(mu) 0
where the matrices U ∈ Rm×m and V ∈ Rp×p have orthonormal columns, the matrix
X ∈ Rn×n is nonsingular, and the vectors sigma and mu have length m + p − n. Notice
that the c, s-pairs sigma and mu are returned in an array sm with two columns, and
that only the last m columns of X are returned when m < n.
A possible fifth output argument returns W = X−1 . The matrices V and W are
not required by any of the routines in this package.
Algorithm:
Calls the Matlab routine gsvd and stores the c, s-pairs in an array with two
columns. The gsvd routine is included in Matlab Version 5.2, and it is based on [1].
Limitations:
The dimensions must satisfy m ≥ n ≥ p or m + p ≥ n ≥ p, which is no restriction
in connection with regularization problems.
See also:
csvd
References:
1. C. F. Van Loan, Computing the CS and the generalized singular value decom-
position, Numer. Math. 46 (1985), 479–491.
58 corner
corner
Purpose:
Find the corner of a discrete L-curve via an adaptive pruning algorithm.
Synopsis:
[k corner,info] = corner (rho,eta,fig)
Description:
Returns the integer k corner so that the corner of the log-log L-curve is located at
(log(rho(k corner)), log(eta(k corner))) .
The vectors rho and eta must contain corresponding values of the residual norm
kA x − bk2 and the solution’s (semi)norm kxk2 or kL xk2 for a sequence of regularized
solutions, ordered such that rho and eta are monotonic and such that the amount of
regularization decreases as their index increases. If a third input argument is present,
then a figure will show the discrete L-curve in log-log scale as well as the corner.
The second output argument describes possible warnings. Any combination of
zeros and ones is possible.
info = 000: No warnings – rho and eta describe a discrete L-curve with a corner.
info = 001: Bad data – some elements of rho and/or eta are Inf, NaN, or zero.
info = 010: Lack of monotonicity – rho and/or eta are not strictly monotonic.
info = 100: Lack of convexity – the L-curve is concave and has no corner.
The warnings described above will also result in text warnings on the command line;
type warning off Corner:warnings to disable all command line warnings.
Examples:
Compute TSVD solutions and find the corner of the discrete L-curve:
[A,b,x] = shaw (32); b = b + 1e-4*randn (size(b));
[U,s,V] = csvd (A); [X,rho,eta] = tsvd (U,s,V,b,1:14);
k corner = corner (rho,eta,1)
Algorithm:
The algorithm uses the two-stage approach from [1] to locate the corner. First a
list of corner candidates is found by successively pruning the points on the discrete
L-curve, and then the “best” corner is found from this list.
See also:
l corner, l curve
References:
1. P. C. Hansen, T. K. Jensen & G. Rodriguez, An adaptive pruning algorithm for
the discrete L-curve criterion, J. Comp. Appl. Math. 198 (2007), 483–492.
csvd 59
csvd
Purpose:
Compute the compact singular value decomposition.
Synopsis:
s = csvd (A)
[U,s,V] = csvd (A)
[U,s,V] = csvd (A,’full’)
Description:
Computes the compact form of the singular value decomposition:
A = U diag(s) VT .
If the matrix A is m × n, then the dimensions of the computed quantities are:

s min(m, n) × 1
U m × min(m, n)
V n × min(m, n)
If a second argument is present, then the full U and V are returned.
Algorithm:
Calls the Matlab routine svd and stores the singular values in a column vector s.
See also:
cgsvd
60 deriv2
deriv2
Purpose:
Test problem: computation of the second derivative.
Synopsis:
[A,b,x] = deriv2 (n,case)
Description:
This is a mildly ill-posed problem; i.e., its singular values decay slowly to zero. It
is a discretization of a first kind Fredholm integral equation (2.1) whose kernel K is
the Green’s function for the second derivative:
½
s(t − 1) , s < t
K(s, t) = .
t(s − 1) , s ≥ t
Both integration intervals are [0, 1], and as right-hand side g and corresponding solu-
tion f one can choose between the following:
case = 1 g(s) = (s3 − s)/6 , f (t) = t
case = 2 g(s) = exp(s) + (1 − e)s − 1 , f (t) = exp(t)

½ 1
½ 1
(4s3 − 3s)/24 , s< 2 t , t< 2
case = 3 g(s) = 1 , f (t) = 1
(−4s3 + 12s2 − 9s + 1)/24 , s≥ 2 1−t , t≥ 2
The first two examples are from [1, p. 315] while the third example is from [2]. The
size of the matrix A is n × n.
References:
1. L. M. Delves & J. L. Mohamed, Computational Methods for Integral Equations,
Cambridge University Press, Cambridge, 1985.
2. A. K. Louis & P. Maass, A mollifier method for linear operator equations of the
first kind, Inverse Problems 6 (1990), 427–440.
0.3
0.25
0.2
case = 2
0.15
0.1 case = 1
0.05 case = 3
0
0 10 20 30 40 50 60 70 80 90 100
discrep 61
discrep
Purpose:
Disrepancy principle criterion for choosing the regularization parameter.
Synopsis:
[x delta,lambda] = discrep (U,s,V,b,delta,x 0)
[x delta,lambda] = discrep (U,sm,X,b,delta,x 0) , where sm = [sigma,mu]
Description:
Least squares minimization with a quadratic inequality constraint:
min kx − x 0k2 subject to kA x − bk2 ≤ delta
min kL (x − x 0)k2 subject to kA x − bk2 ≤ delta
where x 0 is an initial guess of the solution, and delta is a positive constant. The
routine discrep requires either the compact SVD of A stored as U, s, and V, or part
of the GSVD of (A, L) saved as U, sm, and X. The regularization parameter lambda
corresponding to the solution x delta is also returned.
If delta is a vector, then x delta is a matrix such that
x delta = [ x delta(1), x delta(2), . . . ] .
If x 0 is not specified, x 0 = 0 is used.
The “opposite” problem, namely that of minimizing the residual norm subject to
an upper bound on the solution (semi)norm, is treated by routine lsqi.
Examples:
Use discrep to solve the test problem shaw:
[A,b,x] = shaw (n); [U,s,V] = csvd (A); e = 1e-3∗randn (size(b)); b = b + e;
x lambda = discrep (U,s,V,b,norm (e)); plot ([x,x lambda]);
Algorithm:
The algorithm, which uses a Newton method for computing the regularization
parameter lambda (implemented in routine newton), is described in [1]. The starting
value of lambda is set equal to that singular value sk of A, or that generalized singular
value γk of (A, L), for which the corresponding TSVD/TGSVD residual norm kA xk −
bk2 is closest to delta.
See also:
lsqi, l curve, ncp, newton, quasiopt
References:
1. V. A Morozov, Methods for Solving Incorrectly Posed Problems, Springer, New
York, 1984; Chapter 26.
62 dsvd
dsvd
Purpose:
Compute the damped SVD solution.
Synopsis:
[x lambda,rho,eta] = dsvd (U,s,V,b,lambda)
[x lambda,rho,eta] = dsvd (U,sm,X,b,lambda) , where sm = [sigma,mu]
Description:
Computes the damped SVD solution, defined as
x lambda = V (diag(s + lambda))−1 UT b (standard form)

µ ¶−1
diag(sigma + lambda ∗ mu) 0
x lambda = X UT b (general form)
0 In−p
If lambda is a vector, then x lambda is a matrix such that
x lambda = [ x lambda(1), x lambda(2), . . . ] .
The solution norm (standard-form case) or seminorm (general-form case) and the
residual norm are returned in eta and rho.
Algorithm:
Straightforward use of the definitions are used to compute x lambda.
References:
1. M. P. Ekstrom & R. L Rhodes, On the application of eigenvector expansions to
numerical deconvolution, J. Comp. Phys. 14 (1974), 319–340.
fil fac 63
fil fac
Purpose:
Compute the filter factors for some regularization methods.
Synopsis:
f = fil fac (s,reg param,method)
f = fil fac (sm,reg param,method) , where sm = [sigma,mu]
f = fil fac (s,reg param,’ttls’,s1,V1)
Description:
Computes all the filter factors corresponding to the singular values in s and the
regularization parameter reg param, for the following methods:
method = ’dsvd’ : damped SVD or GSVD
method = ’tsvd’ : truncated SVD or GSVD
method = ’Tikh’ : Tikhonov regularizaton
method = ’ttls’ : truncated TLS
If sm = [sigma,mu] is specified, then the filter factors for the corresponding generalized
methods are computed. If method is not specified, ’Tikh’ is default.
If method = ’ttls’ then the singular values s1 and the right singular matrix V1 of (A b)
must also be supplied.
Examples:
Compute the filter factors for Tikhonov regularization corresponding to the regu-
larization parameter lambda = 10−3 :
f = fil fac (s,1e-3);
Algorithm:
The filter factors are computed by means of the definitions from Section 2.7. See
also [1] for more details.
Limitations:
The routine fil fac cannot be used to compute filter factors for the iterative meth-
ods; these filter factors must be computed by the respective routines. Neither can
fil fac be used to compute filter factors for MTSVD or maximum entropy regulariza-
tion. For truncated TLS, the small filter factors are not computed accurately.
References:
1. P. C. Hansen, Rank-Deficient and Discrete Ill-Posed Problems. Numerical As-
pects of Linear Inversion, SIAM, Philadelphia, 1997.
64 foxgood
foxgood
Purpose:
Severely ill-posed test problem.
Synopsis:
[A,b,x] = foxgood (n)
Description:
Discretization of a Fredholm integral equation of the first kind (2.1) with both
integration intervals equal to [0, 1], with kernel K and right-hand side g given by
1 1³ 3
´
K(s, t) = (s2 + t2 ) 2 , g(s) = (1 + s2 ) 2 − s3 ,
3
and with the solution f = t. The problem was first used by Fox & Goodwin.
This is an artificial discrete severely ill-posed problem which does not satisfy the
discrete Picard condition for the smaller singular values. The matrix A is n × n.
References:
1. C. T. H. Baker, The Numerical Treatment of Integral Equations, Clarendon
Press, Oxford, 1977; p. 665.
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 10 20 30 40 50 60 70 80 90 100
gcv 65
gcv
Purpose:
Plot the GCV function and find its minimum.
Synopsis:
[reg min,G,reg param] = gcv (U,s,b,method)
[reg min,G,reg param] = gcv (U,sm,b,method) , where sm = [sigma,mu]
Description:
Plots the generalized cross-validation (GCV) function
kA x − bk22
G=
(trace (Im − A AI ))2
as a function of the regularization parameter reg param, depending on the method

chosen. Here, AI is a matrix which produces the regularized solution x when multi-
plied with the right-hand side b, i.e., x = AI b. The following methods are allowed:
If any output arguments are specified, then the minimum of the GCV function is
identified and the corresponding regularization parameter reg min is returned.
Examples:
Generate a test problem and use gcv to compute the optimal regularization pa-
rameter lambda for Tikhonov regularization in standard form:
[A,b,x] = phillips (n); b = b + 1e-3∗randn (size(b)); [U,s,V] = csvd (A);
lambda = gcv (U,s,b); x lambda = tikhonov (U,s,V,b,lambda);
plot (1:n,x,1:n,x lambda,’o’)
Algorithm:
For Tikhonov regularization and damped SVD/GSVD, 200 logarithmically dis-
tributed regularization parameters are generated, and G is plotted for these values.
Then the minimizer of the GCV function is computed via fmin, using the minimizer
of G as initial guess. For truncated SVD/GSVD/TLS, G is computed and plotted for
all valid values of the discrete regularization parameter.
66 gcv
Limitations:
gcv cannot be used to compute the GCV function for the iterative regularization
methods. Use instead the filter factors and Eq. (2.63). If cgls or lsqr b is used with
reorthogonalization, then G can be approximated by (rho./(m − [1 : k]0 )).^2, where k
is the number of iterations.
Diagnostics:
The number of points on the GCV curve for Tikhonov regularization and damped
SVD/GSVD is determined by the parameter npoints which can easily be changed by
the user to give a finer resolution.
See also:
discrep, l curve, ncp, quasiopt
References:
1. G. Wahba, Spline Models for Observational Data, CBMS-NSF Regional Con-
ference Series in Applied Mathematics, Vol. 59, SIAM, Philadelphia, 1990.
gen form 67
gen form
Purpose:
Transformation of a standard-form solution back to the general-form setting.
Synopsis:
x = gen form (L p,x s,A,b,K,M) (method 1)
x = gen form (L p,x s,x 0) (method 2)
Description:
Transforms the standard-form solution x s back to the required solution x in the
general-form setting. Notice that x s may have more than one column. Two methods
are available, and the distinction between them is controlled by the number of input
parameters to gen form.
In method 1, described in [1], the transformation takes the form
x = L p x s + K M (b − A L p x s) ,
where L p is the pseudoinverse of the matrix L, the matrix K has columns in the
null space of L, and M = To−1 HoT , cf. (2.25). In method 2, described in [2,3,4], the
transformation takes the form
x s = L px s + x 0 ,
where now L p is the A-weighted generalized inverse of L, and x 0 is the component of
x in the null space of L (notice that this component is independent of x s); see (2.32).
Usually, the transformation from general form into standard form is performed by
means of routine std form.
Notice that both gen form and std form are available for pedagogical reasons only—
usually it is more efficient to build the transformations into the solution process, such
as in the iterative methods pcgls, plsqr b, and pnu.
Examples:
Transform a general-form problem into standard form, produce 10 TSVD solu-
tions, and transform back again, using method 1; then compare with the mathemat-
ically identical TGSVD solutions:
[A s,b s,L p,K,M] = std form (A,L,b); [U,s,V] = csvd (A s);
X s = tsvd (U,s,V,b s,1:10); X = gen form (L p,X s,A,b,K,M);
[U1,V1,sm,X1] = gsvd (A,L); XX = tgsvd (U1,sm,X1,b,1:10); norm (X−XX)
Algorithm:
The algorithms are described in details in refs. [1,2,3,4].
See also:
std form
68 gen form
References:
BIT 17 (1977), 134–145.
2. L. Eldén, A weighted pseudoinverse, generalized singular values, and constrained
least squares problems, BIT 22 (1982), 487–502.
3. M. Hanke, Regularization with differential operators. An iterative approach, J.
Numer. Funct. Anal. Optim. 13 (1992), 523–540.
gen hh 69
gen hh
Purpose:
Generate a Householder transformation.
Synopsis:
[x1,beta,v] = gen hh (x)
Description:
Given a vector x, gen hh computes the scalar beta and the vector v determining a
Householder transformation
H = In − beta v vT
such that H x = ± kxk2 e1 , where e1 is the first unit vector. The quantity x1 returned
by gen hh is the first element of H x.
Examples:
The very first step of the bidiagonalization of a matrix A looks as follows:
[A(1,1),beta,v] = gen hh (A(1:m,1)); A(2:m,1) = zeros (m-1,1);
A(1:m,2:n) = app hh (A(1:m,2:n),beta,v);
See also:
app hh
70 get l
get l
Purpose:
Compute discrete derivative operators.
Synopsis:
[L,W] = get l (n,d)
Description:
Computes the discrete approximation L to the derivative operator of order d on a
regular grid with n points, i.e. the matrix L is (n − d) × n. In particular, for d = 1
and d = 2, L has the form
   
1 −1 1 −2 1
 1 −1   1 −2 1 
L=  .. .. 
 and L=  .. .. .. ,

. . . . .
1 −1 1 −2 1
respectively. The matrix L is stored as a sparse matrix.

If required, an orthonormal basis represented by the columns of the matrix W is
also computed. The matrix W is used in the “preconditioned” iterative methods.
Algorithm:
The matrix W is computed by first generating the “trivial” basis vectors (1 1 . . . 1)T ,
(1 2 . . . n)T , etc., and then orthonormalizing these by means of MGS.
gravity 71
gravity
Purpose:
Test problem: one-dimensional gravity surveying model problem.
Synopsis:
[A,b,x] = gravity (n,example,a,b,d)
Description:
Discretization of a one-dimensional model problem in gravity surveying, in which
a mass distribution f (t) is located at depth d, while the vertical component of the
gravity field g(s) is measured at the surface. The resulting problem is a first-kind
Fredholm integral equation with kernel
¡ ¢−3/2
K(s, t) = d d2 + (s − t)2 .
The following three examples are implemented (example = 1 is default):

1. f (t) = sin(π t) + 0.5 sin(2π t),
2. f (t) = piecewise linear function,
3. f (t) = piecewise constant function.
The t integration interval is fixed to [0,1], while the s integration interval [a,b] can be
specified by the user. The default interval is [0,1], leading to a symmetric Toeplitz
matrix. The parameter d is the depth at which the mass distribution is located, and
the default value is d = 0.25. The larger the d, the faster the decay of the singular
values.
Algorithm:
The problem is discretized by the midpoint quadrature rule with n points, leading
to the matrix A and the vector x. Then the right-hand side is computed as b = A ∗ x.
References:
1. P. C. Hansen, Deconvolution and regularization with Toeplitz matrices, Numer-
ical Algorithms 29 (2002), 323–378.
example = 1 example = 2 example = 3

2 2 2
1.5 1.5 1.5
1 1 1
0.5 0.5 0.5
0 0 0
0 50 100 0 50 100 0 50 100
72 heat
heat
Purpose:
Test problem: inverse heat equation.
Synopsis:
[A,b,x] = heat (n,kappa)
Description:
The inverse heat equation used here is a Volterra integral equation of the first kind
with [0, 1] as integration interval. The kernel is K(s, t) = k(s − t) with
µ ¶
t−3/2 1
k(t) = √ exp − .
2 kappa π 4 kappa2 t
Here, the parameter kappa controls the ill-conditioning of the matrix A:
kappa = 5 gives a well-conditioned matrix,
kappa = 1 gives an ill-conditioned matrix.
The default is kappa = 1.
Algorithm:
The integral equation is discretized by means of simple collocation and the mid-
point rule with n points, cf. [1,2]. An exact solution x is constructed, and then the
right-hand side b is produced as b = A x.
References:
1. A. S. Carasso, Determining surface temperatures from interior observations,
SIAM J. Appl. Math. 42 (1982), 558–574.
2. L. Eldén, The numerical solution of a non-characteristic Cauchy problem for
a parabolic equation; in P. Deuflhart & E. Hairer (Eds.), Numerical Treatment
of Inverse Problems in Differential and Integral Equations, Birkhäuser, Boston,
1983.
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 10 20 30 40 50 60 70 80 90 100
i laplace 73
i laplace
Purpose:
Test problem: inverse Laplace transformation.
Synopsis:
[A,b,x] = i laplace (n,example)
Description:
Discretization of the inverse Laplace transformation (a Fredholm integral equation
of the first kind (2.1)) by means of Gauss-Laguerre quadrature. The kernel K is given
by
K(s, t) = exp (−s t),
and both integration intervals are [0, ∞). The following examples are implemented,
where f denotes the solution and g denotes the corresponding right-hand side:
1
example = 1: f (t) = exp (−t/2), g(s) = s+1/2
1 1
example = 2: f (t) = 1 − exp (−t/2), g(s) = s − s+1/2
2
example = 3: f (t) = t2 exp (−t/2), g(s) = (s+1/2)3
½
0, t ≤ 2 exp(−2s)
example = 4: f (t) = , g(s) =
1, t > 2 s
All four examples are from [1]. The size of the matrix A is n × n.
References:
1. J. M. Varah, Pitfalls in the numerical solution of linear ill-posed problems, SIAM
J. Sci. Stat. Comput. 4 (1983), 164–176.
2.5
example = 1
example = 2
2
example = 3
example = 4
1.5
0.5
0
0 10 20 30 40 50 60 70 80 90 100
74 lagrange
lagrange
Purpose:
Plot the Lagrange function for Tikhonov regularization.
Synopsis:
[La,dLa,lambda0] = lagrange (U,s,b,more)
[La,dLa,lambda0] = lagrange (U,sm,b,more) , sm = [sigma,mu]
Description:
Plots the Lagrange function for Tikhonov regularization,
La = kA xλ − bk22 + λ2 kL xλ k22 ,
and its first derivative dLa = dLa/dλ, versus λ. Here, xλ is the Tikhonov regularized
solution. If nargin = 4, then the norms kA xλ − bk2 and kL xλ k2 are also plotted. The
routine returns La, dLa, and an approximate value lambda0 of λ for which dLa has its
minimum.
See also:
l curve, tikhonov
lanc b 75
lanc b
Purpose:
Lanczos bidiagonalization.
Synopsis:
B k = lanc b (A,p,k,reorth)
[U,B k,V] = lanc b (A,p,k,reorth)
Description:
Performs k steps of the Lanczos bidiagonalization process with starting vector p,
producing a lower bidiagonal (k + 1) × k matrix B k,
 
b11
 b21 b22 
 .. 
 . 
B k= b32  , such that AV = UB k .
 .. 
 . bkk 
bk+1,k
The matrix B k is stored as a sparse matrix. The matrices U ∈ Rm×(k+1) and V ∈
Rn×k consist of the left and right Lanczos vectors (which are orthonormal in exact
arithmetic). Reorthogonalization is controlled by means of reorth as follows:
reorth = 1 : reorthogonalization by means of MGS
reorth = 2 : Householder reorthogonalization.
Examples:
Perform 10 steps of Lanczos bidiagonalization without reorthogonalization, then
compute the singular values of the 11 × 10 bidiagonal matrix and compare with the
10 largest singular values of the coefficient matrix A:
B k = lanc b (A,b,10); s k = csvd (B k); s = svd (A); [s k,s(1:10)]
Algorithm:
The algorithm is a straight-forward implementation of the Lanczos bidiagonaliza-
tion process as described in, e.g., [1].
See also:
bidiag, bsvd, lsqr b, plsqr b
References:
1. G. H. Golub & C. F. Van Loan, Matrix Computations, 2. Ed., Johns Hopkins,
Baltimore, 1989; Section 9.3.4.
76 lsolve
lsolve
Purpose:
Utility routine for “preconditioned” iterative methods.
Synopsis:
x = lsolve (L,y,W,T)
Description:
Computes the vector
x = L†A y ,
where L†A is the A-weighted generalized inverse of L. Here, L is a p×n matrix, W holds
a basis for the null space of L, and T is a utility matrix which should be computed
by routine pinit. Alternatively, L is square and dense, and W and T are not needed.
Notice that x and y may be matrices, in which case
x(:, i) = L†A y(:, i) .
Algorithm:
The vector x is computed by means of the following algorithm from [2, §2.3.2],
based on [1]:
x̂ ← L(: , 1: p)−1 y
µ ¶
x̂
x ← − W T(: , 1: p) x̂ .
0
See also:
ltsolve, pcgls, pinit, plsqr b, pnu
References:
lsqi 77
lsqi
Purpose:
Least squares minimization with a quadratic inequality constraint.
Synopsis:
[x alpha,lambda] = lsqi (U,s,V,b,alpha,x 0)
[x alpha,lambda] = lsqi (U,sm,X,b,alpha,x 0) , where sm = [sigma,mu]
Description:
Solves the following least squares minimization problems with a quadratic inequal-
ity constraint:
min kA x − bk2 subject to kx − x 0k2 ≤ alpha
min kA x − bk2 subject to kL (x − x 0)k2 ≤ alpha
where x 0 is an initial guess of the solution, and alpha is a positive constant. The
routine requires either the compact SVD of A saved as U, s, and V, or part of the GSVD
of (A, L) saved as U, sm, and X. The regularization parameter lambda corresponding
to the solution x alpha is also returned.
If alpha is a vector, then x alpha is a matrix such that
x alpha = [x alpha(1), x alpha(2), . . . ] .
If x 0 is not specified, x 0 = 0 is used.
The “opposite” problem, namely that of minimizing the solution (semi)norm sub-
ject to an upper bound on the residual norm, is treated by routine discrep.
Examples:
Generate a “noisy” test problem with a know exact solution x. Then compute two
regularized solutions x1 and x2 with kx1k2 = 1.1 kxk2 and kL x1k2 = 1.1 kL xk2 , where
L is an approximation to the second derivative operator:
[A,b,x] = shaw (32); b = b + 1e-3∗randn (size(b)); [U,s,V] = csvd (A);
L = get l (32,2); [UU,sm,XX] = cgsvd (A,L);
x1 = lsqi (U,s,V,b,1.05∗norm(x));
x2 = lsqi (UU,sm,XX,b,1.05∗norm(L∗x));
plot ([x,x1,x2])
Algorithm:
The algorithm uses Newton iteration with a Hebden rational model (implemented
as routine heb new) to solve the secular equation kL (xλ − x 0)k2 = alpha, cf. [1].
The initial guess of λ is computed as described in [1]. Although it is derived for the
case x 0 = 0, the idea is also applicable for x 0 6= 0 because the initial λ is almost
unaffected by x 0 (since kx0 k2 À kx 0k2 , where x0 is the unregularized solution).
78 lsqi
Diagnostics:
In rare cases the iteration in heb new converges very slowly; then try to increase
the maximum number of allowed iterations in heb new or a slightly different alpha.
See also:
hebnew, discrep
References:
1. T. F. Chan, J. Olkin & D. W. Cooley, Solving quadratically constrained least
squares using black box unconstrained solvers, BIT 32 (1992), 481–495.
lsqr b 79
lsqr b
Purpose:
Solution of least squares problems by the LSQR algorithm based on Lanczos bidi-
agonalization.
Synopsis:
[X,rho,eta,F] = lsqr b (A,b,k,reorth,s)
Description:
Performs k steps of the LSQR Lanczos bidiagonalization algorithm (due to Paige
& Saunders) applied to the system
min kA x − bk2 .
The routine returns all k solutions, stored as columns of the matrix X. The solution
norms and residual norms are returned in eta and rho, respectively.
Reorthogonalization of the Lanczos vectors is controlled by means of reorth as
follows:
If the singular values s of A are also provided, then lsqr b computes the filter
factors associated with each step and stores them columnwise in the array F.
A “preconditioned” version of LSQR for the general-form problem where one min-
imizes kL xk2 instead of kxk2 is available in the function plsqr b.
Examples:
Perform 25 LSQR iterations without reorthogonalization, and plot the correspond-
ing discrete L-curve:
[X,rho,eta] = lsqr b (A,b,25); plot lc (rho,eta,’o’);
Algorithm:
lsqr b is a straightforward implementation of the algorithm described in [1]. The
original algorithm also includes the possibility for adding Tikhonov regularization
with a fixed regularization parameter to the least squares problem, but for simplicity
this feature is not included in Regularization Tools.
See also:
art, bidiag, cgls, lanc b, mr2, plsqr b, rrgmres, splsqr
References:
1. C. C. Paige & M. A. Saunders, LSQR: an algorithm for sparse linear equations
and sparse least squares, ACM Trans. Math. Software 8 (1982), 43–71.
80 ltsolve
ltsolve
Purpose:
Utility routine for “preconditioned” iterative methods.
Synopsis:
x = ltsolve (L,y,W,T)
Description:
Computes the vector
x = (L†A )T y ,
where L†A is the A-weighted generalized inverse of L. Here, L is a p×n matrix, W holds
a basis for the null space of L, and T is a utility matrix which should be computed
by routine pinit. Alternatively, L is square and dense, and W and T are not needed.
If W and T are not specified, then instead the vector
x = (L(: , 1: p)−1 )T y(1: p)
is computed.
Notice that x and y may be matrices, in which case
x(:, i) = (L†A )T y(:, i) or x(: , i) = (L(: , 1: p)−1 )T y(1: p, i) .
Algorithm:
The following algorithm from [2, §2.3.2] (based on [1]) is used. If W and T are
specified, then y is first updated by means of the relation
y ← y(1: p) − T(: , 1: p)T WT y ,
otherwise the update y ← y(1: p) is computed. The latter option is used in the
“preconditioned” iterative methods. Then x is computed as (L(: , 1: p)−1 )T y.
See also:
lsolve, pcgls, pinit, plsqr b, pnu
References:
l corner 81
l corner
Purpose:
Locate the “corner” of the L-curve.
Synopsis:
[reg c,rho c,eta c] = l corner (rho,eta,reg param)
[reg c,rho c,eta c] = l corner (rho,eta,reg param,U,s,b,method,M)
[reg c,rho c,eta c] = l corner (rho,eta,reg param,U,sm,b,method,M) , sm = [sigma,mu]
Description:
l corner locates the corner of the L-curve, given by rho and eta, in log-log scale. It
is assumed that corresponding values of kA x − bk2 , kL xk2 , and the regularization
parameter are stored in the arrays rho, eta, and reg param, respectively (such as the
output from routine l curve). If nargin = 3 then no particular method is assumed, and
if nargin = 2 then it is assumed that reg param = 1:length (rho).
If nargin ≥ 6, then the following regularization methods are allowed:
method = ’mtsvd’ : modified TSVD
method = ’Tikh’ : Tikhonov regularization
method = ’tsvd’ : truncated SVD or GSVD.
If no method is specified, ’Tikh’ is default.
An eighth argument M specifies an upper bound for eta, below which the corner
should be found.
Examples:
l corner is invoked by l curve if any output arguments are specified to the latter
routine. However, l corner can also be use as a stand-alone routine. In the following
example, rho and eta come from the LSQR algorithm with 15 iterations—so the
regularization parameters in reg param are 1:15—and l corner is used to compute the
“corner” of the L-curve associated with rho and eta:
[A,b,x] = shaw (32); b = b + 1e-3∗randn (size(b)); [X,rho,eta] = lsqr b (A,b,15,1);
[reg c,rho c,eta c] = l corner (rho,eta,1:15);
Algorithm:
If method is Tikhonov regularization or DSVD/DGSVD, then the L-curve is dif-
ferentiable, and it is straightforward to compute the curvature and find the point on
the L-curve with maximum curvature, which is then defined as the “corner”. For the
remaining methods, the L-curve is discrete and a different approach is used.
82 l corner
• If the Spline Toolbox is available, then we fit a 2D spline curve to the data
points, compute the point on the spline curve with maximum curvature, and
then define the “corner” of the discrete L-curve as the data point closest to the
point on the spline curve with maximum curvature [3].
• Otherwise we use the adaptive pruning algorithm from [2] implemented in corner.
If the curvature of the L-curve is negative everywhere, then the leftmost point on the
L-curve is taken as the “corner”. This choice is the most natural in connection with
the L-curve criterion for choosing the regularization parameter.
Diagnostics:
If the user decides to always use the adaptive pruning algorithm, even if the Spline
Toolbox is available, then change line 41 in l corner to alwayscorner = 1; (the default
of this logical variable is 0).
For discrete L-curves, if the spline approach is used then the number of data points
in rho and eta must be greater than four, the order of the fitting 2D spline curve. This
algorithm may sometimes mistake a fine-grained local “corner” for the desired corner.
In this case, the user may wish to experiment with the default parameters deg (the
degree of a local smoothing polynomial applied before fitting the 2D spline curve),
q (the half-width of the local smoothing interval in which polynomial smoothing is
applied), and order (the order of the fitting 2D spline curve).
The user may also wish to increase the parameter s thr: for TSVD, TGSVD and
MTSVD, values of rho and eta corresponding the singular values s below s thr are not
used in the search for the “corner”.
See also:
corner, l curve
External routines required:

The following ten routines from the Spline Toolbox [1] are required by l corner, if
the spline approach is used:
fnder, ppbrk, ppcut, ppmak, ppual, sp2pp, sorted, spbrk, spmak, sprpp.
References:
1. C. de Boor, Spline Toolbox, Version 1.1, The Mathworks, MA, 1992.
2. P. C. Hansen, T. K. Jensen & G. Rodriguez, An adaptive pruning algorithm for
3. P. C. Hansen & D. P. O’Leary, The use of the L-curve in the regularization of
discrete ill-posed problems, SIAM J. Sci. Comput. 14 (1993), 1487–1503.
l curve 83
l curve
Purpose:
Plot the L-curve and find its “corner”.
Synopsis:
[reg corner,rho,eta,reg param] = l curve (U,s,b,method)
[reg corner,rho,eta,reg param] = l curve (U,sm,b,method) , where sm = [sigma,mu]
[reg corner,rho,eta,reg param] = l curve (U,s,b,method,L,V)
Description:
Plots the L-shaped curve of eta, the solution norm kxk2 or seminorm kL xk2 ,
versus rho, the residual norm kA x − bk2 , for the following methods:
method = ’Tikh’ : Tikhonov regularization (solid line)
method = ’tsvd’ : truncated SVD or GSVD (o markers)
method = ’dsvd’ : damped SVD or GSVD (dotted line)
method = ’mtsvd’ : modified TSVD (x markers)
The corresponding regularization parameters are returned in reg param. If no method
is specified, ’Tikh’ is default. For other methods, use routine plot lc. Note that ’Tikh’,
’tsvd’ and ’dsvd’ require either U and s (standard-form regularization) or U and sm
(general-form regularization), while ’mtsvd’ requires U and s as well as L and V.
If any output arguments are specified, then the corner of the L-curve is identi-
fied (by means of routine l corner) and the corresponding regularization parameter
reg corner is returned as well as printed on the plot. This is the L-curve criterion
for choosing the regularization parameter. Use l corner if an upper bound on eta is
required when finding the L-curve’s “corner”.
Examples:
Plot the L-curve for Tikhonov regularization and compute the regularization pa-
rameter reg corner corresponding to the “corner”:
reg corner = l curve (U,s,b);
Algorithm:
The algorithm for computing the “corner” of the L-curve in log-log scale is de-
scribed in [2], where the existence of the “corner” is also discussed. A survey of the
L-curve criterion for choosing the regularization parameter, including a discussion of
the limitations of the method, can be found in [1].
Diagnostics:
See the documentation for l corner for diagnostics.
84 l curve
See also:
corner, discrep, gcv, l corner, ncp, plot lc, quasiopt
References:
1. P. C. Hansen, The L-curve and its use in the numerical treatment of inverse
problems; in P. Johnston (Ed.), Computational Inverse Problems in Electrocar-
diography, WIT Press, Southampton, 2001; pp. 119–142.
2. P. C. Hansen & D. P. O’Leary, The use of the L-curve in the regularization of
maxent 85
maxent
Purpose:
Maximum entropy regularization.
Synopsis:
[x lambda,rho,eta] = maxent (A,b,lambda,w,x0)
Description:
Computes the regularized maximum entropy solution x lambda which minimizes
Pn
kA x − bk22 + lambda2 i=1 xi log(wi xi ) ,
Pn
where − i=1 xi log(wi xi ) is the entropy of the vector x, and w = (w1 , . . . , wn )T is a
vector of weights. If no weights are specified, unit weights are used. A starting vector
x0 for the iterative process can be specified, otherwise a constant starting vector is
used. If lambda is a vector, then x lambda is a matrix such that
Examples:
Compare maximum entropy regularization with Tikhonov regularization:
lambda = logspace (-1,-4,12);
X tikh = tikhonov (U,s,V,b,lambda);
X ment = maxent (A,b,lambda);
for i=1:12
dif(i,1) = norm(x−X tikh(:,i));
dif(i,2) = norm(x−X ment(:,i));
dif(i,3) = norm(X tikh−X ment(:,i));
end
loglog (lambda,dif/norm(x))
Algorithm:
Routine maxent uses a nonlinear conjugate gradient algorithm [1, §4.1] with inexact
line search and a step-length control which ensures that all elements of the iteration
vector are positive. For more details, see §2.7.5.
Diagnostics:
Slow convergence may occur for small values of the regularization parameter lambda.
References:
1. R. Fletcher, Practical Methods of Optimization. Vol. 1, Unconstrained Opti-
mization, Wiley, Chichester, 1980.
86 mr2
mr2
Purpose:
Solution of symmetric indefinite problems by the MR-II algorithm.
Synopsis:
[X,rho,eta] = mr2 (A,b,k,reorth)
Description:
MR-II is a variant of the MINRES algorithm for symmetric indefinite linear sys-
tems A x = b, with starting vector A b (instead of b as in MINRES). The function
returns all k iterates, stored as the columns of the matrix X. The solution norm and
residual norm are returned in eta and rho, respectively.
Reorthogonalization is controlled by means of reorth as follows:
A “preconditioned” version of MR-II for the general-form problem where one
minimizes kL xk2 instead of kxk2 is available in the function pmr2.
Examples:
Perform 25 MR-II iterations without reorthogonalization, and plot the correspond-
ing discrete L-curve:
[X,rho,eta] = mr2 (A,b,25); plot lc (rho,eta,’o’);
Algorithm:
mr2 is a straightforward implementation of the MR-II algorithm described in [1].
An analysis of MR-II as a general regularization routine can be found in [2].
See also:
cgls, lsqr, pmr2, rrgmres, splsqr
References:
1. M. Hanke, Conjugate Gradient Methods for Ill-Posed Problems, Longman Sci-
entific and Technical, Essex, 1995.
2. T. K. Jensen & P. C. Hansen, Iterative regularization with minimum-residual
methods, BIT 47 (2007), 103-120.
mtsvd 87
mtsvd
Purpose:
Regularization by means of modified TSVD.
Synopsis:
[x k,rho,eta] = mtsvd (U,s,V,b,k,L)
Description:
The modified TSVD solution is defined as the solution to the problem
min kL xk2 subject to kAk x − bk2 = min ,
where Ak is the truncated SVD of the matrix A. The MTSVD solution is then given
by µ ¶
ξk
xk=V .
ξ0
T
Here, ξk defines the usual TSVD solution ξk = (diag(s(1 : k)))−1 U(:, 1 : k) b, and ξ0
is chosen so as to minimize the seminorm kL xk k2 . This leads to choosing ξ0 as the
solution to the following least squares problem:
min k(L V(:, k + 1 : n)) ξ − L V(:, 1 : k) ξk k2 .
The truncation parameter must satisfy k > n − p, where L ∈ Rp×n The solution and
residual norms are returned in eta and rho.
If k is a vector, then x k is a matrix such that
x k = [ x k(1), x k(2), . . . ] .
Examples:
Compute the MTSVD solutions corresponding to k = 5:12:
X = mtsvd (U,s,V,b,5:12,L);
Algorithm:
The algorithm computes one QR factorization of the matrix L V(:, kmin + 1 : n),
where kmin = min (k). For more details, cf. [1].
See also:
discrep, lsqi, tgsvd, tikhonov, tsvd
References:
1. P. C. Hansen, T. Sekii & H. Shibahashi, The modified truncated-SVD method
for regularization in general form, SIAM J. Sci. Stat. Comput. 13 (1992), 1142-
1150.
88 ncp
ncp
Purpose:
Plot the normalized cumulative periodograms (NCPs) and find the one closest to
a straight line.
Synopsis:
[reg min,dist,reg param] = ncp (U,s,b,method)
[reg min,dist,reg param] = ncp (U,sm,b,method) , where sm = [sigma,mu]
Description:
The normalized cumulative periodogram (NCP) of the residual vector measures
the frequency contents of the residual, and the closer the NCP is to a straight line,
the more the residual vector resembles white noise. This is used to choose the regu-
larization parameter for the following methods:
Examples:
Generate a test problem and use ncp to compute the optimal regularization pa-
rameter lambda for Tikhonov regularization in standard form:
[A,b,x] = phillips (128); b = b + 1e-2∗randn (size(b)); [U,s,V] = csvd (A);
lambda = ncp (U,s,b); x lambda = tikhonov (U,s,V,b,lambda);
plot (1:n,x,1:n,x lambda,’o’)
Algorithm:
The NCP c of a residual vector r is defined as the vector computed by the com-
mands
p = abs(fft(r)).^2; c = cumsum(p(2:q))/sum(p(2:q));
where q = floor(length(r)/2). The function returns the regularization parameter
reg min for which the NCP resembles most a straight line, measure by the 2-norm of
difference between the vectors c and v = (1:q)’/q. The function also returns the
distances dist and the corresponding regularization parameters reg param which can
be plotted by means of plot (reg param,dist).
Limitations:
ncp cannot be used to compute the NCPs for the iterative regularization methods.
Use the above technique and stop when the NCP is “close enough” to a straight line;
see [1] for details.
ncp 89
Diagnostics:
The number of initial NCPs for for Tikhonov regularization and damped SVD/GSVD
is determined by the parameter npoints, while the number of plotted NCPs is deter-
mined by the parameter nNCPs; both can easily be changed by the user.
See also:
discrep, gcv, l curve, quasiopt.
References:
1. P. C. Hansen, M. Kilmer & R. H. Kjeldsen, Exploiting residual information in
the parameter choice for discrete ill-posed problems, BIT 46 (2006), 41–59.
90 nu
nu
Purpose:
Brakhage’s ν-method.
Synopsis:
[X,rho,eta,F] = nu (A,b,k,nu,s)
Description:
Performs k steps of Brakhage’s ν-method for the problem
min kA x − bk2 .
norms and residual norms are returned in eta and rho, respectively.
If nu is not specified, nu = .5 is the default value which gives the Chebychev
method of Nemirovskii and Polyak.
If the singular values s are also provided, nu computes the filter factors associated
with each step and stores them columnwise in the array F.
A “preconditioned” version of the ν-method for the general-form problem where
one minimizes kL xk2 instead of kxk2 in each step is implemented in routine pnu.
Examples:
Perform 50 iterations of the ν-method and plot the corresponding L-curve:
[X,rho,eta] = nu (A,b,50); plot lc (rho,eta,’o’);
Algorithm:
nu is a straightforward implementation of the algorithm described in [1]. The
iteration converges only if kAk2 < 1. Therefore, A and b are scaled before the iteration
begins with a scaling factor given by 0.99/kBk2 , where B is a bidiagonal matrix
obtained from a few steps of the Lanczos bidiagonalization process applied to A with
starting vector b. Hence, kBk2 is a good approximation to kAk2 .
See also:
art, cgls, lsqr b, pnu
References:
1. H. Brakhage, On ill-posed problems and the method of conjugate gradients; in H.
W. Engl & G. W. Groetsch (Eds.), Inverse and Ill-Posed Problems, Academic
Press, New York, 1987.
parallax 91
parallax
Purpose:
Stellar parallax problem with 28 fixed, real observations.
Synopsis:
[A,b] = parallax (n)
Description:
The underlying problem is a Fredholm integral equation of the first kind with
kernel Ã µ ¶2 !
1 1 s−t
K(s, t) = √ exp −
σ 2π 2 σ
with σ = 0.014234, and it is discretized by means of a Galerkin method with n
orthonormal basis functions. The right-hand side b consists of a measured distribution
function of stellar parallaxes from [1], and its length is fixed at 26; i.e., the matrix A is
26×n. The exact solution, which represents the true distribution of stellar parallaxes,
is unknown.
References:
1. W. M. Smart, Stellar Dynamics, Cambridge University Press, Cambridge, 1938;
p. 30.
92 pcgls
pcgls
Purpose:
“Preconditioned” conjugate gradient algorithm applied implicitly to the normal
equations.
Synopsis:
[X,rho,eta,F] = pcgls (A,L,W,b,k,sm)
Description:
Performs k steps of the “preconditioned” conjugate gradient algorithm applied
implicitly to the normal equations AT A x = AT b with “preconditioner” L†A (L†A )T and
with starting vector x0 (2.30) computed by means of Algorithm (2.34). Here, L†A is
the A-weighted generalized inverse of L. Notice that the matrix W holding a basis for
the null space of L must also be specified.
seminorms and the residual norms are returned in eta and rho, respectively.
If the c, s-pairs sm of the GSVD of (A, L) are also provided, then pcgls computes
the filter factors associated with each step and stores them in the array F.
Reorthogonalization of the normal equation residual vectors AT (A X(: , i) − b),
i = 1, . . . , k is controlled by means of reorth as follows:
A simpler version of pcgls for the case L = In is implemented in routine cgls.
Examples:
Choose minimization of the second derivative as side constraint, and perform 20
iterations of the “preconditioned” conjugate gradient method:
[L,W] = get l (n,2); X = pcgls (A,L,W,b,20); mesh (X)
Algorithm:
pcgls is a straightforward implementation of the algorithm described in [1], with
the necessary modifications to the case L 6= In from [2]. The computation of the filter
factors is also described in [2].
See also:
cgls, plsqr b, pmr2, pnu, prrgmres
pcgls 93
References:
1. Å. Björck, Numerical Methods for Least Squares Problems, SIAM, Philadelphia,
1996.
94 phillips
phillips
Purpose:
Phillips’s “famous” test problem.
Synopsis:
[A,b,x] = phillips (n)
Description:
Discretization of the “famous” Fredholm integral equation of the first kind (2.1)
deviced by D. L. Phillips [1]. Define the function
½ ¡ ¢
1 + cos π3x , |x| < 3
φ(x) = .
0, |x| ≥ 3
Then the kernel K, the solution f , and the right-hand side g are given by:
K(s, t) = φ(s − t)
f (t) = φ(t)
µ ³ π s ´¶ µ ¶
1 9 π |s|
g(s) = (6 − |s|) 1 + cos + sin .
2 3 2π 3
Both integration intervals are [−6, 6]. The size of the matrix A is n × n.
Limitations:
The order n must be a multiple of 4.
References:
1. D. L. Phillips, A technique for the numerical solution of certain integral equa-
tions of the first kind, J. ACM 9 (1962), 84–97.
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 10 20 30 40 50 60 70 80 90 100
picard 95
picard
Purpose:
Visual inspection of the Picard condition.
Synopsis:
xi = picard (U,s,b,d)
xi = picard (U,sm,b,d) , where sm = [sigma,mu]
Description:
picard plots the singular values s, the absolute values of the Fourier coefficients,
T
|UT b|, and a (possible smoothed) curve of the solution coefficients |(U b)./s|. If the
c, s-values sm = [sigma,mu] are specified, where γ = sigma./mu are the generalized
T
singular values, then the routine plots γ, |UT b|, and (smoothed) |(U b)./γ|.
The smoothing is a geometric mean over 2d + 1 points. If nargin = 3 then d = 0
(i.e., no smoothing).
The quantities plotted by picard are useful for visually inspecting whether the
discrete Picard condition is satisfied for the given problem: for the large singular
values s, or the large generalized singular values γ, the Fourier coefficients |UT b|
should decay at least as fast as the s and γ, respectively.
Examples:
Generate a test problem, add noise to the right-hand side, and use picard to check
the discrete Picard condition visually:
[A,b,x] = shaw (32); [U,s,V] = csvd (A); b = b + 1e-3∗randn (size(b));
picard (U,s,b);
See also:
l curve
References:
1. P. C. Hansen, The discrete Picard condition for discrete ill-posed problems, BIT
30 (1990), 658–672.
96 pinit
pinit
Purpose:
Utility initialization-procedure for the “preconditioned” iterative regularization
methods.
Synopsis:
T = pinit (W,A)
[T,x 0] = pinit (W,A,b)
Description:
pinit is a utility routine used for initialization inside the iterative regularization
methods. Given a matrix W whose columns span the null space of L, pinit computes
a matrix T which is needed in the routines for treating general-form regularization
problems.
If b is also specified then x 0, the component of the regularized solution in the null
space of L, is also computed.
Algorithm:
T and x 0 are computed by the following procedure:
S ← (A W )† , T←SA
x 0 ← WSb .
The Matlab command pinv is used to compute the pseudoinverse (A W )† .
See also:
get l
References:
Numer. Funct. Anal. Optim. 13 (1992), 523-540.
plot lc 97
plot lc
Purpose:
Plot the L-curve.
Synopsis:
plot lc (rho,eta,marker,ps,reg param)
Description:
Plots the L-curve, i.e., the L-shaped curve of the solution norm
eta = kxk2 if ps = 1
eta = kL xk2 if ps = 2
versus the residual norm rho = kA x − bk2 . If ps is not specified, the value ps = 1 is
assumed.
The text string marker is used as marker for the L-curve. If marker is not specified,
the marker ’–’ is used, giving a solid line.
If a fifth argument reg param is present, holding the regularization parameters
corresponding to rho and eta, then some points on the L-curve are identified by their
corresponding regularization parameter.
Diagnostics:
The default number of identified points on the L-curve is 10. It is controlled by
the parameter np inside plot lc.
See also:
l curve
98 plsqr b
plsqr b
Purpose:
“Preconditioned” version of the LSQR Lanczos bidiagonalization algorithm.
Synopsis:
[X,rho,eta,F] = plsqr b (A,L,W,b,k,reorth,sm)
Description:
Performs k steps of the “preconditioned” LSQR Lanczos bidiagonalization algo-
rithm applied to the system min kA x − bk2 with “preconditioner” L†A (L†A )T and with
starting vector x0 (2.30). Here, L†A is the A-weighted generalized inverse of L. Notice
that the matrix W holding a basis for the null space of L must also be specified.
seminorms and the residual norms are returned in eta and rho, respectively.
Reorthogonalization of the Lanczos vectors is controlled by means of reorth
If the c, s-values sm of the GSVD of (A, L) are also provided, then plsqr b computes
the filter factors associated with each step and stores them columnwise in the array F.
A simpler version of plsqr b for the case L = In is implemented in routine lsqr b.
Examples:
Choose minimization of the second derivative as side constraint, perform 20 itera-
tions of the “preconditioned” LSQR algorithm with no reorthogonalization including
computation of the filter factors, and display the filter factors:
[L,W] = get l (n,2); sm = gsvd (A,L);
[X,rho,eta,F] = plsqr b (A,L,W,b,25,0,s); mesh (F), axis (’ij’)
Algorithm:
plsqr b is an implementation of the LSQR algorithm [2] with the necessary modi-
fications for “preconditioning” described in [2].
See also:
lsqr b,pcgls, pmr2, pnu, prrgmres
References:
2. C. C. Paige & M. A. Saunders, LSQR: an algorithm for sparse linear equations
pmr2 99
pmr2
Purpose:
Preconditioned MR-II algorithm for symmetric indefinite problems.
Synopsis:
[X,rho,eta] = pmr2 (A,b,k,reorth)
Description:
This function applies smoothing-norm preconditioning to the MR-II method, which
is a variant the MINRES algorithm for symmetric indefinite linear systems A x = b,
with starting vector A b (instead of b as in MINRES). The function returns all k it-
erates, stored as the columns of the matrix X. The solution norm and residual norm
are returned in eta and rho, respectively.
The preconditioned version seeks to minimize kL xk2 instead of kxk2 . The precon-
ditioner uses two matrices: the matrix L that defines the smoothing norm, and the
matrix N whose columns span the null space of L. It is assumed that L is p × n with
p < n.
Reorthogonalization is controlled by means of reorth as follows:
Examples:
Perform 25 preconditioned MR-II iterations without reorthogonalization, using
the second derivative smoothing norm kL xk2 :
[L,N] = get l (n,2); [X,rho,eta] = pmr2 (A,L,N,b,25);
Algorithm:
The underlying MR-II algorithm is a straightforward implementation of the al-
gorithm described in [1], while the smoothing-norm preconditioning used here is de-
scribed in [2].
See also:
mr2, prrgmres, rrgmres
References:
1. M. Hanke, Conjugate Gradient Methods for Ill-Posed Problems, Longman Sci-
entific and Technical, Essex, 1995.
2. P. C. Hansen and T. K. Jensen, Smoothing-norm preconditioning for regularizing
minimum-residual methods, SIAM J. Matrix Anal. Appl. 29 (2006), 1–14.
100 pnu
pnu
Purpose:
“Preconditioned” version of Brakhage’s ν-method.
Synopsis:
[X,rho,eta,F] = pnu (A,L,W,b,k,nu,sm)
Description:
Performs k steps of the “preconditioned” version of Brakhage’s ν-method applied
to the system min kA x − bk2 with “preconditioner” L†A (L†A )T and with starting vector
x0 (2.30). Here, L†A is the A-weighted generalized inverse of L. Notice that the matrix
W holding a basis for the null space of L must also be specified. The routine returns
all k solutions, stored as columns of the matrix X. The solution seminorms and the
residual norms are returned in eta and rho, respectively.
If nu is not specified, nu = .5 is the default value which gives the “preconditioned”
version of the Chebychev method of Nemirovskii and Polyak.
If the c, s-values sm of the GSVD of (A, L) are also provided, then pnu computes
the filter factors associated with each step and stores them columnwise in the array F.
A simpler version of pnu for the case L = In is implemented in routine nu.
Examples:
Perform 50 iterations of the “preconditioned” ν-method with ν = 0.2 and plot the
corresponding L-curve:
[X,rho,eta] = pnu (A,L,W,b,50,.2); plot lc (rho,eta,’o’);
Algorithm:
pnu is a straightforward implementation of the algorithm described in [1], with
the necessary modifications for “preconditioning” from [2]. The iteration converges
only if kA L†A k2 < 1. Therefore, A and b are scaled with a scaling factor given by
0.99/kBk2 , where B is a bidiagonal matrix obtained from a few steps of the Lanczos
bidiagonalization process applied to A L†A with starting vector b.
See also:
cgls, nu, pcgls
References:
1. H. Brakhage, On ill-posed problems and the method of conjugate gradients; in H.
W. Engl & G. W. Groetsch, Inverse and Ill-Posed Problems, Academic Press,
New York, 1987.
prrgmres 101
prrgmres
Purpose:
Preconditioned range-restricted GMRES algorithm for square systems.
Synopsis:
[X,rho,eta] = prrgmres (A,L,N,b,k)
Description:
The function applies smoothing-norm preconditioning to the RRGMRES algo-
rithm (range-restricted GMRES), which is a variant of the GMRES algorithm for
square linear systems A x = b, with starting vector A b (instead of b as in GMRES).
The function returns all k iterates, stored as the columns of the matrix X. The solution
norm and residual norm are returned in eta and rho, respectively.
For symmetric matrices, use the function pmr2 instead.
The preconditioned version seeks to minimize kL xk2 instead of kxk2 . The precon-
ditioner uses two matrices: the matrix L that defines the smoothing norm, and the
matrix N whose columns span the null space of L. It is assumed that L is p × n with
p < n.
Examples:
Perform 25 preconditioned RRGMRES iterations, using the second derivative
smoothing norm kL xk2 :
[L,N] = get l (n,2); [X,rho,eta] = prrgmres (A,L,N,b,25);
Algorithm:
The underlying RRGMRES algorithm was originally developed in [1] for incon-
sistent systems, and the updating formula for the RRGMRES iterate is given in [2].
The smoothing-norm preconditioning used here is described in [3].
See also:
mr2, pmr2, rrgmres
References:
1. D. Calvetti, B. Lewis & L. Reichel, GMRES-type methods for inconsistent sys-
tems, Lin. Alg. Appl. 316 (2000), 157–169.
2. P. C. Hansen, Regularization Tools Version 4.0 for Matlab 7.3. Numer. Algo.
(2007).
3. P. C. Hansen and T. K. Jensen, Smoothing-norm preconditioning for regularizing
minimum-residual methods, SIAM J. Matrix Anal. Appl. 29 (2006), 1–14.
102 quasiopt
quasiopt
Purpose:
Quasi-optimality criterion for choosing the regularization parameter.
Synopsis:
[reg min,Q,reg param] = quasiopt (U,s,b,method)
[reg min,Q,reg param] = quasiopt (U,sm,b,method) , where sm = [sigma,mu]
Description:
Plots the quasi-optimality function Q, cf. (2.64), for the following methods:
method = ’Tikh’ : Tikhonov regularization (solid line)
method = ’tsvd’ : truncated SVD or GSVD (o markers)
method = ’dsvd’ : damped SVD or GSVD (dotted line)
If no method is specified, ’Tikh’ is default. Returns Q and the corresponding regular-
ization parameters in reg param.
If any output arguments are specified, then the approximate minimum of Q is
identified and the corresponding regularization parameter reg min is returned.
Examples:
Use the quasi-optimality criterion to find the optimal regularization parameter for
Tikhonov regularization in general form:
[A,b,x] = shaw (n); b = b + 1e-3*randn (size(b)); [U,s,V] = svd (A);
lambda opt = quasiopt (U,s,b);
Algorithm:
For Tikhonov regularization and damped SVD/GSVD, 200 logarithmically dis-
tributed regularization parameters are generated, and Q is plotted for these values.
Then the minimizer of the quasi-optimality function is computed via fmin, using the
minimizer of Q as initial guess. For truncated SVD and GSVD, the quasi-optimality
function is discrete and simply becomes Q = (U0 ∗ b)./s and Q = (U0 ∗ b)./(sigma./mu),
respectively, and the minimum is easily found.
Limitations:
The method is not implemented for MTSVD; for this and other discrete regular-
ization methods use the relations Q(k) = kxk − xk−1 k2 and Q(k) = kL (xk − xk−1 )k2 .
See also:
gcv, discrep, l curve, ncp
regudemo 103
regudemo
Purpose:
Tutorial introduction to Regularization Tools.
Synopsis:
regudemo
Description:
This script contains all the Matlab commands used in the Tutorial in Section 3,
with appropriate pause statements added.
Diagnostics:
The user may want to experiment with other seeds for the random number gen-
erator in the beginning of the script, as well as other noise levels.
104 regutm
regutm
Purpose:
Generates random test matrices for regularization methods.
Synopsis:
[A,U,V] = regutm (m,n,s)
Description:
Generates a random m × n matrix A such that A AT and AT A are oscillating.
Hence, in the SVD of A,
A = U diag(s) VT ,
the number of sign changes in the column vectors U(: , i) and V(: , i) is exactly i − 1.
The third argument s specifies the singular values of A. If not present, then s =
logspace(0,round(log10(eps)),n).
Algorithm:
If m = n then U and V are computed as the left and right singular matrices of
a random bidiagonal n × n matrix with nonnegative elements. If m 6= n then V is
computed as above, and U is computed analogously via a random m × m bidiagonal
matrix. Finally, A is computed as A = U diag(s) VT .
See also:
The test problem routines.
References:
1. P. C. Hansen, Test matrices for regularization methods, SIAM J. Sci. Comput.
16 (1995), 506–512.
rrgmres 105
rrgmres
Purpose:
Range-restricted GMRES algorithm for square systems.
Synopsis:
[X,rho,eta] = rrgmres (A,b,k)
Description:
Range-restricted GMRES is a variant of the GMRES algorithm for square linear
systems A x = b, with starting vector A b (instead of b as in GMRES). The function
returns all k iterates, stored as the columns of the matrix X. The solution norm and
residual norm are returned in eta and rho, respectively.
For symmetric matrices, use the function mr2 instead.
A “preconditioned” version of RRGMRES for the general-form problem where one
minimizes kL xk2 instead of kxk2 is available in the function prrgmres.
Examples:
Perform 25 RRGMRES iterations, and plot the corresponding discrete L-curve:
[X,rho,eta] = rrgmres (A,b,25); plot lc (rho,eta,’o’);
Algorithm:
The RRGMRES algorithm was originally developed in [1] for inconsistent systems,
while the underlying GRMES algorithm is from [4]. An analysis of RRGMRES as a
regularization routine can be found in [3]. The updating formula for the RRGMRES
iterate is given in [2].
See also:
cgls, lsqr, mr2, prrgmres, splsqr
References:
1. D. Calvetti, B. Lewis & L. Reichel, GMRES-type methods for inconsistent sys-
tems, Lin. Alg. Appl. 316 (2000), 157–169.
2. P. C. Hansen, Regularization Tools Version 4.0 for Matlab 7.3. Numer. Algo.
(2007).
3. T. K. Jensen & P. C. Hansen, Iterative regularization with minimum-residual
methods, BIT 47 (2007), 103-120.
4. Y. Saad & M. H. Schultz, GMRES: A generalized minimal residual algorithm
for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (1986),
856–869.
106 shaw
shaw
Purpose:
Test problem: one-dimensional image restoration model.
Synopsis:
[A,b,x] = shaw (n)
Description:
Discretization of a Fredholm integral equation of the first kind (2.1) with [−π/2, π/2]
as both integration intervals. The integral equation is a one-dimensional model of an
image reconstruction problem from [1]. The kernel K and the solution f are given by
µ ¶2
2 sin(u)
K(s, t) = (cos(s) + cos(t))
u
u = π (sin(s) + sin(t))
¡ ¢ ¡ ¢
f (t) = a1 exp −c1 (t − t1 )2 + a2 exp −c2 (t − t2 )2 .
The kernel K is the point spread function for an infinitely long slit. The parameters
a1 , a2 , etc., are constants that determine the shape of the solution f ; in this imple-
mentation we use a1 = 2, a2 = 1, c1 = 6, c2 = 2, t1 = .8, t2 = −.5 giving an f with
two different “humps”.
The kernel and the solution are discretized by simple collocation with n points to
produce A and x. Then the discrete right-hand side is produced as b = A x.
Limitations:
The order n must be even.
Reference:
1. C. B. Shaw, Jr., Improvements of the resolution of an instrument by numerical
solution of an integral equation, J. Math. Anal. Appl. 37 (1972), 83–112.
2.5
1.5
0.5
0
0 10 20 30 40 50 60 70 80 90 100
spikes 107
spikes
Purpose:
Test problem with a “spiky” solution.
Synopsis:
[A,b,x] = spikes (n,t max)
Description:
Artificially generated discrete ill-posed problem. The size of the matrix A is n × n.
The solution x consists of a unit step at t = .5, and a pulse train of spikes of decreasing
amplitude at t = .5, 1.5, 2.5, . . . The parameter t max controls the length of the pulse
train (i.e., the time unit of x is t max/n); t max is optional, and its default is 5 giving
five spikes.
25
20
15
10
0
0 10 20 30 40 50 60 70 80 90 100
108 spleval
spleval
Purpose:
Evaluation of a spline or spline curve.
Synopsis:
points = spleval (f)
Description:
Computes points on the given spline or spline curve f between its extreme breaks.
The representation f for the spline is assumed to be consistent with that used in the
Spline Toolbox [1]. Routine spleval is used in l corner only.
Algorithm:
Original routine fnplt from de Boor’s Spline Toolbox; modified by P. C. Hansen to
skip the plot. The number of points is set by npoints inside spleval; currently, npoints
= 300.
References:
1. C. de Boor, Spline Toolbox, Version 1.1, The Mathworks, MA, 1992.
splsqr 109
splsqr
Purpose:
Subspace preconditioned LSQR (SP-LSQR) for discrete ill-posed problems.
Synopsis:
x = splsqr (A,b,lambda,Vsp,maxit,tol,reorth)
Description:
Subspace preconditioned LSQR for solving the Tikhonov problem in general form
min{kA x − bk2 + λ2 kxk2 }
with a preconditioner based on the subspace defined by the (preferably orthonormal)

columns of the matrix Vsp. The output x holds all the solution iterates as columns,
and the last iterate x(:,end) is the best approximation to xλ .
The parameter maxit is the maximum allowed number of iterations (default value
is maxit = 300), while tol is used a stopping criterion for the norm of the least squares
residual relative to the norm of the right-hand side (default value is tol = 1e-12). A
seventh input parameter reorth 6= 0 enforces MGS reorthogonalization of the Lanczos
vectors.
Examples:
Use the first 6 DCT basis vectors as the preconditioning subspace:
[A,b,x] = shaw (32); b = b + 1e-2*randn (size(b)); lambda = 2e-2;
Vsp = dct (eye(32))’; Vsp = Vsp(:,1:6);
xsp = splsqr (A,b,lambda,Vsp); plot ([x,xsp(:,end)])
Algorithm:
The algorithm is based on the two-level Schur complement algorithm [1], but
implemented via LSQR such that the normal equations are avoided [2].
Limitations:
This is a model implementation of SP-LSQR. In a real implementation the House-
holder transformations should use LAPACK routines, only the final iterate should be
returned, and reorthogonalization is not used. Also, if Vsp represents a fast transfor-
mation (such as the DCT) then explicit storage of Vsp should be avoided. See the
reference for details.
See also:
cgls, lsqr
110 splsqr
References:
1. M. Hanke & C. R. Vogel, Two-level preconditioners for regularized inverse prob-
lems I: Theory, Numer. Math. 83 (1999), 385-402.
2. M. Jacobsen, P. C. Hansen & M. A. Saunders, Subspace preconditioned LSQR
for discrete ill-posed problems, BIT 43 (2003), 975–989.
3. M. E. Kilmer, P. C. Hansen & M. I. Español, A projection-based approach to
general-form Tikhonov regularization, SIAM J. Sci. Comput. 29 (2007), 315–
330.
std form 111
std form
Purpose:
Transform a general-form regularization problem into one in standard form.
Synopsis:
[A s,b s,L p,K,M] = std form (A,L,b) (method 1)
[A s,b s,L p,x 0] = std form (A,L,b,W) (method 2)
Description:
Transforms a regularization problem in general form
© ª
min kA x − bk22 + λ2 kL xk22
into one in standard form

© ª
min kA s xs − b sk22 + λ2 kxs k22 .
Two methods are available, and the distinction between them is controlled by the
number of input parameters to std form.
In method 1, described in [1], A s, L p, K, M, and b s are computed by Eqs.
(2.22)–(2.24). In particular, L p is the pseudoinverse of L, K = Ko has columns in
the null space of L, and M = To−1 HoT . Then the transformation back to the original
general-form setting is given by
x = L p x s + K M (b − A L p x s) .
In method 2, described in [2,3,4], a matrix W holding a basis for the null space
of L is also required, and A s, L p, b s, and x 0 are computed by Eqs. (2.29)–(2.31).
Here L p is the A-weighted generalized inverse of L, and x 0 is the component of the
solution in the null space of L (which, in this method, is independent of xs ). For
method 2, the tranformation back to the general-form setting is simply
x s = L px s + x 0 .
The transformations back to general form can be performed by means of routine

gen form.
Notice that both gen form and std form are available for pedagogical reasons only—
usually it is more efficient to build the transformations into the solution process, such
as in the iterative methods pcgls, plsqr b, and pnu.
Examples:
Transform a general-form problem into standard form, produce 10 TSVD solu-
tions, and transform back again, using method 1; then compare with the mathemat-
ically identical TGSVD solutions:
112 std form
[A s,b s,L p,K,M] = std form (A,L,b); [U,s,V] = csvd (A s);

X s = tsvd (U,s,V,b s,1:10); X = gen form (L p,X s,A,b,K,M);
[U1,sm,X1] = cgsvd (A,L); XX = tgsvd (U1,sm,X1,b,1:10); norm (X-XX)
Algorithm:
The formulas used in method 1 are:
µ ¶
T Rp
L = (Kp Ko ) (QR factorization), L p = Kp Rp−T
0
µ ¶
To
A Ko = (Ho Hq ) (QR factorization), A s = HqT A L p
0
K = Ko , M = To−1 HoT , b s = HqT b .
The formulas used in method 2 are the following (see also pinit):
[T, x 0] = pinit (W, A, b) ;

µµ ¶ ¶
Ip
Lp= − W T(: , 1: p) L(: , 1: p)−1
0
A s = AL p .
See also:
gen form
References:
BIT 17 (1977), 134–145.
2. L. Eldén, A weighted pseudoinverse, generalized singular values, and constrained
tgsvd 113
tgsvd
Purpose:
Compute the truncated GSVD solution.
Synopsis:
[x k,rho,eta] = tgsvd (U,sm,X,b,k) , where sm = [sigma,mu]
Description:
Computes the truncated GSVD solution, defined as
p
X Xn
U(:, i)T b
xk= X(:, i) + U(:, i)T b X(:, i) .
sigma(i) i=p+1
i=p−k+1
x k = [ x k(1), x k(2), . . . ] .
The solution seminorm and the residual norm are returned in eta and rho.
Examples:
Compute the first 10 TGSVD solutions to an inverse Laplace transform model-
problem:
[A,b,x] = ilaplace (n,2); L = get l (n,1); b = b + 1e-4∗randn (size(b));
[U,sm,X] = cgsvd (A,L); X tgsvd = tgsvd (U,sm,X,b,1:10);
Algorithm:
Straightforward use of the above definition is used to compute x k. For motivations
for the TGSVD solution, cf. [1].
See also:
mtsvd, tsvd
References:
1. P. C. Hansen, Regularization, GSVD and truncated GSVD, BIT 29 (1989), 491–
504.
114 tikhonov
tikhonov
Purpose:
Compute the Tikhonov regularized solution.
Synopsis:
[x lambda,rho,eta] = tikhonov (U,s,V,b,lambda,x 0)
[x lambda,rho,eta] = tikhonov (U,sm,X,b,lambda,x 0) , where sm = [sigma,mu]
Description:
Compute the solution by means of Tikhonov regularization. If the SVD of A is
used, i.e. if U, s, and V are given as input, then Tikhonov regularization in standard
form is used, and x lambda solves the problem
© ª
min kA x − bk22 + lambda2 kx − x 0k22 .
On the other hand, if the GSVD of (A, L) is used, i.e. if U, sm, and X are given as
input, then Tikhonov regularization in general form is used, and x lambda then solves
the problem © ª
min kA x − bk22 + lambda2 kL (x − x 0)k22 .
If x 0 is not specified, then x 0 = 0 is used. If lambda is a vector, then x lambda is a
matrix such that
The solution norm (standard-form case) or seminorm (general-form case) and the
residual norm are returned in eta and rho.
Examples:
Solve a discrete ill-posed problem for three values of the regularization parameter,
namely 10−1 , 10−3 , and 10−5 :
x lambda = tikhonov (U,s,V,b,[1e-1,1e-3,1e-5]); plot (x lambda)
Algorithm:
x lambda is computed by straightforward use of the formulas from Section 2.4.
See also:
discrep, lsqi, mtsvd
References:
1. A. N. Tikhonov & V. Y. Arsenin, Solutions of Ill-Posed Problems, Winston &
Sons, Washington, D.C., 1977.
tomo 115
tomo
Purpose:
Test problem: two-dimensional tomography problem.
Synopsis:
[A,b,x] = tomo (N,f)
Description:
This function creates a simple two-dimensional tomography test problem. A 2-D
domain [0, N] × [0, N] is divided into N2 cells of unit size, and a total of round(f ∗ N2 )
rays in random directions penetrate this domain. The default value is f = 1.
Once a solution x reg has been computed, it can be visualized by means of
imagesc (reshape (x reg,N,N)).
The exact solution, reshape (x,N,N), is identical to the exact image in the function
blur.
Algorithm:
Each cell is assigned a value (stored in the vector x), and for each ray the corre-
sponding element in the right-hand side b is the line integral along the ray, i.e.,
X
xcellj · `cellj ,
cells ∈ ray
where `cellj is the length of the ray in the jth cell. The matrix A is sparse, and each
row (corresponding to a ray) holds the value `cellj in the jth position. Hence b = A∗x.
Limitations:
The integer N should not be too small; we recommend N ≥ 16.
References:
1. P. C. Hansen, Discrete Inverse Problems. Insight and Algorithm., manuscript,
IMM, 2007.
116 tsvd
tsvd
Purpose:
Compute the truncated SVD solution.
Synopsis:
[x k,rho,eta] = tsvd (U,s,V,b,k)
Description:
Computes the truncated SVD solution, defined as
Xk
U(:, i)T b
xk= V(:, i) .
i=1
s(i)
x k = [ x k(1), x k(2), . . . ] .
The solution and residual norms are returned in eta and rho.
Examples:
Compute the TSVD solutions for the truncation parameters 3, 6, 9, 12, and 15:
[A,b,x] = shaw (n); b = b + 1e-3∗rand (size(b)); [U,s,V] = csvd (A);
X = tsvd (U,s,V,b,[3,6,9,12,15]); plot ([x,X])
Algorithm:
Straightforward use of the definition is used to compute x k.
See also:
mtsvd, tgsvd
ttls 117
ttls
Purpose:
Compute the truncated TLS solution.
Synopsis:
[x k,rho,eta] = ttls (V1,k,s1)
Description:
Computes the truncated TLS solution, given by
x k = −V1(1 : n, k + 1 : n + 1) V1(n + 1, k + 1 : n + 1))† ,
where V1 is the right singular matrix in the SVD of the compound matrix
(A b) = U1 diag(s1) V1T .
x k = [ x k(1), x k(2), . . . ] .
If k is not specified, the default value is k = n and x k is then the standard TLS
solution.
The solution norms kx kk2 and the corresponding residual norms k(A b)−(Ã b̃)k kF ,
where (Ã b̃)k = U1 (: , 1 : k) diag(s1(1 : k)) V1(: , 1 : k)T , are returned in eta and rho,
respectively. The singular values s1 of (A b) are required to compute rho.
Examples:
Compute the truncated TLS solutions for k = n − 5, n − 10, and n − 15:
[U1,s1,V1] = csvd ([A,b],’full’); X = ttls (V1,n−[5,10,15]);
Algorithm:
x k is computed by means of the definition, using the following relation for the
pseudoinverse of a row vector vT :
(vT )† = kvk−2
2 v .
The norms are given by

q
kx kk2 = kV1(n + 1, k + 1 : n + 1)k−2
2 −1
k(A b) − (Ã b̃)k kF = ks1(k + 1 : n + 1)k2 .
References:
1. R. D. Fierro, G. H. Golub, P. C. Hansen & D. P. O’Leary, Regularization by
truncated total least squares, SIAM J. Sci. Comput. 18 (1997), 1223–1241.
118 ursell
ursell
Purpose:
Test problem: integral equation with no square integrable solution.
Synopsis:
[A,b] = ursell (n)
Description:
Discretization of a Fredholm integral equation of the first kind (2.1) from [1] with
kernel K and right-hand side g given by
1
K(s, t) = , g(s) = 1 ,
s+t+1
where both integration intervals are [0, 1]. This integral equation does not satisfy the
Picard condition and has no square integrable solution (hence, no x is produced). The
size of the matrix A is n × n.
Examples:
Generate A and b and check the discrete Picard condition:
[A,b] = ursell (16); [U,s,V] = csvd (A); picard (U,s,b);
References:
1. F. Ursell, Introduction to the theory of linear integral equations, Chapter 1 in L.
M. Delves & J. Walsh (Eds.), Numerical Solution of Integral Equations, Claren-
don Press, Oxford, 1974.
wing 119
wing
Purpose:
Test problem with a discontinuous solution.
Synopsis:
[A,b,x] = wing (n,t1,t2)
Description:
Discretization of a Fredholm integral equation of the first kind (2.1) from [1, p. 109]
with both integration intervals equal to [0, 1], with kernel K and right-hand side g
given by
exp(− s t12 ) − exp(− s t22 )

K(s, t) = t exp(− s t2 ) , g(s) = ,
2s
and with the solution f given by
½
1, for t1 < t < t2
f (t) =
0, elsewhere.
Here, t1 and t2 are constants satisfying 0 < t1 < t2 < 1. If they are not specified, the
values t1 = 1/3 and t2 = 2/3 are used. The size of the matrix A is n × n.
References:
1. G. M. Wing & J. D. Zahrt, A Primer on Integral Equations of the First Kind,
SIAM, Philadelphia, 1991; p. 109.
0.1
0.08
0.06
0.04
0.02
0
0 10 20 30 40 50 60 70 80 90 100
120 wing
Bibliography
[1] R. C. Allen, Jr., W. R. Boland, V. Faber & G. M. Wing, Singular values and
condition numbers of Galerkin matrices arising from linear integral equations of
the first kind, J. Math. Anal. Appl. 109 (1985), 564–590.
[2] R. S. Anderssen & P. M. Prenter, A formal comparison of methods proposed for

the numerical solution of first kind integral equations, J. Austral. Math. Soc.
(Series B) 22 (1981), 488–500.
[3] C. T. H. Baker, Expansion methods, Chapter 7 in [19].
[4] C. T. H. Baker, The Numerical Treatment of Integral Equations, Clarendon Press,

Oxford, 1977.
[5] C. T. H. Baker & G. F. Miller (Eds.), Treatment of Integral Equations by Nu-

merical Methods, Academic Press, New York, 1982.
[6] D. M. Bates, M. J. Lindstrom, G. Wahba & B. S. Yandell, GCVPACK – routines

for generalized cross validation, Commun. Statist.-Simula. 16 (1987), 263–297.
[7] M. Bertero, C. De Mol & E. R. Pike, Linear inverse problems with discrete data:
II. Stability and regularization, Inverse Problems 4 (1988), 573–594.
[8] Å. Björck, A bidiagonalization algorithm for solving large and sparse ill-posed
systems of linear equations, BIT 28 (1988), 659–670.
[9] Å. Björck, Numerical Methods for Least Squares Problems, SIAM, Philadelphia,
1996.
[10] Å. Björck & L. Eldén, Methods in numerical algebra for ill-posed problems, Report
LiTH-MAT-R33-1979, Dept. of Mathematics, Linköping University, 1979.
[11] H. Brakhage, On ill-posed problems and the method of conjugate gradients; in

H. W. Engl & C. W. Groetsch (Eds.), Inverse and Ill-Posed Problems, Academic
Press, Boston, 1987.
[12] D. Calvetti, B. Lewis & L. Reichel, GMRES-type methods for inconsistent sys-
tems, Lin. Alg. Appl. 316 (2000), 157–169.
122 BIBLIOGRAPHY
[13] T. F. Chan & P. C. Hansen, Some applications of the rank revealing QR factor-
ization, SIAM J. Sci. Stat. Comput. 13 (1992), 727–741.
[14] J. A. Cochran, The Analysis of Linear Integral Equations, McGraw-Hill, New

York, 1972.
[15] D. Colton & R. Kress, Integral Equation Methods for Scattering Theory, John
Wiley, New York, 1983.
[16] I. J. D. Craig & J. C. Brown, Inverse Problems in Astronomy, Adam Hilger,

Bristol, 1986.
[17] J. J. M. Cuppen, A Numerical Solution of the Inverse Problem of Electrocardio-

graphy, Ph. D. Thesis, Dept. of Mathematics, Univ. of Amsterdam, 1983.
[18] L. M. Delves & J. L. Mohamed, Computational Methods for Integral Equations,

Cambridge University Press, Cambridge, 1985.
[19] L. M. Delves & J. Walsh (Eds.), Numerical Solution of Integral Equations, Claren-
don Press, Oxford, 1974.
[20] J. B. Drake, ARIES: a computer program for the solution of first kind integral
equations with noisy data, Report K/CSD/TM-43, Dept. of Computer Science,
Oak Ridge National Laboratory, October 1983.
[21] M. P. Ekstrom & R. L. Rhodes, On the application of eigenvector expansions to

numerical deconvolution, J. Comp. Phys. 14 (1974), 319-340.
[22] L. Eldén, Algorithms for regularization of ill-conditioned least-squares problems,

BIT 17 (1977), 134–145.
[23] L. Eldén, A program for interactive regularization, Report LiTH-MAT-R-79-25,

Dept. of Mathematics, Linköping University, Sweden, 1979.
[24] L. Eldén, A weighted pseudoinverse, generalized singular values, and constrained

[25] L. Eldén, A note on the computation of the generalized cross-validation function

for ill-conditioned least squares problems, BIT 24 (1984), 467–472.
[26] H. W. Engl & J. Gfrerer, A posteriori parameter choice for general regularization
methods for solving linear ill-posed problems, App. Numerical Math. 4 (1988),
395–417.
[27] R. D. Fierro & J. R. Bunch, Collinearity and total least squares, SIAM J. Matrix
Anal. Appl., 15 (1994), pp. 1167–1181.
BIBLIOGRAPHY 123
[28] R. D. Fierro, G. H. Golub, P. C. Hansen & D. P. O’Leary, Regularization by

truncated total least squares, SIAM J. Sci. Comput. 18 (1997), 1223–1241.
[29] R. Fletcher, Practical Optimization Methods. Vol. 1, Unconstrained Optimiza-
tion, Wiley, Chichester, 1980.
[30] G. H. Golub & C. F. Van Loan, Matrix Computations, 3. Ed., Johns Hopkins,
Baltimore, 1996.
[31] G. H. Golub & U. von Matt, Quadratically constrained least squares and quadratic
problems, Numer. Math. 59 (1991), 561–580.
[32] C. W. Groetsch, The Theory of Tikhonov Regularization for Fredholm Equations
of the First Kind, Pitman, Boston, 1984.
[33] C. W. Groetsch, Inverse Problems in the Mathematical Sciences, Vieweg Verlag,
Wiesbaden, 1993.
[34] C. W. Groetsch & C. R. Vogel, Asymptotic theory of filtering for linear operator
equations with discrete noisy data, Math. Comp. 49 (1987), 499–506.
[35] J. Hadamard, Lectures on Cauchy’s Problem in Linear Partial Differential Equa-
tions, Yale University Press, New Haven, 1923.
[36] M. Hanke, Regularization with differential operators. An iterative approach, J.
[37] M. Hanke, Iterative solution of underdetermined linear systems by transformation
to standard form; in Numerical Mathematics in Theory and Practice, Dept. of
Mathematics, University of West Bohemia, Plzeň, pp. 55–63 (1993).
[38] M. Hanke, Conjugate Gradient Methods for Ill-Posed Problems, Longman Scien-
tific and Technical, Essex, 1995.
[39] M. Hanke & P. C. Hansen, Regularization Methods for Large-Scale Problems,
Surv. Math. Ind. 3 (1993), 253–315.
[40] P. C. Hansen, The truncated SVD as a method for regularization, BIT 27 (1987),
543–553.
[41] P. C. Hansen, Computation of the singular value expansion, Computing 40
(1988), 185–199.
[42] P. C. Hansen, Perturbation bounds for discrete Tikhonov regularization, Inverse
Problems 5 (1989), L41–L44.
[43] P. C. Hansen, Regularization, GSVD and truncated GSVD, BIT 29 (1989), 491–
504.
124 BIBLIOGRAPHY
[44] P. C. Hansen, Truncated SVD solutions to discrete ill-posed problems with ill-de-
termined numerical rank, SIAM J. Sci. Stat. Comput. 11 (1990), 503–518.
[45] P. C. Hansen, Relations between SVD and GSVD of discrete regularization prob-
lems in standard and general form, Lin. Alg. Appl. 141 (1990), 165–176.
[46] P. C. Hansen, The discrete Picard condition for discrete ill-posed problems, BIT
30 (1990), 658–672.
[47] P. C. Hansen, Analysis of discrete ill-posed problems by means of the L-curve,

SIAM Review 34 (1992), 561–580.
[48] P. C. Hansen, Numerical tools for analysis and solution of Fredholm integral
equations of the first kind, Inverse Problems 8 (1992), 849–872.
[49] P. C. Hansen, Rank-Deficient and Discrete Ill-Posed Problems. Numerical As-

[50] P. C. Hansen, Regularization Tools Version 3.0 for Matlab 5.2, Numer. Algo. 20
(1999), 195–196.
[51] P. C. Hansen, Deconvolution and regularization with Toeplitz matrices, Numer.

Algo. 29 (2002), 323–378.
[52] P. C. Hansen & T. K. Jensen, Smoothing-norm preconditioning for regularizing

minimum-norm methods, SIAM J. Matrix Anal. Appl. 29 (2006), 1–14.
[53] P. C. Hansen, T. K. Jensen & G. Rodriguez, An adaptive pruning algorithm for

[54] P. C. Hansen, M. Kilmer & R. H. Kjeldsen, Exploiting residual information in

the parameter choice for discrete ill-posed problems, BIT 46 (2006), 41–59.
[55] P. C. Hansen & D. P. O’Leary, The use of the L-curve in the regularization of
[56] P. C. Hansen, T. Sekii & H. Shibahashi, The modified truncated SVD method for
regularization in general form, SIAM J. Sci. Stat. Comput. 13 (1992), 1142-1150.
[57] B. Hofmann, Regularization for Applied Inverse and Ill-Posed Problems, Teubner,
Leipzig, 1986.
[58] M. Jacobsen, P. C. Hansen & M. A. Saunders, Subspace preconditioned LSQR

for discrete ill-posed problems, BIT 43 (2003), 975–989.
[59] T. K. Jensen & P. C. Hansen, Iterative regularization with minimum-residual

methods, BIT 47 (2007), 103–120.
BIBLIOGRAPHY 125
[60] R. Kress, Linear Integral Equations, Springer, Berlin, 1989.
[61] C. L. Lawson & R. J. Hanson, Solving Least Squares Problems, Prentice-Hall,

Englewood Cliffs, 1974.
[62] P. Linz, Uncertainty in the solution of linear operator equations, BIT 24 (1984),
92–101.
[63] Matlab Reference Guide, The MathWorks, Mass., 1996.
[64] K. Miller, Least squares methods for ill-posed problems with a prescribed bound,
SIAM J. Math. Anal. 1 (1970), 52–74.
[65] V. A. Morozov, Methods for Solving Incorrectly Posed Problems, Springer Verlag,
New York, 1984.
[66] F. Natterer, The Mathematics of Computerized Tomography, John Wiley, New

York, 1986.
[67] F. Natterer and F. Wübbeling, Mathematical Methods in Image Reconstruction,

[68] F. Natterer, Numerical treatment of ill-posed problems, in [75].
[69] D. P. O’Leary & J. A. Simmons, A bidiagonalization-regularization procedure for

large scale discretizations of ill-posed problems, SIAM J. Sci. Stat. Comput. 2
(1981), 474–489.
[70] C. C. Paige & M. A. Saunders, LSQR: an algorithm for sparse linear equations
[71] R. L. Parker, Understanding inverse theory, Ann. Rev. Earth Planet Sci. 5 (1977),
35–64.
[72] D. L. Phillips, A technique for the numerical solution of certain integral equations
of the first kind, J. ACM 9 (1962), 84–97.
[73] F. Santosa, Y.-H. Pao, W. W. Symes & C. Holland (Eds.), Inverse Problems of
Acoustic and Elastic Waves, SIAM, Philadelphia, 1984.
[74] C. Ray Smith & W. T. Grandy, Jr. (Eds.), Maximum-Entropy and Bayesian
Methods in Inverse Problems, Reidel, Boston, 1985.
[75] G. Talenti (Ed.), Inverse Problems, Lecture Notes in Mathematics 1225, Springer
Verlag, Berlin, 1986.
[76] H. J. J. te Riele, A program for solving first kind Fredholm integral equations by
means of regularization, Computer Physics Comm. 36 (1985), 423–432.
126 BIBLIOGRAPHY
[77] A. N. Tikhonov, Solution of incorrectly formulated problems and the regulariza-

tion method, Dokl. Akad. Nauk. SSSR 151 (1963), 501–504 = Soviet Math. Dokl.
4 (1963), 1035–1038.
[78] A. N. Tikhonov & V. Y. Arsenin, Solutions of Ill-Posed Problems, Winston &
Sons, Washington, D.C., 1977.
[79] A. N. Tikhonov & A. V. Goncharsky, Ill-Posed Problems in the Natural Sciences,
MIR Publishers, Moscow, 1987.
[80] A. van der Sluis, The convergence behavior of conjugate gradients and Ritz values
in various circumstances; in R. Beauwens & P. de Groen (Eds.), Iterative Methods
in Linear Algebra, North-Holland, Amsterdam, 1992.
[81] A. van der Sluis & H. A. van der Vorst, SIRT- and CG-type methods for iterative
solution of sparse linear least-squares problems, Lin. Alg. Appl. 130 (1990), 257–
302.
[82] J. M. Varah, On the numerical solution of ill-conditioned linear systems with
applications to ill-posed problems, SIAM J. Numer. Anal. 10 (1973), 257–267.
[83] J. M. Varah, A practical examination of some numerical methods for linear dis-
crete ill-posed problems, SIAM Rev. 21 (1979), 100–111.
[84] J. M. Varah, Pitfalls in the numerical solution of linear ill-posed problems, SIAM
J. Sci. Stat. Comput. 4 (1983), 164–176.
[85] C. R. Vogel, Optimal choice of a truncation level for the truncated SVD solution
of linear first kind integral equations when data are noisy, SIAM J. Numer. Anal.
23 (1986), 109–117.
[86] C. R. Vogel, Solving ill-conditioned linear systems using the conjugate gradi-
ent method, Report, Dept. of Mathematical Sciences, Montana State University,
1987.
[87] G. Wahba, Spline Models for Observational Data, CBMS-NSF Regional Confer-
ence Series in Applied Mathematics, Vol. 59, SIAM, Philadelphia, 1990.
[88] G. M. Wing, Condition numbers of matrices arising from the numerical solution
of linear integral equations of the first kind, J. Integral Equations 9 (Suppl.)
(1985), 191–204.
[89] G. M. Wing & J. D. Zahrt, A Primer on Integral Equations of the First Kind,
[90] H. Zha & P. C. Hansen, Regularization and the general Gauss-Markov linear
model, Math. Comp. 55 (1990), 613–624.

RTV 4 Manual - Regu Tools

Uploaded by

Copyright:

Available Formats

RTV 4 Manual - Regu Tools

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

RTV 4 Manual - Regu Tools

Uploaded by

Copyright:

Available Formats

Regularization Tools

A Matlab Package for

Version 4.1 for Matlab 7.3

Per Christian Hansen

The software described in this report was originally published in

The current version is published in Numer. Algo. 46 (2007), pp. 189–194,

Changes from Earlier Versions 3

2 Discrete Ill-Posed Problems and their Regularization 9

3 Regularization Tools Tutorial 39

4 Regularization Tools Reference 47

• Replaced gsvd by cgsvd which has a different sequence of output arguments.

• Deleted the function mgs.

• Changed the storage format of bidiagonal matrices to sparse, instead of a dense

• Removed the obsolete function bsvd.

• Added the blur test problem.

• Added a priori guess x 0 as input to tikhonov.

• Corrected get l such that the sign of L*x is correct.

• Added MGS reorthogonalization of the normal equation residual vectors in the

• Added the method ’ttls’ to the function fil fac.

• More precise computation of the regularization parameter in gcv, lcurve, and

• Changed heb new and newton to work with λ instead of λ2 .

• Added legend to lagrange and picard.

• Corrected the routines to work for complex problems.

Ill-posed problems—and regularization methods for computing stabilized solutions to

with coefficient matrix A and right-hand side b given by

Here, the right-hand side b is generated by adding a small perturbation to an exact

min kA x − bk2 subject to kxk2 ≤ α .

1. the condition number of the matrix A is large

2. replacing A by a well-conditioned matrix derived from A does not necessarily

3. care must taken when imposing additional constraints.

The purpose of numerical regularization theory is to provide efficient and numerically

Acknowledgement. I am indebted to the many users of Regularization Tools

In this chapter we give a brief introduction to discrete ill-posed problems, we discuss

2.1. Discrete Ill-Posed Problems

and due to the Riemann-Lebesgue lemma it follows that ∆g → 0 as p → ∞ [32, p. 2].

and linear least-squares problems

min kA x − bk2 , A ∈ Rm×n , m>n. (2.3)

2.2. Regularization Methods

where the regularization parameter λ controls the weight given to minimization of

2.3. SVD and Generalized SVD

2.3.1. The Singular Value Decomposition

where U = (u1 , . . . , un ) and V = (v1 , . . . , vn ) are matrices with orthonormal columns,

2.3.2. The Generalized Singular Value Decomposition

where the columns of U ∈ Rm×n and V ∈ Rp×p are orthonormal, X ∈ Rn×n is

and they are normalized such that

error in the perturbed solution xλ satisfies

2.4. The Discrete Picard Condition and Filter Factors

The discrete Picard condition. The unperturbed right-hand side b̄ in a discrete

2.5. The L-Curve

Figure 2.1: The generic form of the L-curve.

yields the following expression for the error in xreg :

For a given fixed right-hand side b = b̄ + e, there is obviously an optimal reg-

In the Regularization Tools package, routine l curve produces a log-log plot of

2.6. Transformation to Standard Form

2.6.1. Transformation for Direct Methods

respectively). First, compute a QR factorization of LT ,

Ā = HqT A L† = HqT A Kp Rp−T , b̄ = HqT b , (2.24)

x = L† x̄ + Ko To−1 HoT (b − A L† x̄) . (2.25)

The SVD of Ā is related to the GSVD of (A, L) as follows: let Ā = Ū Σ̄ V̄ T denote