Computing Symmetric Rank-Revealing Decompositions via Triangular Factorization.
We present a family of algorithms for computing symmetric rank-revealing VSV decompositions based on triangular factorization of the matrix. The VSV decomposition consists of a middle symmetric matrix that reveals the numerical rank in having three blocks with small norm, plus an orthogonal matrix whose columns span approximations to the numerical range and null space. We show that for semidefinite matrices the VSV decomposition should be computed via the ULV decomposition, while for indefinite matrices it must be computed via a URV-like decomposition that involves hypernormal rotations.
Introduction
. Rank-revealing decompositions of general dense matrices are
widely used in signal processing and other applications where accurate and reliable
computation of the numerical rank, as well as the numerical range and null space,
are required. The singular value decomposition (SVD) is certainly a decomposition
that reveals the numerical rank, but what we have in mind here are the RRQR and
(i.e., URV and ULV) decompositions which can be computed and, in particular,
updated more eciently than the SVD. See, e.g., [7, xx2.7.5{2.7.7], [20, x2.2] and [33,
Chapter 5] for details and references to theory, algorithms, and applications.
The key to the eciency of RRQR and UTV algorithms is that they consist of an
initial triangular factorization which can be tailored to the particular matrix, followed
by a rank-revealing post-processing step. If the matrix is mn with m  n and with
numerical rank k, then the initial triangular factorization requires O(mn 2 )
ops, while
the rank-revealing step only requires O((n k)n 2 )
ops if k  n, and O(kn 2 )
ops if
n. The updating can always be done in O(n 2 )
ops, when implemented properly.
We refer to the original papers [9], [10], [16], [18], [19], [23], [31], [32] for details about
the algorithms.
For structured matrices (e.g., Hankel and Toeplitz matrices), the initial triangular
factorization in the RRQR and UTV algorithms has the same complexity as the
rank-revealing step, namely, O(mn)
ops; see [7, x8.4.2] for signal processing aspects.
However, accurate principal singular values and vectors can also be computed by
means of Lanczos methods in the same complexity, O(mn)
ops [13]. Hence the
advantage of a rank-revealing decomposition depends on the matrix structure and
the numerical rank of the matrix.
Rank-revealing decompositions of general sparse matrices are also in use, e.g., in
optimization and geometric design [27]. For sparse matrices, the initial pivoted triangular
factorization can exploit the sparsity of A. However, the UTV post-processors
may produce a severe amount of ll, while the ll in the RRQR post-processor is
P. Y. Yalamov was supported by a fellowship from the Danish Rectors' Conference and by Grants
MM-707/97 and I-702/97 from the National Scientic Research Fund of the Bulgarian Ministry of
Education and Science.
y Department of Mathematical Modelling, Technical University of Denmark, Building 321, DK-
2800 Lyngby, Denmark (pch@imm.dtu.dk).
z Center of Applied Mathematics and Informatics, University of Rousse, 7017 Rousse, Bulgaria
(yalamov@ami.ru.acad.bg).
P. C. HANSEN AND P. Y. YALAMOV
restricted to lie in the columns that are permuted to the right of the triangular factor
[7, Thm. 6.7.1]. An alternative sparse URL decomposition A = U RL, where U is
orthogonal and R and L are upper and lower triangular, respectively, was proposed
in [26]. This decomposition can be computed with less ll, at the expense of working
with only one orthogonal matrix.
Numerically rank-decient symmetric matrices also arise in many applications,
notably in signal processing and in optimization algorithms (such as those based on
interior point and continuation methods). In both areas, fast computation and efcient
updating are key issues, and sparsity is also an issue in some optimization
problems. Symmetric rank-revealing decompositions enable us to compute symmetric
rank-decient matrix approximations (obtained by neglecting blocks in the rank-
revealing decomposition with small norm). This is important, e.g., in rank-reduction
algorithms in signal processing where one wants to compute rank-decient symmetric
semidenite matrices. In addition, utilization of symmetry leads to faster algorithms,
compared to algorithms for nonsymmetric matrices.
In spite of this, very little work has been done on symmetric rank-revealing decom-
positions. Luk and Qiao [24] introduced the term VSV decomposition and proposed
an algorithm for symmetric indenite Toeplitz matrices, while Baker and DeGroat [2]
presented an algorithm for symmetric semi-denite matrices.
The purpose of this paper is to put the work in [2] and [24] into a broader perspective
by surveying possible rank-revealing VSV decompositions and algorithms,
including the underlying theory. Our emphasis is on algorithms which, in addition
to revealing the numerical rank, provide accurate estimates of the numerical range
and null space. We build our algorithms on existing methods for computing rank-
revealing decompositions of triangular matrices, based on orthogonal transformations.
Our symmetric decompositions and algorithms inherit the properties of these underlying
algorithms which are well understood today.
We emphasize that the goal of this paper is not to present detailed implementations
of our VSV algorithms, but rather to set the stage for such implementations.
The papers [4] and [28] clearly demonstrate that careful implementations of ecient
and robust mathematical software for numerically rank-decient problems requires a
major amount of research which is outside the scope of the present paper.
Our paper is organized as follows. After brie
y surveying general rank-revealing
decompositions in x2, we dene and analyze the rank-revealing VSV decomposition of
a symmetric matrix in x3. Numerical algorithms for computing VSV decompositions
of symmetric semi-denite and indenite matrices are presented in x4, and we conclude
with some numerical examples in x5.
2. General Rank-Revealing Decompositions. In this paper we restrict our
attention to real square n  n matrices. The singular value decomposition (SVD) of
a square matrix is given by
where u i and v i are the columns of the orthogonal matrices U and V , and
with
. The numerical rank k of A, with respect to the threshold  , is the number of
singular values greater than or equal to  , i.e.,  k   >  k+1 [20, x3.1].
The RRQR, URV, and ULV decompositions are given by
Here, Q, UR , UL , VR , and VL are orthogonal matrices,  is a permutation matrix, T
and R are upper triangular matrices, and L is a lower triangular matrix. Moreover,
if we partition the triangular matrices as
then the numerical rank k of A is revealed in the triangular matrices in the sense that
are k  k and
F
The rst k columns of the left matrices Q, UR , and UL span approximations to the
numerical range of A, dened as spanfu g, and the last n k columns of
the right matrices VR and VL span approximations to the numerical null-space of A,
dened as spanfv g. See, e.g., [20, x3.1] for details.
Precise denitions of RRQR decompositions and algorithms are given by Chandrasekaran
and Ipsen [11], Gu and Eisenstat [19] and Hong and Pan [23], and associated
large-scale implementations are available in Fortran [4]. Denitions of UTV
decompositions and algorithms are given by Stewart [31], [32]. Matlab software for
both RRQR and UTV decompositions is available in the UTV Tools package [17].
3. Symmetric Rank-Revealing Decompositions. For a symmetric n  n
matrix A, we need rank-revealing decompositions that inherit the symmetry of the
original matrix. In particular this is true for the eigenvalue decomposition (EVD)
are the right singular vectors, while
Corresponding to the UTV decompositions, Luk and Qiao [24] dened the following
VSV decomposition
where VS is an orthogonal matrix, and S is a symmetric matrix with partitioning
in which S 11 is k  k. We say that the VSV decomposition is rank-revealing if
4 P. C. HANSEN AND P. Y. YALAMOV
This denition is very similar to the denition used by Luk and Qiao, except that they
use ktriu(S 22 )k 2
F instead of kS 22 k 2
F , where \triu" denotes the upper triangular part.
Our choice is motivated by the fact that kS 22 k 2
n as kS 12
Given the VSV decomposition in (3.2), the rst k columns of VS and the last
columns of VS provide approximate basis vectors for the numerical range and
null space, respectively. Moreover, given the ill-conditioned problem b, we can
compute a stabilized \truncated VSV solution" x k by neglecting the three blocks in
consists of the rst k columns
of VS . We return to the computation of x k in x4.4.
Instead of working directly with the matrix S, it is more convenient to work
with a symmetric decomposition of S and, in particular, of S 11 . The form of this
decomposition depends on both the matrix A (semi-denite or indenite) and the
rank-revealing algorithm. Hence, we postpone a discussion of the particular form of
S to the presentation of the algorithms. Instead, we summarize the approximation
properties of the VSV decomposition.
Theorem 3.1. Let the VSV decompositions of A be given by (3.2), and partition
the matrix S as in (3.3) where k is the numerical rank. Then the singular values
of diag(S 11 ; S 22 ) are related to those of A as
Moreover, the angle  between the subspaces spanned by the rst k columns of V and
VS , dened by sin
bounded as
sin
Proof. The bound (3.4) follows from the standard perturbation bound for singular
values:

where we use that the singular values of the symmetric \perturbation matrix" appear
in pairs. To prove the upper bound in (3.5), we partition
columns. Moreover, we write
k. If we insert these partitionings as well as
(3.1) and (3.2) into the product AV S;0 then we obtain
Multiplying from the left with V T
k we get
from which we obtain
Taking norms in this expression and inserting sin
we get
sin    1
which immediately leads to the upper bound in (3.5). To prove the lower bound, we
use that
Taking norms and using sin
k+1 , we obtain the left bound in (3.5).
We conclude that if there is a well-dened gap between  k and kS 22 k 2 , and if the
norm kS 12 k 2 of the o-diagonal block is suciently small, then the numerical rank k
is indeed revealed in S, and the rst k columns of VS span an approximation to the
singular subspace spanfv g. The following theorem shows that a well-dened
gap is also important for the perturbation bounds.
Theorem 3.2. Let e
S , and let  denote the angle between
the subspaces spanned by the rst k columns of VS and e
sin   4
g.
Proof. The bound follows from Corollary 3.2 in [14].
We see that a small upper bound is guaranteed when kAk 2 as well as  and
k+1 are somewhat smaller than  k .
4. Algorithms for Symmetric Rank-Revealing Decompositions. Similar
to general rank-revealing algorithms, the symmetric algorithms consist of an initial
triangular factorization and a rank-revealing post-processing step. The purpose of the
latter step is to ensure that the largest k singular values are revealed in the leading
submatrix S 11 and that the corresponding singular subspace is approximated by the
span of the rst k columns of VS .
For a semi-denite matrix A, our initial factorization is the symmetrically pivoted
where P is the permutation matrix, and C is the upper triangular (or trapezoidal)
factor. The numerical properties of this algorithm are discussed by Higham
in [22]. If A is a symmetric semi-denite Toeplitz matrix, then there is good evidence
(although no strict proof) that the Cholesky factor can be computed eciently and
reliably without the need for pivoting by means of the standard Schur algorithm [30].
When A is indenite, then it would be convenient to work with an initial factorization
of the form P T
C where C is again triangular
and
diag(1).
Unfortunately such factorizations are not guaranteed to exist. Therefore our initial
factorization is the symmetrically pivoted LDL T factorization
where P is the permutation matrix, L is a unit lower triangular matrix, and D is a
block diagonal matrix with 11 and 22 blocks on the diagonal. The state-of-the-art
in LDL T algorithms is described in [1], where it is pointed out that special care must
be taken in the implementation to avoid large entries in L when A is ill conditioned.
Alternatively, one could use the factorization
G
6 P. C. HANSEN AND P. Y. YALAMOV

Table
The four post-processing rank-revealing steps for a symmetric semi-denite matrix.
Post-proc. Decomposition Symmetric matrix
R T
22 R 22
RRQR
22 T 22
22 L 21 L T
22 L 22
described in [29], where G is block triangular. If A is a symmetric indenite Toeplitz
matrix, then the currently most reliable approach to computing the LDL T factorization
seems to be via orthogonal transformation to a Cauchy matrix [21].
The reason why we need the post-processing step is that the initial factorization
may not reveal the numerical rank of A| there is no guarantee that small eigenvalues
of A manifest themselves in small diagonal elements of C or in small eigenvalues
of D. In particular, since
2 , we obtain
2showing that a small n may not be revealed in D when L is ill conditioned.
4.1. Algorithms for Semi-Denite Matrices. For symmetric semi-denite
matrices there is a simple relationship between the SVDs of A and C.
Theorem 4.1. The right singular vectors of P T AP are also the right singular
vectors of C, and
Proof. The result follows from inserting the SVD of C into P T
Hence, once we have computed the initial pivoted Cholesky factorization (4.1),
we can proceed by computing a rank-revealing decomposition of C, and this can be
done in several ways. Let E denotes the exchange matrix consisting of the columns
of the identity matrix in reverse order, and write P T AP as
Then we can compute a URV or RRQR decomposition of C, a ULV decomposition of
ECE, or an RRQR decomposition of (ECE) T , as shown in the left part of Table 4.1.
The approach using the URV decomposition of C was suggested in [2]. Table 4.1 also
shows the particular forms of the resulting symmetric matrix S, as derived from the
following relations:
R
The rst, third and fourth approaches lead to a symmetric matrix S that reveals the
numerical rank of A by having both an o-diagonal block S 12 and a bottom right
block S 22 with small norm. The second approach does not produce blocks S 12 and
S 22 with small norm; instead (since T 11 is well conditioned) this algorithm provides a
permutation P that is guaranteed to produce a well-conditioned leading
The remaining three algorithms yield approximate bases for the range and null
spaces of A, due to Theorem 3.1. It is well known that among the rank-revealing
decompositions, the ULV decomposition can be expected to provide the most accurate
bases for the right singular subspaces, in the form of the columns of VL ; see, e.g., [32]
and [15]). Therefore, the algorithm that computes the ULV decomposition of ECE is
to be preferred. We remark that the matrix UL in the ULV decomposition need not
be computed.
In terms of the blocks S 12 and S 22 , the ULV-based algorithm is the only algorithm
that guarantees small norms of both the o-diagonal block S
22 and the bottom
right block S
22 L 22 , because the norms of both L 12 and L 22 are guaranteed
to be small. From Theorem 4.1 and the denition of the ULV decomposition we have
and therefore kS 12 k 2 ' kS 22 k 2 '  k+1 .
For a sparse matrix the situation is dierent, because the UTV post-processors
may produce severe ll, while the RRQR post-processor produces only ll in the n k
rightmost columns of T . For example, if A is the upper bidiagonal matrix
in which B p is an upper bidiagonal p  p matrix of all ones, and e p is the pth column
of the identity matrix, then URV with threshold produces a full k  k
upper triangular R 11 , while RRQR with the same threshold produces a k  k upper
bidiagonal T 11 . Hence, for sparsity reasons, the UTV approaches may not be suited
for computing the VSV decomposition, depending on the sparsity pattern of A.
An alternative is to use the algorithm based on RRQR decomposition of the
transposed and permuted Cholesky factor (ECE) and we note that the
permutation matrix  is not needed. In terms of the matrix S, only the bottom
right submatrix of S is guaranteed to have a norm of the order  k+1 , because of the
relations kS
and kS 22 k
In practice the situation can be better, because the RRQR-algorithm|when
applied to the matrix E C T E |may produce an o-diagonal block T 12 whose norm
is smaller than what is guaranteed (namely, of the order  1=2
The reason is that the
initial Cholesky factor C often has a trailing (n
whose norm is close to  1=2
k+1 , which may produce a norm kS 12 k 2 close to  k+1 . From
the partitionings
22 En k En k C T
and the fact that the RRQR post-processor leaves column norms unchanged and may
permute the leading n k columns of E C T E to the back, we see that the norm of
the resulting o-diagonal block T 12 in the RRQR decomposition can be bounded by
Our numerical examples in x5 illustrate this.
However, we stress that in the RRQR approach we can only guarantee that kS 12 k 2
is of the order  1=2
, and this point is illustrated by the matrix
8 P. C. HANSEN AND P. Y. YALAMOV
1. Compute the eigenvalue decomposition
2. Write  as
jj 1=2 .
3. Compute an orthogonal W such that C
is lower triangular.
Fig. 4.1. Interim processor for symmetric indenite matrices.
where K is the \infamous" Kahan matrix [7, p. 105] that is left unchanged by QR
factorization with ordinary column pivoting, yet its numerical rank is
factorization with symmetric pivoting computes the Cholesky factor
and when we apply RRQR to E C T E we obtain an upper triangular matrix T in which
only the (n; n)-element is small, while kT 12
n .
4.2. Algorithms for Indenite Matrices. No matter which factorization is
used for an indenite matrix, such as (4.2) or (4.3), there is no simple relationship
between the singular values of A and the matrix factors. Hence the four \intrinsic"
decompositions from Table 4.1 do not apply here, and the diculty is to develop a
new factorization from which the numerical rank can be determined.
All rank-revealing algorithms currently in use maintain the triangular form of the
matrix in consideration, but when we apply the algorithms to the matrix L in the
LDL T factorization (4.2) we destroy the block diagonal form of D. We can avoid
this diculty by inserting an additional interim stage between the initial LDL T factorization
and the rank-revealing post-processor, in which the middle block-diagonal
matrix D is replaced by the signature
matrix
diag(1). At the same time, L is
replaced by the product of an orthogonal matrix and a triangular matrix. The interim
processor, which is summarized in Fig. 4.1, thus computes the factorization
where W is orthogonal and C is upper triangular
The interim processor is simple to implement and requires at most O(n 2 ) opera-
tions, because W and W are block diagonal matrices with the same block structure
as D. For each 1  1 block d ii in D the corresponding 1  1 blocks in W , jj 1=2 , and
W are equal to 1, jd ii j 1=2 , and 1, respectively. For each 2  2 block in D we compute
the eigenvalue decomposition
d ii d i;i+1
d i;i+1 d i+1;i+1
then the corresponding 22 block in W is W ii , and the associated 22 block in W is a
Givens rotation chosen such that C stays triangular. If A is sparse, then some ll may
be introduced in C by the interim processor, but since the Givens transformations are
applied to nonoverlapping 22 blocks, ll introduced in the treatment of a particular
block does not spread during the processing of the other blocks. The same type of
processor can also be applied to the
G
G T factorization (4.3) in order to
turn the block triangular matrix G into triangular form.
Future developments of rank-revealing algorithms for more general matrices than
the triangular ones may render the interim processor super
uous. It may also be
possible to compute the factorization (4.5) directly.
We shall now explore the possibilities for using triangular rank-revealing post-processors
similar to the ones for semi-denite matrices, but modied such that they
yield a decomposition of C in which the leftmost matrix U is hypernormal with respect
to the signature
matrices
and b
i.e., we require U
Hypernormal matrices
and the corresponding transformations are introduced in [8] in connection with up-
and downdating of symmetric indenite matrices. Here we use them to maintain the
triangular form of the matrix C.
The following theorem shows that a small singular value of A is guaranteed to be
revealed in the triangular matrix C.
Theorem 4.2. If n (C) denotes the smallest singular value of C in the interim
factorization (4.5), then
Proof. We have  1
, from which the result follows.
Unfortunately, there is no guarantee that n (C) does not underestimate  1=2
dramatically, neither does it ensure that the size of n is revealed in S. We illustrate
this with a small 5  5 numerical example from [1] where A is given by
with
and . The singular values of A are
such that A has full rank with respect to the threshold . The corresponding
matrix C has singular values
Thus,  5 (C) is not a good approximation of  1=2
5 , and if we base the rank decision on
5 (C) and the threshold  wrongly conclude that A is numerically
rank decient.
The conclusion is that for indenite matrices, a well conditioned C ensures that A
is well conditioned, but we cannot rely solely on C for determination of the numerical
rank of A. This rules out the use of RRQR factorization of C and ECE. The following
theorem (which expands on results in [24]) shows how to proceed instead.
Theorem 4.3. Let wn be an eigenvector of C
C corresponding to the eigenvalue
n that is smallest in absolute value, and let ~
wn be an approximation to wn .
Moreover, choose the orthogonal matrix b
V such that b
the last column of
the identity matrix, and partition the matrix
P. C. HANSEN AND P. Y. YALAMOV
such that S 11 is (n 1)  (n 1). Then
ks
wn wn k 2
and
wn wn
Proof. Let
C and consider rst the quantity
A ~
wn
s 22
Next, write ~
A ~
Au
wn
A u:
Combining these two results we obtain
s 22 n
A n I) u
and taking norms we get
ks
Both ks 12 k 2
2 and js 22 n j are lower bounds for the left-hand side, while kuk 2 is
bounded above by tan . Combining this with the bound k ^
obtain the two bounds in the theorem.
The above theorem shows that in order for n to reveal itself in S, we must
compute an approximate null vector of C
apply Givens rotations to this vector
to transform it into e n , and accumulate these rotations from the right into C. At the
same time, we should apply hypernormal rotations from the left in order to keep C
upper triangular. Theorem 4.3 ensures that if ~
wn is close enough to wn then ks 12 k 2
is small and s 22 approximates n . We note that hypernormal transformations can be
numerically unstable, and in our implementations we use the same stabilizations as
in the stabilized hyperbolic rotations [7, x3.3.4].
Once this step has been performed, we de
ate the problem and apply the same
technique to the (n 1)(n 1) submatrix S
are
the leading submatrices of the updated factors. This is precisely the algorithm from
[24]. When the process stops (because all the small singular values of A are revealed)
we have computed the URV-like decomposition
R such that U T
R
and the middle rank-revealing matrix is given by
R Tb
where
1 is k  k.
The condition estimator used in the URV-like post-processor must be modied,
compared to the standard URV algorithm, because we must now estimate the smallest

Table
Summary of approaches for symmetric indenite matrices. Note that the RRQR approaches
do not reveal the numerical rank, and that the ULV-like approach is impractical.
Post-proc. Decomposition Comments to decomposition
R
R
R 11
RRQR
singular value of the matrix C
C. In our implementation we use one step of inverse
iteration applied to C
C, with starting vector from the condition estimator of the
ordinary URV algorithm applied to C.
Finally consider a ULV-like approach applied to ECE. Again we must compute
an approximate null vector of C
C and transform it into the form e n by means of
an orthogonal transformation. This transformation is applied from the right to ECE,
and a hypernormal transformation from the left is then required to resotre the lower
triangular form of
To de
ate this factorization, note that the leading (n 1)  (n 1) block of
L is given by
L 11
and
This shows that we cannot merely work on the block L 11 ; also
the 1  (n 1) block ' T
21 is needed. Hence, after the de
ation step we must work with
trapezoidal matices instead of triangular matrices. This fact renders the ULV-like
approach impractical.
To summarize, for symmetric indenite matrices only the approach using the
URV-like post-processor leads to a practical algorithm for revealing the numerical
rank of A. Moreover, a well conditioned C signals a well conditioed A, but C cannot
reveal A's numerical rank. Our analysis is summarized in Table 4.2, and the URV-
based algorithm is summarized in Fig. 4.2 (following the presentations from [17]),
where  is the rank-decision tolerance for A.
4.3. Updating the VSV Decomposition. One of the advantages of the rank-
revealing VSV decomposition over the EVD and SVD is that it can be updated
eciently when A is modied by a rank-one change v v T . From the relation
we see that the updating of A amounts to updating the rank-revealing matrix S by
the rank-one matrix ww T with
v, i.e., S . This can be done in
while the EVD/SVD updating requires O(n 3 ) operations.
Consider rst the semi-denite case, and let M denote one of the triangular
12 P. C. HANSEN AND P. Y. YALAMOV
1. Let k n and compute an initial factorization P T
2. Apply the interim processor to compute P T
3. Condition estimation: let e  k estimate  k C(1: k; 1:
and let w k estimate the corresponding right singular vector.
4. If e  k >  1=2 then exit.
5. Revealment: determine an orthogonal Q k such that Q T
6. update C(1: k; 1:
7. update C(1: k; 1:
the hypernormal
is chosen such that the updated C is
8. De
ation: let k k 1.
9. Go to step 3.
Fig. 4.2. The URV-based VSV algorithm for symmetric indenite matrices.
matrices R, L, or T T from the algorithms in Table 4.1. Then
and we see that the VSV updating is identical to standard updating of a triangular
RRQR or UTV factor, which can be done stably and eciently by means of Givens
transformation as described in [5], [31] and [32].
Next we consider the indenite case (4.9), where the updating takes the form
showing that the VSV updating now involves hypernormal rotations. Hence, the up-dating
is computationally similar to UTV downdating, whose stable implementation
is discussed in [3] and [25]. Downdating the VSV decomposition will, in both cases,
also involve hypernormal rotations.
4.4. Computation of Truncated VSV Solutions. Here we brie
y consider
the computation of the truncated VSV solution which we dene as
consists of the rst k columns of VS . Both URV-based decompositions are
simple to use. For the ULV-based decomposition we have S
and we can safely neglect the term L T
whose norm is at most of the order  k+1 .
Finally, for the RRQR-based decomposition we can use the following theorem.
Theorem 4.4. If
T is the triangular QR factor of (T
then
Alternatively, if the columns of the matrix
form an orthonormal basis for the null space of (T
I W 1
(4.

Table
Numerical results for the rank-revealing VSV algorithms.
Post-processor
(semi-def.)
ULV mean 2:7
(semi-def.)
RRQR mean 1:5
(semi-def.)
URV-like mean
(indef.)
Proof. If
T is a QR factorization then S
T and S 1
T whici is (4.11). The same relation leads to S 1
denotes the pseudoinverse. In [6] used that
which, combined with the relation (I W W T immediately leads
to (4.12).
The rst relation (4.11) in Theorem 4.4 can be used when k  n, while the
second relation (4.12) is more useful when k  n. Note that W can be computed by
orthonormalization of the columns of the matrix
I
This approach is particularly useful for sparse matrices because we only introduce ll
when working with the \skinny" n  (n
5. Numerical Examples. The purpose of this section is to illustrate the theory
derived in the previous sections by means of some test problems. Although robust-
ness, eciency and
op counts are important practical issues, they are also tightly
connected to the particular implementation of the rank-revealing post-processor, and
not the subject of this paper.
All our experiments were done in Matlab, and we use the implementations of the
ULV, URV, and RRQR algorithms from the UTV Tools package [17]. The condition
estimation in all three implementations is the Cline-Conn-Van Loan (CCVL)
estimator [12]. The modied URV algorithm used for symmetric indenite matrices is
based on the URV algorithm from [17], augmented with stabilized hypernormal rotations
when needed, and with a condition estimator consisting of the CCVL algorithm
followed by one step of inverse iteration applied to the matrix C
C.
Numerical results for all the rank-revealing algorithms are shown in Table 5.1,
where we present mean and maximum values of the norms of various submatrices
associated with the VSV decompositions. In particular, X
or T 12 , and X 22 denotes either R 22 , L 22 , or T 22 . The results are computed on the
basis of randomly generated test matrices of size 64, 128, and 256 (100 matrices of
each size), each with n 4 eigenvalues geometrically distributed between 1 and 10 4 ,
14 P. C. HANSEN AND P. Y. YALAMOV

Table
Numerical results with improved singular vector estimates.
URV mean 1:8
(semi-def.)
URV-like mean 8:6
(indef.)
and the remaining four eigenvalues given by
the numerical rank with respect to the threshold
The test matrices were produced by generating random orthogonal matrices and
multiplying them to diagonal matrices with the desired eigenvalues. For the indenite
matrices the signs of the eigenvalues were chosen to alternate.

Table

5.1 illustrates the superiority of the ULV-based algorithm for semi-denite
matrices, for which the norm kS 12 k 2 of the o-diagonal block in S is always much
smaller than the norm kS 22 k 2 of the bottom right submatrix. This is due to the fact
that the ULV algorithm produces a lower triangular matrix L whose o-diagonal block
L 21 has a very small norm (and we emphasize that the size of this norm depends on
the condition estimator). The second best algorithm for semi-denite matrices is the
one based on the RRQR algorithm, for which kS 12 k 2 and kS 22 k 2 are of the same size.
Note that it is the latter algorithm which we recommend for sparse matrices. The
URV-based algorithm for semi-denite matrices produces results that are consistently
less satisfactory than the other two algorithms. All these results are consistent with
our theory.
For the indenite matrices, only the URV-like algorithm can be used, and the
results in Table 5.1 show that this algorithm also behaves as expected from the theory.
In order to judge the backward stability of this algorithm, which uses hypernormal
rotations, we also computed the backward error kA
k 2 for all three hundred
test problems. The largest residual norm was 1:9  10 11 , and the average is
We conclude that we loose a few digits of accuracy due to the use of the hypernormal
rotations.
It is well known that the norm of the o-diagonal block in the triangular URV
factor depends on the quality of the condition estimator | the better the singular
vector estimate, the smaller the norm. Hence, it is interesting to see how much the
norms of the o-diagonal blocks in R and S decrease if we improve the singular vector
estimates by means of one step of inverse iteration (at the expense of additional
ops). In the semi-denite case we now apply an inverse iteration step to
the CCVL estimate, and in the indenite case we use two steps of inverse iteration
applied to C
C instead of one. The results are shown in Table 5.2 for the same
matrices as in Table 5.1. As expected, the norms of the o-diagonal blocks are now
smaller, at the expense of more work. The average backward errors kA VS S V T
did not change in this experiment.
6. Conclusion. We have dened and analyzed a class of rank-revealing VSV
decompositions for symmetric matrices, and proposed algorithms for computing these
decomposition. For semi-denite matrices, the ULV-based algorithm is the method
of choice for dense matrices, while the RRQR-based algorithm is better suited for
sparse matrices because it preserves sparsity better. For indenite matrices, only the
URV-based algorithm is guaranteed to work.



--R

Accurate Symmetric Inde
A correlation-based subspace tracking algorithm
An algorithm and a stability theory for downdating the ULV decomposition







On rank-revealing QR factorizations
Generalizing the LINPACK condition estima- tor

Perturbation analysis for two-sided (or complete) orthogonal decompositions
Bounding the subspaces from rank revealing two-sided orthogonal decompositions

Matlab templates for rank- revealing UTV decompositions
Rank and null space calculations using matrix decomposition without column interchanges


Transformation techniques for Toeplitz and Toeplitz-plus-Hankel matrices II
Analysis of the Cholesky decomposition of a semi-de nite matrix
The rank revealing QR decomposition and SVD
A symmetric rank-revealing Toeplitz matrix decomposition

A Sparse URL Rather Than a URV Factorization
Sparse multifrontal rank revealing QR factorization


Cholesky factorization of semi-de nite Toeplitz matrices
An updating algorithm for subspace tracking
Updating a rank-revealing ULV decomposition
Matrix Algorithms Vol.
--TR
rank-revealing decompositions;hypernormal rotations;matrix approximation;symmetric matrices