行政院國家科學委員會專題研究計畫 成果報告
應用在正交矩陣上之二重移位 QR 算則收斂之探討
計畫類別： 個別型計畫
計畫編號： NSC93-2115-M-004-002-
執行期間： 93 年 08 月 01 日至 94 年 07 月 31 日
執行單位： 國立政治大學應用數學學系
計畫主持人： 王太林
報告類型： 精簡報告
處理方式： 本計畫可公開查詢
中 華 民 國 94 年 12 月 23 日
ON CONVERGENCE OF THE DOUBLE-SHIFT QR ALGORITHM FOR REAL ORTHOGONAL HESSENBERG MATRICES
TAI-LIN WANGy
Abstract. In applying the double-shift QR algorithm to compute the eigenvalues of real or-thogonal Hessenberg matrices, a unimodular shift strategy with a distinct monotonicity property and a cubic asymptotic rate of convergence is devised and tested, which is more e¢ cient than the conventional Francis shift in producing swift convergence and, in extreme circumstances, signi…cantly accelerates the de‡ation process, as numerical experiments demonstrate.
Key words. QR algorithm, shift strategy, convergence, orthogonal Hessenberg matrices AMS sub ject classi…cation. 65F15, 15A18
1. Introduction. In applying the double-shift QR algorithm [7][10] to compute the eigenvalues of real orthogonal Hessenberg matrices, we seek to customize the conventional Francis shift strategy to make the QR algorithm perform more e¢ ciently. Speci…cally, a unimodular shift is proposed and analyzed, with which the QR iteration has more rapid numerical convergence and still enjoys a unique monotonicity property and a cubic asymptotic rate; the distinctive advantage of using this customized shift is that, for an iterating matrix with elements under extreme conditions, the number of iterations required for de‡ation can be signi…cantly reduced. Furthermore, since the explicit invariant form under the double-shift QR can be fully determined in this case, preventive measures are devised in advance to avoid the unusual situations of slow- and no-convergence through the inclusion of an auxiliary shift. For orthogonal matrices, improvements in practical convergence of QR are demonstrated through numerical testing.
2. Notation and Background. In this paper, A will represent a real (upper) Hessenberg matrix of order n, with positive elements k := ek+1Aek, 1 k < n, on
the subdiagonal and zero entries below it, where ek is the k th column of the identity
matrix, and the superscript denotes the transposition of a given vector or matrix. The Euclidean norm jj jj is used to measure the length of a row or column vector. In case A is orthogonal, we use the letter U to respresent the matrix. The same structure and similar expressions are also extended to bA and A(k), which are de…ned later. A Hessenberg matrix is said to be unreduced if none of its subdiagonal elements is zero, and with no loss of generality in convergence analysis we may assume that the Hessenberg matrices considered in this paper have only positive subdiagonal elements, and hence are unreduced.
2.1. Double-shift QR iteration. Given a real Hessenberg matrix A, consider the orthogonal-triangular factorization of
(A 1I)(A 2I) 2(A) =: QR; (2.1) where Q is orthogonal, R is upper triangular with nonnegative diagonal elements, and
1, 2 are the eigenvalues of the trailing principal 2-by-2 submatrix of A. Observe
Research supported by the National Science Council of the Republic of China under grant 93-2115-M-004-002
yDepartment of Mathematical Sciences, National Chengchi University, Taipei, Taiwan. E-mail:
wang@math.nccu.edu.tw
that 2is a quadratic polynomial with real coe¢ cients, though 1, 2may appear as a complex conjugate pair. With Q we de…ne bA, the double QR transform of A, by
b
A := Q AQ:
It can be shown [8][13] that bA is a unique, real Hessenberg matrix with positive subdiagonal elements if 1, 2 are not eigenvalues of A, and it is easy to check that
( bA 1I)( bA 2I) 2( bA) = RQ: (2.2) The real Hessenberg QR algorithm (HQR) [10] iterates the transformation A ! bA (usually called a double QR step), and a sequence fA(k)g of unreduced, orthogonally similar Hessenberg matrices is produced if the (k)2 (A(k)) are nonsingular. It was
shown [4][14] that, if A is normal, then
bn 2bn 1 n 2 n 1 (2.3) for each double step; hence in the QR iteration, (k)n 2 (k)n 1form a bounded decreasing sequence for which the limit is essentially 0, if A(k)is not invariant under double QR; the asymptotic rate of (k)n 2 (k)n 1 ! 0 is very fast (see [17] or §3.4 for an analysis); consequently, as (k)n 2 (k)n 1becomes negligibly small, one or two (approximate) eigen-values of A are obtained after de‡ating A(k). The procedure continues on the de‡ated
submatrix until all the eigenvalues are found. See [10] for more details.
In this paper the HQR is employed on an orthogonal matrix U , and the speci…c structure of U is exploited to make the iterative process converge more e¢ ciently by customizing the conventional shift polynomial 2 while maintaining the merits, namely, with a monotonicity property similar to (2.3) and a cubic rate of convergence. We begin with a parametric representation of U .
2.2. Parametrization of the orthogonal matrix. It is known that each or-thogonal Hessenberg matrix A =: U of order n with positive subdiagonal elements f jgn 1j=1 can be uniquely expressed in a parametric form [1][9]
U = 2 6 6 6 6 6 6 4 0 1 0 1 2 0 1 2 n 2 n 1 0 1 2 n 1 n 1 1 2 1 2 n 2 n 1 1 2 n 1 n 2 . .. ... ... . .. n 2 n 1 n 2 n 1 n n 1 n 1 n 3 7 7 7 7 7 7 5 ; (2.4) where 1; 2; : : : ; n are usually called the Schur parameters of U , 0:= 1, j nj = 1,
and
j jj2+ 2j= 1, j 2 R, j > 0, 1 j < n. (2.5)
Observe from (2.4) that these parameters can be determined successively from the top row and the subdiagonal of U [9][15]:
1= e1U e1;
j = e1U ej= 1 2 j 1; 2 j n:
For unitary eigenvalue problems, this parametrization plays a fundamental role in simplifying intricate relations and deriving useful properties, because of the conve-nience of using n parameters to fully represent a structured matrix with more than n2=2 elements; see, for example, a recent paper by Bohnhorst, Bunse-Gerstner, and Faßbender [3], and the numerous references cited therein.
To …nd the eigenvalues of U , there is no loss of generality in assuming that n is even and n= 1; (2.6)
hence the eigenvalues are all in complex conjugate pairs, which further implies that, in the double QR iteration,
(k)
2j 1; j = 1; 2; 3; : : : ; n=2; are bounded away from 0: (2.7)
This is because existence of the real eigenvalues (+1 and/or 1) of U can be detected beforehand by inspecting the parity of n and the sign of n (for more details, see [1]
or [14]), and removal of these eigenvalues from U can be carried out through QR steps with shifts 1.
3. Convergence of the QR Iteration. In applying the double-shift QR to an orthogonal matrix U with the form (2.4), there are circumstances in which the use of the Francis shifts 1; 2 in (2.1) may not be appropriate; the shifts could be
fairly small or even zero in extreme cases, since they are the eigenvalues of a 2-by-2 submatrix of U and hence with magnitude always less than unity (because n 2> 0) [16]. If 1 = 2 = 0 (corresponding to n 2 = n 1 = 1), then bU = U and the sequence U(k) is invariant under double QR because
2(U ) is orthogonal. If 1; 2
are very small, then 2(U ) is still close to being orthogonal and bU would not be quite di¤erent from U , resulting in very slow decreases of (k)n 2 (k)n 1 in the early stages of the QR iteration; hence it may take numerous steps before (k)n 2 (k)n 1 ultimately gets into the asymptotic regime, where the convergence is cubic and swift. (See the data presented in Experiment 2 of Section 4.) Besides, 1; 2 could be real numbers with di¤erent magnitude, which apparently are not close to any conjugate pair of eigenvalues of U .
3.1. Construction of a unimodular-shift strategy. Based upon the fact that the eigenvalues of U are all situated on the unit circle and by (2.6) all appear in complex conjugate pairs, we should choose instead a unimodular conjugate pair f ; g as shifts, and it is only natural (by anticipating that n 2 0 and n 2 1
in (2.4)) to take the eigenvalues of the 2-by-2 orthogonal matrix
n 1 n 1 n 1 n 1 ;
which are n 1 i n 1, as and . This matrix has the characteristic polynomial 2( ) = 2+ 2 n 1 + 1 (3.1)
which depends only on n 1 and is used in place of 2 in (2.1) when A =: U is
orthogonal. For brevity, we shall call 2the unimodular shift (polynomial).
3.2. Monotonicity property of 2. Under the assumption that A is normal, a basic property of the QR transformation A ! bA : p(A) = QR; bA = Q AQ with any shift polynomial p is
jjenp( bA)jj jjenp(A)jj; (3.2)
which states that the last row of p(A) is decreasing after each QR step. To see this, observe that
jjenp( bA)jj = jjenRQjj = jjenRjj
jjRenjj = jjQRenjj = jjp(A)enjj (3.3)
= jjenp(A)jj;
in which properties of the unitary invariance of jj jj, triangular structure of R, and normality of p(A) have been used.
Applying the parametric form (2.4) of U and relation (2.5) to 2de…ned by (3.1), we obtain en 2(U ) = en 2 n 2 n 1+ en 1(1 n 2) n 1 n 1 (3.4) +en(1 n 2) 2n 1 and jjen 2(U )jj = p 2(1 n 2) n 1: (3.5)
For each double QR step U ! bU , it is easy to check that, ifbn 16= n 1, then
jjenb2( bU )jj < jjen 2( bU )jj; (3.6)
where b2( ) = 2+ 2bn 1 + 1. Summing up these calculations, we have
p
2(1 bn 2) bn 1
= jjenb2( bU )jj; from (3.5), jjen 2( bU )jj; by (3.6),
jjen 2(U )jj; by (3.2) with p = 2and A = U ,
= p2(1 n 2) n 1; from (3.5),
(3.7)
for each double step, which implies that q
2(1 (k)n 2) (k)n 1is monotonically decreas-ing in the QR iteration. If the limit of this sequence is zero, then clearly
(k) n 2=
q
1 (k)2n 2! 0
is implied, because (k)n 1is bounded below from 0 by (2.7), and a pair of complex conjugate eigenvalues is obtained after de‡ation. This is the ordinary case and it happens almost all the time. In the following subsection we consider the (unstable) situation in which the limit of the declining sequence is not zero.
3.3. Asymptotic analysis in the extreme case. Now suppose q
2(1 (k)n 2) (k)n 1! > 0 (3.8) as k ! 1; then from the above analysis outlined in (3.7) and (3.3) we know
jjenR(k)jj jjR(k)enjj !
and hence
R(k)en ! en;
which implies, through the QR factorization of
(k)
2 (U(k)) =: Q(k)R(k);
that asymptotically, the last column of (k)2 (U(k)) is orthogonal to all its preceding
n 1 columns and has a length of , or equivalently,
(k) 2 (U(k))
(k)
2 (U(k))en! 2en: (3.9)
This indicates that and en are approached by a singular value and vector of (k)
2 (U(k)). We seek to extract from (3.9) all the basic asymptotic relations among
the entries of U(k) in terms of the parameters (k)
j . The inner product of ej, j < n,
with (3.9) gives ej (k)2 (U(k)) (k)
2 (U(k))en ! 0 which, after switching the order of (k)
2 (U(k)) and (k)
2 (U(k)), is equivalent to
ej (k)2 (U(k)) (k)2 (U(k)) en ! 0 for j < n: (3.10)
The trick used here is that, due to the (upper) Hessenberg structure of U(k), the inner
product of ej (k)2 (U(k)) and en (k)2 (U(k)) in (3.10) is much simpler to compute, when j is close to n. After somewhat lengthy but elementary calculations using (2.4), (3.1) to compute ej (k)2 (U(k)) and applying (3.4) to e
n (k)
2 (U(k)), we reach the following
conclusions from (3.10):
for j = n 1, after rearrangements and factorizations,
(k) n 3 1 + (k)n 2 3 (k)n 2 ! (k) n 1! 0 (3.11)
because (1 (k)n 2) (k)n 1is bounded away from 0, by (3.8); for j = n 2, after collecting terms and simplifying with (3.11),
(k)
n 4! 1 and hence (k)
n 4! 0; (3.12)
because (k)2n 3 (k)n 2 (k)n 1is bounded away from 0, by (2.7) and (3.8).
From the asymptotic condition (3.12), it follows that a 4-by-4 submatrix is ultimately decoupled from U(k) in the lower right corner with
E4U(k)E4! U4 as k ! 1; (3.13)
where
E4:= [en 3; en 2; en 1; en];
and U4 is a 4-by-4 orthogonal matrix with Schur parameters fajg4j=1 satisfying the
conditions
a1
1 + a2
3 a2
= a3; (3.14)
with 0 jajj < 1 for 1 j 3, and a4= 1, because of (3.11) and (2.6). And there
is no need to pursue further information from (3.10) for j < n 2 because of the decoupling.
Some comments are given on the asymptotic behavior of (3.13) and the limit form U4 de…ned by (3.14):
1. Through the derivation of (3.11) and (3.12) from (3.8), it is obvious that the invariance ofp2(1 n 2) n 1under double QR implies that of
n 3
1 + n 2
3 n 2
= n 1 and n 4= 1 ( n 4= 0);
hence no unreduced U with n > 4 can be invariant because n 4 > 0. If n = 4, then clearly U4 is invariant, since 2(U4) gives a multiple of an
orthogonal matrix with 0 < =p2(1 a2)(1 a23) < 2, and this is the only
invariant form with 2. Also note that given the real numbers (0 j j < 1), (0 < < 2), and the quadratic polynomial
(z) = z2+ 2 z + 1;
there exist exactly 4 distinct numbers fzjg4j=1, on the unit circle and in two
conjugate pairs, such that j (zj)j = for 1 j 4; more precisely, if x1and
x2are the real parts of these two conjugate pairs, then 1
2(x1+ x2) =
2jxj+ j = ; j = 1; 2:
2. The highly unusual situation of (3.13) might exist only if four of the eigen-values j of U (in two conjugate pairs) satisfy the asymptotic condition
j (k)2 ( j)j ! ;
a direct consequence of (3.9), which is equivalent to jjen
(k)
2 (U(k))jj ! : (3.15)
This extreme case is unstable [4][2], in the sense that any small perturbations (e.g., rounding errors in numerical computation) in the limit space of (3.15) will translate into perturbations of (k)n 1 in the shift (k)2 which, due to the strict property of (3.6), may well cause
q
2(1 (k)n 2) (k)n 1 to drive past the value and eventually down to 0. So we would not observe this kind of decoupling of U(k)to occur in practical problems.
3. However, for an iterate U(k) with
E4U(k)E4 U4; (3.16)
it usually takes numerous steps in the QR iteration for q
2(1 (k)n 2) (k)n 1& 0
before the asymptotic regime is established, though convergence is ultimately cubic (which will be shown in the next subsection). Condition (3.16) therefore provides the very speci…c circumstances under which 2 does not function e¤ectively; fortunately, the situation of (3.16) is predictable and measurable through (3.11) and (3.12) before each QR step U(k) ! U(k+1) is carried out. To hasten convergence, an auxiliary shift can be used to break up the closeness in (3.16) and give the iteration a new starting U . See the formula devised later in (3.32). This special treatment is implemented and tested in Section 4.
3.4. Rate of convergence. In this subsection we …rst present, for normal ma-trices, a simple and illustrative analysis on the asymptotic rate of (k)n 2 (k)n 1 ! 0 in the double-step QR iteration with shift 2, using an elementary technique applied to the single-step QR by Wilkinson [18]. We then demonstrate that, with the unimod-ular shift 2 described by (3.1),
(k)
n 2 ! 0 in the QR iteration U(k) still enjoys a
cubic asymptotic rate, the same as that with the Francis shift 2. Note that a general
treatment for a generalized problem (including the convergence of multishift QR and LR) has been given by Watkins and Elsner [17] using subspace iteration techniques. Let
A =: A11 A12 A21 A22
; (3.17)
where A222 R2 2, and let
(A 1I)(A 2I) =: B11 B12 B21 B22 = Q11 Q12 Q21 Q22 R11 R12 O R22 := QR; (3.18)
where the partitionings are conformal with that of A. Premultiplying this equation by Q , we have Q11 Q21 Q12 Q22 B11 B12 B21 B22 = R11 R12 O R22
and, by equating respectively the (2; 1) and (2; 2) block entries on each side,
Q12B11+ Q22B21= O (3.19)
Q12B12+ Q22B22= R22 (3.20)
where, from (3.17) and (3.18),
B11= (A11 1I)(A11 2I) + A12A212 R(n 2) (n 2) (3.21)
B22= (A22 1I)(A22 2I) + A21A122 R2 2: (3.22)
In the following it is illustrated that, if A is normal Hessenberg, 1; 2are the eigenval-ues of A22, and n 2= " for some arbitrarily small number " > 0, then bn 2bn 1=
O("3) after one double-QR step. We use the Frobenius norm jjCjjF :=rP i;jjc
ijj2 to
measure the magnitude of a matrix C = [cij]. Since A is normal and (unreduced)
Hessenberg, we know
jjA12jjF= jjA21jjF = n 2= " (3.23)
and, in (3.22) with 2(A22) = (A22 1I)(A22 2I) = O (by the Cayley-Hamilton
theorem),
jjB22jjF "2: (3.24)
From the partitions de…ned by (3.17) and (3.18) we have B21= A21(A11 2I) + (A22 1I)A21;
and hence we know
jjB12jjF= jjB21jjF= O(") (3.25)
because of (3.23) and the fact that jjAjjF is always bounded. Furthermore, since all
the eigenvalues of A are simple, we see from (3.21) that, if " is su¢ ciently small and (3.23), (3.24) hold, then B11 is “bounded away” from being singular (see Lemma 6.4
and Theorem 6.5 of [17]) and
jjB111jjF (3.26)
for some …nite number . Next, let us examine the partitioned block of Q. From (3.19) we have
Q12= Q22B21B111
and therefore
jjQ12jjF= jjQ21jjF= O(") (3.27)
because of (3.25), (3.26) and the fact that Q is orthogonal and
jjQ22jjF jjQjjF=pn: (3.28)
Hence from (3.20) we infer that
jjR22jjF= O("2) (3.29)
by (3.24), (3.25), (3.27), and (3.28). Now let the partition of b A =: " b A11 Ab12 b A21 Ab22 #
be conformal with that of A and equate the (2; 1) block entry on each side of (2.2): b
A21( bA11 2I) + ( bA22 1I) bA21= R22Q212 R2 (n 2):
Further equating the (2; n 2) entry on each side of this matrix equation (observing that bA21= bn 2e2en 2 and that e2Ab22e1= bn 1) gives
bn 2bn 1= e2R22Q21en 2 jjR22Q21jjF= O("3);
because of (3.27) and (3.29). If (k)n 1are bounded away from zero (e.g., the eigenvalues sought are in complex conjugate pairs, as they are in the orthogonal case), then
jj bA21jjF= bn 2= O("3)
and this demonstrates that convergence of the double-shift QR is of cubic order. In reaching this conclusion, note that condition (3.29) is crucial, which in turn (cf. (3.20)) hinges on condition (3.24) that jjB22jjF is of order "2, which is derived from (3.22),
where we take advantage of the double shift 2 (i.e., 2(A22) = O) and normality of
A (cf. (3.23)). Note also that in the orthogonal case, cubic convergence of the QR process is still preserved if instead the unimodular shift 2 is used. This is because in (3.22) (with the A’s replaced by the U ’s)
B22= 2(U22) + U21U12 = 2(U22) 2(U22) + U21U12 = (1 n 2)( n 1U22+ I2) + U21U12; and we have jjB22jjF 2 p 2(1 n 2) + jjU21jj2F ( p 2 + 1)"2; since j n 1j 1, jjU22jjF jjI2jjF=p2, jjU21jjF = n 2= ", and
1 n 2= 1
p
1 "2 1
2"
2
for su¢ ciently small "1. Therefore, jjB22jjF is still of order "2 as before, and hence
cubic convergence of the iteration.
3.5. Su¢ cient conditions for convergence. We summarize our analyses in the previous subsections and conclude with a theorem on the convergence of double-step QR for orthogonal matrices; a more e¢ cient shift strategy based on the assump-tions of this theorem for such speci…cally structured matrices is proposed.
Theorem 3.1. Let the double-step QR algorithm with the unimodular shift poly-nomial 2 be applied to an unreduced orthogonal Hessenberg matrix U of order n 4 expressed in the Schur parametric form (2.4) under the assumption (2.6). If there is a number > 0 such that either
(j) n 3 1 + (j)n 2 3 (j)n 2 ! (j) n 1 (3.30) or (j) n 4 for n > 4 (3.31)
1Note that eventually
n 2 = j n 2j =
q 1 2
n 2 ! 1 if n 2 = " ! 0, since n 2 is
becoming the last Schur parameter of the de‡ated U of order n 2, which is +1 by (2.6). 9
holds for some subsequence with index j in the QR iteration U(k), then (k)
n 2! 0 as k ! 1
and two eigenvalues of U in a complex conjugate pair are being approximated; the rate of convergence is cubic.
Proof. If (k)n 2 9 0 as k ! 1, then from the previous analysis, (3.8) must be true, which implies both (3.11) and (3.12), a contradiction to either (3.30) or (3.31), whichever is supposed to hold in the QR iteration. That the rate of convergence is cubic has been proved in the last subsection.
Note that in the forecast of slow- or no-convergence (cf. (3.16) and (3.13)), condition (3.11) is of primary concern to us and more informative than (3.12), as all the subdiagonal elements of U(k)are inherently becoming smaller in the iterative process. In practice we may choose a “suitable” value for and examine, before each QR step, if condition (3.30) (with j := k) in the above theorem is satis…ed; and the auxiliary shift
e2( ) = 2+ 2 + 1 (3.32) is employed if it is not, to give the iteration a di¤erent starting U . (Another choice is e2( ) = 2 2 + 1.) The aim of using such an ad hoc shift is mainly to accelerate numerical convergence, when condition (3.16) is close enough and 2may not operate e¢ ciently.
It is clear that this ad hoc shift e2 can be viewed as the limit case of 2 with n 1 ! 1. Note also that e2 cannot be ine¤ective under condition (3.16) without
the parameters (k)n 1, (k)n 2, and (k)n 3 all being very close to unity (cf. (3.11)), in other words, the subdiagonal elements (k)n 1, (k)n 2, and (k)n 3 of U(k) all becoming vanishingly small in addition to (k)n 4 0 (cf. (3.12)), an extreme situation in which two conjugate pairs of eigenvalues are clustered near the point 1, and in this case, e2 (with double zeros at 1) should indeed be considered as a most desired shift polynomial to deal with these clustered eigenvalues.
4. Numerical Experiments. We demonstrate in this section that, for orthog-onal matrices, shift 2 provides more e¢ cient numerical convergence than shift 2 does, especially before the establishment of the asymptotic regime, though both were shown in §3.5 to have cubic asymptotic rates; the latter is incorporated with double QR in software packages. The testing was focused on determining the numbers of iterations required for (k)j 2 (k)j 1, j = n; n 1; : : : ; 3, to become negligible [10]. The computations were done in double-precision Fortran on an IBM compatible PC-80486 with unit roundo¤ " 10 19.
The following contractions are used to describe how the shift scheme in each case
is devised or modi…ed in carrying out the double QR iteration: HQR the subroutine in EISPACK [12] with shift 2 and
an exceptional shift [10]2introduced at steps k = 11
and k = 21 to combat the extreme situations;
HQR0 the HQR with shift 2but skipping the exceptional shift; HQRu the HQR with shift 2 instead and a check on condition
(3.30) in Theorem 3.1 with j := k and = 10 12 before using the auxiliary shift e2.
We also use the data from HQR0 to accentuate, for orthogonal matrices in particular, the signi…cance of the incorporated exceptional shift.
Given an orthogonal Hessenberg matrix U of size n (with n= 1 and n even), let
itmax= the maximum number of double-step QR iterations required to get one conjugate pair of eigenvalues :
To test the e¤ectiveness of the various shift strategies, we measure and list in the following tables the numerical averages of itmax over 10; 000 sample matrices of size n, each of which is constructed by a set of “randomly selected”real Schur parameters f jgnj=1with certain prescribed constraints as stated in each experiment.
Experiment 1. Randomly selected f jgn 1j=1 with n = 1.
itmax HQR0 HQR HQRu n = 4 5.06 5.02 4.11 n = 10 5.88 5.77 5.16 n = 20 6.48 6.30 5.81 n = 30 6.70 6.61 6.18
Observe the slight di¤erences in numbers between the columns HQR0 and HQR, which indicate that, even for randomly selected orthogonal matrices, the exceptional shift does have a role to play; in other words, there is still a small percentage (around 1-3%) of matrices for which the itmax exceeds at least 10 iterations. On average, HQRu is only about 7-18% more e¢ cient (in terms of itmax) than HQR for matrices of size 30, as the data listed in Experiment 1 show. However, under “speci…c” circumstances, the di¤erences among the three are very substantial, and HQRu works signi…cantly better than HQR, as Experiments 2 and 3 demonstrate.
Experiment 2. Randomly selected f jgn 1j=1 with j n 2j 10 7, j n 1j 10 7,
and n= 1. itmax HQR0 HQR HQRu n = 4 23.6 15.4 5.44 n = 10 22.6 16.1 5.67 n = 20 22.8 16.3 6.10 n = 30 22.9 16.4 6.34
2Later versions of both LAPACK and Matlab have made further modi…cations over the
excep-tional shift in the double-step QR [6][5][11]. 11
Sluggish convergence shown in the columns HQR0 and HQR of this table was pre-dicted and explained in the opening paragraph of Section 3.
Moreover, the conventional shift 2 also becomes ine¤ective for those U ’s with parameters under the conditions
j n 3 n 2 n 1j " and j n 4 1j ";
where " is small [14]; orthogonal matrices of this special type that cause HQR to fail are not hard to come by, even with the exceptional shift, as the next table indicates. Those “extreme matrices” for which convergence (to a pair of eigenvalues) does not occur within the iteration limit 30n set in HQR are not counted in the itmax averages. Instead, the number of such matrices (out of 10; 000 for each n) is placed in parentheses after the (itmax) average.
Experiment 3. Randomly selected f jgn 5j=1 [ f n 3; n 2g with n 4 =
q 1 (10 7)2 (for n > 4), n 1= n 3 n 2, and n= 1. itmax HQR0 HQR HQRu n = 4 49.0 (1349) 16.0 (16) 6.18 n = 10 40.5 (358) 16.0 (8) 6.30 n = 20 32.6 (69) 15.6 (2) 6.66 n = 30 30.4 (35) 15.4 (1) 6.93
The unimodular shift 2 does have its Achilles’heel, as described by (3.16). But HQRu takes the precautionary measure by checking this extreme condition before each QR sweep; if (3.16) is close enough, for example,
(k) n 3 1 + (k)n 2 3 (k)n 2 ! (k) n 1 < = 10 12;
then the auxiliary shift e2is applied to break up the unusual con…guration. Guarded
with this measure, HQRu is again doing better than HQR, as the data listed in the following table indicate.
Experiment 4. Randomly selected f jgn 5j=1 [ f n 3; n 2g with n 4 =
q 1 (10 7)2 (for n > 4), n 1= n 3(1 + n 2)=(3 n 2), and n= 1. itmax HQR0 HQR HQRu n = 4 7.76 7.76 4.72 n = 10 7.82 7.80 4.98 n = 20 7.89 7.93 5.62 n = 30 7.99 8.04 6.01
It is interesting to observe that, under the speci…c conditions set in this experiment, HQR0 in general works slightly better than HQR as n becomes larger.
5. Concluding Remarks. For orthogonal matrices, it appears just natural to use unimodular shifts in the QR iteration to approximate the eigenvalues. In theory we have derived the distinctive convergence properties that it possesses, and in practice
we have demonstrated its superiority in numerical convergence through both general and speci…c examples.
Though rarely encountered in practical problems, extreme examples of orthogonal matrices of order n up to 30, for which the subroutine HQR in EISPACK fails to reveal the eigenvalues (within the iteration limit 30n), can easily be found, as in Experiment 3 of Section 4, if a su¢ cient number of sample matrices, say 10; 000, is taken.
Through a better understanding of the convergence properties with the orthogonal matrices, a more e¢ cient shift strategy can be speci…cally devised in the QR iteration for such structured matrices; still, more investigations and experiments are needed for further improvements.
REFERENCES
[1] G. S. Ammar, W. B. Gragg, and L. Reichel, On the eigenproblem for orthogonal matrices, in Proceedings of the 25th IEEE conference on Decision and Control, Athens, Greece, 1986, pp. 1963–1966.
[2] S. Batterson, Convergence of the Francis shifted QR algorithm on normal matrices, Linear Algebra Appl., 207 (1994), pp. 181–195.
[3] B. Bohnhorst, A. Bunse-Gerstner, and H. Faßbender, On the perturbation theory for unitary eigenvalue problems, SIAM J. Matrix Anal. Appl., 21 (2000), pp. 809–824. [4] H. J. Buurema, A geometric proof of convergence for the QR method, doctoral dissertation,
Rijksuniversiteit Te Groningen, Groningen, The Netherlands, 1970.
[5] D. Day, How the QR algorithm fails to converge and how to …x it, Technical Report 96-0913J, Sandia National Laboratory, Albuquerque, NM, April 1996.
[6] J. W. Demmel, Applied Numerical Linear Algebra, SIAM, Philadelphia, 1997.
[7] J. G. F. Francis, The QR transformation: a unitary analogue to the LR transformation, parts I and II, Comput. J., 4 (1961–1962), pp. 265–271, 332–345.
[8] G. H. Golub and C. F. Van Loan, Matrix Computations, third edition, Johns Hopkins Univ. Press, Baltimore, 1996.
[9] W. B. Gragg, The QR algorithm for unitary Hessenberg matrices, J. Comput. Appl. Math., 16 (1986), pp. 1–8.
[10] R. S. Martin, G. Peters, and J. H. Wilkinson, The QR algorithm for real Hessenberg matrices, Numer. Math., 14 (1970), pp. 219–231.
[11] C. Moler, The QR algorithm, Cleve’s Corner, MATLAB News & Notes, The MathWorks Inc., Natick, MA, Summer 1995.
[12] B. T. Smith, J. M. Boyle, B. S. Garbow, Y. Ikebe, V. C. Klema, and C. B. Moler, Matrix Eigensystem Routines–EISPACK Guide, Springer-Verlag, Berlin, 1974.
[13] G. W. Stewart, Matrix Algorithms Volume II : Eigensystems, SIAM, Philadelphia, 2001. [14] T.-L. Wang, Invariant normal matrices under the double-shift QR iteration, preprint
submit-ted to Linear Algebra Appl. for publication, 2003.
[15] T.-L. Wang and W. B. Gragg, Convergence of the shifted QR algorithm for unitary Hessen-berg matrices, Math. Comp., 71 (2002), pp. 1473–1496.
[16] T.-L. Wang and W. B. Gragg, Convergence of the unitary QR algorithm with a unimodular Wilkinson shift, Math. Comp., 72 (2003), pp. 375–385
[17] D.S. Watkins and L. Elsner, Convergence of algorithms of decomposition type for the eigen-value problem, Linear Algebra Appl., 143 (1991), pp. 19–47.
[18] J. H. Wilkinson, The Algebraic Eigenvalue Problem, Clarendon Press, Oxford, 1965.
