More Accurate Bidiagonal Reduction for Computing the Singular Value Decomposition.
Bidiagonal reduction is the preliminary stage for the fastest stable algorithms for computing the singular value decomposition (SVD) now available. However, the best-known error bounds on bidiagonal reduction methods on any matrix are of the form  \[ A are orthogonal, $\varepsilon_M$ is machine precision, and f(m,n) is a modestly growing function of the dimensions of A.A preprocessing technique analyzed by Higham [Linear Algebra Appl., 309 (2000), pp. 153--174] uses orthogonal factorization with column pivoting to obtain the factorization  \[ A=Q \left( \begin{array}{c} C^T \\ 0 \end{array} \right) P^T, \]  where Q is orthogonal, C is lower triangular, and P is permutation matrix. Bidiagonal reduction is applied to the resulting matrix C.To do that reduction, a new Givens-based bidiagonalization algorithm is proposed that produces a bidiagonal matrix B that satisfies $C bounded componentwise and $\delta C$ satisfies a columnwise bound (based upon the growth of the lower right corner of C) with U and V orthogonal to nearly working precision. Once we have that reduction, there is a good menu of algorithms that obtain the singular values of the bidiagonal matrix B to relative accuracy, thus obtaining an SVD of C that can be much more accurate than that obtained from standard bidiagonal reduction procedures. The additional operations required over the standard bidiagonal reduction algorithm of Golub and Kahan [J. Soc. Indust. Appl. Math. Ser. B Numer. Anal., 2 (1965), pp. 205--224] are those for using Givens rotations instead of Householder transformations to compute the matrix V, and 2n3/3 flops to compute column norms.
Introduction
. We consider the problem of reducing an m \Theta n matrix A to
bidiagonal form. That is, we find orthogonal matrices U n\Thetan such
that
To denote B in (1.1) we use the shorthand
or use the MATLAB-like form
We will also use MATLAB notation for submatrices. Thus A(i: j; k: ') denotes the
submatrix of A consisting of rows i through j and columns k through '. Likewise,
denotes all of columns k through ' and A(i: all of rows i through
j.
For a matrix
denote the i th singular value of X , for ng. We also let X y be
the Moore-Penrose psuedoinverse of X and we let J(i; be a Givens rotation
through an angle ' ij applied to columns i and j.
Department of Computer Science and Engineering, The Pennsylvania State University, University
Park, PA 16802-6106, e-mail: barlow@cse.psu.edu, URL: http://trantor.cse.psu.edu/~barlow.
The research of Jesse L. Barlow was supported by the National Science Foundation under grants
no. CCR-9201612 and CCR-9424435. Part of this work was done while the author was visiting the
University of Manchester, Department of Mathematics, Manchester M13 9PL UK
The reduction (1.1) is usually done as a preliminary stage for computing the
singular value decomposition of A. There are now a number of very good algorithms
for computing the singular value decomposition of bidiagonal matrices. We know that
the "zero-shift" Q-R algorithm [13], bisection [4], and the dqds algorithm [15] can
compute all of the singular values of B to relative accuracy. We also know that it
is not reasonable to expect any algorithm to compute all of the singular values of a
matrix to relative accuracy unless that matrix has an acyclic graph [11] or is totally
sign compound [12].
Thus, it is not surprising that no algorithm can be expected to produce the
bidiagonal form of a general matrix to relative accuracy in fixed precision arithmetic.
The Jacobi algorithm is a more accurate method for finding the singular values of a
general matrix than any algorithm that requires bidiagonal reduction. Unfortunately,
the Jacobi algorithm is usually slower. For simplicity, assume that
reduction followed by the Q-R algorithm can produce that SVD in about 20n 3 flops.
One Jacobi sweep requires about 7n 3 flops. Thus, for Jacobi to be competitive, it
must converge in about three sweeps, and that rarely happens.
In this paper, we present a bidiagonal reduction method that will often preserve
more of the accuracy in the singular value decomposition.
The reduction is computed in two stages. In the first stage, using a Householder
factorization method of Cox and Higham [10], we reduce A to a lower triangular
n\Thetan . In floating point arithmetic with machine unit " M , the first stage
reduction satisfies
an );
where
Here f(n) is a modestly sized functions and ae A is a growth factor given in [10]. A
similar reduction is recommended by Demmel and Veseli'c [14] before using the Jacobi
method. Thus the difference in our algorithm is the reduction of C.
In the second stage, we apply a new bidiagonal reduction algorithm to C. This
algorithm produces a bidiagonal matrix B such that for some n \Theta n matrices ~
and some orthogonal matrices U and V , we have
where s is the smallest integer such that kC(:
M kCkF . The growth
ae V is bounded provided that the s \Theta s principal submatrix of the corresponding
Krylov matrix is nonsingular. If that submatrix is singular, the standard backward
error bounds from [17] apply. This is not as good a bound as the Jacobi method
achieves [14, 27, 28], but can be much better than that achieved by the standard algo-
rithm. Moreover, the algorithm can be implemented in less than 2 8
flops than the Lawson-Hanson-Chan SVD [26, pp.107-120],[9]. Using fast versions of
Givens rotations, that additional overhead may be reduced to about 8=27n 3
flops.
Our procedure for bidiagonal reduction has some important differences from the
Golub-Kahan Householder transformation based procedure [17]. They are
Givens transformations are used in the construction of the right orthogonal
matrix V . (Clearly, 2 \Theta 2 Householder transformations could also be used.)
ffl The matrix A is preprocessed by Householder factorization with maximal
column pivoting using a new procedure due to Cox and Higham [10].
ffl The first column of V is not e 1 , the first column of the identity matrix.
ffl The computation of the matrices U and V are interleaved in a different manner
to preserve accuracy in the small columns.
In the next section, we give our algorithm for producing the bidiagonal form. In
x4 we prove the bounds (1.2)-(1.3). In x5, we give some tests and a conclusion.
2. Reduction to triangular form. Before giving the bidiagonal reduction pro-
cedure, we preprocess A by reducing it to lower triangular form using a Householder
transformation based procedure due to Cox and Higham [10]. It is based upon the
row and column pivoting procedure of Powell and Reid [30], but uses a simpler form
of pivoting.
The procedure is as follows.
1. Reorder the rows of A so that
2. Using the maximal column pivoting algorithm of Businger and Golub [8],
factor A into
n\Thetan is
lower triangular.
This particular Householder factorization algorithm has very strong numerical
stability properties. If we let C be the computed lower triangular factor, then Cox
and Higham [10] (based on the analysis of a more complicated pivoting strategy by
Powell and Reid [30]) showed that for some orthgonal matrix U 0 ,
and
(2.
ae A is a growth factor bounded by p
2.
The column oriented error bound (2.2) holds for standard Householder factorization
[31]. The second rowwise error bound (2.3) can only be shown for algorithms that do
some kind of row and column permutations for stability.
Similar results for Givens based algorithms are given by Barlow [3] and Barlow and
Handy [5]. The Givens based algorithm in [3] could be substituted for the Cox-Higham
algorithm. Gulliksen [21] has given a new framework for orthogonal transformations
with this property. Cox and Higham demonstrate that Householder's original version
of Householder transformation must be used for these bounds to hold. The bound
does not hold if Parlett's [29] version of the Householder transformation is used.
The columns of the matrix C satisfy
In fact, any reduction of A satisfying the property (2.4) will be a suitable preprocessing
step and will lead to the results given here.
We now give algorithms for computing the bidiagonal reduction of C.
3. Bidiagonal Reduction Algorithms.
3.1. The Golub-Kahan Bidiagonal Reduction. The usual bidiagonal reduction
algorithm is that given by Golub and Kahan [17, pp.208-210,Theorem 1], see
also Golub and Reinsch [19, pp.404-405]. It is given below for the square matrix C.
Since C is lower triangular, it is exactly rank deficient only if it has zero diagonals.
Minor and obvious modifications to the bidiagonal reduction algorithms in this section
are necessary if C has zero diagonals.
Algorithm 3.1 (Golub-Kahan Bidiagonal Reduction).
1. Find an orthogonal transformation U 1 such that
2. for
(a) Find an orthogonal transformation V k such that
(b) Find an orthogonal transformation U k such that
3.
% The bidiagonal reduction of C is given by
Golub and Kahan [17] used Householder transformations to describe their algo-
rithm. In that case, it requires 4n 3 +O(n 2 ) flops if the matrices U and V are kept in
factored form. If U and V are accumulated, it requires 8n 3 +O(n 2 ) flops.
3.2. A Givens Based Bidiagonal Reduction Algorithm. Below is our new
algorithm for the bidiagonal reduction of C.
Algorithm 3.2 (New Procedure for Bidiagonal Reduction).
We now present a Givens-based bidiagonal reduction procedure for an n \Theta n matrix
C satisfying (2.4). In x4, we show that this new algorithm will achieve error bounds
of the form (1.5)-(1.8).
1. Determine the smallest integer s such that
22 is an empty matrix and we use Algorithm 3.1, complete
steps 2-4.
2. Compute the vector z given by
22
be the product of Givens rotations
such that
Compute
Let U 1 be an orthogonal transformation such that
3. for
(a) Let V k be the product of Givens rotations
that satisfies
(b) Find an orthogonal transformation U k such that
C(k: n; k: n) / U T
4.
% The bidiagonal reduction of C is given by
If we use standard Givens rotations, steps 2-4 of this algorithm require 10n 3
flops. The use of fast Givens rotations as described in Hammarling [22], Gentleman
[16], Bareiss [2], Barlow [3], or Anda and Park [1] would produce an algorithm
with approximately the same as the Golub-Kahan Householder-
based procedure.
Step one requires the minimum length solution of the least squares problem
11 is already upper triangular, a procedure in Lawson and Hanson [26, pp.77-
83] would allow us to compute the orthogonal factorization of the above matrix in
flops. The maximum complexity of this step is for
this step never costs more than 8

Table

3.1 summarizes the complexity of the two bidiagonal reduction algorithms.

Table
Complexity of Bidiagonal Reduction Algorithms
Compute 3.1 Algorithm 3.2(SG) Algorithm 3.2(FG)
Givens rotations
fast Givens rotations
4. Growth Factors and Error Bounds.
4.1. Bounding the columns of C. To show the advantages of Algorithm 3.2
for bidiagonal reduction, we consider the effects of orthogonal transformations from
the right used to form V in Algorithm 3.1.
Let U be the matrix from Algorithm 3.1, then consider
By orthogonal equivalence,
If we let
~
then
has the form
where
and F (k)
12 is zero except for the last row. Therefore, ~
in effect, zeros out lower
k \Theta (n \Gamma k) block of F .
The following lemma bounds the effect of a large class of orthogonal transformations
from the right.
Lemma 4.1. Let F n\Thetan be partitioned according to
where V 11 is nonsingular. Let
Then
~
22
where
~
Proof: Matching blocks in (4.1) and (4.2) yields
Using the fact that V 11 is nonsingular, block Gaussian elimination yields (4.3).
The result of Lemma 4.1 generalizes to the case given next.
Lemma 4.2. Let F 2 ! m\Thetan be partitioned according to
22 0
and V has the form in (4.1). Let G = FV be partitioned
Then
where
~
Proof: We note that
22 V (1)
If we partition V according to (4.1) then
First we note that V 11 is nonsingular if and only if V (i)
are nonsingular.
Evaluating ~
22 leads to
~
\GammaV (2)
Thus k ~
. Now we have that
22 F 23
where the (1,3) and (2,3) blocks of C are unaffected. We then apply Lemma 4.1 to
F (1)
F (1)
!/
to obtain (4.5).
The following corollary relates the growth of the (2; 3) block to growth in Gaussian
elimination.
Corollary 4.3. Let F , G and V be as in Lemma 4.2. Then
V is the growth factor for k steps of Gaussian elimination on V in the Euclidean
norm.
Proof: If we let
where L is unit lower triangular and R is upper triangular then
Using the fact that V we have that ~
22 solves
'' ~
~
22
22
and that ~
V 22 are just the result of performing k steps of Gaussian elimination
on V . Thus
~
Taking norms yields
22
V is the growth factor from k steps of Gaussian elimination on V .
Unfortunately, the bound on growth for Gaussian elimination for a given row
ordering is no better for orthogonal matrices than it is for all matrices. The following
result is proven by Barlow and Zha [7].
Proposition 4.4. Let n\Thetan be nonsingular and have the P-L-R factoriza-
tion
by Gaussian elimination with the row ordering P where L is lower triangular and R
is upper triangular. For each have the growth factor ae (k)
. Let
X have the factorization
where V is orthogonal and Y is upper triangular. Then Gaussian elimination with
partial on V obtains the P-L-R factorization
where
and for each k, ae (k)
X .
Suppose that C 2 ! n\Thetan is a lower triangular matrix satisfying
j. Note that row pivoting procedure used in the
previous section assures us that a small block C 22 may be isolated.
n\Thetan be given by
Then X has blocks satisfying
where
22 C 21 =j;
(4.
22 C 22
Note that 2. We can prove some results above the Krylov
matrix associated with X .
We now give the following lemma about Krylov matrices for X .
Lemma 4.5. Let n\Thetan be given by (4.7)-(4.10) and the matrix X 11 in (4.7)
be nonsingular. Let y be such that
s y
Assume that y is chosen so that kzk
Proof: The proof is just a verification of the formulae. First consider
have that
The first term of both rows satisfy the lemma, and second term may be bounded by
kp (1)
kp (1)
The induction step is simply,
where
For both of the reccurances (4.11) and (4.12), we can bound the norms by
where
Simple substitution shows that i k - k(1
We can now prove the following lemma.
Lemma 4.6. Let z n\Thetan and n\Thetan be as in Lemma 4.5 and let
n\Thetan be the Krylov matrix
Let K be partitioned in the form
s K 11 K 12
where s is as defined in (4.7). If K 11 is nonsingular, then
Proof: First we note that if two lower triangular matrices are given by
s I 0
s I 0
then
s I 0
above be given by
Then
s K 11 K 12
~
K 22
Here
and
s:
Thus,
leading to the bound
If we let
s K 11 K 12
K 22
If we reconstruct this factorization, then
Taking norms yields
which is the desired result.
It is well known that matrix V from bidiagonal reduction is the orthgonal factor
of the Krylov matrix K [18, pp.472-473]. That allows us to find the following bound
of the growth factor for bidiagonal reduction. A caveat to all of our results is that the
s \Theta s matrix K 11 be nonsingular.
Proposition 4.7. Let V be as in Lemma 4.2. Let C and s be as in (3.1), let X
be as in (4.7)-(4.10), and let K be as in Lemma 4.6. Assume that K 11 is nonsingular.
Then
where
Proof: From Proposition 4.4, since the matrix K has the orthogonal factorization
where V is the set of Lanczos vector for the bidiagonal reduction of C with the first
we have that
where
and
By a classical equivalence,
Since jX
22
Corollary 4.8. Let C satisfy (4.6) for some value of s. Let z be defined by
s y
for some y 6= 0. Let Algorithm 3.2 be applied to C and let U
be the orthogonal transformations generated by that algorithm. Define
~
where
Let C (k) be defined by
and let ae V be defined by (4.13)-(4.15). Then
Proof: Since ~
k is product of Givens rotations in the standard order, we directly
apply the results in Lemma 4.2.
First, let
Then
Since ~
taking advantage of rotations that
commute, we can write
~
where
Y
Y
Y
Y
Thus, -
have the structure in (4.1). Using the terminology of Lemma 4.2, we
have
(2)and that
22
Combining the results of Lemma 4.2 and Proposition 4.7 obtains
Note that the ability to factor ~
2 is a feature of a Givens rotation
based algorithm for bidiagonal reduction.
4.2. Error Bounds and Implications. The error bounds for this paper are
stated in two theorems. The first one is proven in x7, the second is a consequence of
the first.
Theorem 4.9. Let C 2 ! n\Thetan and let n\Thetan be
the bidiagonal matrix computed by Algorithm 3.2 in floating point arithmetic with
machine be the contents of C after k passes thorough
the main loop of Algorithm 3.2. Then there exist U; modestly growing
functions
where
and for
Theorem 4.10. Let C 2 ! n\Thetan satisfy (2.4), and let
n\Thetan be the bidiagonal matrix computed by Algorithm
3.2 in floating point arithmetic with machine unit " M . Let C r be the
contents of C after k passes thorough the main loop of Algorithm 3.2. Let s be the
smallest integer such that
and let ae V be defined by (4.13). Then there exist U; satisfying (4.16) and a
modestly growing function g 5 (\Delta) such that B and C satisfy (4.17), ffiB is as in Theorem
4.9 and
Proof: The proof of this result is simply a matter of bounding
for each value of k.
For k - s, we have that
For k ? s, orthogonal equivalence yields
2:
An application of (4.21) yields (4.22).
The standard error bounds on bidiagonal reduction are of the form (4.17) but
where ffiC only satisfies a bound of the form
For the Golub-Kahan procedure, this bound is probably as good as we can expect.
These lead to error bounds on the singular values of the form
This is satisfactory for large singular values, but of little use for small ones. The following
4 \Theta 4 example illustrates the difference between standard bidiagonal reduction
and the approach advocated here.
Example 4.1. Let A be the 4 \Theta 4 matrix
are small parameters. Using the MATLAB value
we chose That yields a matrix that has two singular values
clustered at 1, and two distinct singular values smaller than ffl.
To the digits displayed,
If we perform Algorithm 3.1 on A without the preprocessing in x2, we obtain
The use of the bisection routine in [4] obtain the singular values
The computed singular vector matrices are
The invariant subspaces for the double singular value at 1 and for singular values 3
and 4 are correct. However, the individual singular vectors for singular values 3 and
4 are wrong.
Algorithm 3.2 used after the reduction in x2 obtains
The computed singular values to the number of digits displayed are
Moreover, these corresponded to those computed by the Jacobi method to about 15
significant digits. The computed singular vector matrices were
\Gamma0:408248 \Gamma0:408248 \Gamma0:553113 0:600611
\Gamma0:408248 \Gamma0:408248 \Gamma0:243588 \Gamma0:779315
\Gamma0:408248 \Gamma0:408248 0:7967 0:178704C C A
Our version of the Jacobi method (coded by P.A. Yoon) obtains the slightly different
singular vector matrices
However, singular vectors for oe 3 and oe 4 are essentially the same, and subspace for
the clustered singular value near 1 is also essentially the same.
Quite recently, there have been a number of papers on the singular values and
vectors of matrices under structured perturbations. We give two results below that
are relevant to singular values for the perturbation given in Theorem 4.9.
The first is due to Kahan [13].
Lemma 4.11. Let
~
n\Thetan , and let - 1. If- ~
~
Lemma 4.11 says the forward errors in the matrix B make only a small relative
difference in the singular values.
The following result shows that the non-orthogonality of U and V in Theorem
4.9 causes only a small relative change in the singular values. This theorem is given
in [24, pp.423-424,problem 18].
Lemma 4.12. Let A 2 ! m\Thetan and let n\Thetap where p - n. Then for
Standard bounds on eigenvalue perturbation [18, Chapter 8] and taking square
roots leads to
The errors in B and non-orthogonality of the transformations U and V , make
only small relative changes in the singular values. The important effect is that of the
error ffiC in Theorems 4.9 and 4.10.
To characterize the effect of the error ffiC in (4.17), we use a generalization of
results that have been published by Barlow and Demmel [4] , Demmel and Veseli'c
[14], and Gu [20]. This version was proven by Barlow and Slapni-car [6] in a work in
preparation.
Consider the family of symmetric matrices
and the associated family of eigenproblems
(4.
Let (- i (i); x i (i)) be the ith eigenpair of (4.25) and define S(ffi) be the set of indices
given by
The set S(ffi) is the set of eigenvalues for which relative error bounds can be found.
The next theorem gives such a bound. Its proof follows that of Theorem 4 in [4, p.773].
Lemma 4.13. Let (- i (i); x i (i)) be the ith eigenpair of the Hermitian matrix in
(4.24). Let S(ffi) be defined by (4.26). If i 2 S(ffi) then
where
jx
jx
Proof: First assume that - i (i) is simple at the point i. Then from standard
eigenvalue perturbation theory for sufficiently small - we have
x
x
x
x
Thus we have
jx
If - i (i) is simple for all i 2 [0; ffi], then the bound (4.29) follows by integrating from
0 to ffi. In Kato[25, Theorem II.6.1,p.139], it is shown that the eigenvalues of H(i) in
are real analytic, even when they are multiple. Moreover, Kato [25, p.143] goes
on to point out that there are only a finite number of i where - i (i) is multiple, so
that - i (i) is continuous and piecewise analytic throughout the interval [0; ffi]. Thus
we can obtain (4.27)-(4.28) by integrating over each of the intervals in which - i (i) is
analytic.
The proof of the following proposition was given by Slapni-car [6].
Proposition 4.14. Let
have the singular
value decomposition
n\Thetan are orthogonal and
Then for each we have that
where
Proof: We prove this lemma is by noting that oe i
of H(i) where
A T (i) 0
Its associated eigenvector is
The matrix EH in (4.24) is given by
Thus from Lemma 4.13, oe i (i) satisfies (4.30) where
A u T
which is the formula for - A
Proposition 4.14 can be used to bound the error in the singular values caused by
both stages of the bidiagonal reduction. For the stage in x2, we have bounds of the
form (1.3)-(1.4). Thus for
From (2.3), the value of ae A is a growth factor bounded by
y is the Moore-Penrose pseudoinverse of A.
Now consider the error bound in Theorem 4.9. In context, we now consider u
v i to be the ith left and right singular vectors of C rather than A as in the preceding
paragraph. If we let
where d i is defined in (1.8). Then for
i is bounded by
A matrix satisfying (2.4) will often have small singular values that have modest values
of - C
. As pointed out by Demmel and Veseli'c [14], only an algorithm with columnwise
backward error bounds can take advantage of that fact.
The error bounds on computed singular vectors of C are also better for a bidi-
agonal reduction satisfying Theorem 4.9. For such bounds see the papers by Veseli'c
and Demmel [14] or Barlow and Slapni-car [6].
5. Numerical Tests. We performed two sets of numerical tests on the bidiag-
onal reduction algorithms. Our tests compared the singular values to those obtained
by the Jacobi method described by Demmel and Veseli'c [14].
The two test sets are as follows
Example 5.1 (Test Set 1). We use the set
was the Cholesky factor of the Hilbert matrix of dimension k. That is, R k
was the upper triangular matrix with positive diagonals that satisfied
where H k is a k \Theta k matrix whose (i; j) entry is
Example 5.2 (Test Set 2). The set
These are just the transposes of the matrices in Example 5.1. The reduction in x2
produces an upper triangular matrix from that. Thus the resulting bidiagonal matrices
tended to be more graded.
Example 5.3 (Test Set 3). We computed the SVD of L 35 from Example 5.2
using the MATLAB SVD routine to obtain
We then constructed 50 matrices of the form
where F was 50\Theta35 matrix generated by the MATLAB function randn (for generating
matrices with normally distributed entries).
Each matrix in both test sets were reduced to a matrix C using the algorithm in
x2. Three separate routines were used to find the SVD of C.
Plots for Bidiagonal Reduction for the SVD
for
Factor
Hilbert
Matrices
FromtoFig. 5.1. Relative Error Plots from Example 5.1
Algorithm J The Jacobi method described in [14].
ffl Algorihtm G The bidiagonalization method of Algorithm 3.2 followed by
the bisection routine of Demmel and Kahan [13].
ffl Algorihtm H The bidiagonalization method of Algorithm 3.1 followed by
the same bisection routine.
For each matrix, we calculated the two ratios
1-i-n
joe G
oe J
1-i-n
oe J
where oe J
i are the ith singular values as calculated by Algorithms J, G, and
H respectively. Thus we were trying to measure how well the SVDs calculated from
the two bidiagonal reduction algorithms agreed with that from the Jacobi algorithm.
The graphs of these values for the three test sets are given in Figures 5.1, 5.2, and 5.3.
A logarithmic plot of singular values of the matrix R 90 from Example 5.1 is given in

Figure

5.4. Clearly the singular values of L 90 in Example 5.2 span the same range.
All three Examples had singular values that spanned a wide range.
As can be seen from the figures, the maximum relative error for any singular value
from either test set for Algorithm G or H is about 10 \Gamma13 or approximately 10 3 times
Plots for Bidiagonal Reduction for the SVD
for
Factor
Hilbert
Matrices
FromtoFig. 5.2. Relative Error Plots from Example 5.2
machine precision in MATLAB. There seems to be no measurable difference between
Algorithms G and H in the quality of the singular values produced.
The error analysis produced in the x4 indicates that we should expect better
accuracy from Algorithm G, but does not explain why the singular values produced
from both algorithms are so accurate. We suspect that matrix C produced from
the reduction in x2 will almost always have good bidiagonal reduction from either
algorithm, but we know of no analysis to explain why this happens for Algorithm H.
6. Conclusion. We have presented a new bidiagonal reduction algorithm that
have four difference from the standard algorithm. Although it cost from 8
to 2 8
flops than the Golub-Kahan algorithm , our analysis shows that the
new reduction gives a better guarantee of accurate singular values.
Our numerical tests seemed to indicate that we performed the Cox-Higham
Householder factorization routine as a pre-processor, both the new algorithm and the
Golub-Kahan routine produced singular values that were even more accurate than
any known theory predicts. Thus we strongly recommend this preprocessing step
whenever it is feasible. The matrix C resulting from this reduction is usually highly
graded and we suspect that there is still much to understand about the behavior of
bidiagonal redcution on graded matrices.
The extra cost of the new bidiagonal reduction method is far less than that of
any implementation of the Jacobi method to date. It gives a reasonable guarantee of
relative accuracy singular values larger than " 3=2
M and tests confirm that behaves well
Plots for Bidiagonal Reduction for the SVD
fortrials
for
Very
Small
Singular
values
Fig. 5.3. Relative Error Plots from Example 5.3
for singular values smaller than " 3=2
M .



--R

Fast plane rotations with dynamic scaling.
Numerical solution of the weighted linear least squares problem by G- transformations
Stability analysis of the G-algorithm and a note on its application to sparse least squares problems
Computing accurate eigensystems of scaled diagonally dominant matrices.
The direct solution of weighted and equality constrained least squares problems.
Optimal perturbation bounds for the Hermitian eigenvalue prob- lem
Growth factors in Gaussian elimination
Linear least squares solutions by Householder transformations.
An improved algorithm for computing the singular value decomposition.
Stability of Householder QR factorization for weighted least squares problems.
On computing accurate singular values and eigenvalues of matrices with acyclic graphs.
Computing the Singular Value Decomposition with High Relative Accuracy
Accurate singular values of bidiagonal matrices.
Jacobi's method is more accurate than QR.
Accurate singular values and differential qd algorithms.
Least squares computations by Givens rotations without square roots.
Calculating the singular values and pseudoinverse of a matrix.
Matrix Computations
Singular value decomposition and least squares solutions.
Studies in Numerical Linear Algebra.
Backward error analysis for the constrained and weighted linear least squares problem when using the weighted QR factorization.
A note on the modifications to the Givens plane rotation.
Accuracy and Stability of Numerical Algorithms.
Matrix Analysis.
A Short Introduction to Perturbation Theory for Linear Operators.
Solving Least Squares Problems.
Fast accurate eigensystem computation by Jacobi methods.
Fast accurate eigenvalue methods for graded positive definite matrices Numer.
Analysis of algorithms for reflectors in bisectors.
On applying Householder's method to linear least squares prob- lems
The Algebraic Eigenvalue Problem.
--TR
singular values;bidiagonal form;orthogonal reduction;accuracy