Choosing Regularization Parameters in Iterative Methods for Ill-Posed Problems.
Numerical solution of ill-posed problems is often accomplished by discretization (projection onto a finite dimensional subspace) followed by regularization. If the discrete problem has high dimension, though, typically we compute an approximate solution by projecting the discrete problem onto an even smaller dimensional space, via iterative methods based on Krylov subspaces. In this work we present a common framework for efficient algorithms that regularize after this second projection rather than before it. We show that determining regularization parameters based on the final projected problem rather than on the original discretization has firmer justification and often involves less computational expense. We prove some results on the approximate equivalence of this approach to other forms of regularization, and we present numerical examples.
Introduction
. Linear, discrete ill-posed problems of the form
(1)
or
min
x
equivalently, A
(2)
arise, for example, from the discretization of first-kind Fredholm integral equations
and occur in a variety of applications. We shall assume that the full-rank matrix A
in (2) and in (1). In discrete ill-posed problems, A is ill-conditioned
and there is no gap in the singular value spectrum. Typically, the right
hand side b contains noise due to measurement and/or approximation error. This
noise, in combination with the ill-conditioning of A, means that the exact solution of
(1) or (2) has little relationship to the noise-free solution and is worthless. Instead, we
use a regularization method to determine a solution that approximates the noise-free
solution. Regularization methods replace the original operator by a better-conditioned
but related one in order to diminish the effects of noise in the data and produce a
regularized solution to the original problem. Sometimes this regularized problem is
too large to solve exactly. In that case, we typically compute an approximate solution
by projection onto an even smaller dimensional space, perhaps via iterative methods
based on Krylov subspaces.
The conditioning of the new problem is controlled by one or more regularization
parameters specific to the method. A large regularization parameter yields a new well-conditioned
problem, but its solution may be far from the noise-free solution since
the new operator is a poor approximation to A. A small regularization parameter
generally yields a solution very close to the noise-contaminated exact solution of (1)
or (2), and hence its distance from the noise-free solution also can be large. Thus,
This work was supported by the National Science Foundation under Grants CCR 95-03126 and
CCR-97-32022 and by the Army Research Office, MURI Grant DAAG55-97-1-0013.
y Dept. of Computer and Electrical Engineering, Northeastern University, Boston, MA 02115
z Dept. of Computer Science and Institute for Advanced Computer Studies, University of Mary-
land, College Park, MD 20742 (oleary@cs.umd.edu).
a key issue in regularization methods is to choose a regularization parameter that
balances the error due to noise with the error due to regularization.
A wise choice of regularization parameter is obviously crucial to obtaining useful
approximate solutions to ill-posed problems. For problems small enough that a rank-
revealing factorization or singular value decomposition of A can be computed, there
are well-studied techniques for computing a good regularization parameter. These
techniques include the Discrepancy Principle [8], generalized cross-validation (GCV)
[9], and the L-curve [15]. For larger problems treated by iterative methods, though,
the parameter choice is much less understood. If regularization is applied to the
projected problem that is generated by the iterative method, then there are essentially
two regularization parameters: one for the standard regularization algorithms, such
as Tikhonov or truncated SVD, and one controlling the number of iterations taken.
One subtle issue is that the standard regularization parameter that is correct for the
discretized problem may not be the optimal one for the lower-dimensional problem
actually solved by the iteration, and this observation leads to the research discussed
in this paper. At first glance, there can appear to be a lot of work associated with
the selection of a good regularization parameter, and many algorithms proposed in
the literature are needlessly complicated. But by regularizing after projection by the
iterative method, so that we are regularizing the lower dimensional problem that is
actually being solved, much of this difficulty vanishes.
The purpose of this paper is to present parameter selection techniques designed
to reduce the regularization work for iterative methods such as Krylov subspace tech-
niques. Our paper is organized as follows. In x2, we will give an overview of the
regularization methods we will be considering, and we follow up in x3 by surveying
some methods for choosing the corresponding regularization parameters. In x4, we
show how parameter selection techniques for the original problem can be applied instead
to a projected problem obtained from an iterative method, greatly reducing the
cost without much degradation in the solution. We give experimental results in x5
and conclusions and future work in x6.
2. Regularization background. In the following we shall assume that
e, where b true denotes the unperturbed data vector and e denotes zero-mean
white noise. We will also assume that b true satisfies the discrete Picard condition;
that is, the spectral coefficients of b true decay faster, on average, than the singular
values.
Under these assumptions, it is easy to see why the exact solution to (1) or (2) is
hopelessly contaminated by noise. Let "
denote the singular value decomposition
of A, where the columns of "
U and "
V are the singular vectors, and the singular values
are ordered as oe 1  oe . Then the solution to (1) or (2) is given by
"
"
As a consequence of the white noise assumption, j"u
i ej is roughly constant for all i,
while the discrete Picard condition guarantees that j"u
than oe i does . The matrix A is ill-conditioned, so small singular values magnify the
corresponding coefficients "
e in the second sum, and it is this large contribution of
noise from the approximate null space of A that renders the exact solution x defined
in (3) worthless. The following regularization methods try in different ways to lessen
the contribution of noise to the solution. For further information on these methods,
see, for example, [17].
2.1. Tikhonov regularization. One of the most common methods of regularization
is Tikhonov regularization [34]. In this method, the problem (1) or (2) is
replaced with the problem of solving
where L denotes a matrix, often chosen to be the identity matrix I or a discrete
derivative operator, and  is a positive scalar regularization parameter. For ease in
notation, we will assume that Solving (4) is equivalent to solving
In analogy with (3) we have
In this solution, the contributions from noise components " u
e for values of oe are
much smaller than they are in (3), and thus x  can be closer to the noise-free solution
than x is. If  is too large, however, A   I is very far from the original operator
A   A, and x  is very far from x true , the solution to (2) when Conversely, if
is too small, the singular values of the new operator A   A are close to those of
A   A; thus x  x, so small singular values again greatly magnify noise components.
2.2. Truncated SVD. In the truncated SVD method of regularization, the regularized
solution is chosen simply by truncating the expansion in (3) as
"
Here the regularization parameter is ', the number of terms to be dropped from the
sum. Observe that if ' is small, very few terms are dropped from the sum, so x '
resembles x in that the effects of noise are large. If ' is too large, however, important
information could be lost; such is the case if "
An alternative, yet related, approach to TSVD is an approach introduced by Rust
[31] where the truncation strategy is based on the value of each spectral coefficient
itself. The strategy is to include in the sum (3) only those terms corresponding to
a spectral coefficient "
b whose magnitude is greater than or equal to some tolerance
ae, which can be regarded as the regularization parameter.
2.3. Projection and iterative methods. Solving (5) or (7) can be impractical
if n is large, but fortunately, regularization can be achieved through projection onto
a subspace; see, for example, [7]. The truncated SVD is an example of one such
projection: the solution is constrained to lie in the subspace spanned by the singular
vectors corresponding to the largest values. Other projections can be
more economical. In general, we constrain our regularized solution to lie in some
k-dimensional subspace of C n , spanned by the columns of an n \Theta k matrix Q (k) . For
example, we choose x (k)
min
kAQ
or equivalently
The idea is that with an appropriately chosen subspace, the operator (Q (k) )   A   AQ (k)
will be better conditioned than the original operator and hence that x (k)
reg will approximate
x true well on that subspace.
This projection is often achieved through the use of iterative methods such as conjugate
gradients, GMRES, QMR, and other Krylov subspace methods. The matrix
contains orthonormal columns generated via a Lanczos tridiagonalization or
bidiagonalization process [27, 1]. In this case, Q (k) is a basis for some k-dimensional
Krylov subspace (i.e., the subspace K k (c; K) spanned by the vectors c;
for some matrix K and vector c). The regularized solutions x (k)
reg are generated iteratively
as the subspaces are built. Krylov subspace algorithms such as CG, CGLS,
GMRES, and LSQR tend to produce, at early iterations, solutions that resemble x true
in the subspace spanned by (right) singular vectors of A corresponding to the largest
singular values. At later iterations, however, these methods start to reconstruct increasing
amounts of noise into the solution. This is due to the fact that for large k,
the operator (Q (k) )   A   AQ (k) approaches the ill-conditioned operator A   A. There-
fore, the choice of the regularization parameter k, the stopping point for the iteration
and the dimension of the subspace, is very important. 1
2.4. Hybrid methods: projection plus regularization. Another important
family of regularization methods, often referred to as hybrid methods [17], was introduced
by O'Leary and Simmons [27]. These methods combine a projection method
with a direct regularization method such as TSVD or Tikhonov regularization. The
problem is projected onto a particular subspace of dimension k, but typically the
restricted operator in (9) is still ill-conditioned. Therefore, a regularization method
is applied to the projected problem. Since the dimension k is usually small relative
to n, regularization of the restricted problem is much less expensive. Yet, with an
appropriately chosen subspace, the end results can be very similar to those achieved
by applying the same direct regularization technique to the original problem. We will
become more precise about how "similar" the solutions are in x4.5. Because the projected
problems are usually generated iteratively by a Lanczos method, this approach
is useful when A is sparse or structured in such a way that matrix-vector products
can be handled efficiently with minimal storage.
3. Existing parameter selection methods. In this section, we discuss a sampling
of the parameter selection techniques that have been proposed in the literature.
They differ in the amount of a priori information required as well as in the decision
criteria.
3.1. The Discrepancy Principle. If some extra information is available - for
example, an estimate of the variance of the noise vector e - then the regularization
parameter can be chosen rather easily. Morozov's Discrepancy Principle [25] says that
if ffi is the expected value of kek 2 , then the regularization parameter should be chosen
so that the norm of the residual corresponding to the regularized solution x reg is
that is,
Usually, small values of the regularization parameter correspond to a closer solution to the noisy
equation, but despite this, we will call k, rather than 1=k, the regularization parameter.
x
l
||Fig. 1. Example of a typical L-curve. This particular L-curve corresponds to applying Tikhonov
regularization to the problem in Example 2
predetermined real number. Note that as
Other methods based on knowledge of the variance are given, for example, in [12, 5].
3.2. Generalized Cross-Validation. The Generalized Cross-Validation (GCV)
parameter selection method does not depend on a priori knowledge about the noise
variance. This idea of Golub, Heath, and Wahba [9] is to find the parameter  that
minimizes the GCV functional
denotes the matrix that maps the right hand side b onto the regularized
solution x . In Tikhonov regularization, for example, A ]
is
GCV chooses a regularization parameter that is not too dependent on any one data
measurement [11, 12.1.3].
3.3. The L-Curve. One way to visualize the tradeoff between regularization
error and error due to noise is to plot the norm of the regularized solution versus
the corresponding residual norm for each of a set of regularization parameter values.
The result is the L-curve, introduced by Lawson and popularized by Hansen [15].

Figure

1 for a typical example. As the regularization parameter increases, noise
is damped, so that the norm of the solution decreases while the residual increases.
Intuitively, the best regularization parameter should lie on the corner of the L-curve,
since for values higher than this, the residual increases without reducing the norm
of the solution much, while for values smaller than this, the norm of the solution
increases rapidly without much decrease in residual. In practice, only a few points
on the L-curve are computed and the corner is located by approximate methods,
estimating the point of maximum curvature [19].
Like GCV, this method of determining a regularization parameter does not depend
on specific knowledge about the noise vector.
3.4. Disadvantages of these parameter choice algorithms. The appropriate
choice of regularization parameter - especially for projection algorithms - is a
difficult problem, and each method has severe flaws.
Basic cost Added Cost
Disc. GCV L-curve
Rust's TSVD O(mn 2 ) O(m log m) O(m log m) O(m log m)
Projection

Table
Summary of additional flops needed to compute the regularization parameter for each four
regularization methods with various parameter selection techniques. Notation:
q is the cost of multiplication of a vector by A.
p is the number of discrete parameters that must be tried;
k is the dimension of the projection.
m and n are problem dimensions.
The Discrepancy Principle is convergent as the noise goes to zero, but it relies on
knowing information that is often unavailable or incorrectly estimated. Even with a
correct estimate of the variance, the solutions tend to be oversmoothed [20, pg. 96]
(see also the discussion in x6.1 of [15]).
One noted difficulty with GCV is that G can have a very flat minimum, making
it difficult to determine the optimal  numerically [35].
The L-curve is usually more tractable numerically, but its limiting properties are
nonideal. The solution estimates fail to converge to the true solution as n !1 [36]
or as the error norm goes to zero [6]. All methods that assume no knowledge of the
error norm - including GCV - have this latter property [6].
For further discussion and references about parameter choice methods, see [5, 17].
The cost of these methods is tabulated in Table 1.
3.5. Previous work on parameter choice for hybrid methods. At first
glance, it appears that for Tikhonov regularization, multiple systems of the form
(5) must be solved in order to evaluate candidate values of  for the Discrepancy
Principle or the L-curve. Techniques have been suggested in the literature for solving
these systems using projection methods.
Chan and Ng [4], for example, note that the systems involve the closely related
matrices matrices suggest solving the systems simultaneously
using a Galerkin projection method on a sequence of "seed" systems. Although
this is economical in storage, it can be unnecessarily expensive in time because they
do not exploit the fact that for each fixed k, the Krylov subspace K k
the same for all values of .
Frommer and Maass [8] propose two algorithms for approximating the  that
satisfies the Discrepancy Principle (10). The first is a "truncated cg" approach in
which they use conjugate gradients to solve k systems of the form (5), truncating
the iterative process early for large  and using previous solutions as starting guesses
for later problems. Like Chan and Ng, this algorithm does not exploit any of the
redundancy in generating the Krylov-subspaces for each  i . The second method
they propose, however, does exploit the redundancy so that the CG iterates for all k
systems can be updated simultaneously with no extra matrix-vector products. They
stop their "shifted cg" algorithm when kAx for one of their  values.
Thus the number of matrix-vector products required is twice the number of iterations
for this particular system to converge. We note that while the algorithms we propose
in x4 for finding a good value of  are based on the same key observation regarding
the Krylov subspace, our methods will usually require less work than the shifted cg
algorithm.
Calvetti, Golub, and Reichel [3] compute upper and lower bounds on the L-curve
generated by the matrices C() using a Lanczos bidiagonalization process. From this,
they approximate the best parameter for Tikhonov regularization without projection.
In x4, we choose instead to approximate the best parameter for Tikhonov regularization
on the projected problem, since this is the approximation to the continuous
problem that is actually being used.
Kaufman and Neumaier [21] suggest an envelope guided conjugate gradient approach
for the Tikhonov L-curve problem. Their method is more complicated than the
methods we propose because they maintain nonnegativity constraints on the variables.
Substantial work has also been done on TSVD regularization of the projected
problems. Bjorck, Grimme, and van Dooren [2] use GCV to determine the truncation
point for the projected SVD. Their emphasis is on stable ways to maintain an accurate
factorization when many iterations are needed, and they use full reorthogonalization
and implicit restart strategies. O'Leary and Simmons [27] take a somewhat different
viewpoint that the problem should be preconditioned appropriately so that a massive
number of iterations is unnecessary. That viewpoint is echoed in this current work,
so we implicitly assume that the problem has been left-preconditioned or "filtered"
[27]. For example, in place of (4), we solve
min
x
2for a square preconditioner M . See [14, 26, 24, 23] for preconditioners appropriate for
certain types of ill-posed problems. Note that we could alternately have considered
right preconditioning, which amounts to solving, in the Tikhonov case,
min
y
A
I
for y then setting . Note that either left or right preconditioning
effectively changes the balance between the two terms in the minimization.
4. Regularizing the projected problem. In this section we develop nine approaches
to regularization using Krylov methods. Many Krylov methods have been
proposed; for ease of exposition we focus on just two of these: the LSQR algorithm
of Paige and Saunders [29] and the GMRES algorithm of Saad and Schultz [33].
The LSQR algorithm of Paige and Saunders [29] iteratively computes the bidiag-
onalization introduced by Golub and Kahan [10]. Given a vector b, the algorithm is
as follows [29, Alg. Bidiag 1]:
Compute a scalar fi 1 and a vector u 1 of length one so that fi 1
Similarly, determine ff 1 and v 1 so that ff 1
For
where the non-negative scalars ff i+1 and fi i+1 are chosen
so that u i+1 and v i+1 have length one.
End for
The vectors u are called the left and right Lanczos vectors respectively. The
algorithm can be rewritten in matrix form by first defining the matrices
. ff k
denoting the ith unit vector, the following relations can be established:
A T U
where the subscript on I denotes the dimension of the identity.
Now suppose we want to solve
min
where S denotes the k-dimensional subspace spanned by the first k Lanczos vectors
. The solution we seek is of the form x vector y (k) of length k.
to be the corresponding residual. From the relations above,
observe that in exact arithmetic
Since U k+1 has, in exact arithmetic, orthonormal columns, we have
Therefore, the projected problem we wish to solve is
min
y
Solving this minimization problem is equivalent to solving the normal equations involving
the bidiagonal matrix:
Typically k is small, so reorthogonalization to combat the effects of inexact arithmetic
might or might not be necessary. The matrix B k may be ill-conditioned because some
of its singular values approximate some of the small singular values of A. Therefore
solving the projected problem might not yield a good solution y (k) . However, we
can use any of the methods of Section 3 to regularize this projected problem; we
discuss options in detail below. As alluded to in x4, the idea is to generate y (k)
reg , the
regularized solution to (18), and then to compute a regularized solution to (16) as
reg .
If we used the algorithm GMRES instead of LSQR, we would derive similar
relations. Here, though, the U and V matrices are identical and the B matrix is
upper Hessenberg rather than bidiagonal. Conjugate gradients would yield similar
relationships.
For cost comparisons for these methods, see Tables 1 and 2. Storage comparisons
are given in Tables 3 and 4.
4.1. Regularization by projection. As mentioned earlier, if we terminate the
iteration after k steps, we have projected the solution onto a k dimensional subspace
and this has a regularizing effect that is sometimes sufficient. Determining the best
value of k can be accomplished, for instance, by one of our three methods of parameter
choice:
1. Discrepancy Principle.
In this case, we stop the iteration for the smallest value of k for which kr k k
ffi . Both LSQR and GMRES have recurrence relations for determining kr k k
using scalar computations, without computing either r k or x k [29, 32].
2. GCV.
For the projected problems (see x4.1) defined by either LSQR or GMRES,
the operator AA ] is given by
U
is the pseudo-inverse of the matrix B k . Thus from (11), the GCV
functional is [17]
We note that there are in fact two distinct definitions for B y
and hence two
definitions for the denominator in G(k); for small enough k, the two are
comparable, and the definition we use here is less expensive to calculate [18,
x7.4].
3. L-Curve.
To determine the L-curve associated with LSQR or GMRES, estimates of
are needed for several values of k. Using either algorithm,
we can compute kr k k 2 with only a few scalar calculations. Paige and Saunders
give a similar method for computing kx k k 2 [29], but, with GMRES, the cost
for computing In using this method or GCV, one must go a
few iterations beyond the optimal k in order to verify the optimum [19].
4.2. Regularization by projection plus TSVD. If projection alone does not
regularize, then we can compute the TSVD regularized solution to the projected
problem (19). We need the SVD of the . This requires O(k 3 )
operations, but can also be computed from the SVD of B k\Gamma1 in O(k 2 ) operations [13].
Clearly, we still need to use some type of parameter selection technique to find a
good value of '(k). First, notice that it is easy to compute the norms of the residual
and the solution resulting from neglecting the ' smallest singular values. If  jk is the
component of e 1 in the direction of the j-th left singular vector of B k , and if fl j is
the j-th singular value (ordered largest to smallest), then the residual and solution
2-norms are
and fi 1@ k\Gamma'(k) X
Using this fact, we can use any of our three sample methods:
1. Discrepancy Principle.
Let r (k)
denote the quantity b \Gamma Ax (k)
and note that by (13) and orthonor-
mality, kr (k)
k 2 is equal to the first quantity in (20). Therefore, we choose
'(k) to be the largest value for which kr (k)
if such a value exists.
2. GCV.
Another alternative for choosing '(k) is to use GCV to compute '(k) for
the projected problem. The GCV functional for the kth projected problem
is obtained by substituting B k for A and B
for A ] , and substituting the
expression of the residual in (20) for the numerator in (11):
3. L-Curve.
We now have many L-curves, one for each value of k. The coordinate values
in (20) form the discrete L-curve for a given k, from which the desired value
of '(k) can be chosen without forming the approximate solutions or residuals.
As k increases, the value '(k) chosen by the Discrepancy Principle will be monotonically
nondecreasing.
4.3. Regularization by projection plus Rust's TSVD. As in standard TSVD,
to use Rust's version of TSVD for regularization of the projected problem requires
that we compute the SVD of the . Using the previous notation,
Rust's strategy is to set
y
ae
ae
ik
where q (k)
are the right singular vectors of B k and I (k)
aeg. We
focus on three ways to determine ae:
1. Discrepancy Principle.
Using the notation from the previous section, the norm of the regularized solution
is given by fi 1 (
ae
ik
According to the discrepancy
principle, we must choose ae so that the residual is less than  ffi . In practice,
this would require that the residual be evaluated by sorting the values j ik j
and adding terms in that order until the residual norm is less than  ffi .
2. GCV.
Let us denote by card(I (k)
ae ) the cardinality of the set I (k)
ae . From (11), it is
easy to show that the GCV functional corresponding to the projected problem
for this regularization technique is given by
ae
ik
ae
In practice, for each k we first sort the values j ik smallest
to largest. Then we define k discrete values ae j to be equal to these values
with ae 1 being the smallest. We set ae that because the values of
are the sorted magnitudes of the SVD expansion coefficients,
we have
(j
Finally, we take the regularization parameter to be the ae j for which G k (ae j )
is a minimum.
3. L-Curve.
As with standard TSVD, we now have one L-curve for each value of k. For
fixed k, if we define the ae as we did for GCV above and we
reorder the fl i in the same way that the j ik j were reordered when sorted,
then we have
When these solution and residual norms are plotted against each other as
functions of ae, the value of ae j corresponding to the corner is selected as the
regularization parameter.
4.4. Regularization by projection plus Tikhonov. Finally, let us consider
using Tikhonov regularization to regularize the projected problem (18) for some integer
k. Thus, for a given regularization parameter , we would like to solve
min y
or, equivalently,
min y
I
The solution y
to either formulation satisfies
Using (13) and (15), we see that y (k)
also satisfies
A   AV k
A   b:
Therefore,
y
A
I
Using x
, we have
Thus as k ! n, the backprojected regularized solution x
approaches the solution
to (4).
We need to address how to choose a suitable value of .
1. Discrepancy Principle.
Note that in exact arithmetic, we have
r
Hence kB k y (k)
. Therefore, to use the Discrepancy Principle
requires we choose  so that kr (k)
discrete trial values  j .
For a given k, we take  to be the largest value  j for which kr (k)
it exists; if not, we increase k and test again.
2. GCV.
Let us define (B k ) y
to be the operator mapping the right hand side of the
projected problem onto the regularized solution of the projected problem:
Given the SVD of B k as above, the denominator in the GCV functional
defined for the projected problem (refer to (11)) is@ k
The numerator is simply kr (k)
2 . For values of k  n, it is feasible to compute
the singular values of B k .
3. L-Curve.
The L-curve is comprised of the points (kB k y (k)
using
(25) and the orthonormality of the columns of V k , we see these points are
precisely (kr (k)
discrete values of ,  the
quantities kr (k)
k 2 can be obtained by updating their respective
estimates at the (k \Gamma 1)st iteration. 2
4.5. Correspondence between Direct Regularization and Projection
Plus Regularization. In this section, we argue why the projection plus regularization
approaches can be expected to yield regularized solutions nearly equivalent to
the direct regularization counterpart. The following theorem establishes the desired
result for the case of Tikhonov vs. projection plus Tikhonov.
Theorem 4.1. Fix  ? 0 and define x
to be the kth iterate of conjugate
gradients applied to the Tikhonov problem
Let y (k)
be the exact solution to the regularized projected problem
are derived from the original problem A   A = A   b, and set z (k)
Then z
Proof: By the discussion at the beginning of x4.4 and equations (23) and (24), it
follows that y (k)
solves
A   b:
Now the columns of V k are the Lanczos vectors with respect to the matrix A   A and
right-hand side A   b. But these are the same as the Lanczos vectors generated with
respect to the matrix A   I and right-hand side A   b. Therefore V k y (k)
is precisely
the kth iterate of conjugate gradients applied to pg. 495].
Hence z (k)
. 2
2 The technical details of the approach are found in [28, pp. 197-198], from which we obtain
. The implementation details for estimating kx (k)
k and kr (k)
were
taken from the Paige and Saunders algorithm at http://www.netlib.org/linalg/lsqr.
Projection plus - Disc. GCV L-curve

Table
Summary of flops for projection plus inner regularization with various parameter selection
techniques, in addition to the O(qk) flops required for projection itself. Here k is the number of
iterations (ie. the size of the projection) taken and p is the number of discrete parameters that must
be tried.
Let us turn to the case of TSVD regularization applied to the original problem
vs. the projection plus TSVD approach. Direct computation convinces us that the
two methods compute the same regularized solution if and arithmetic is exact.
An approximate result holds in exact arithmetic when we take k iterations, with
n. Let the singular value decomposition of B k be denoted by
and define the s \Theta j matrix W s;j as
I
Then the regularized solution obtained from the TSVD regularization of the projected
problem is
reg
denotes the leading j \Theta j principle submatrix of \Gamma k . If k is taken to be
enough larger than j so that V k Q k W k;j  "
U T and
the leading principle submatrix of \Sigma, then we expect x (k)
reg to be a
good approximation to x ' . This is made more precise in the following theorem.
Theorem 4.2. Let k ? j such that
contain the first j columns of "
U respectively. Let
Then
reg
kbk:
Proof: Using the representations x
2 )b, we obtain
reg
and the conclusion follows from bounding each term. 2
Note that typically oe j AE oe n so that 1=oe j is not too large. For some results
relating to the value of k necessary for the hypothesis of the theorem to hold, the
interested reader is referred to theory of the Kaniel-Paige and Saad [30, x12.4].
Basic cost Added Cost
Disc. GCV L-curve
TSVD O("q) O(1) O(m) O(m)
Rust's TSVD O("q) O(m) O(m) O(m)
Projection O(kn) O(1) O(k) O(k)

Table
Summary of additional storage for each of four regularization methods under each of three
parameter selection techniques. The original matrix is m \Theta n with q nonzeros, p is the number of
discrete parameters that must be tried, k iterations are used in projection, and the factorizations are
assumed to take " q storage.
Projection plus - Disc. GCV L-curve
Rust's TSVD O(k) O(k

Table
Summary of storage, not including storage for the matrix, for projection plus inner regularization
approach, various parameter selection techniques. Here p denotes the number of discrete
parameters tried. Each of these regularization methods also requires us to save the basis V or else
regenerate it in order to reconstruct x.
5. Numerical results. In this section, we present two numerical examples. All
experiments were carried out using Matlab and Hansen's Regularization Tools [16],
with IEEE double precision floating point arithmetic. Since the exact, noise-free
solutions were known in both examples, we evaluated the methods using the two-
norm difference between the regularized solutions and the exact solutions. In both
examples when we applied Rust's method to the original problem, the ae i were taken
to be the magnitudes of the spectral coefficients of b sorted in increasing order.
5.1. Example 1. The 200\Theta200 matrix A and true solution x true for this example
were generated using the function baart in Hansen's Regularization Toolbox. We
generated true and then computed the noisy vector b as b + e, where e was
generated using the Matlab randn function and was scaled so that the noise level,
kb truek
. The condition number of A was on the order of 10 19 .
Many values of  were tested: log displays the
values of the regularization parameters chosen when the three parameter selection
techniques were applied together with one of the four regularization methods on the
original problem. Since 5:3761E\Gamma4, we set  ffi that defines the discrepancy
principle as the very close approximation 5:5E\Gamma4.
The last column in the table gives the value of the parameter that yielded a
regularized solution with the minimum relative error when compared against the true
solution. The relative error values for regularized solutions corresponding to the
parameters in Table 5 are given in Table 6. Note that using GCV to determine a
regularization parameter for Rust's TSVD resulted in an extremely noisy solution
with huge error.
The corners of the L-curves for the Tikhonov, projection, and TSVD methods
were determined using Hansen's lcorner function, with the modification that points
corresponding to solution norms greater than 10 6 for the TSVD methods were not
Rust's TSVD ae 1:223E\Gamma4 9:645E\Gamma7 1:223E\Gamma4 1:259E\Gamma4 or 1:223E\Gamma4
Projection

Table
Example 1: parameter values selected for each method.
Disc. GCV L-curve optimal
Rust's TSVD .1213 7E+14 .1213 .1213
Projection .1134 .1207 .1134 .1134

Table
Example 1: comparison of kx true for each of 4 regularization methods on
the original problem, where the regularization method was chosen using methods indicated.
considered (otherwise, a false corner resulted).
Next, we projected using LSQR and then regularized the projected problem with
one of the three regularization methods considered. For each of the three methods,
we computed regularization parameters for the projected problem using Discrepancy,
GCV, and L-curve, then computed the corresponding regularized solutions; the parameters
that were selected in each case at iterations 10 and 40 are given in Tables 7
and 9 respectively. As before, the lcorner routine was used to determine the corners
of the respective L-curves.
Comparing Table 6 and 8, we observe that computing the regularized solution
via projection plus Tikhonov for projection size of 10 using either the Discrepancy
Principle or the L-curve to find the regularization parameter gives results as good as if
those techniques had been used with Tikhonov on the original problem to determine
a regularized solution. Similar statements can be made for projection plus TSVD
and projection plus Rust's TSVD. We should also note that for Tikhonov, with and
without projection, none of the errors in the tables is optimal; that is, no parameter
selection techniques ever gave the parameter for which the error was minimal.
5.2. Example 2. The 255 \Theta 255 matrix A for this example was a symmetric
Toeplitz matrix with bandwidth 16 and exponential decay across the band. 3 The
true solution vector x true is displayed as the top picture in Figure 2. We generated
true and then computed the noisy vector b as b + e, where e was generated
using the Matlab randn function and was scaled so that the noise level, kek
kb truek
, was
. The vector b is shown in the bottom of Figure 2. The condition number of A
was
We generated our discrete  i using log \Gamma1. The norm of the
noise vector was 7:16E\Gamma2, so we took the value of  ffi that defines the discrepancy
principle to be 8:00E\Gamma2.
In this example, it took 61 iterations for LSQR to reach a minimum relative error
of 9:48E\Gamma2, and several more iterations were needed for the L-curve method to
3 It was generated using the Matlab command
Rust's TSVD ae(k) 1:679E\Gamma4 1:773E\Gamma4 1:679E\Gamma5 1:679E\Gamma5

Table
Example 1, iteration 10: regularization parameters selected for projection plus Tikhonov,
TSVD, and Rust's TSVD.
Disc. GCV L-curve optimal
Rust's TSVD .1213 .1663 .1213 .1213

Table
Example 1, iteration 10: comparison of kx true projection plus Tikhonov,
TSVD, and Rust's TSVD.
Disc. GCV L-curve optimal
Rust's TSVD ae(k) 9:201E\Gamma5 1:225E\Gamma4 9:201E\Gamma5 9:201E\Gamma5

Table
Example 1, iteration 40: regularization parameters selected for projection plus Tikhonov,
TSVD, and Rust's TSVD.
Disc. GCV L-curve optimal
Rust's TSVD .1162 .1162 .1162 .1162

Table
Example 1, iteration 40: comparison of kx true projection plus Tikhonov,
TSVD, and Rust's TSVD.
50 100 150 200 250
-22610exact solution
50 100 150 200 250
-226
Fig. 2. Example 2: Top: exact solution. Bottom: noisy right hand side b.
Rust's TSVD ae 2:183E\Gamma2 2:586E\Gamma6 1:477E\Gamma2 1:527E\Gamma2
Projection

Table
Example 2: parameter values selected for each method. The projection was performed on a left
preconditioned system.
Disc. GCV L-curve optimal
Rust's TSVD
Projection

Table
Example 2: comparison of kx true for each of 4 regularization methods on
the original problem.
estimate a stopping parameter. Likewise, the dimension k of the projected problem
had to be around 60 to obtain good results with the projection-plus-regularization ap-
proaches, and much larger than 60 for the L-curve applied to the projected, Tikhonov
regularized problem to give a good estimate of the corner with respect to the Tikhonov
regularized original problem. Therefore, for the projection based techniques, we chose
to work with a left preconditioned system (refer to the discussion at the end of x 3.5).
Our preconditioner was chosen as in [22] where the parameter defining the preconditioner
was taken to be
The values of the regularization parameters chosen when the three parameter
selection techniques were applied together with one of the four regularization methods
on the original problem are given in Table 11. The last column in the table gives the
value of the parameter that gave a regularized solution with the minimum relative
error over the range of discrete values tested, with respect to the true solution. The
relative errors that resulted from computing solutions according to the parameters in

Table

11 are in Table 12. We note that GCV with TSVD and Rust's TSVD were
ineffective.
The corners of the L-curves for the Tikhonov, projection, and TSVD methods
were determined using Hansen's lcorner function, with the modification that points
corresponding to the largest solution norms for the TSVD methods were not considered
(otherwise, a false corner was detected by the lcorner routine).
Next, we projected using LSQR (note that since the matrix and preconditioner
were symmetric, we could have used MINRES as in [22]) and then regularized the
projected problem with one of the three methods considered. For each of the three
methods, we computed regularization parameters for the projected problem using Dis-
crepancy, GCV, and L-curve, then computed the corresponding regularized solutions;
the parameters that were selected in each case at iterations 15 and 25 are given in Tables
13 and 15, respectively. The relative errors of the regularized solutions generated
accordingly are given in Tables 14 and 16.
Again, we used the lcorner routine to determine the corners of the respective
L-curves, except in the case of Rust's TSVD method. In the latter case, there was
Rust's TSVD ae(k) 3:558E\Gamma2 3:558E\Gamma2 3:558E\Gamma2 3:558E\Gamma2

Table
Example 2, iteration 15: regularization parameters selected for projection plus Tikhonov,
TSVD, and Rust's TSVD.
Disc. GCV L-curve optimal
Rust's TSVD

Table
Example 2, iteration 15: comparison of kx true projection plus Tikhonov,
TSVD, and Rust's TSVD.
always a very sharp corner that could be picked out visually.
Comparing Table 11 with Tables 13 and 15, we see that the parameter chosen
by applying the L-curve method to projected-plus-Tikhonov problem was the same
parameter chosen by applying the L-curve to the original problem. Moreover, a comparison
of Table 12 with Tables 14 and 16 shows that relative errors of the regularized
solutions computed accordingly are comparable to applying Tikhonov to the original
problem with that same parameter. Similar results are shown for the other cases,
with the exception that the discrepancy principle did not work well for the projection-
plus-TSVD problems, and GCV was not effective for the projected problems when
6. Conclusions. In this work we have given methods for determining the regularization
parameter and regularized solution to the original problem based on regularizing
a projected problem. The proposed approach of applying regularization
and parameter selection techniques to a projected problem is economical in time
and storage. We presented results that in fact the regularized solution obtained by
backprojecting the TSVD or Tikhonov solution to the projected problem is almost
equivalent to applying TSVD or Tikhonov to the original problem, where "almost"
depends on the size of k. The examples indicate the practicality of the method, and
illustrate that our regularized solutions are usually as good as those computed using
the original system and can be computed in a fraction of the time, using a fraction of
the storage. We note that similar approaches are valid using other Krylov subspace
methods for computing the projected problem.
In this work, we did not address potential problems from loss of orthogonality
as the iterations progress. In this discussion, we did, however, assume that either k
was naturally very small compared to n or that preconditioning had been applied to
enforce this condition. Possibly for this reason, we found that for modest k, round-off
did not appear to degrade either the LSQR estimates of the residual and solution
norms or the computed regularized solution in the following sense: the regularization
parameters chosen via the projection-regularization and the corresponding regularized
solutions were comparable to those chosen and generated for the original discretized
problem.
For the Tikhonov approach in this paper, we have assumed that the regularization
Disc. GCV L-curve optimal
Rust's TSVD ae(k) 4:828E\Gamma2 7:806E\Gamma3 4:828E\Gamma2 4:828E\Gamma2

Table
Example 2, iteration 25: regularization parameters selected for projection plus Tikhonov,
TSVD, and Rust's TSVD.
Disc. GCV L-curve optimal
Rust's TSVD

Table
Example 2, iteration 25: comparison of kx true projection plus Tikhonov,
TSVD, and Rust's TSVD.
operator L was the identity or was related to the preconditioning operator; this allowed
us to efficiently compute kr (k)
k and kx (k)
k for multiple values of  efficiently for each k.
If L is not the identity but is invertible, we can first implicitly transform the problem
to "standard form" [17]. With
Lx, we can solve the equivalent system
min
Then the projection plus regularization schemes may be applied to this transformed
problem. Clearly the projection based schemes will be useful as long as solving systems
involving L can be done efficiently.



--R



Estimation of the L-curve via Lanczos bidiagonal- ization
Galerkin projection method for solving multiple linear systems
The 'minimum reconstruction error' choice of regularization pa- rameters: Some more efficient methods and their application to deconvolution problems
Using the L-curve for determining optimal regularization pa- rameters
Equivalence of regularization and truncated iteration in the solution of ill-posed image reconstruction problems
Fast CG-based methods for Tikhonov-Phillips regularization
Generalized cross-validation as a method for choosing a good ridge parameter
Calculating the singular values and pseudo-inverse of a matrix
Matrix Computations
Theory of Tikhonov Regularization for Fredholm equations of the First Kind
A stable and fast algorithm for updating the singular value decom- position
Preconditioned iterative regularization for ill-posed problems
Analysis of discrete ill-posed problems by means of the L-curve
a Matlab package for analysis and solution of discrete ill-posed problems


The use of the L-curve in the regularization of discrete ill-posed problems
Regularization for Applied Inverse and Ill-Posed Problems
Regularization of ill-posed problems by envelope guided conjugate gradients
Symmetric Cauchy-like preconditioners for the regularized solution of 1-d ill-posed problems

Pivoted Cauchy-like preconditioners for regularized solution of ill-posed problems
On the solution of functional equations by the method of regularization
Iterative image restoration using approximate inverse preconditioning
A bidiagonalization-regularization procedure for large scale discretization of ill-posed problems
Algorithm 583

The Symmetric Eigenvalue Problem
Truncating the singular value decomposition for ill-posed problems
Iterative Methods for Sparse Linear Systems
GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems
Solutions of Ill-Posed Problems
Pitfalls in the numerical solution of linear ill-posed problems

--TR

--CTR
Angelika Bunse-Gerstner , Valia Guerra-Ones , Humberto Madrid de La Vega, An improved preconditioned LSQR for discrete ill-posed problems, Mathematics and Computers in Simulation, v.73 n.1, p.65-75, 6 November 2006
E. T. F. Santos , A. Bassrei, L- and -curve approaches for the selection of regularization parameter in geophysical diffraction tomography, Computers & Geosciences, v.33 n.5, p.618-629, May, 2007
G. Landi, The Lagrange method for the regularization of discrete ill-posed problems, Computational Optimization and Applications, v.39 n.3, p.347-368, April 2008
Alexander B. Konovalov , Vitaly V. Vlasov , Olga V. Kravtsenyuk , Vladimir V. Lyubimov, Space-varying iterative restoration of diffuse optical tomograms reconstructed by the photon average trajectories method, EURASIP Journal on Applied Signal Processing, v.2007 n.1, p.18-18, 1 January 2007
projection;l-curve;truncated singular value decomposition;regularization;discrepancy principle;iterative methods;ill-posed problems;tikhonov;krylov subspace