Computing smallest singular triplets with implicitly restarted Lanczos bidiagonalization.
A matrix-free algorithm, IRLANB, for the efficient computation of the smallest singular triplets of large and possibly sparse matrices is described. Key characteristics of the approach are its use of Lanczos bidiagonalization, implicit restarting, and harmonic Ritz values. The algorithm also uses a deflation strategy that can be applied directly on Lanczos bidiagonalization. A refinement postprocessing phase is applied to the converged singular vectors. The computational costs of the above techniques are kept small as they make direct use of the bidiagonal form obtained in the course of the Lanczos factorization. Several numerical experiments with the method are presented that illustrate its effectiveness and indicate that it performs well compared to existing codes.
Introduction
Consider the singular value decomposition (SVD) of a matrix
loss of generality,
m. Denote its singular triplets by
In this paper we are interested in
computing few, say k, of the smallest singular triplets of a general large sparse
matrix. This problem arises in several important applications including image
and signal processing [38], control [6] and matrix pseudospectra [37].
The computation of few extremal singular triplets of large sparse matrices
has been the focus of many research e#orts, see [3,9,22,21,30,33,35] as well as
Preprint submitted to Elsevier Science
[2,10,29,36,15] and numerous references therein. Recent needs in applications
such as the ones mentioned earlier, however, have motivated research oriented
towards the development of algorithms for the computation of the smallest
singular triplets, a problem that is acknowledged to challenge the capabilities
of current state-of-the-art software, e.g. see [1,5,7,11,13,14,24].
It is common practice to approximate singular values by computing the eigen-values
of equivalent hermitian eigenproblems. Furthermore, since computing
the smallest eigenvalues values of a matrix is equivalent to computing the
largest eigenvalues values of its inverse, significant work has been done on
"shift-and-invert" techniques. For example, this approach was adopted in the
MATLAB (version svds routine, that is based on ARPACK ([22]); the latter,
implements one of the most successful theoretical frameworks for the e#ective
implicitly restarted Arnoldi technique, based on seminal work of Sorensen,
Lehoucq and collaborators. However, as the size of the matrices increases, this
approach becomes too expensive in terms of storage and computational costs,
as it requires the factorization of and solution with large sparse, possibly indefinite
matrices. Developments that attempt to remedy this problem concern
inexact inverse iteration and inexact inverse Lanczos methods (see for example
[19] and [2, Sec. 11.2]). An alternative approach that avoids such solves and
is frequently e#ective is based on the use of harmonic Ritz values [26,33].
In this paper we propose and investigate an algorithm, we call IRLANB, that is
based on Lanczos bidiagonalization (LBD), a method for computing singular
values originally due to Golub and Kahan [8]. This is a matrix-free method
for the computation of the singular triplets, thus the only operations with A
are matrix vector multiplications with it and its hermitian adjoint A # . We
enhance the LBD algorithm with state-of-the-art technology for the e#ective
computation of few small singular triplets of large and possibly sparse ma-
trices. These improvements are described in the paper, whose structure is as
follows. In Section 2 we review Lanczos bidiagonalization and describe its limitations
when deployed to compute the smallest singular triplets. In Section 3
we show how to incorporate implicit restarts, introduced in [35], that permit
Lanczos bidiagonalization to maintain limited storage and computational requirements
per restart. In Section 4 we study the use of Ritz and harmonic
Ritz values as implicit shifts. In Section 5 we show how to apply the orthogonal
deflation transformation proposed in [34] in the context of Lanczos bidiago-
nalization to also make it more e#ective when singular values are clustered
singular. In Section 6 we show how to use refinement, originally proposed for
eigenvectors in [16], to enhance the computation of singular pairs. In Section
7, we describe the overall structure of IRLANB. Finally, in Section 8 we describe
numerical experiments that illustrate the behavior of IRLANB in various cases
and compares its performance with related methods.
1.1 Definitions and Equivalent Symmetric Eigenproblems
The following well known connections (see e.g. [10, Sec. 8.6]) between the SVD
and the eigendecompositions of the following hermitian matrices
are fundamental to our discussion:
| {z }
n-m
Partitioning U as columns and setting
Y is an orthonormal eigenbasis for the augmented matrix C and
| {z }
n-m
The above equivalences provide a convenient framework when seeking few singular
values of large matrices because they permit the computation of singular
triplets using symmetric eigensolvers as black box. The problem has been studied
in the literature and there exist several software packages for its solution.
Nevertheless, when seeking few small singular triplets, as we do in this paper,
several complications arise that must be addressed [28], [33].
In particular, since we are interested in the smallest singular values of A, equivalent
targets are the smallest eigenvalues of either AA # or A # A, or interior
eigenvalues of C (in the latter two, excluding spurious zeroes). Observe that,
while squaring the singular values of A will induce an increase of the separation
of the largest ones, it will also cause a corresponding clustering of the smallest
ones; this can cause problems to symmetric eigensolvers [27]. Furthermore, if A
is ill-conditioned, and we denote by #(A) its condition number with respect to
the 2-norm, the squaring of the condition number,
is likely to cause significant loss of accuracy for small singular values. Note
Algorithm Lanczos bidiagonalization
Input: A # C m-n , starting vector p 0 # C m and scalar k
Output: Bidiagonal matrix B k # R (k+1)-k and orthogonal bases
1.
2. for
3. r
4.
5.
7.
8.
9. end

Table
Lanczos bidiagonalization (cf. [2,8,21]) The real scalars # i , # i are the diagonal and
subdiagonal elements respectively, of the bidiagonal matrix B k .
that for nonsquare matrices the above analysis holds if we refer instead to the
"e#ective denotes the pseudoiverse of A (see
[13], [14]. If, on the basis of relation (3), we select instead to recover the singular
triplets of A from the eigenvalues of the augmented matrix C, we have
to approximate interior eigenvalues. Unfortunately, such a computation also
challenges the performance of symmetric eigensolvers, e.g. their convergence
behavior becomes irregular [27]. Furthermore, since each singular value corresponds
to an eigenvalue pair, -# i , symmetric eigensolvers tend to take twice
the number of iterations. An additional di#culty stems from the increased
length (m + n) of the basis vectors and corresponding increase in the storage
requirements, from which approximations to the singular values are drawn.
We next describe Lanczos bidiagonalization (LBD) that holds a central role in
our framework. LBD was originally proposed by Golub and Kahan (cf. [8] and
[10, Sec. 9.3.3]) as a process for transforming a matrix A # C m-n to upper
bidiagonal In line with the bidiagonalization algorithms
presented elsewhere in the literature, we will consider a version of the process
that transforms A to lower bidiagonal form. In fact, our discussion owes a lot
to the work of Larsen in [21]. After k < m (successful) steps, LBD produces
two blocks of Lanczos vectors
whose columns are orthonormal bases for the Krylov subspaces K k+1 (AA # ,
respectively 1 and satisfy the following relations:
where the matrix B k # R (k+1)-k has real elements and is lower bidiagonal:
. # k
. (6)

Table

1 provides an algorithmic outline of LBD. Following the execution of
LBD, the singular values of B k could be used as approximations to the singular
values of A. If we premultiply both sides of (5) with A and use (4) we obtain
However, from the LBD algorithm (cf. lines 6-8 of Table 1) we can also write
Matrix
real symmetric and tridiagonal, therefore, in
exact arithmetic, relation (8) is a symmetric Lanczos factorization and hence
LBD is equivalent to symmetric Lanczos iteration on AA # .
It is also known that there is an equivalence between LBD applied on A and
Lanczos applied on the augmented matrix C [10, Sec. 9.3.2]. In particular,
consider the starting vector
| {z }
Following standard terminology, Km (A, r) # span{r, Ar, ., A m-1 r}.
After 2k steps of Lanczos with starting vector q 1 the following relation holds:
where q
. # k
After an odd-even permutation of rows and columns of (9), we obtain a Lanczos
factorization that contains both LBD factorizations (4) and (5):
We next disuss some of the di#culties of the LBD algorithm. An important
di#culty with LBD, typical of Lanczos type algorithms, is the loss of orthogonality
among the basis vectors in V k and U k+1 [25]. The application of
reorthogonalization schemes can remedy the problem, though this is at an
extra computational cost. A compromise is to use partial reorthogonalization
schemes that dynamically update the level of orthogonality among the bases
vectors at each step. Recent work of Larsen has produced MATLAB codes that
implement partial reorthogonalization in the context of LBD; see [21,31,32].
In order to obtain acceptable approximations to the smallest singular triplets,
even with sophisticated schemes for partial reorthogonalization, convergence is
slow and the bases U k+1 , V k often need to become so long that computational
and storage costs become overwhelming.
As we show in the next sections, to counter these problems, we adopt implicit
restarting mechanisms to LBD that maintain computational and memory requirements
constant at each step. Furthermore, we combine implicit restarting
with harmonic Ritz values for the approximation of the smallest singular
triplets.
Algorithm Bulgechasing
Input: Tridiagonal matrix T
Output: Updated upper bidiagonal matrix
l
1. Set
2. for
3. Determine sin(#) such that
-s c
4. Apply to B l the Givens rotation from the right: B
5.
sin(#) such that
-s c
z
7. Apply to B l the Givens rotation from the left B
8. if i < l - 1 then
9.
10. end
11. end

Table
Bulgechasing algorithm (Golub-Kahan SVD step [10, Sec. 8.6.2]).
Implicitly Restarted LBD
Implicit restarting, proposed by D. Sorensen in [35] for the Arnoldi and Lanczos
iterations, through its practical implementation in ARPACK [22], is widely
acknowledged to be one of the most successful frameworks for solving very
large eigenproblems. In this section we describe how to apply this framework
in the case of LBD. Implicit restarting in the context of LBD was first studied
by Bj-orck, Grimme and Van Dooren in [4] and later Larsen combined it with
partial reorthogonalization in [?].
In Section 2 we established that LBD is equivalent to Lanczos applied on AA # ,
according to factorization (7). Therefore, after steps of LBD we can
apply p implicitly shifted QR steps on matrix T
l , which is real symmetric
and tridiagonal. Alternatively, we can apply Golub-Kahan SVD steps ([10,
Sec. 8.6.2]) directly on the bidiagonal matrix B l in order to enhance stability
[27]. The implicitly shifted QR step is applied directly on an upper bidiagonal
matrix by means of bulgechasing; cf. Table 2. The first Givens rotation (line
creates a "bulge" (i.e a nonzero element in the subdiagonal) and the trailing
Givens rotations "chase" the bulge out of the matrix in order to restore its upper
bidiagonal form. Since we work with a lower bidiagonal matrix, the update
can be written as
are
orthogonal matrices that implement Givens rotations. Therefore, by updating
the bases V l and U l+1 we can recover the bidiagonalization
l ,
This updated LBD factorization is
what we would have obtained after l steps of LBD with the special starting
vector
using shift -. If the previous procedure is repeated for p-1 shifts - 2 , - 3 , ., - p
we obtain a bidiagonalization that corresponds to the starting vector
Y
and therefore we can apply polynomial filtering with implicit restarts of LBD
as an equivalent to implicitly restarted Lanczos on AA # .
We also showed in relation (10) that LBD is equivalent to Lanczos applied on
the augmented matrix C. It is thus natural to ask whether implicitly restarted
LBD can be equivalent to implicitly restarted Lanczos on C? As is shown in
the following proposition that is stated assuming exact arithmetic, the answer
is negative.
Proposition 3.1 It is not possible, in general, to apply implicit QR steps on
the Lanczos factorization (9) of the augmented matrix C, and obtain a Lanczos
factorization that can be computed by LBD.
Proof Implicit restarts essentially perform polynomial filtering on the starting
vector u 1 . After p implicit QR steps on factorization (9), the updated Lanczos
factorization can be written as
with starting vector q
is a non-trivial polynomial of
the augmented matrix C of degree p. Observe now that the powers of C have
the following special structure
If we define the polynomials # containing strictly odd and even powers
respectively such that then for the polynomial #(C)
it holds that
Since for the starting vector it holds that q # we have that
Observe now that according to (2) it holds that
and thus #A # e (AA # )u 1
we have used the SVD of A. Since V is orthonormal, if we denote by #
diag[# 1 , ., #m ], it follows that
Notice that U # u 1 cannot be zero since U is orthonormal and has full rank.
Furthermore, for a general matrix with m distinct nonzero singular values,
the above norm would be zero only if # e
m. Since the
degree of # e is p < m, however, this can only happen if # e is identically zero.
Therefore, in general, the updated vector q
cannot have the special structure
and thus the updated Lanczos factorization on the augmented
matrix C cannot be equivalent to an LBD factorization. #
We next consider shift selection for the implicitly restarted LBD. In particular,
we examine two strategies: i) exact Ritz values and ii) exact harmonic Ritz
values.
4.1 Ritz Values
Using relation (5) and premultiplying with U #
l we see that after l = k+p steps
of LBD, the following relationship holds:
l AA # U l+1 =U #
l
l Av l+1 e # l+1 .
Applying relation (4) and considering only the first l columns of each side
it follows that U #
l AA # U
l , where -
l # R l-l denotes the square lower
subdiagonal matrix that we obtain by omitting the last row of B l . Therefore,
the squares of the singular values of the matrix -
B l are Ritz values of the
hermitian matrix AA # and therefore provide approximations to the singular
values of A. Our exact Ritz values strategy is to pick as implicit shifts the
largest p of the squared singular values of -
l . It is worth noting that since
our target is to compute singular values of -
l and not eigenvalues of -
l , we
do not expect loss of precision due to squared conditioning. Furthermore, by
not approximating squared singular values we do not aggravate any existing
clustering of the smallest singular values of A.
4.2 LBD and Harmonic Ritz Values
Ritz values readily provide a straightforward shift strategy. It is often the case,
however, that the smallest singular values of A are clustered. This is a situation
that can significantly slow down the convergence of implicitly restarted
Lanczos. In order to secure satisfactory convergence rates we can try to approximate
the smallest singular values of A by computing the largest Ritz
values of (AA # ) -1 . In the remainder of this section we would be assuming
that A is of full rank. In line with the matrix-free approach aspired to in this
paper, however, we prefer to avoid explicit computations with (AA # ) -1 . This
becomes possible using the concept of harmonic Ritz values [26]:
Definition 4.1 A value -
C is a harmonic Ritz value of a matrix A #
C m-m with respect to some linear subspace W k if -
is a Ritz value of A -1
with respect to W k .
Returning to the Lanczos factorization (7), since we are interested in the Ritz
values of (AA # ) -1 we could compute harmonic Ritz values of AA # . We do
this by means of oblique projection and the corresponding Petrov-Galerkin
condition. Our presentation in the remainder of this section owes a lot to the
discussion of Sleijpen and van der Vorst in [33] regarding harmonic Ritz values
(the reader can also refer to [14] for a relevant discussion). In particular, if the
search space U l+1 is of dimension l and the test space is W
then the corresponding Petrov-Galerkin condition becomes
# l+1 is a harmonic Ritz value of AA # . Furthermore, if U l+1 and W l+1 are
bases that span the subspaces U l+1 and W l+1 respectively, then the harmonic
Ritz values of AA # are the eigenvalues of matrix -
It is clear now how to compute the shifts for the implicit restart. At each restart
we compute the harmonic Ritz values and use as shifts the p largest ones. It
is worth noting that we are actually using an "exact shift" strategy with harmonic
rather than ordinary Ritz values. As we show next, the harmonic Ritz
values can be computed from the eigenvalues of a symmetric rank-one modification
of a symmetric tridiagonal matrix. In particular, expanding relation
(11) with the help of LBD relations (4), (5), and (7) it follows that
and
l
l e l+1 e #
Notice next that C
l+1 where, similarly to the previous section,
is the square lower bidiagonal matrix -
Therefore, -
B l+1 is the Cholesky factor of C l+1 . Furthermore, the term D l+1
can be written as
Note also that e l+1 e # l+1 is idempotent. Since -
B l+1 is the Cholesky factor of
C l+1 , the eigenvalues of C
l+1 D l+1 are also eigenvalues of
therefore, if we set
The above relation shows that the harmonic Ritz values are real and non-negative
and can be computed from the eigenvalues of a symmetric rank-one
modification of a symmetric tridiagonal matrix. If needed for large l,
therefore, one could deploy fast algorithms that exploit this special structure.
It is worth noting that the harmonic Ritz values in [33] were derived from the
eigenvalues a rank one update of an Arnoldi matrix. Furthermore, the term
l+1 above is identical to
since u l+1 is orthonormal and u # l+1 Av . But from the LBD algorithm
and thus follows an alternative way of writing equation (13) above:
so the terms follow naturally after another LBD step. Furthermore, if -
# S# is the singular value decomposition of -
thus we have to compute the eigenvalues of a symmetric
rank-one modification of a diagonal matrix (see [10] sec. 8.6).
One important issue in the design of implicitly restarted Arnoldi algorithms is
the implementation of e#cient deflation techniques that enhance convergence
and stability and provide an e#ective way to compute multiple and clustered
eigenvalues. This is so as to let the methods become an e#ective alternative
to block methods. It is also worth noting that implicit restarting Arnoldi has
also been combined with block methods to deal with the computation of few
selected eigenpairs and singular triplets in an algorithm recently proposed by
Baglama, Calvetti and Reichel [1]. We thus need to consider how to implement
deflation in the context of implicitly restarted LBD. Our scheme builds
upon results presented in [2,23,34]. As in [23] we employ "locking", that de-couples
converged approximate singular values and singular subspaces and
"purging", that removes unwanted but converged singular pairs. In this section
we describe the modification and application of the "orthogonal deflating
transformation" (ODT for short), a scheme originally proposed by Sorensen
in [34] in the context of implicitly restarted Arnoldi for eigenvalues. We show
that the transformation can be applied directly on the bidiagonal matrix that
results from implicitly restarted LBD. The deflation scheme enables the efficient
stable and e#cient locking of approximate singular values that have
converged with relative accuracies that may be much larger than the machine
epsilon.
The ODT is based upon a special unitary matrix, say Q, that is built, as
shown in [34] to satisfy Qe suitably chosen unit norm vector
for the construction of Q. Furthermore, Q has the form
where R is upper triangular, its first column is zero and R # It may also
be written as
where L is lower triangular and
Assuming now that such a
Q can be built, the following lemma shows how to apply the ODT in the case
of implicitly restarted LBD.
Lemma 5.1 Let (#, y L , y R ) be an approximate singular triplet of A # C m-n
computed from the bidiagonal matrix B resulting after k steps of LBD. Let
also be the unitary
matrices produced for ODT from the vectors y L and y R respectively. Then the
updated matrix -
lower bidiagonal and has the special form
is the approximate singular value and -
B is also lower
bidiagonal.
Proof Using the same notation as above, the following relations hold for
Similarly, the following relations also hold for QR :
R , LR e
We will prove that -
upper Hessenberg as well as lower trian-
gular, and therefore lower bidiagonal. In particular,
since Q # L y Therefore,
R RR
since y #
relations (19). Matrix L # L BRR +#e 1 e # 1 is upper Hessenberg
because L # L and RR are upper triangular and B is lower bidiagonal, thus -
B is
upper Hessenberg. Furthermore,
R
R
since Q # L y of (15). Therefore,
(R
R
since y # R LR
because of (20). Since both R # L and LR are lower tri-
angular, B would be lower bidiagonal while the rank-one update would not
modify the lower triangular form, therefore -
B is also lower triangular. #
It is worth noting that the observations concerning the numerical stability
of ODT discussed in [34] carry over to the present case. In particular, note
that matrices Q L , QR are built from y L and y R respectively, therefore, some of
their implicit properties are not exactly satisfied in finite precision arithmetic.
Therefore, in order for -
BQR to be numerically upper Hessenberg, special
care must be taken so that #g L (y # L B)RR # 2 would remain small in practice.
If we write y #
that #g L (y # L B)RR
1 denotes the first component of
y L . Unfortunately, for small values of # L
1 the above factor could be large and
a rescaling strategy, such as the one described in [34], must be applied. On
the other hand, R #
R , if we apply the aforementioned rescaling strategy, the
norm #g R #
# Ris kept small and therefore -
B would be numerically lower
triangular since #e 1 z # LR # 2 will be small for small z.
6 Refined Singular Vector Approximations
It is often the case when computing eigenvalues that a Ritz vector may exhibit
poor converge even though the corresponding Ritz value has converged. Jia
proposed in [16] a refined Ritz vector strategy. The key is to approximate the
eigenvector by means of a refined Ritz vector designed to minimize the norm
of the residual over the subspace involved. We first outline the refinement
process in the case of Lanczos from which follows naturally its application to
LBD.
Let us assume that we have performed steps of Lanczos, - # is the
approximation to an eigenvalue and -
l z is the corresponding refined Ritz
vector that is extracted from a Krylov subspace K l (A, v 1 ). Moreover, let
be the corresponding Lanczos factorization, where V l # C m-k is the basis of
the Krylov subspace and -
(l+1)-l is the augmented tridiagonal matrix.
We seek to find an approximation -
u that minimizes the norm of the residual.
Therefore, the following relations hold:
min
z#C l #AV l z - #V l z#
z#C l #V l+1
z#C l #V l+1 ( -
z#C l #( -
I),
since the norm of the residual is minimized when z is the right singular vector
associated with the smallest singular value # min ( -
I). This singular value
is called the refined residual. Jia reports that the angle between the refined
Ritz vector - z and the exact eigenvector is better than the corresponding
angle of the standard Ritz vector. Furthermore, notice that we can use the
Rayleigh quotient in an attempt to obtain an improved eigenvalue,
since # may be more accurate than - #; cf. [36, Sec. 4.3].
In Section 3 we established that implicitly restarted Lanczos on C cannot be
equivalent to implicitly restarted LBD while implicitly restarted Lanczos on
remains equivalent to implicitly restarted LBD. However, in Section 2 we
saw that the LBD decompositions are equivalent to Lanczos decompositions
on either AA # or the augmented matrix C (with a starting vector of special
structure). Therefore, we can compute the refined residual (and vector) using
either AA # or C.
Decomposition (7) suggests that if -
# min is the current approximation to the
smallest singular value of A, the refined residual and refined singular vector
can be retrieved by computing the smallest singular value and right singular
vector of
min
I. (21)
In the case of the augmented matrix C, according to the Lanczos decomposition
(10) we obtain the refined residual and refined singular vectors by
computing the smallest singular value and right singular vector of the matrix
l 0
I. (22)
We next have to decide which refined residual to compute, the one from AA #
or from C? Since (21) involves the tridiagonal matrix BB # one might expect
stability problems in contrast to (22). Furthermore, the refined residual for
approximations only to the left singular vector so that to obtain
approximations to the right singular vector we would need to use the relation
or also work with the refined residual of A # A. It is thus
preferable to use the augmented matrix C which also facilitates the concurrent
approximation of both the left and right singular vectors of A. For more details,
see also the discussion in [36, Sec. 4.3].
Algorithm IRLANB
Input matrix A # C m-n , k, p, eignum. Starting vector u 1 . Set
Output eignum of the smallest singular triplets
1. Compute bases U l+1 and V l and bidiagonal B l using LBD
2. Repeat
3. if (shifts == Ritz) then
4. Compute the eigenvalues # i ,
l
5. elseif (shifts == Harmonic )
6. Compute the eigenvalues # i ,
7. end
8. Perform p implicit QR steps using bulgechasing on B with the p largest
eigenvalues # i as shifts and update the LBD factorization: AV
9. Compute the approximation -
10. Compute the refined residual of -
11. if -
# min has converged then
12. Compute the left and right refined singular vectors of -
13. Compute QL and QR matrices using ODT and perform deflation
14. Discard the first column of U l+1 , V l and the first row and column of B l
16. end
17. Reorthogonalize u
against all previous (even converged) basis vectors
k to length using LBD
19. Until convergence of all eignum singular values

Table
IRLANB: A method to compute few of the smallest singular triplets of large sparse
matrices.
7 IRLANB: Implicitly Restarted Harmonic Lanczos Bidiagonaliza-
tion
Based on the previous discussion, we will next construct an algorithm, we call
IRLANB and depict in Table 3, for the computation of few of the smallest singular
triplets of large sparse matrices. We will first proceed with its algorithmic
description and will then highlight some of its important fine points and its
implementation.
Parameter is the maximum dimension of the bidiagonalization,
where p is the number of implicitly shifted QR steps applied on B l . Parameter
eignum determines the number of smallest singular values that we seek. The
first step of IRLANB constructs an LBD factorization of length l. For this
purpose we have used the function lanbpro from Larsen's PROPACK [20] (see
also [21]) which is a set of MATLAB codes for the symmetric eigenvalue and
SVD problems based on Lanczos and Lanczos bidiagonalization with partial
reorthogonalization. As described in Section 4, if we select to shift with Ritz
values, we prefer, in terms of stability, to compute singular values of B l rather
than eigenvalues of B l B #
l . If, instead, we select to shift by means of harmonic
Ritz values, we could use the eigenvalues of (13).
The next step is to compute the 2-norm of the refined residual according to
either one of the strategies described in Section 6. If convergence has not taken
place, we proceed to the reorthogonalization steps (line 17) and repeat the pro-
cess. As soon as the current approximation to # min satisfies the convergence
criterion we compute the corresponding left and right refined singular vectors
and proceed with the deflation procedure. We compute the orthogonal matrices
using ODT, as described in Section 5. Purging is accomplished
by discarding the first column of the bases U
k as well as the first
row and column of B
k . As a result, we obtain an LBD factorization of length
(k-1) while the deflated factorization no longer contains the targeted singular
values. However, in subsequent restarts we reorthogonalize the updated vectors
k against all previous vectors, even purged ones, since roundo#
may introduce components towards the directions of converged vectors. Note
that since we are computing a small number of singular triplets, the extra cost
incurred is low. Computational practice indicates that this limited reorthog-
onalization su#ces to maintain orthogonality among basis vectors that may
have been degraded by the implicit restart.
8 Numerical Experiments
In this section we present numerical experiments designed to illustrate the
numerical and computational performance of IRLANB. All codes were written
in MATLAB 6.1 and ran on a 866 MHz Pentium III equipped with 1GB of
RAM and 512 Kb of cache memory running Windows 2000 Server. We also
illustrate the performance of IRLANB vs. two recent methods for which MATLAB
codes are publicly available and which are matrix-free, so as to permit the
solution of very large sparse problems in computational environments such as
the above. These methods were:
IRBLSVDS-IRBLEIGS: Code due to Baglama, Calvetti, and Reichel and based
on implicitly restarted block Lanczos ([1]) designed to compute one or more
eigenvalues and/or singular values 2 .
JDQZ: Code based on Jacobi-Davidson QZ method due to Fokkema, Sleijpen
and van de Vorst and implemented in MATLAB ([7]) 3 .
Note that if asked to compute few of the smallest singular values of sparse
matrices, the MATLAB 6 built-in function svds, that is based on a compiled
implementation of ARPACK (eigs), applies shift-and-invert and requires an LU
decomposition of the augmented matrix C. Therefore, we do not include svds
in our experiments. It is also worth noting that in [1], IRBLSVDS-IRBLEIGS was
compared to methods selected based on criteria similar to the ones described
herein.
8.1 Ritz and Harmonic Ritz Shift Strategies
The first set of experiments is designed to illustrate the convergence behavior
of Ritz values vis-a-vis that of harmonic Ritz values, when used as shifts in
the implicitly restarted LBD algorithm. We constructed a sequence of diagonal
matrices A s # R n-n , increasing clustering of
their smallest singular values. In MATLAB notation:
The test space dimension was while at each restart we performed
steps. We used a random starting vector normalized to
have unit length and convergence tolerance tol=1e-8. Figure 1 illustrates the
true relative error for Ritz shifts as well as for harmonic Ritz shifts. It is
evident that as the clustering of the smallest singular values of A s increases,
harmonic Ritz values either converge significantly faster (they require fewer
restarts) or, with the same backward error (2-norm of residual) used for the
Ritz values, produce results with better forward relative error. Therefore, in
case of severe clustering of the smallest singular values we use harmonic Ritz
values, which, as shown in Section 4.2, can be computed at relatively small
cost.
Number of restarts
Relative
error
Ritz
Harmonic Ritz
Number of restarts
Relative
error
Ritz
Harmonic Ritz
Number of restarts
Relative
error
Ritz
Harmonic Ritz
Number of restarts
Relative
error
Ritz
Harmonic Ritz
Fig. 1. Experiments with the diagonal matrices (23). Starting from the left top
corner and moving clockwise we depict relative errors for Convergence
tolerance was set to tol=1e-8. Solid lines correspond to standard Ritz shifts, dashed
lines to harmonic Ritz shifts.
Number of restarts
Abs.
of
relative
error
50 100 150 200 250
Number of restarts
Abs.
of
relative
error
Fig. 2. Experiments with the diagonal matrices (24), with increasing condition numbers
8.2 Experiments with Ill-Conditioned Matrices
We next investigate the behavior of IRLANB with harmonic Ritz values and
ill-conditioned matrices. We constructed a sequence of diagonal matrices A s #
2 Available at http://hypatia.math.uri.edu/#jbaglama.
3 Available at http://www.math.ruu.nl/people/sleijpen/JD software/JDQZ.html.
Number of restarts
of
residual
Number of restarts
of
residual
Fig. 3. Experiments with IRLANB on grcar(1000). Up: Convergence tolerance
tol=1e-6. Down: Convergence tolerance tol=1e-10.
R n-n , increasing condition numbers:
We used the same starting vector and parameters as in the previous examples
illustrates the absolute value of the
relative error achieved by IRLANB. In all cases IRLANB computed successfully
the smallest singular value. In particular, in the case of modest condition
numbers (left plot) we observed a smooth convergence behavior. However, as
the condition numbers of the matrices A s deteriorates, the behavior of the
smooth error behavior vanishes especially towards the end of the restarts.
Still, IRLANB converges with a relative error of order 10 -5 in the worst case
8.3 Computing Few Singular Values
We next illustrate the ability of IRLANB to quickly detect a few additional
singular values that lie near the smallest one, once the latter has converged.
We continue using IRLANB with harmonic Ritz values.
l
Fig. 4. Experiments with IRLANB on dw 2048. Left: Convergence tolerance tol=1e-8.
Right: Convergence tolerance tol=1e-12.
We first experiment with matrix grcar of dimension included
in MATLAB's function gallery. Our target is to compute its 10 smallest singular
values. The length of LBD was while we used
implicit shifts per step. Figure 3 illustrates the norms of the residual for each
iteration. The dashed lines represent the convergence criterion that was set
equal to normest(A)-tol, where normest(A) is an estimation of the norm of
A which we approximate by #A# 2 #B l # 2 (computed before the first restart).
We conducted two experiments. In the first case (top of figure 3) we used
convergence tolerance equal to tol=1e-6 while in the second case we used
tol=1e-10. The plots on the right of figure 3 are detailed versions of the plots
on the left. We immediately notice that IRLANB can contine computing singular
values at subsequent restarts. This behavior is even more pronounced
when we employ a stricter convergence tolerance (tol=1e-10). Notice that
after the smallest singular value has been approximated (at restart number
98), then in the subsequent 9 restarts each of the remaining singular values is
approximated. Observe that the ratio among the largest and smallest singular
values computed is #N-9
Obviously, deflation has helped to deal
e#ectively with this level of clustering.
We next experiment with matrix dw 2048 from Matrix Market 4 , the 10 smallest
singular values of which are not as clustered as in the previous case
experimented with convergence
tolerances tol=1e-8 and tol=1e-12. Figure 4 illustrates the results. As
in the experiment with grcar we observe that IRLANB rapidly approximates
the remaining singular values once convergence for the smallest one has been
achieved. Because of the decreased clustering of the smallest singular values,
however, convergence is not as fast as in the previous case.
Matrix IRLANB IRBLEIGS-IRBLSVDS JDQZ
jpwh 991 4.9 14.9 13.2 34.3 23.8 136.5
well 1850 22.2 26.8 31.7 39.1 FAILED

Table
Runtimes (in seconds) for matrices jwph 991 and well 1850 (in seconds) in order
to compute one or two of the smallest singular triplets.
8.4 Comparisons With Related Methods
In this section we provide numerical experiments that illustrate the behavior of
IRLANB vis-a-vis the methods selected above, namely IRBLEIGS-IRLBSVDS and
JDQZ. In both algorithms, the singular values are obtained via the augmented
matrix C. The first two examples are with two matrices, namely jwph 991
obtained from Matrix Market. We used the algorithms under consideration to
compute one as well as two of the smallest singular triplets. The convergence
tolerance was set to tol = 1e-6. The minimum search space dimension (k for
IRLANB, jmin for JDQZ, and BLSZ for IRBLEIGS-IRBLSVDS) was set to 3. The
maximum search space dimensions used were set as k+p=jmax=NBLS-BLSZ=15;
cf. the help pages of IRBLEIGS-IRBLSVDS and JDQZ for detailed explanation
regarding the input parameters. Table 4 illustrates the corresponding runtimes
and indicates that IRLANB competes well with modern available methods.
Our last experiment originates from the computation of pseudospectra of large
matrices. Since the #-pseudospectrum of a matrix can be defined as the locus
of points z of the complex plane that satisfy the inequality # min (zI -A) #, it
becomes of critical importance to use fast algorithms to estimate the smallest
singular value (see [37] for a comprehensive survey). We experiment with a
family of matrices, studied in [39], that originate from specific bidiagonal ones
to which we add random sparse entries. In MATLAB notation, the matrices are
defined as
where N is the size of the matrix. Specifically, we seek # min
dimensions 200000. The parameters
for IRLANB and JDQZ were: Minimum dimension of search space k=jmin=15
and maximum dimension of search space k+p=jmax=30. Convergence tolerance
was set to tol=1e-10. The corresponding parameters for IRBLEIGS-IRBLSVDS
were BLSZ=3, NBLS=10 and tol=1e-6. The maximum number of restarts was
set to MAXIT=1000. Table 5 illustrates the runtimes and convergence results.
100000 0.3719 163 0.3719 183 0.3719 429
200000 0.3737 332 0.3737 376 *
50000 1.2373e-4 182 1.2373e-4 243 1.8391e-1 2270
100000 2.4861e-4 474 2.4861e-4 561 1.8754 11000
150000 5.9600e-5 705 5.9600e-5 850 *
200000 6.2756e-5 943 6.2756e-5 1210 *

Table
Run times (in seconds) and approximations of # min for the family of random matrices
(25). The star ("*") indicates that the method ran out of memory (1 GByte).
We observe that for the shift z = 3.5, all three methods return similar results
(up to 4 digits); however, IRLANB is significantly faster. Furthermore, when
we observe that JDQZ ran out of memory for both shifts.
Finally, we note that the -
# min computed by JDQZ for the shift
di#erent than the result of the other two methods that are in agreement in 9
to digits.
9 Conclusions
In this paper we described the design of IRLANB, an implicitly restarted Lanczos
bidiagonalization algorithm for the computation of a few of the smallest
singular values of a matrix. We investigated Ritz as well as harmonic Ritz values
as shifts in the implicit QR steps and demonstrated the superiority of the
latter in the case of clustered smallest singular values. We showed how to e#-
ciently compute the harmonic Ritz values, only at a very small additional cost
compared to Ritz values. Furthermore, we demonstrated that IRLANB with
harmonic Ritz values can successfully compute the smallest singular value of
matrices with very large condition numbers. We proved that the orthogonal
deflation transformation can be applied directly on Lanczos bidiagonalization.
Numerical experiments demonstrate that this deflation scheme can e#ciently
compute clustered singular values. Finally, we demonstrated the application
of refined residuals and vectors in the case of Lanczos bidiagonalization. The
computation of the smallest singular values is a di#cult and computationally
challenging problem. We believe that the above framework will prove to be
very helpful in future investigations as well as in practical computations.

Acknowledgments
The first two authors wish to thank the Bodossaki Foundation for its financial
support. A first version of this work was developed in the context of the
first author's Diploma thesis ([18]) and presented during the International
Workshop on Parallel Matrix Algorithms and Applications (PMAA'02) held
in Neuch-atel. We thank Errikos Kontoghiorghes for his hospitality during our
stay in Neuch-atel. We would also like to thank Valeria Simoncini, for bringing
to our attention [21] and PROPACK; Henk van der Vorst for helpful discussions
during HERCMA'01 in Athens and his suggestion that we consider using
harmonic Ritz values; Gerard Sleijpen and Michiel Hochstenbach for their
useful comments concerning stability issues; and Andreas Stathopoulos for
many insightful discussions regarding all aspects of this work. The IRLANB
code is available from the authors upon request.



--R

IRBL: An Implicitly Restarted Block Lanczos Method for large-scale Hermitian eigenproblems
Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide.
Large scale singular value decomposition.

Computing the field of values and pseudospectra using the Lanczos method with continuation.
An algorithm for computing the distance to uncontrollability.

Calculating the singular values and pseudo-inverse of a matrix
Eigenvalue Comptutation in the 20th Century.
Matrix Computations.
Parallel computation of spectral portrait of large matrices by Davidson type methods.
The Test Matrix Toolbox for MATLAB (version 3.0).
A Jacobi-Davidson type SVD method
Harmonic and refined extraction methods for the singular value problem
Computational methods for large eigenvalue problems.
Refined Iterative algorithms based on Arnoldi's process for large unsymmetric eigenproblems.
An implicitly restarted refined bidiagonalization lanczos method for computing a partial singular value decomposition.
Parallel iterative methods for the computation of the smallest singular values of large sparse matrices and applications.
An inexact inverse iteration for large sparse eigenvalue problems.
PROPACK: A software package for the symmetric eigenvalue problem and singular value problems on Lanczos and Lanczos bidiagonalization with partial reorthogonalization
Lanczos bidiagonalization with partial reorthogonalization.
Arpack User's Guide: Solution of Large-Scale Eigenvalue Problems With Implicitly Restarted Arnoldi Methods
Deflation techniques for an implicitly restarted arnoldi iteration.
Parallel computation of the pseudospectrum of large matrices.
The Computation of Eigenvalues and Eigenvectors of Very Large Sparse Matrices.
Approximate solutions and eigenvalue bounds from Krylov subspaces.
The Symmetric Eigenvalue Problem.
Rational Krylov algorithms for nonsymmetric eigenvalue problems II.
Numerical Methods for Large Eigenvalue Problems.
The trace minimization method for the symmetric generalized eigenvalue problem.
The Lanczos algorithm with partial reorthogonalization.

A Jacobi-Davidson iteration method for linear eigenvalue problems
Deflation for implicitly restarted arnoldi methods.
Implicit application of polynomial filters in a k-step Arnoldi method
Matrix Algorithms.
Computation of pseudospectra.


--TR
Implicit application of polynomial filters in a k-step Arnoldi method
An algorithm for computing the distance to uncontrollability
Deflation Techniques for an Implicitly Restarted Arnoldi Iteration
Matrix computations (3rd ed.)
The symmetric eigenvalue problem
Jacobi--Davidson Style QR and QZ Algorithms for the Reduction of Matrix Pencils
A Jacobi--Davidson Iteration Method for Linear Eigenvalue Problems
Implicitly Restarted GMRES and Arnoldi Methods for Nonsymmetric Systems of Equations
LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares
Templates for the solution of algebraic eigenvalue problems
Low-Rank Matrix Approximation Using the Lanczos Bidiagonalization Process with Applications
Eigenvalue computation in the 20th century
The trace minimization method for the symmetric generalized eigenvalue problem
Matrix algorithms
A Jacobi--Davidson Type SVD Method
Large-Scale Computation of Pseudospectra Using ARPACK and Eigs
Methods for Large Scale Total Least Squares Problems
Parallel computation of pseudospectra of large sparse matrices
Lanczos Algorithms for Large Symmetric Eigenvalue Computations, Vol. 1
An Implicitly Restarted Refined Bidiagonalization Lanczos Method for Computing a Partial Singular Value Decomposition

--CTR
C. Bekas , E. Kokiopoulou , E. Gallopoulos, The design of a distributed MATLAB-based environment for computing pseudospectra, Future Generation Computer Systems, v.21 n.6, p.930-941, June 2005
deflation;implicit restarting;harmonic Ritz values;pseudospectrum;refined singular vectors;lanczos bidiagonalization