4v: '1':
‘ 17.9‘ J. '-

Vsé‘. v
‘7‘ 15?.” J A
3‘. 1 V6,“ ‘
- . v

’ "V.
...

l , '

Hz - .

"EL , ‘ l
,5 ‘y ,

1 ‘

 

a "K“ I
.37" -"-_

 

\lllllllllllllllllll fillzllllll

W

293 01399

 

l

 

\\\\l\\

m ls
l

 

This is to certify that the

dissertation entitled

GUASI-LAGUERRE ITERATJON AND ITS
APPLI CATION w SOLVIIV6 SYMMETRIC

TRIDIAGONAL ElGENVALUE PROBLEMS

presented by

MING JIM

has been accepted towards fulfillment
of the requirements for

PH.D degreein MATHEMATICS

4s %’

Majo4proessof

 

Date JULY 204 L995

 

s:

 

LlBRARY
Michigan State
University

 

 

 

PLACE ll RETURN BOX to remove this chockoulfmm your record.
TO AVOID FINES Mam on or baton date duo.

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU I. An Afflnnafivo ActlonlEqual Opponunlty Institution
cum

m1

——-___

 

QUASI-LAGUERRE ITERATION AND ITS
APPLICATION IN SOLVING SYMMETRIC
TRIDIAGONAL EIGENVALUE PROBLEMS

Ming J in

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements

for the degree of

DOCTOR OF PHILOSOPHY

Department of Mathematics

1995

ABSTRACT

QUASI-LAGUERRE ITERATION AND ITS
APPLICATION IN SOLVING SYMMETRIC
TRIDIAGONAL EIGENVALUE PROBLEMS

By

Ming J in

In this work, the quasi-Laguerre iteration is established in the spirit
of Laguerre’s iteration for solving polynomial f with all real zeros. The
new iteration, which maintains the monotonicity and global conver-
gence of the Laguerre iteration, no longer needs to evaluate f". The
ultimate convergence rate is \/§ + 1.

Based on the quasi-Laguerre iteration, an algorithm is developed
and tested for solving the eigenvalue problem of symmetric tridiagonal
matrices. In general, it improves the speed of Laguerre’s iteration sub-
stantially and the accuracy is at least as good as any standard algorithm
for the symmetric tridiagonal eigenproblem. More preciously, the al-
gorithm, like the Laguerre iteration, enjoys the flexibility in evaluating
partial spectrum.

Although the quasi-Laguerre iteration is derived for solving polyno-
mials with only real zeros, the iteration itself can be extended to solve
polynomials with complex zeros and the local convergence can still be

guaranteed with the same convergence rate \/§ + 1 as the real case.

To my parents

ii

ACKNOWLEDGMENTS

I would like to express my deepest gratitude to my thesis advisor,
Professor Tien- Yien Li, for his guidance through the years I have spent
at Michigan State University. His encouragement, helpful directions
and lively discussions with me make this work possible.

I would like to thank my dissertation committee members Professor
Chichia Chin, Professor Richard Hill, Professor William T. Sledd and
Professor Zhengfang Zhou for their valuable suggestions and their time.

I would also like to thank Dr. Zhonggang Zeng for his precious advice

and collaboration.

iii

Contents

0 Introduction 1
1 Derivation of the quasi-Laguerre iteration 5
2 The formulation and convergence 13
2.1 The formulation ....................... 13
2.2 The convergence properties ................. 15
3 On symmetric tridiagonal eigenvalue problems 26
3.1 Evaluation of the logarithmic derivative f’ / f of the de—
terminant .......................... 26
3.2 The split-merge process ................... 29
3.3 Cluster and deflation .................... 32
3.4 Partial spectrum ...................... 33
3.5 Practical considerations in implementation ........ 34
3.5.1 Stopping criteria .................. 34
3.5.2 Initial points for the quasi-Laguerre iteration . . . 35
3.5.3 An alternative approach for clusters of eigenvalues 39
3.6 Numerical tests ....................... 42
3.6.1 Testing matrices .................. 45
3.6.2 Speed test in evaluating eigenvalues without com-
puting eigenvectors ................. 46
3.6.3 Speed test in evaluating eigenpairs ........ 49
3.6.4 Accuracy test .................... 52
4 The complex quasi-Laguerre iteration 54
4.1 The formulation ....................... 54
4.2 Convergence properties ................... 56

iv

List of Tables

1 Accelerated quasi-Laguerre iteration ............ 44
2 Execution time for evaluating eigenpairs ......... 51
3 Accuracy test of various algorithms ............ 53

List of Figures

OOxlovmrhOOMl—A

Algorithm DETEVL ..................... 28
Split and merge processes ................. 31
Cluster of eigenvalues .................... 32
The sign pattern of —f’/f ................. 36
Algorithm INIPTS for the initial points .......... 38
Algorithm ACCQ-LAG for acceleration .......... 43
Execution time in evaluating eigenvalues only (I) . . . . 47
Execution time in evaluating eigenvalues only (II) . . . . 48

vi

0 INTRODUCTION 1
0 Introduction

The globally and monotonically convergent Laguerre’s iteration

111(17): :1: + n (1)

(4(5)) i \ (n ‘1)[(" ‘ 1) (4(5) )2 ’ " (ffl) )]

has been successfully used to find the eigenvalues of an n x n symmetric

 

 

 

 

tridiagonal matrix T with nonzero sub-diagonal entries by approximat—

ing the zeros, always real and simple, of its characteristic polynomial
f(/\) = det[T — AI]. (2)

Remarkable numerical results, in terms of both speed and accuracy, on
a substantial variety of matrices have been obtained [8]. The algorithm
employs the split-merge process, similar to Cuppen’s divide and con-
quer strategy, which provides an excellent set of starting values making
the algorithm naturally parallel. Impressive speedups of the parallel
version of the algorithm were reported in [13]. For 2,000 x 2,000 ran-
dom matrices, the algorithm can reach a speedup of 52 on an N -cube
with 64 nodes.

Laguerre’s iteration, or modified Laguerre’s iteration when clusters
occur, converges ultimately with a very fast cubic convergence rate.
Nonetheless, the most important advantage of Laguerre’s iteration in
solving (2), in contrast to Newton’s iteration, is its monotone conver-
gence. However, comparing with Newton’s iteration, a major disadvan-
tage of Laguerre’s iteration is the evaluation of f”, which is relatively
time consuming. In this work, a new algorithm is designed in the spirit
of Laguerre’s iteration for finding zeros of (2). The algorithm, without

evaluating f”, maintains the monotonicity of Laguerre’s iteration and

0 INTRODUCTION 2

converges globally with ultimate convergence rate \/2 + 1 in finding
zeros of any polynomial with all real zeros.
The formula (1) of Laguerre’s iteration can be derived in diverse

ways. The best one seems to be to answer the following question [7]:

QUESTION 1: Among all polynomials p(:1:) of degree n with 72 real
zeros and with p(:1:0) = f(a:0) aé 0,p'(:z:0) = f’(:r0) and p”(a:0) =
f”(a:0) at a specified real 1130, which one has a zero closest to 11:0?

and where?

In general, L+(a:0) in (1) gives the closest zero from the right of all
those polynomials and L_(:1:0) gives the closest one from the left.
To avoid the evaluation of f” , we revise the above optimization prob-

lem as follows:

QUESTION 2: Given two specified reals a < b, among all polyno-
mials p(:r) of degree n with n real zeros, none of which lie in [a, b],
and with

11(3)._.f_'(_‘1_) d fl9=fl (3)

p<a>‘f<a) a“ p<b> f(b)’

which one has a zero closest to a from the right or from the left?

and where?

The optimization problem for this class of polynomials can be easily
solved and the solution has a closed form. To convert the solution of
this optimization problem to an iteration scheme, which we call the
quasi-Laguerre iteration, with the assumption that x(k‘1) and 2:“) are

available, we may

(2) let a = $(k—1), b = 17““) and let 513““) be the closest zero to 11:“)

from the right when sew-1) < 116(k) and a zero of f is bigger than
(k)
a:

3

0 INTRODUCTION 3

(“l let a = 113(k), b = (Cw-1) and let 13““) be the closest zero to 22““)

from the left when 550-1) > 27““) and a zero of f is less than 17““).

In this manner, the sequence generated will converge, in fact, mono-
tonically to a zero of f.

The optimization problem in Question 2 is formulated in a more
general form in Chapter 1 to account for multiple zeros. The answer
to this problem leads to the quasi-Laguerre iteration. In Section 2.1,
the iteration is formulated. The ultimate rate of convergence of our
iteration so obtained is \/2 + 1 which Will be shown in Section 2.2.

One of our main goals is to use the newly derived quasi-Laguerre
iteration for approximating eigenvalues of symmetric tridiagonal ma-
trices. The details of this implementation are described in Chapter 3,
including some practical considerations in implementation (Section 3.5)
and comprehensive numerical results on diverse types of matrices (Sec-
tion 3.6). Apparently, our algorithm improves the speed of Laguerre’s
iteration considerably in general. Like Laguerre’s iteration, our algo-
rithm is inherently parallel and has a great capacity for vectorization.
An intensive experiment in this regard will be reported in a future work.

Although the quasi-Laguerre iteration is derived for solving polyno-
mials with only real zeros, the iteration itself can be extended to solve
polynomials with complex zeros in the same way the Laguerre iteration
is extended. In Chapter 4, we will first define and formulate the com-
plex quasi-Laguerre iteration. Then in Section 4.2, we will show that
the local convergence of our iteration is guaranteed for simple zeros
and the ultimate convergence rate is \/2 + 1, similar to the real case.
Comparing with Newton’s and Laguerre’s iteration, the most attractive
aspect of our quasi-Laguerre iteration in complex case is its fast con-

vergence as well as the elimination of evaluating the second derivative,

0 INTROD UCTION

even though it no longer enjoys the monotone or global convergence.

1 DERIVATION OF THE QUASI—LAGUERRE ITERATION 5
1 Derivation of the quasi-Laguerre iteration

Let f be a polynomial of degree n with all of its zeros being real. For
a < b with f(a)f(b) 75 0, let .7: be the class of polynomials p(a:) of

degree n that satisfy the following conditions:

(i) all zeros of p(a:) are real;

(ii) none of the zeros of p(a:) are in [a,b];

(M) p<a>p<b> at 0 and

 

 

 

 

We pose an optimization problem more general than Question 2:

QUESTION 3: For 19 E .7, consider its m-th (m < n) zero to the
right (or to the left) of a. Which polynomial in f has the closest

one to a from the right (or from the left)? and where?

To answer this question, let z1,22, - - - , 2,, be all the real zeros of a

given polynomial p(a:) in .7: such that
Z1SZ2S“°SZk<a<b<Zk+1S"'Szn- (5)

It follows from (4) that

W)— " 1 _ p'<b)_ " _
P(a)—.-=Zia—z.-_ mpg, ,—q<b>. (6)

ForlSiSnJet

 

 

 

 

b—Zi b—a
+

a—z, a—zi

 

 

then, X,- > 0 for all i since no 2,: falls between a and b. Moreover,

 

iXi=i(l+:—::-i)=n+(b—a)q(a)Er(a,b), ()
zznl 1:711 7
z 31: :- 2; (1+ 5:) = n+ (a — b)q(b) = r(b,a),

' 1 1

S
N

1 DERIVAT ION OF THE QUASI—LAGUERRE ITERATION 6

and from (5),
XkZXk—IZ"'ZX1>1>XnZXn—1Z"'2Xk+1>0- (8)
For convenience, we rewrite (8) as
W12W22~2Wk>1>Wk+1z-o-ZW. (9)

with

X _ f l < k
m = k+l l OI’ _ (10)
Xn+k+1_[ for l > k.

From (5), the m-th zero of p(a:) to the left of a is zk_m+1 which, from
(10), corresponds to Wm in (9), and minimizing (a —— zk_m+1) becomes
maximizing Wm. Similarly, the m-th zero of p(:c) to the right of a is
zk+m and minimizing (zk+m — a) becomes minimizing Wn_m+1. As a
convention, if m > k then the m-th zero of p to the left of a will be
taken as the (n — m + 1)-st zero to the right of a and if m > n — k then
the m-th zero of p to the right of a is the (n — m +1)-st zero to the left
of a.

The optimization problem in Question 3 can now be converted to

the following optimizations:

(P1)
Max Wm
subject to Z W,- = r(a, b)
1 b (11)
E W: -’ 74( 3a)

1 DERIVATION OF THE QUASI—LAGUERRE ITERATION 7

(P2)
M2” I/Vn—m+1
subject to Z W,- = r(a, b)
i=1
1.21% = r(b, a)

W12W2.>_"'2Wn>0-

The solution of (P1) gives the m-th closest zero to a from the left
while the solution of (P2) gives the m—th closest zero to a from the

right.
We first consider the optimization problem (P1).
Lemma 1.1 r(a,b) > 0,r(b,a) > 0 and r(a,b)r(b,a) 2 n2.

Proof. By (8), we have r(a, b) > 0 and r(b, a) > 0. From the Cauchy-

Schwartz inequality,

r(a, b)r(b, a) (because of (11: ). D

 

 

 

Lemma 1.2 The point W = (W1, - - - ,Wn) satisfying (11) with W1:

 

W2=~~=Wm=~~=W—n_1 2 Wu > 0 and
_ b
Wm Z T(a7 ) (12)
n
is feasible to problem (P1).
Proof. From (11), Wm and V7,, must satisfy
(n—1)Wm+Wn = r(a,b). (13)
—1 1
TL. + __ = r(b,a). (14)

 

 

Wm Wn

1 DERIVATION OF THE QUASI—LAGUERRE ITERATION 8

Eliminating Wm in (13) and (14) yields
r(b,a)Wi — [r(a,b)r(b,a) — n2 + 277.] W, + r(a,b) = 0. (15)

By Lemma 1.1, it can be easily verified that both solutions WE,” and
WE?) of (15) are real. In addition, they are positive since W£1)+W[,2) > 0
and Wian) > 0. Let W?) 3 W53). Since
:(a, b): r,(a b)r( n

n 0)

 

 

WSW?) :-

s(

so, W[,1)< M. Substituting Wnl) into (13),

 

102) (because of Lemma 1. 1 ),

 

 

W“) = 1 (r(a, b) — W‘”) >

"I ”—1 n — n

It follows that Wm > W41 and so, W = (an, --,W,(n),

feasible. D

n )9

Now, for any W 2 (W1, - -- ,W n) satisfying (11), we have
(Wj—Wm)2 " (W; W) " ( W2)

jg" Wj g1 WJ H]; WJ

= r,(a b) — 2nWm + r(b, a)W,3,. (16)

 

 

K}.

On the other hand,
2(Wm—Wj) = (n-m)Wm- Z W,-

j>m j>m
= (n — m)Wm — r(a,b) + Z W]-

JSm

2 nVVm — r(a, b).

From (12), nWm — r(a,b) Z 0 for a certain feasible point, so we may

assume, in general, Wm satisfies

nWm — r(a,b) Z 0

1 DERIVATION OF THE QUASI—LAGUERRE ITERATION 9

since we are maximizing Wm. It follows that

 

[nWm _ ,(a, 5)]? g [3;(Wm — Wj)[2 = [g (Es/$1 . [HZ-Hz
s z: (W: X“) 2%. W].

(by Cauchy-Schwartz inequality)
3 [r(a,b) -— 2nWm + r(b, a)W3,] (r(a, b) — me)
(because of (16) ).

Or,
mr(b, a)W3, — [7°(a,b)r(b,a) — n2 + 2mm] Wm + mr(a, b) g 0 (17)
and the equality holds iff W, = Wm for all j < m and W; = Wn for all

l>m.

Similarly, when solving the optimization problem (P2), one can show

that, with the assumption
TU), a’)Wn—m+l T" n S O

which we may assume without loss of generality, Wn_m+1 satisfies the

same inequality, that is,

mr(b, a)W3_m+1 — [r(a, b)r(b, a) — n2 + 2mn]Wn_m+1+ mr(a, b) S 0

(18)
and the equality holds iff W, = Wn_m+1 for all j > n — m + 1 and
l/V1=W1foralll<n—m+1.

Now, taking equalities in (17) and (18) yields the quadratic equation,

mr(b, a)Y2 — [r(a, b)r(b, a) — n2 + 2mm] Y + mr(a, b) = 0.

1 DERIVATION OF THE QUASI—LAGUERRE ITERATION 10

Its solutions Ymi =

 

r(a,b)r(b, a) —n2+2mn:l:\/[r(a,b)r(b, a) —n2] [r(a,b)r(b, a) — (n—2m)2]
2mr(b,a)

 

are both real by Lemma 1.1. Furthermore, since
rm, + Y ,_ = r(a,b)r(b, a) — n2 + 2mm 2 2mm > 0

and

r(a,b)
r(b,a)
so, Ymi > 0. On the other hand, it follows from (17) and (18) that

Ym+-Y _=

 

>0,

)fm- S Wm _<_. Ym+ (19)

and
Ym— S Wn—m+1 S Ym+- (20)
Accordingly, W 2 (W1,- - - ,Wn) with W1 2 2 Wm = Ym+ and
Wm+1 = = W" is the solution of the maximization problem (P1),

and, from (11),
mYm+ + (n — m)Wn = r(a, b),

m n—m

— —:— = b .
Ym+ + Wn T( ’0’)
Solving for W”, one can easily see that Wu is actually I/(n_m)_. So,
_ Ym ' < m
WJ = + J _ (21)
3/(n—m)— .7 > m°

Similarly, the solution of the minimization problem (P2) is W 2 (W1,

--,Wn) with

(22)

Y(n—m)+ j< n—m'l'l
Ym_ jZn—m+1.

1 DERIVATION OF THE QUASI-LAGUERRE ITERATION 11

Recall that Wj’s are defined by
b — a

a—uj

 

Wj=1+ j=1,---,n

where uj’s are zeros of a certain polynomial in 3:. Thus any polynomial

satisfying (21) must be in the following form
P103) = C(13 - urn—WU? — ”On—m (23)

where C is a constant, and

 

 

b—a b—a
1 —— = m ' 1 t1 m- = — —— 24
+ a _ Um_ + or equlva en y u a Ym+ _ 1 ( )
and
b — b _
1 + a = I/(n_m)_ or equivalently v1 = a — a .
a _ U1 }/(n—m)— _ 1

And any polynomial satisfying (22) must be in the form
192(17) 2 C(1: — um+)m(.r — v2)"_’" (25)

where C is a constant, and

 

 

b - a b — a
1 + —— = K”. or u'val tl m = — —— 26
a—um+ eq1 eny 11+ a Ym——1 ( )
and
b — a b — a
1 + = Yn_m or e uivalentl v = — .
a "' ’02 ( )+ q y 2 a nn—m)+ _ 1

From (24) and (26),

b — um_ b — m

——= m+>0 and —u—+-=Ym_>0,

a — Um_ (I — Um+
so both um+ and um_ are outside the interval (a, b). Similarly, both 121
and 112 are outside the interval (a, b) since Y(,,_m)i > 0. Consequently,

both 101(1)) and 192(3)) are in .73.

1 DERIVATION OF THE QUASI—LAGUERRE ITERATION 12

Moreover, from (24) and (26)
b — a
Ym; — 1

 

’llmi = a —
— 2mr(b,a)(b—a) _
r(a,b)r(b,a)—n2+2mn—2mr(b,a)2}:\/[r(a,b)r(b,a)—n2][r(a,b)r(b,a)—(n—2m)2]
2min — (b — 604(5)]

 

 

 

 

= a + —(b — a)R — 2mq(b) :l: \/R[(b — a)2R + 4m(n — m)] (27)
where

_fW0 _fw)m, =naM%20)_ a
gm)—fmy am—fw, ciR ( ,_a ) q(mo)

The solutions of the optimization problems in Question 3 can now

be described as follows:

(i) If f has at least m (counting multiplicities) zeros to the left of a,

then the zero um_ in (27) of the polynomial
mb)=ce—us)%x—mvv:

in (23) is the closest m-th zero to the left of a among all polyno-

mials in .77 . And it is clear that
zm_§um__<_---_<_u2_Su1_<a<b (28)
where 2",. is the m-th zero of f to the left of a.
(ii) If f has at least m (counting multiplicities) zeros to the right of b,

then the zero um+ in (27) of the polynomial

m

192(13) = C(IB - um+)m($ - v2)"-
in (25) is the closest m-th zero to the right of b among all poly-

nomials in f . It is also clear that

a < b < “1+ _<_ “2+ S _<_ Um+ S Zm+ (29)

where zm+ is the m-th zero of f to the right of b.

2 THE FORMULATION AND CONVERGENCE 13
2 The formulation and convergence

2.1 The formulation

Let f be a polynomial of degree n with all its zeros being real. Let 2M
be a zero of f with multiplicity M Z 1. To approximate 2M from its
left, we shall use um+ in (27) with m g M to generate a monotonically
increasing sequence {$151} which converges to 2M as k —-) 00. Similarly,
to approximate 2M from its right, a monotonically decreasing sequence
{13%).} converging to 2M, as k —-) 00, can be generated by using um- in
(27).

To be more precise, suppose 2M is on the right hand side of two
starting points 1:121 < 159,11 and none of the zeros of f lie between 2:121

and 2M. Then, for k 2 1, we may let

 

 

 

22251” = wlfil”+
2m [ — (2)21 — Eflmath]

1.. (30)
_($(n-)l-— —$1n+1 )_)R 2mq(a: $171+) W¢R( xm+ _xm+ 1)_)2R+4m(n m)]
where (k_1) (k)

(k——l) _ f($m+ ) ,(k) __ f'($m+)
Q($m+ )_ (k 1) 7 q(LLm+) _ (k)
f(xm+ ) f($m+)
and (k 1) (k)
9013me )_9($m I k— 1 k
R = n (k) (k_1) + — 9(1 (12+ )chl’inil-
a: —:1:
m+ m+
In other wo d l b d ' (k 1) (k) “(+1)
r s, we rep ace a, an um+ 1n (27) by xm+ ,aa:m+ and at",
respectively.
Remark 2.1 Interchanging 13:54:” and 11,51, (30) can also be written
as
(k+1) (k)

$m+ Z $m++

2 THE FORMULATION AND CONVERGENCE 14

k—l L k—l
2m [n — (12., 1 — 12111.12 ’1]

 

.(31)

 

—(a:£,f;l)—a:§: +)R— 21nq( a:,,,+1) )+\/R([.1:m+1) —:1:m+)2R+4m(n— m)]

From what have been derived in the last Chapter, it is easy to see

that
k k-l-l
{Bin-)+— < min-l- ) S 2(M

and no zeros of f fall between 2:11:11) and 2M. So the sequence {13%.}

satisfies

0 1 k
332111 < $(n)+ < < 33an < S 2M-

Similarly, when two starting points (13912. < 2:121 are available on

the right hand side of zM with no zeros of f lie 1n (2M,:cm_), then for

k 2 1, let a = ESL, b = .’L‘[,€_1) and um- = 13%?) in (27). Namely,

(1+1)_ (1)

 

 

 

 

xm— — aim—4'
2m [nm — ((5731) — mg1)q(x)£fi_lm))]

, (32)
—(.1:£n:1)—a:m_)R— 2mq(:r 1))—((\/R[(:1:,,,_1) m)_)—-2R+4m(n m)]
where (1: k 1

(k) _ f’(~’L‘m)—) (k—l) _ f’(113(n: ))
Q(~’L‘m_) — (k) a q(xm- ) — —TT)—
f($m—) f(5cm—- )
and (k) (k 1)
q($m—) — q($m: ) k k—l
R :71 ( (k—l) (k) ) '7 q(4r(nZ-)q(x(n— ))'
.rm_ — 13m-
Again it is clear that
Z)» M < A: +) < xiii).
and none of the zeros of f lie between 2M and 2:11:11). Moreover, the

sequence {mffil} satisfies

zMg-~<x§,'§)_<-~<a:$,1,)_<x§3)_.

2 THE FORMULATION AND CONVERGENCE 15

We call the sequences {xiii} defined by (30) (or (31) ) and (32) the

quasi-Laguerre iterations.

2.2 The convergence properties
Theorem 2.2 Both sequences {Elsi} converge to zM monotonically.

Proof. We will prove the theorem for the sequence {x151 }‘,§°=1. Similar
arguments hold for the sequence {xifiifiip We shall write yk for 13:71:).
when there is no ambiguity. Obviously, by the way it is constructed,
{yk} is a decreasing sequence bounded below by 2M. Let 2(1), - - - , 2("’M)
be the rest of the zeros of f besides 2M. Recall that given yk_1 and yk,

yk+1 satisfies

 

 

m n — m M ”—M 1
— + = + ——p
111 — yk+1 31k — ”U yk — ZM 1:1 111 —,z
m n — m M "—M 1

 

 

Elk—1"" yk+1 111—1- 11 _ y1¢_1— zM ,-___1 yk_1 — 2(1)
for a certain real number 11 ¢ (yk, yk_1). Subtracting the second from

the first and dividing the result by (311-1 —- yk) yields,

 

 

m + n — m
(311 - yk+l)(yk—1 - yk+1) (31k - ’U)(yk—1 - ’0)
M "‘M l

 

 

(yk - 2M)(y1—1 - 2M) + 1:21 (we - 2(1))(111—1 - 200'

Each term in the above equation is positive and 2M is the closest zero

 

 

 

 

to yk, so
m
00 < as k —-) 00
(31k — yk+1)(yk—1 - 9H1) ( )
M "‘M 1
< + . .
(111 — mum/1-1 — 2M) g (31k — 212))(121— M)
M "‘M 1

 

 

(yk — ZM)(1Jk-—1 — ZM)

2 THE FORMULATION AND CONVERGENCE 16

Therefore yk ——> 2M. 1]

The rate of convergence of both sequences {xiii} is given in the

following theorem.

Theorem 2.3 (i) When m < M, the convergence of {xiii} to 2M is

linear and the convergence ratio

 

 

 

(1+1 _Z
. mi 11/!
w_11m
11—100 $(3t-ZM
is the only real solution of
:1: —a: —a: :0
n—M + M

in (0,1).

(ii) When m = M, the convergence of {$153,} to 2M is super-linear

with convergence rate V5 + 1.

Proof. Again we shall prove the theorem for the sequence {mlfilfiil

and write yk for 2:55). when there is no ambiguity. Similar arguments

hold for the sequence {xifiiflir
Recall that

yk+1 = Elk—
2m[n - (vi—1 - yk)q(yk_1)]
(we-1 — 110113 + 2mg(yk—1)+ \/R[(3/k—1 — 11021? + 4mm - 171)]

 

 

(33)

 

where

_ f’(3/k—1) _
9(91—1)- f(y1s_1)’ Q(yk)— f(yk)

 

 

and

 

_ (1(3/kl-Q(yk—1)
R—n( sis—y). )—q(yk)q(yk-1).

2 THE FORMULATION AND CONVERGENCE

 

 

 

 

 

17
Let z“) ,z(" M be the rest of the zeros of f besides 2M. Then
f'(yk_1) M
_ : : + _ ,
q(yk 1) f(yk—1) yk—1 — ZM Qty]: 1)
f,(3/k M . "’M
and

R = n (MEL—133?”) — 1(ys)q(y1_1)

: MW “ M) _ MQ(yk—1) _ M62011)
(yk—ZM)(yk-1—2M) yk — 2M

— Q(y1)Q(y1_1)

yk—l-ZM
R ' h R EM 1
+n W1t E . . .
1 1 121 (3111—1 — 2(1))(yk — 3(1))

Since yk —1 2M as k —) 00, so

 

 

 

 

 

1 1
k- 2 )—* < oo,
Q<y 1) "1:2:1 yk—1_ZQ:HZ;Z:1MZM—Z(——Tl
o( ) "EM 1 o
91: = . —)
i=1 yk— 2(2)
and
R EM 1 R .211 1
= ' - ‘2 * ‘3 . < oo.
1 (yk—l—z<l>)(yk—z<z>) (zM—zw)?
Let 6k = Ali-Ta"— V]; Z 1 then yk—l-yk _
Elk—1 ZM

yk-I‘ZM — 1 — 6],. Obviously,
0<ek<1and0<1—ek<1. From (33),

yk+1 — ZM _ 1_
yk - 2M

2min—(yk—l_yk)q(yk—l)l . (34)
(vie-2M)[(yk—l—y1~)R+2mQ(yk—1)+x/R[(yk—1-yk)2R+4m(n-m)l}

Now,

 

(a) (yk—l — yk)9(3/k—1) = (91—1 - yk)(yk_1_ ZM + Q(yk—1))

= M(1—€k)+Q(€k)(l—€k)(yk_1—ZM) = M(1—5k)+0(yk—1_ZM)

2 THE FORMULATION AND CONVERGENCE 18

(b) (yk-ZM)(yk—1"ZM)R = M(n_M) — (311—1 - 2M)MQ(yk—1)
—(3/k — ZMIMQWU — Q(yk)Q(yk—1)(yk-ZM)(yk—1-ZM)
+n(yk - ZM)(3/k—1 — ZMlRl = MOI - M) + 0(y1_1— ZM)

(C) (yk — 211)(;1/1_1- y1)R = (1 - 600/1- - ZM)(3/1~—1 - 2M)R
= (1 — ek)M(n — A!) + 0(y1._1 — 2M) (because of (b) above)

(d) (M — 3M)2R = 61(3/1 - ZM)(yk—1 - 2M)R
= eleI(n — M) + 0(yk_1 — zM) (because of (b) above)

(e) (yk — 2M)?R [(yk_1 — yk)2R + 4m(n — m)]
= [0/1 - 2M)(yk—1 - .111)ng + 47M" — m)(yk - ZMVR
=(1"61-)21V12(n—1V1)2 + 4m(n—m)€kM(n-M) + 0(yk—1-ZM)
(because of (c) and (d) above).

Substituting (a), (c), (e) into (34) yields

Zl/k+1 - 2M _
= m *1“

Zlnln-(l—CklMl‘l'OWk—l"ZMI
M(n-M)(1—€k)+2mM€k+M(n-M)\/(1-€k)2+Wék+0(yk—1—2M)
= C(61)+ 0(y1_1 - 2M) (35)

 

where
2m[n — (1 — e)M]
M(n—M)(1—€) + 2mMe + M(n—M)\/(1—¢E)2 + (WWW

 

 

G(e) =1—

We claim that 0(6) can be rewritten as
2(M — m)
M(1+e+[/(1—6)2+4M"%Ej[—4n%e

 

 

G(6) = (36)

 

To see this, let A = \/(1— 6)2 + We and

A = m[n— (1—6)M]

2 THE FORMULATION AND CONVERGENCE 19

B = (n—M)(1—e)+2me+(n—M)A
C = .M—m
D = (1+€)+A.

Then, C(e) =1 — 2% and

AD 2 m(n—.M)(1+e)+m1Me(1+e)+[m(n—lW)+mMe]A
BC = (M—m)(n—M)(1--e)+2(M—m)me+(M—m)(n—M)A.
So,
AD+BC = mil162+(A/IQ—Mn+2mn+Mm—2m2)e
+M(n—M)+[M(n—M)+mM€]A.
Meanwhile,
MBD = Al(n—M)(1—62)+M(n—M)(1—e)2+4m(n—m)e+2mMe
+2mM62+[M(n—M)(1+e)+M(n—M)(1—€)+2mMe]A
= 2mMe2+(2M2—2Mn+4mn+2Mm—4m2)e
+2M(n—M)+[2M(n—M)+2mMe]A

 

 

 

: 2(AD+BC).
Thus,
1 A A C C AD+BC C M C
aMll‘G<€)1-§—§+5‘5- BD ‘B-T'E
Therefore,
G() 2C 2(M—m)
e = —= .
MD M<1+€+\/(1—€)2+%7%?_——$‘)l€)

In this form, it can be easily shown that C(e) is monotonically decreas-

ing. To prove this, let

 

m(n — m)

h(€)=1+€+\/(1—€)2+4’7€ where 7=M(n—M)'

 

2 THE FORMULATION AND CONVERGENCE 20

 

 

 

 

For 6 6 (0,1),
— 1 2
h’(e) 2 1+ 6 + 7 > 0,
\/(1 — (E)2 + 476

so, from (36), C(e) is monotonically decreasing.

Now,

(i) For m < .M, let

M — m m(n — m)
(3 A4 €(0,) an 7 M(n—M)>0

Since C(e) is monotonically decreasing for e E [0,1] and 6k 6 (0,1) for
all k 2 1, we have

 

fl
0< =G1<G.<GO= <1.
1+¢1 () (s) () 5
Hence, when k is sufficiently large,
G1 G0 +1
0<-——(2-—)-<ek+1=G(ek)+O(yk_1—ZM)<~i~%———<1.

Consequently, the convergence of {ykhiozl to 2M is linear.

Moreover, (35) implies that
|€k+1 - G(€k)| = |O(y1_1 — ZM)| S C(yk—l - ZM) for k 2 1 (37)

for some positive number C independent of 11. Since yk ——> zM, there

exists an integer N > 1 such that
G(1) — C(yN_.1— 2M) > 0 and C(0) + C(yN—l — 2M) < 1. (38)
Let
01 2: G(1) — C(yN_1— 2M), d1: G(0) + C(yNfll —- 2M),

Ck+1 = C(dk) " C(yN+k—1 — 2M), dk+1 = C(Ck) + C(yN+k—1 — 2M)»
for k 2 1. Then,

2 THE FORMULATION AND CONVERGENCE 21

(i) {ck} is monotonically increasing and {dk} is monotonically de-

creasing, namely,

0<C1<C2<C3<°"<Ck<"'<dk<"'<d3<d2<d1<1.

(ll) Ck < €N+k < (if, for k 2 1.

This can be shown by mathematical induction. For k = 1 ,

0 < C1 2 C(1) — C(yN_1— ZM) (because of (38) )
< C(eN) — C(yN_1— 2M) (0<€N<1 and G is decreasing)
< €N+1 < C(eN) + C(yN_1— 2M) (because of (37) )

< 6(0) + C(yN_1— 2M) 2 d1 < 1. (because of (38) )
So, both (i) and (ii) hold. Suppose the assertions hold for k 2 1, that
is,

0<c1<-~-<c1._1 <ck<eN+k<dk<d1_1<---<d1<1. (39)

Then it follows from the monotonicity of both C and the sequence {yk}
that

Ck+1: G(dk) — C(yN+k—l — 2M) > G(dk-1) — C(yN+k—2 - ZM) = Ck
dk+1 = C(01) + C(yN+k—1 — 2M) < G(Ck—1)+ C(yN+k—2 - 2M) = (1k-
Furthermore,

Ck+1 = G(dk)-C(yN+k—1-ZM)
< C(6N+k)—C(yN+k_1—2M) (because of (39) )
< eN+k+1 < G(€N+k) + C(yN+k_1 — 2M) (because of (37))
< C(ck) + C(yN+k_1 — 2M) = dk+1. (because of (39) )

Hence,

0<cl<~~<ck<ck+1<eN+k+1<dk+1<dk<~~<d1<1.

2 THE FORMULATION AND CONVERGENCE 22

We now let ck —) c and dk ——> d. Then C(c) = d and C(11) 2
Since C(e) is monotonically decreasing, G2 = G o G is monotonically
increasing and can have at most one fixed point. It follows that c = d

and consequently C(w) = w, where

 

 

 

 

 

 

 

 

 

 

 

yk—
w— — lim 6),: lim
k—ioo k—mo yk_1 — ZM
Namely,
2(M — m)
w :
All (1+w+\/(1 —w)2+%fi%w)
or, w is a solution to the equation
g(x)En;£lA;mx3—x2—E+MA_J-m=0. (40)
This equation has exactly one real solution in (0,1), since
M — m
n—M—m M—m n—M—m M-—-m
= —1—1 = —— —
9(1) n—M + M [n—M 1l+i M ll
_ —m + —m _ —mn < 0
_ n — M M _ M(n — M) ’
and

g1($)_3"__A_LMfl$2 _2$—1<0 forxe(01)

7?,

(ii) For m = M, from (36), C(61) = 0 for k 2 1; accordingly, from
(35),

61,11: 0(yk_1-— 2M) —) 0, as k —-) 00.

Thus the convergence of {ykfiozl to 2M is super-linear. Further, from
the derivation of yk+1 described in Chapter 1, yk_1, yk, and yk+1 satisfy

the following two equations

M n — M M
—— + = 1011) = 621111) +
311 - yk+1 yk - v 91 - 2M

 

 

2 THE FORMULATION AND CONVERGENCE 23

 

 

 

 

 

 

 

M n — .M 11/!
+ = (1011—1) = Q(yk—1)+
Elk—1 - yk+1 Elk-1 — ”U Elk—1 — 2M
for a certain number v. From the above two equations, we have
n — M M Ms 1 — ZM
= 12(1).) — ( + l (41)
21/1 - v (311 — 210(11): — 11m)
11 — M M yr — ZM
ilk—1 — U (Elk—1 — ZM)(3/k—1 — yk+1)

Subtracting (42) from (41) and dividing the result by (yk_1 — yk) yields,

n—M =R1— M(yk+1-ZM)(yk-1+yk—yk+1—ZM)

(yk—v)(1/k—1 -’U) (yk—ZM)(yk—yk+1)(yk—1 -ZM)(yk—1—yk+1)'
(43)

Substituting (41) and (42) into (43) and multiplying both sides by
(n — 114), we have

_ M(yk+1 - 2M) _ M(yk+1 - ZM)
(Q(yk) (yk _ ZMlUlk — yk+1)) (ka-l) (yk—l — ZM)(1/k—1— yk+1))

_ ,,_ __ M(n—MXZ/kH—ZM)(yk-1+yk—yk+1-ZM)
_ ( MRI (.111 - 2M)(yk - yk+1)(y1_1— 2M)(yk—1 — 111.11) '

 

 

 

 

 

Rewrite the above equation as follows:

yk+1 — 2M
(Elk - ZM)(3/k-yk+1)(yk—1—ZM)

 

'13 =(n-M)R1— Q(y1)Q(;1/1_1) (44)

 

 

 

 

where
yk—l + yk — yk1- ZM
E E M(n — M) + — MQ(yk—1)(yk—1 — ZM)
yk—l — yk+1
-MQ(yk)(yk _ 2111) Me — yk+1 + M2 yk+1— ZM
yk—l — yk+1 yk—l — yk+1
6/:
= M(n — M)(1+ -,———) — M12<yk_1)(yk_1— 2M)
‘— €k€k+1
e 1 — e 6 e
—MQ(y1)(yk _ ZM) k( (+1) + M2__£_’°i_
1 - €k€k+l 1 — €k€k+1

—> M(n—M) as k—)oo.

2 THE FORMULATION AND CONVERGENCE 24

It follows from (44) that
yk+1- ZM = (.91 — yk+l> (n - MIRI - Q(3/k)Q(3/k-1)

 

 

 

(31k _ ZM)2(yk—-1 " 2M) y), — ZM E
_ _ (n — M)R1— Q(yk)Q(yk—1)
— (1 €k+1) E
(n — Mm - a?

 

—)F

<00 as k-—>oo.

M(n—M)
Iiere,

1 11—111 1 1 "—M 1 2
=_ . _ __._____ > .
F M 1:21 (mu—£50»2 M(n—M) (1:21 2(M—2(2)) —O

 

 

If F = 0, then 2“) = 2(1), for i : 1,---,n — Ill. That is, f(:1:) has only
two distinct zeros: 2(1) and 2M. Hence 1:13,). = 2M and no convergence
discussion is necessary. Consequently, only f (:13) which has more than
two distinct zeros is considered. Then, F > 0. By the following lemma,

the convergence rate of {yk}f_’__1 to 2M is indeed \/2 + 1. El

Lemma 2.4 Suppose a sequence {ykfiil converges to 2. Let e;c be
lyk - 2|- If

61,, = Mkezek_1 and 3320 M1, = M“

where 0 < M” < 00, then the convergence rate of {yk} is x/2+ 1, i.e.,

 

€k+1 __ * 1/\/§
11—100 6,2544 — (M ) .

Proof. We assume first M * = 1, that is, lim),_,OO Mk 2 1. Let
dkzlnek+1—(\/2+1)lnek, 111721.
From ek+1 = Mkeiek_1,

1n ek+1 = 2111 e), +lnek_1+1an.

2 THE FORMULATION AND CONVERGENCE 25
So,

(1;, = 2lnek +lnek_1 +1an — (\/2+1)lnek
(1 — \/2)lne1c +lnek_1+lan

= (1 - \/2)[lnek — (Vi+ 1) 11161-1] +111 Mk
(1 — \/2)dk—1+lan.

Since0>1—\/2>—1andlanIk—>0ask——>oo,wehave
dk——)0 as k—>oo.

That is,
€k+1

ex/‘§+1——)1 as k——>oo.
k

 

Now, if [11* 79 1, let

Mk

e[.=\/M*e/c and M,',=M*.

 

Then,

M
830” : VM*ek+1= ME<VM*ek)2(VM*ek—l)

M12(e’1.)263._1 and

lim M]t = 1.
k-—>oo
Therefore, it follows from the above derivation,

el+1 :
k—>oo (ems/1+1

 

which implies

1. 6k+1 _ (VM’ilfl'H

_ 11: 1/\/2
111120 ell/5+1 '— M T (M ) ' U

 

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 26
3 On symmetric tridiagonal eigenvalue problems

In this chapter, we shall use our quasi-Laguerre iteration derived in pre-
vious chapters to approximate all eigenvalues of symmetric tridiagonal

matrices.
3.1 Evaluation of the logarithmic derivative f’ / f of the de-
terminant

Let T be a symmetric tridiagonal matrix of the form

(01 H1 l
31 a2 52 O

T=(fi._1,a.—,a)= If ,. . (45)

O fin—2 an—l flit-1
\ 1811—1 an f

 

 

We may assume, without loss of generality, that T is unreduced; that
is, O, 74 0, j = 1, - - - ,n — 1. For an unreduced T, the characteristic
polynomial

f(/\) s det(T —— )1) (46)

has only real and simple zeros ( [14], p300). In order to use our quasi-
Laguerre iteration developed in previous Chapters for finding zeros of
f (A), or the eigenvalues of T, it is necessary to evaluate f and f’ effi—
ciently with satisfactory accuracy in the first place.

It is well known that the characteristic polynomial f ()1) :— det(T —
A1) of T = [fit—11011351] and its derivative with respect to A can be

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 27

evaluated by three-term recurrences ( [14], p423):

= 1, = — A
100 P1 C11 2 . (47)
Pi = (011‘- A)/)i—1 — (31-1P1—2, l = 2131'“ an
I = 0, I = _1
pi) pl 1 2 1 - (48)
pi : ((15 — A)pi—l _ pi—l — fii—lpi—21 2 = 2131- ‘ : an

and

f(/\) = pm f'()\) = pi;-
However, these recurrences may suffer from a severe underflow/over-
flow problem and require constant testing and scaling. The following
modified recurrence equations avoid the above problem subtly by com-
puting the quotient q(A )2 f(A ) / f ( ) directly [8]. After all, only q(A)
is really needed in the formulae (30), (31) and (32). Let

 

m 102-
Si: 1 771: __1
101—1 P1
£1 = a1 _/\1 (49)
2
Si : Qi—A—?]:11,Z—2,3,‘ ana

1

770:0, 771:5,

2 (50)
771:3; (01i_)\)7li—1+1-(?:_11)77i—2]1 i:2939”°3n

2

 

and
f’(/\) __
_ _ 77”.
f (A)
To prevent the algorithm from breaking down when 5,- : 0 for some

1 S i g 11, an extra check is provided:

0 If {1 = 0 (i.e., 011: A), set £1: 252;

312—152.

oIf£,:0,i>1,set§,-= {_1

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS

 

Algorithm DET EV L:

Input: T=[)’31—1.011,fii], A
f' A
Output: _f()\) 21],, (where f(A) Edet(T—AI))
and 11(A) = neg_count
(=the number of eigenvalues of T less than A)
Begin DETEVL
£1 = 01 — /\
If 51: 0 then 51: (e2
1
z 0, : —, _ , t : 0
770 7]] £1 neg COUH
If £1 < 0 then neg-count =1

For i=2:n

2
€1=(C¥1-A)-(%:—11)

(3’2“) 2
If£,=0then§,= E: E

If f,- < 0 then neg_count = neg-count +1

2
771: "El: [(011 - (1)771—1 +1 — (21:11) 771—2]

End for
End DETEVL

 

 

 

 

Figure 1: Algorithm DETEVL

where e is the machine precision. A determinant evaluation subroutine
DETEVL (see Figure 1) is formulated according to the recurrences (49)
and (50). When 5,, i = 1,---,n, are known, the Sturm sequence is
available ([9], p47). Thus, as a by-product, DETEVL also evaluates the

number of eigenvalues of T which are less than A. This number is

denoted by 11(A).

Let A1<A2<---<A,,bethezerosoff(A)andA1<A2<-~<A,,

be the zeros of the numerical approximation f (A), it was shown in [8]

28

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 29

that
111(1) — 11s 2111,2212) + (1)-+11) + 1115. (51)

3.2 The split-merge process

Let
A1<A2<~~<An

be the zeros of f in (46). Using the quasi-Laguerre iteration to approx-
imate any /\,’, i = 1, 2, - - - ,n, it is essential to provide a pair of starting
points 23(0) and 17(1), being either (1(0) < pm < A,- or A,- < 112(1) < 113(0),

with no other Aj’s lying between 17(0), 13(1) and A). For this purpose, we

._ To 0
T_(O T1) (52)

(011 51 ) (aim _fik 18k+1 )
161 °'- °'- _ 18k+1 ' °

° ' file—1 fln—l
\ (Bk—1 air—18k} \ 312—1 an f

split the matrix T into

where

 

 

 

 

(53)
Obviously, the eigenvalues of T consist of eigenvalues of To and T1.
Without loss of generality, we may assume B, > 0, for all i = 1, 2, - - - ,
n — 1, since in (47)—(50), Bis always appear in their square form. The
following interlacing property for this rank-one tearing is important to

our algorithm.

Theorem 3.1 Let A1 < A2 < < An and A1 __<_ A2 _<_ S An be

eigenvalues of T and T respectively. Then

Sis/)1sissAQs-~~si\nsx\.<x\s+s

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 30

with the convention An+1 = An + 231,.

Proof. See [4, Theorem 8.6.2, p462].

The eigenvalues of T will be used critically to approximate the eigen-
values of T by our quasi-Laguerre iteration. We shall call this proce-
dure, splitting T into T 0 and T1 of T and using eigenvalues of T, consist-
ing of eigenvalues of T0 and T1, to approximate eigenvalues of T, the
split-merge process, similar to Cuppen’s divide-and-conquer strategy
[2], of course.

To evaluate certain eigenvalue A,-, i = 1, 2, - - - , n, of T, both A,- and
A,“ are suitable candidates for the first starting point of our algorithm,
since by Theorem 3.1,

311 S A1 S (9+1-

For the second starting point, let

;\i + 31241
2

I
and evaluate 6%? by the subroutine DETEVL given in the last section.

C
As we mentioned before, the Sturm sequence at c, which decides the

position of 0 relative to A,-, is a by-product of this evaluation. When c >

A,-, we may let 2:10) = 3241, (129) = c and when 0 < A,-, let 2):?) = A,- and

512$) = c. Then A,- < :19) < 11:10) or .1739) < 33$) < A,. Respectively,

equations (31) and (32) with m = 1 are used to generate 11$“) for

k 2 1, namely,
1'1““) = 121%

 

 

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 31

00
And sequences {13%)} are obtained which converge monotonically,

either increasingly or decreasingly, to A,- with ultimately convergence
rate \/2 + 1 by Theorem 2.3.

The eigenvalues of T in (52) consist of eigenvalues of To and T1 in
(53). To find eigenvalues of To and T1, the split-merge process described
above may be applied to T0 and T1 again. Indeed, the splitting process
should be applied to T recursively (See Figure 2) until 2 x 2 and 1 x 1

matrices are reached.

 

 

/\ 0
SP“! /\ /\ Merge
1’ Tm Tm T011) T011 T100 T1 1 T1 T111

Figure 2: Split and merge processes

After T is well split into a tree structure as shown in Figure 2, the
merging process in the reverse direction from 2 X 2 and 1 x 1 matrices

can be started. More specifically, let T, be split into T 00 and T01. Let
T00 0

‘1’, - . - A3,, be eigenvalues of To = 0 T in ascending order. Then
01

the quasi-Laguerre iteration is applied to the polynomial equation

A

f,,(A) .2 det[T, — AI] = 0

from either 1)”) or Am to obtain the corresponding eigenvalue Ag, 1 =

1,2, - - - ,m, by the merging process described above. This process is
continued until To and T1 are merged into T. That is, in the final step
all the eigenvalues of T are obtained by applying the quasi-Laguerre
iteration to f (A) = det(T — AI) from eigenvalues of T0 and T1.

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 32

3.3 Cluster and deflation

By Theorem 3.1, A,- E (A,,A,-+1) for each i = 1,---,n with the con-
vention An+1 = An + 2%. If A,“ — A,- is less than the error tolerance,
then either A,- or A,“ can be accepted as A,. In general, if T has a
cluster of r + 1 very close eigenvalues, for instance, Aj+, — A,- is less
than the error tolerance for certain 1 S j S n — r, then r eigenvalues
Aj, Aj+1, - - - ,Aj+,_1 of T can be obtained free of computations. They
can be set to any one of Aj, - - - ,Aj+,.

Since the matrix T in (45) is unreduced, its eigenvalues are all simple.
Therefore using m = 1 in the quasi-Laguerre iterations given in (31)
or (32) seems appropriate in all cases to obtain ultimate super-linear

convergence with convergence rate x/2+ 1. However, in some occasions,

there may exist a group of r > 1 eigenvalues of T, say,
Ai+1 < /\i+2 < < /\i+r

which are relatively close to each other, comparing to the distance from

(0) (1)

the two starting points, say 23+ < 23+ < A,-+1 as shown in Figure 3.

Some of them may even be numerically indistinguishable.

x2?) x1!)

A l LLLLJJII >x

A-i A1+1 11+] A'ivo-r

 

Figure 3: converging to relatively closer eigenvalues

Numerical evidence shows that to reach A,“ from 11318 )(= (9+1) and
333,1), it may take many steps of the quasi-Laguerre iteration with m = 1
before the super-linear convergence shows. In this situation, the quasi-
Laguerre iteration with m = r in (31) may be used to speed up the

convergence. Exhibited in (29), the sequences {rig} in (31) with

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 33

different m’s satisfy

k— - k 1.-
11:1 1) <.1:[k+)<a:§+)<~-<xl+) <A,-+,.

Accordingly, rig.) is relatively closer to Ai+1. Therefore if the number r

of those relatively closer zeros can be estimated, then instead of using

m = 1 as usual, we may let xii) = 51:53,). To estimate r in our algorithm,

we put A,- in the same group with A,-+1 if

[Aj — 21,42] < 0-01[:\i+2 — A141].
It might happen that Ai+1 < 12$”)

Sturm sequence given at 1:25); in this case, the number r may be re-

, revealed by the information of the

duced. Therefore an overestimate of the number r causes essentially
no harm because our algorithm can dynamically reduce r, while an

underestimate may result in slow convergence.

3.4 Partial spectrum

From the strong interlacing property given in Theorem 3.1 one can

easily obtain the following:

Proposition 3.2 If [a,b] contains k eigenvalues of T, then [a,b] con-
tains at least k—l and at most k+1 eigenvalues of T. More precisely, let
11(33) be the number of eigenvalues of T which are less than :1: E R and

As+1,- - - ,A,+k be all the eigenvalues ofT' in [a,b]. Then s—1 3 m(a) S s
ands+k~1g m(b) g s+k.

In our algorithm, either A,- or A,“ of T is used as the first starting
point of our quasi-Laguerre iteration to obtain A,- of T. To find eigenval-
ues of T in a given interval [a, b], the eigenvalues of T in [a,b] are known,

A

say A3+1,- - - ,A,+k. By evaluating m(a) and m(b), the actual member of

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 34

eigenvalues of T in [a, b] is a = 11(b)— 11(a), so that An(a)+1,- - - , A,,(a)+,,
are the eigenvalues of T in [a, b]. By Proposition 3.2, s — 1 g m(a) S s
and s+k—1_<_ 11(b) S s+k, so, a can either be k—l, k, or k+1. Thus,
at most k + 2 values are needed to be considered as the first starting

points to evaluate these a eigenvalues of T. Let

A
A A

A
A

Asza,

A

3+1 = /\s+l) )‘s+2 = )‘3-1—2) ' ' '1 )‘s+k = /\s+k1 /\s+k+1 = b-

>4»

Then 0 values among them can serve as the first staring points which

will lead to all a eigenvalues of T in [a, b].

In some occasions, the eigenvalues of T of interest are identified in
magnitude, say the largest 20%, then they can be evaluated by using
the largest 20% eigenvalues of T as starting points without computing
the other 80% eigenvalues of both T and T. i

3.5 Practical considerations in implementation

3.5.1 Stopping criteria

The following stopping criterion was suggested by Kahan [7] and used

for our implementation:
[m(Hl) _ 2301)]? S 013(k) _ 3411-1)] _ [m(kH) _ $01)”, (54)

where T is the error tolerance. This criterion is based on the following

observations. Let
5,301+” _ m(k)

m(k) _ m(k-l)

 

9k:

 

 

then as {rm ”2:, converges to A when k —-) oo, q;c is normally decreasing.
Thus

< E ]$(k+2+i) __ $(k+1+i)|

i=0

EWMHH) _ x(k+l+i))
i=0

i)‘ _ $(k+1)| :

 

 

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 35

 

 

oo k+1 __ k
_<_ [$(k+1)_ $01)] 2 q: : (Iklml ) :1:( )]
i=1 1 — qk
lx(k+1) _ $(k)]2
: (1(1) _ $11—11] _ (x1141) _ m(kll'

From (51), an obvious error tolerance at A, can be chosen as

55 .
T = 7 m1X<IaI+1111))+Ix“+1>)e.

3.5.2 Initial points for the quasi-Laguerre iteration

At every stage of the merging process, we evaluate the eigenvalues A1 <
- - - < Am of an m x m submatrix, given m initial values A1 3 - . - 3 Am

and an upper bound that interlace those 111 eigenvalues:
11311sighs-SESAME“

(see Theorem 3.1). To evaluate an eigenvalue A,- by the quasi-Laguerre
iteration, two initial points are needed on the same side of A,- without
any other eigenvalues between them. However, for i = 1, - - - ,m, the
interlacing property
311 S A1 S 3\i+1

provides only one initial point from each side. The simplest way of
obtaining the other initial point, preferably in the interval (A1, A5“), is
to bisect the interval by choosing c = ’\1_+é\1_-1-_1_ and evaluate €532 along
with the Sturm sequence at c by the algorithm DETEVL as described in
Section 3.2. If c < A,- then A,- and c are a suitable pair of initial points
for approximating A,-. Otherwise, the pair c and A,“ will be used.

While the above approach is the simplest, it may not be the most effi-
cient. Several improvements are made in our practical implementation.
First, since f’ / f is evaluated at every initial value A,- (j = 1, - - - ,m) ,

one step of Newton’s iteration can be performed without extra cost and

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 36

critical point critical point

1
1
1
1
1
1
1
1
1
1

b

 
     

 
   

 

 

---— —----¢

---- ----a
:2
X
v

critical point critical point

Figure 4: The two initial points should be between one of the nearby critical point of
f and the zero. This can be identified by the sign pattern of —f’/f.

the result can be considered as a candidate, besides c, for 512(1). Secondly
and more importantly, from our computing experience, the high order
convergence of the quasi-Laguerre iteration occurs only when there is
no critical point of f (i.e. zero of f’) between 11(0) and A,- (see Figure 4).
In other words, if 33(0) is to the left (resp. right) of A,-, then it is de—
sirable that — f’ (17(0)) / f (1(0)) is positive (resp. negative). If there is
one or more critical point in [A,-, Ai+1], then bisection or one step N ew-
ton’s iteration is used repeatedly until the above requirement at 22(0)
is satisfied. Based on these observations, our procedure for determin-
ing a pair of initial points among different choices to achieve the best
possible efficiency in practice is summarized below and the algorithm
INIPTS is given in Figure 5. For notational convenience, we assume the
eigenvalue A,- is in the interval [a, b] and there are no other eigenvalues

in this interval.
START.

I I
1. If —%§)2 < 0 and —%(£)l > 0, then neither a nor b should be

used as an initial point (there are critical points of f in both

(a,A,-) and (A,,b)). Set c = 23—9. Evaluate —%3, the Sturm

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 37

sequence at c is available in the process. Reset [a, b] to be [c, b]
or [a, c] in accordance with c < A,- or c > A,- respectively. Go

to START.

I
2. If —if-(%1 < 0 (there is a critical point between a and A;, so a
I

should not be used as an initial point ) and —%g)1 < 0, then

I
b can be one of the initial point. Set 0 = max {b — fl glib}

f(b)’ 2
and evaluate —If,—((*CE)2. If c E (/\,-, b), then b and c will be used

as a pair of initial points; go to END. Otherwise, reset [a,b]

to be [c, b]; go to START.

I I
3. If —If(%2 > 0 and _ngZ} > 0, then a can be an initial point but

b can not (there is a critical point of f between b and A.) Set

0 = min{a — 1%3’97442} and evaluate —J}%§)l. If c E (a,/\,-),

then a and c are a good pair of initial points; go to END.

Otherwise, reset [a, b] to be [a, c]; go to START.

’ ' b . . . .
4. If — f(:) > 0 and —%b_)2 < 0, then there 18 no critical pomt of

f in [a, b]. Thus either a or b can be used as one of the initial

points. Set
- _flel a_+_b - M £92
mi“ f(a)’ 2 i If ]f(a) S m]

 

 

max {b — %§)21 9-3-9} , otherwise.

and evaluate —%CE)2. Choose {a, c} or {0, b} to be the initial

points depending on c < A,- or c > A,- respectively; go to END.

END.

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS

38

 

 

Algorithm IN IPTS

Input: subscript i, initial end points 3,331,“.
Output: starting iterates r(o),x(l) and m(x(0)),n(:z:(1)),— “flon— (Wm)
Begin INIPTS
(a, b): (;\,-,;\,-+1); Evaluate —Lé3], m(a), —fT’((—b%, m(b) by DETEVL;
(1111) If —7{J) < 0 and —-(U > 0 then
c: 42L; Evaluate _ff_’((2), m(c ) by DETEVL;
(a,b)={ (c,b), if c<)\,- (i.e. m(c) <z') ; Goto (NI);
(a,,c) ifC>/\,' (i..e K(c)2i)
else if if} < 0 and ——(l < 0 then
c=max{b— 'b, 4119]» )Evaluate ——((—%, I€( c) by DETEVL;
If c> A; then (r(o),a:(l))- _ (b, c) ; Goto (t) ;
else (a, b): (c, b); Goto (it) ;
else if —%3)2 >0 and —fT(%l>0 then
c: min {a— %,—3fi}f;( Evaluate —+],( K:( )by DETEVL;
If (c< x\,) then (m(o),x(1))- — (a,c); Goto (it) ;
else (a,b) = (a,c); Goto (It) ;

C:{ min{a—-f7'£%l, 92%},
—f: , 3%}, otherwise

         

<flb§l
f(b
><l 1!
Evaluate —%(§)z, m(c) by DETEVL;
set (3(0),$(1)) : ((1,6), if C < /\,‘

(b, C), if C > Al:

End if ;

(13) End INIPTS

1W) “13(0) °

 

 

Figure 5: Algorithm INIPTS

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 39

3.5.3 An alternative approach for clusters of eigenvalues

The multiplicity estimation and the quasi-Laguerre iterations (31) and
(32) with m > 1 can successfully overcome the slow convergence when
a cluster of eigenvalues occur. An alternative approach without multi-
plicity estimation as well as dynamical adjustments of the multiplicities
can also work very efficiently in practice:

By Theorem 2.3, when the quasi-Laguerre iteration for simple zero
(m = 1) is used to approximate an M -fold zero :1:* (M _>_ 2), or a cluster
of M simple zeros, of a polynomial f of degree n, it converges linearly

with a ratio qu which is the only solution of

a: — :1: — :1: =
n — M M
in (0, 1). The following theorem gives the properties of qn, M for M _>_ 2

ananM+1.

 

 

0

Theorem 3.3 (i) For M fired, qu increases as n increases (n> M).

(ii) For a given n, qu Z qmg if 2 S M g n — 2.

(iii) When M = n — 1, qn,n_1 increases as n increases (n 2 3).
Proof. (2) For a fixed M, write

n — M — 1 3 2 M — 1

:1: —a: —:1:+

n—M M

 

 

912(33)

Then when n1 > n2,

9mm — g..<x> = (

SO,qn1’M > gn2,M.

1 1

TLQ—M nl—M

 

 

>r3>0 for xE(0,1).

(ii) For a given n, write

cc—a: —:c+ .

n—M M

 

 

9M0”)

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 40

When2SM_<_-3,

 

 

 

 

1 1 1 1
9M($)_92($) : (N—Z—R—M)x3+2—M
>1+ 1 _1_ 1
- 2 n—2 M n—M
n Tl

 

= m_M(n—M) 20 for x€(0,1).

50, gn,M 2 (111,2-
If g < M S n — 2, then 2 _<_ n — M < 721. It follows from above,
qn,n_M Z qmg. On the other hand,

gM(1‘) — gn_M(113) = (n 3M — fi) (1 — 333) > 0 for £1? 6 (0,1).

 

Then gn,M > qn,n—M- Hence: qn,M > qn,n-—M 2 qn,2-

(iii) If M = n — 1, write

n—M—l3 2 M—1 2
9,,(17) n—M :1: —:1: —a:+ M ——a: —.7:+1—n_1.

 

 

 

Then when n1 > n2,

1 1

 

 

gn1(a:)— gn2(:1:) = > 0 for :13 6 (0,1).

n2—1 n1—1

SO, gnlml—l > gngmg—l' D

For n 2 3 and M Z 2, by Theorem 3.3 (iii),
(111,114 2 (13,2 When M = n -1,
and by Theorem 3.3 (i), (ii),
qn,M 2 (171,2 2 (13,2 when M S n — 2-

In any case, q... E q3,2, the zero of 1:2 + x — 0.5 = 0 in (0, 1), is the

smallest linear ratio for multiple zeros. Let 33(1),a:(2),- - - ,x("), - -- be

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 41

the monotone sequence of the quasi-Laguerre iterates. When {513(k)}
converges super-linearly, then

|$<k+n _ mm]

9k: ‘0 as k—)oo,

wk) _ mac—1)] ’

 

so, qk should be smaller than q... when k is sufficiently large. Therefore,
when qk > q..., linear convergence seems apparent. In this case, for

sufficiently large k,

 

 

 

 

(Ic+l)_ =1:
:1: a: ,.. ...
I Q; qu or ]$(k+1)—$ I z q].”L‘(k)—.’C I where q E (171M,
I13“) — 117*]
and
_ |:r(k+1)—:r(k)| _ |x(k)—a:*|—|x(k+1)—x*| ~ (l-q)I:B(k)—a:*| ~
9k — |x(k)_x(k—l)| “‘ ]$(k—1)_x*|_|x(k)_x*| "’ (1_q)lx(k—1)_x¢| N q.

In practice, qk’s can be very close to q after several iterations. If qj = q

forj Z k then

0+1)

Ix -$*|=q|:v(j)-~’v*|, for jzk

and

__1_($(k+1) _ aw).

*:,(k)
a: :L +1-q

So, at the k—th iterate :13“), if linear convergence is revealed, namely

qk > q..., then instead of performing the regular quasi-Laguerre iteration

130*“) : (Elk) + 6k

where (5;. is the correction m(k“) — 56““) from 1:“) calculated by the quasi-
Laguerre iteration with m = 1, one can accelerate the iteration by
setting

5k

m(w) z m(n +

 

1—%°
The accelerated iteration can be summarized as follows: Assume we

are evaluating the eigenvalue X" E [a, b] and initial points 23(0), 17(1) have

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 42

been obtained by the process described in Section 3.5.2 along with data

I 13(0) I (1)
_LLFJ and 475%}

(#)

For k = 1, 2, - - - do (while the stopping criteria is not satisfied),

1. Evaluate 6;, = 517““) — :13“) using the quasi-Laguerre formula

with m =1.

. If q (Le—i 5—3: < q,., then most likely the target eigenvalue x\* is

neither a multiple zero nor a member of any cluster. Go to #
and continue the loop.

Otherwise,

(a) Let c = mm + 11:15:; and check for a possible overshot by
evaluating the Sturm sequence. If c overshoots, then re—
set c = r(k) + w for l = 8,4,2,1 until overshooting
disappears.

(b) set 113““) = c and go to #.

End for

The algorithm ACCQ-LAG in Figure 6 summarizes the fundamental

features of the acceleration process described in this section. Table 1

gives the comparison between the original quasi-Laguerre iteration and

the accelerated quasi-Laguerre iteration. The Wilkinson matrix W9?) is
used here with the starting points mm) = 11.5 and x“) = L_(2:(0)) =
11.27 targeting A23 -_- 11.0+5x 10-15 of W93, while A22 = 11.0—5x10"15.

3.6 Numerical tests

Our algorithm is implemented and tested on SPARC stations with

IEEE floating point standard. The machine precision is e z 2.2 x 10—16.

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS

43

 

Algorithm ACCQ-LAG

Input: subscript i, starting iterates 55(0), mm with Ii($(0)),
m(rm), —%§50)1)2 and —%(§Tl;))- evaluated by INIPTS.
output: the i-th eigenvalue A,- of T.
Define q... = q3,2: the positive solution of x2 + :1: — 0.5 = 0
(the smallest linear ratio for multiple zeros);
Begin ACCQ-LAG

For k: 1,2,-~-
U+, 1f 517(k) < Ai _
c = (iii are from quasi-Lag.formulae);
21-, if 1:0“) > A;

(Ska-1““: q= 3%.: t01=2-5€[ma><j(lfijl+|fij+1l)]+lz(")|€ ;
If (Icikl < tol) or (k > 2 and 17517:;me < tol) Goto (it) ;

If (q > q...) then
I .

c: m(k)+ {253; evaluate — f(cc) and m(c) by DETEVL;

If (c is on the same side of A; as r(k),a:(k‘l)) then
Goto (fit);

else, For I: 8,4,2,1,

l
c = xa) + 311.011.).
l—q ’
evaluate the Sturm sequence n(c) ;

If (c is on the same side of A,- as m(k),a:(k'l))

Goto (ti);
End for
End if
(it) 12”“) = c ;
End for

 

(11) End ACCQ-LAG

 

 

Figure 6: Algorithm ACCQ-LAG

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS

Table 1:

QOQNQO‘AOJNI—

“(A)““QMNNNNNMNNNHHV-‘Hh—I—‘D—HHH
#0383"OCDWNC3U‘AOONI-‘CCDWQC5U‘ACJNHO

513““) (original)

1 1 .500000000000000
1 1 .270700712327294
11 .149428542073966
11.0 63735435125402
1 1.0 26235830973979
11 .0 10578521008206
11.00 4261938886133
11.00 1715645778917
11.000 690622605815
11.000 277993423158
1 1 .000 111899071095
11.0000 45041905516
11.0000 18130368451
11.00000 7297872838
11.00000 2937554150
11.00000 1182429996
11.000000 475954003
11.000000 191581922
11 .0000000 77115924
1 1 .0000000 31040850
11.0000000 12494622
11.00000000 5029360
11.00000000 2024429
11.000000000 814877
1 1 .000000000 328006
1 1 .000000000 132029
1 1.0000000000 53145
1 1 .0000000000 21393
1 1.00000000000 8612
11.00000000000 3467
1 1.00000000000 1396
l 1.000000000000 565
11.000000000000 233
l 1 .000000000000 107

I(k+l)_ (1‘)]
qk: |$(k)_$( (k— 1)]

 

0. 706618
0.4 37604
0.4 17533
0.40 3427
0.40 3112
0.4025 55
0.4025 56
0.40252 7
0.40252 5
0.40252 3
0.402522
0.402522
0.402522
0.402522
0.402522
0.402522
0.402522
0.402522
0.402522
0.402522
0.402522
0.402522
0.402522
0.402520
0.40252 1
0.40252 6
0.40252 1
0.402 477
0.40 1680
0. 399604
0. 381937
0. 305425

Acceleration of the quasi-Laguerre iteration for W55

the table, the digits before a space are correct.

 

3:0“) (accelerated)

1 1. 500000000000000
11. 270700712327294
11.0 13287361699255
11 .000 124299930121
11 .0000 61856842477
11.000000 198368395
11.0000000 99023529
1 1.0000000000 12077
11.00000000000 6040
1 1.00000000000 2304
1 1 .000000000000 156
1 1.0000000000000 64

44

For each number in

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 45

3.6. 1 Testing matrices

There are 12 types of matrices employed in testing our algorithm. In the
following description of these matrix types, 0;, i = 1, - - - ,n, denote the

diagonal entries and 6;, i = 1, - - - ,n — 1, are the sub-diagonal entries.

Matrices with known eigenvalues

1. Toeplitz matrices [b, a, b]. Exact eigenvalues:

{a + 2b cos 'gkfillsksn ([5], Example 7.4, p137).

2. alza—b,a;=a for i=2,-~,n—1,an=a+b. fljzb,
j=1,---,n—1. Exact eigenvalues:

{a + 2b cos Wllgg ([5], Example 7.6, p138).
{ a for odd i
(15::

, fl,- = 1. Exact eigenvalues :
b for even i

 

 

2 "+1 (plus {a} if n is odd)
1§kSn/2
([5], Example 7.8 and 7.9, p139).

4. a; = 0, fl; = \/z(—n——i_). Exact eigenvalues:
{—n + 2k — lbs/,9. ([5], Example 7.10, p140).

5. oz.- = -[(2i—1)(n - 1) — 2(i—1)2], fl,- = i(n — i - 2). Exact
eigenvalues: {—k(k — liilsksn ([5], Example 7.11, p141).

{a+b :t \/(a—b)2+16cos2 k—”

Wilkinson and random matrices ,

6. Wilkinson matrices W; 3,- = 1,
n/2—i+1 forevennandlgign/2
i—n/2 for evennand n/2<i_<_ n
(n—1)/2—i+1 foroddnandlgig(n+1)/2
i-(n+1)/2 foroddnand(n+1)/2<i_<_n

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 46

([14], pp308-309). Most of the eigenvalues are in pairs, con-

sisting of two numerically indistinguishable eigenvalues.

7. Random matrices. a;’s and fli’s are random numbers in [0,1].

LAPACK testing matrices (generated by the LAPACK test

matrix generator [1])

8. Matrices with eigenvalues evenly distributed between its

smallest and largest eigenvalues.

9. Matrices with geometrically distributed eigenvalues. Namely,

eigenvalues can be written as {qkhsksn for some q 6 (0,1).
10. Matrices with an eigenvalue 1 and the remaining eigenvalues
in (—e,e).
11. Matrices with eigenvalues evenly distributed in the interval

(0, 1] except one eigenvalue with very small magnitude.

12. Matrices with an eigenvalue 1 and the rest of the eigenvalues

are evenly distributed in a small interval [10‘12 - 5, 10‘12 + 5].

3.6.2 Speed test in evaluating eigenvalues without computing eigenvec-

tors
We compare the performance of the following codes for evaluating eigen-
values of an n x n matrix:

(1) Q-LAG: our split-merge algorithm using the quasi—Laguerre itera-

tion. Storage requirement: 9n;

(2) LAG: split-merge algorithm using Laguerre’s iteration developed
by T.Y. Li and Z. Zeng [8]. Storage requirement: 9n;

(3) B/M: bisection/multisection subroutine DSTEBZ in LAPACK.

Storage requirement: 12n;

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 47

(a) matrices of type 1

 

 

 

 

 

 

 

 

 

 

 

 

 

11011..
19 ‘ \
3.5
a
£3 3
.23
s
32 1\%
.§
§ Q-LAG
3 1
g Rob—‘11

 

.40
o
o

i

200 300 100 A 500

Matrix size

(c) matrices of type 3

3

N
U'l

N

y.-

Execution time relative to Q—LAG

 

no
G
O

400

Matrix size

(e) matrices of type 5

L.»
U"

 

 

 

 

 

 

 

0 BM
5 3
c5 /
92.5» 1
o
g 2
9 LAG
o .
gm
5 1 Q-LAG
0.5 RFQ
iii
loo zoo no 100 500

 

 

 

 

 

 

Matrix size

(b) matrices of type 2

3.5

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

g 3
c5
2 2.5
.3
E 2
9
0 1.5
5
g 1
30.5
loo
Matrixslze
((1) matrices of type 4

(9 F/ \I W
3
9, 3
E
5 LAG
; 2 1‘ I
.5.
§ Q-LAG
§ 1
.11 W1

ioo zoo 300 A A Am A A A son

Matrixsize
(f) matrices of type 6

6
o /BIM
5
6 5 l/
9
g ‘ z/
E
E 3//
a
C 2
8 1 Pi'v ..... -—
111 Q-LAG

110 TA 210 AA A 3110 A A A mo A A 510

 

 

 

 

 

 

Matn'x size

Figure 7: Relative execution time in evaluating all eigenvalues without computing
eigenvectors. B / M: DSTEBZ; LAG: Laguerre’s Iteration; Q-LAG: the quasi-Laguerre

Iteration; RFQR: DSTERF.

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS

Execution time relative to Q-LAG

3

N
U"

Execution time relative to Q-LAG

3

2.5

2

Execution time relative to Q-LAG

qOO

(a) matrices of type 7

 

200 300 400

Matrix size

(c) matrices of type 9

 

Matrix size

(e) matrices of type 11

 

Matrix size

Execution time relative to Q-LAG Execution time relative to Q-LAG

Execution time relative to Q-LAG

(b) matrices of type 8

3

N
U!

N

1.5

5 a

200 300

Matrix size

(d) matrices of type 10

48

 

400 500

 

 

 

 

 

 

 

 

] /\
/ no“
. .......... \ .......................... 449..
0.1.10
L W

 

 

 

 

200

300

Matrix size

100 500

(f) matrices of type 12

 

 

\

 

 

 

 

Q-LAG‘

 

 

 

 

 

200 300

Matrix size

400 500

Figure 8: Relative execution time in evaluating all eigenvalues without computing

eigenvectors. B / M: DSTEBZ; LAG: Laguerre’s Iteration; Q-LAG: the quasi-Laguerre

Itera

tion; RFQR: DSTERF.

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 49

(4) RFQR: root-free QR routine DSTERF in LAPACK, as recom-
mended in LAPACK for evaluating eigenvalues only. Storage re-

quirement: Zn.

The result of this test for matrices of type 1 through type 6 is given
in Figure 7 and the result for matrices of types 7 through types 12 in
Figure 8. Matrices of type 9 to type 12 involve tiny eigenvalues or clus-
ters and are used more in stability tests (Please see Section 3.6.4). The
speed comparison in this category may not be relevant. In particular,
matrices of types 10 and 12 have big and dense clusters, and intensive
deflations make Q—LAG as well as LAG outscore RF QR in speed by a
wide margin.

For all other matrices (type 1 through type 8), Q-LAG is much faster
than LAG (especially for types 1, 4 and 5) and LAG is much faster than
B / M. RFQR is apparently the fastest code in almost all cases, except
for matrices of type 6 and type 7 where Q-LAG is about 70% and 40%
(n = 500) faster respectively. Overall, Q-LAG performs very similarly
to RFQR and is much more competitive than B / M. Taking into account
Q-LAG’s flexibility in evaluating partial spectrum and the potential in
parallel computing, Q-LAG is by far the most efficient algorithm among
all the algorithms being tested here.

3.6.3 Speed test in evaluating eigenpairs

We compare the performance of the following codes for evaluating the

eigenpairs of an n X n matrix:

(1) Q-LAG: our split-merge algorithm using the quasi-Laguerre iter-
ation plus the inverse iteration subroutine DSTEIN in LAPACK.

Storage requirement: n2 + O(n);

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 50

(2) B / M: DSTEBZ plus DSTEIN. Storage requirement: n2 + O(n);

(3) INV: DSTEIN with eigenvalues computed by Q—LAG (for analysis

purposes);

(4) D&C: divide and conquer subroutine DSTEDC in LAPACK. Stor-

age requirement: 3n2 + O(n);

(5) QR: QR routine DSTEQR in LAPACK, as recommended in LA-
PACK for evaluating eigenpairs only. Storage requirement: n2 +

O(n).

The storage requirement for D&C is relatively higher than those
which use the inverse iteration code. The results of evaluations of all
eigenpairs and 1 / 3 of largest eigenvalues with corresponding eigenvec-
tors are listed in columns A and B respectively in Table 2. In comparing
the execution time in evaluating all eigenpairs, our algorithm leads QR
substantially. For matrices of types 6, 9 and 10 where a number of
clusters exist, our algorithm is slower than D&C due to drastic defla-
tions of clusters. But this sort of matrices scarcely occurs in practice,
they are often used for testing the robustness of an algorithm. For all
other types of matrices, our algorithm leads D&C for large matrices
(n 2 300), while for small matrices, D&C is faster since it takes much
longer time for INV to compute eigenvectors in our algorithm.

Only Q—LAG, LAG and B / M can be used efficiently to evaluate 1 / 3
of the largest eigenvalues and their corresponding eigenvectors, where
D&C and QR must evaluate all the eigenpairs to obtain the largest 1 / 3.
This makes our algorithm more flexible and efficient. In comparing
with B / M in this category, our algorithm leads clearly for all types of

matrices.

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 51

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Type n = 100 n = 200 n = 300 n = 400 n = 500
A B A B A B A B A B

Q-LAG 1.15 0.42 4.31 1.49 9.95 3.42 18.1 6.49 31.1 11.2

B/M 1.69 0.60 6.84 2.42 14.9 5.17 27.7 9.57 43.9 15.8

1 INV 0.90 0.32 3.49 1.16 8.14 2.76 15.2 5.44 26.4 9.60
D&C 0.54 — 3.18 — 9.36 — 21.0 — 40.4 —

QR 1 89 — 13 8 44.0 — 102 -—- 196 —
Q-LAG 1.18 0.43 4.44 1.57 10.3 3.60 19.0 6.55 31.9 11.6

B/M 1.73 0.60 6.50 2.36 14.9 5.22 27.4 9.60 44.1 16.1

2 INV 0.92 0.31 3.50 1.22 8.17 2.84 15.6 5.31 26.4 9.58
D&C 0.76 —— 4.62 — 14.0 — 31.3 — 59.6 —

QR 1.85 — 13.9 — 45.0 — 104. — 198. —
Q-LAG 1.26 0.42 4.74 1.61 10.7 3.67 19.8 6.79 33.0 11.8

B/M 1.67 0.62 6.60 2.34 15.1 5.14 27.0 0.68 44.8 16.1

3 INV 0.93 0.31 3.56 1.21 8.21 2.81 15.8 5.47 26.8 9.78
D&C 0.77 — 4.92 —— 14.8 — 33.6 — 63.6 —

QR 1.83 — 13.1 — 40.8 — 93.4 — 179. —
Q-LAG 1.13 0.41 4.29 1.60 9.57 3.20 16.9 5.69 26.5 8.95

B/M 1.69 0.62 6.69 2.33 14.7 5.15 26.0 9.14 40.4 14.1

4 INV 0.92 0.32 3.54 1.20 7.83 2.63 14.0 4.77 21.9 7.34
D&C 0.51 — 3.18 — 9.32 — 21.6 — 39.8 ——

QR 1.78 — 13.9 — 45.2 — 105. — 203. —
Q-LAG 1.16 0.43 4.57 1.67 10.8 4.35 21.5 10.2 40.2 22.1

B/M 1.67 0.58 6.56 2.38 15.6 5.93 29.8 13.2 52.7 26.6

5 INV 0.90 0.32 3.66 1.31 8.89 3.60 18.1 8.93 35.0 20.2
D&C 0.55 — 3.23 —— 9.31 -— 21.2 — 39.1 —

QR 1.74 — 13.7 — 44.3 — 103. — 197. —
Q-LAG 1.10 0.37 4.02 1.33 8.73 2.86 15.2 5.07 23.5 7.76

B/M 1.48 0.50 5.23 1.77 11.7 4.11 20.4 6.93 31.9 10.8

6 INV 0.93 0.32 3.59 1.21 7.87 2.63 14.0 4.71 22.0 7.31
D&C 0.29 — 0.79 — 1.37 —— 2.05 — 2.77 —-

QR 1.73 — 12.5 — 41.8 — 98.8 — 183. ——
Q-LAG 1.14 0.39 4.33 1.49 9.43 3.25 16.6 5.63 25.8 8.67

B/M 1.70 0.67 6.68 2.34 14.9 5.20 26.1 9.19 40.9 14.2

7 INV 0.91 0.31 3.66 1.25 8.15 2.78 14.5 4.90 22.9 7.65
D&C 0.70 — 2.13 — 4.09 — 7.29 — 10.5 —

QR 2.01 — 16.3 — 53.7 —— 123. — 234. —-—
Q-LAG 1.24 0.44 4.79 1.69 10.6 3.58 18.8 6.36 29.1 9.85

B/M 1.71 0.66 6.67 2.41 14.9 5.14 26.0 9.08 40.5 14.1

8 INV 0.91 0.31 3.60 1.25 8.08 2.71 14.4 4.87 22.3 7.62
D&C 0.83 — 4.95 — 14.7 — 32.8 —— 59.6 —

QR 1.86 — 14.4 —- 45.0 — 104. — 200. —
Q-LAG 2.21 0.54 14.0 2.40 41.5 6.34 92.5 12.7 181. 24.3

B/M 2.48 0.65 14.9 3.01 43.6 7.51 96.4 15.0 187. 27.7

9 INV 2.00 0.44 13.2 2.13 40.1 5.77 90.0 11.8 177. 22.7
D&C 0.36 — 1.69 —— 4.19 — 8.95 — 15.3 —

QR 1.49 —— 11.2 — 35.3 — 79.1 — 157. —
Q-LAG 2.71 0.51 0.24 0.12 2.82 1.16 132. 18.4 28.4 8.68

B/M 2.57 0.52 16.8 2.76 52.4 7.83 120. 17.3 229. 30.6

10 INV 2.68 0.48 0.22 0.09 2.55 0.98 131. 18.1 27.7 8.38
D&C 0.55 - O.49 -— 0.75 — 1.38 —— 1.69 —

QR 1.56 — 11.8 —- 37.0 —— 87.4 —— 172. —
Q-LAG 1.24 0.42 4.78 1.66 10.6 3.59 18.5 6.36 28.8 9.72

B/M 1.72 0.65 6.60 2.32 14.8 5.09 26.0 9.05 40.2 14.1

11 INV 0.92 0.31 3.63 1.25 8.12 2.74 14.2 5.86 22.4 7.52
D&C 0.81 -— 4.90 —- 14.5 - 32.6 — 59.5 —

QR 1.85 — 14.7 — 45.4 — 105. — 205. —
Q-LAG 0.13 0.13 0.07 0.07 1.06 0.98 7.45 6.81 0.18 0.17

B/M 2.51 0.54 16.9 2.74 53.0 8.07 120. 17.2 229. 31.1

12 INV 0.12 0.11 0.07 0.04 1.04 0.93 7.42 6.75 0.17 0.10
D&C 0.23 -— 0.52 —- 0.79 — 1.51 — 1.75 —

QR 1.38 —— 9.36 — 31.1 — 67.9 — 138. —

 

 

 

 

 

 

 

 

Table 2: Speed test on all matrices (for types 1, 2, 3, we choose a = 2,b = 1). The execution
time in seconds in evaluating all eigenpairs is listed in column A. The execution time in seconds in
evaluating only 1/3 of the largest eigenpairs is listed in column B (Only Q-LAG and B/M have this

feature).

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 52

3.6.4 Accuracy test

For those matrices with known eigenvalues (type 1 to type 5), the ac-

curacy of a method can be determined by

X — A.-
direct error: D = max I—‘——'
' “Till
where 5‘,- is the approximation of the exact eigenvalue A,- of T and I] ~ ”1

is the 11 norm. For matrices of other types, the exact eigenvalues are
unavailable, the accuracy of a method is determined by the residual
error and the orthogonality error in the computed solution. For T =

QDQT with computed solution QDQT, these errors are computed by

residual error: 72 = max l|(TQ|;|QD)eill2

orthogonality error: (9 = max ”(QTQ - [)62'H2

 

If the residual and orthogonality error are small, then the error matrix
T — QDQT is small [11].

The results are listed in Table 3. It appears that our algorithm Q-
LAG achieves the smallest direct error on all matrices of the first 5
types. The direct error of our algorithm as well as B / M is independent
of the matrix size, whereas RFQR and QR seem to have larger direct
error when the matrix size becomes larger. For nearly all matrices,
the eigenvectors generated by D&C suffer the least loss of orthogonal-
ity, whereas our algorithm and B / M enjoy the smallest residual errors,
except on those occasions where clusters of eigenvalues exist and the

inverse iteration fails to compute orthogonal vectors.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

3 ON SYMMETRIC TRIDIAGONAL EIGENVALUE PROBLEMS 53
Matrix n = 100 n = 200 n = 300 n z 400 n = 500
Type 72 o D R o v 72 o v 72 o 1) 12 o v

Q-LAG 0.84 76.1 0.50 1.27 189. 1.00 1.41 152. 1.00 1.89 186. 1.00 2.16 158. 1.00
B/M 0.84 75.7 1.00 1.27 154. 1.00 1.49 143. 1.00 1.85 145. 1.00 1.96 173. 1.00
1 D&C 3.60 11.5 1.50 3.52 13.9 2.00 3.12 16.6 2.00 3.53 21.6 2.00 3.38 20.3 2.00
RFQR 2.97 35.4 2.50 3.72 140. 4.00 2.97 121. 3.00 3.39 161. 3.00 6.27 146. 6.50
QR 6.01 17.4 2.00 8.05 24.5 2.75 10.2 29.6 4.25 11.2 35.0 5.50 12.8 38.7 3.00
Q-LAG 0.71 66.7 1.00 1.20 156. 1.00 1.50 144. 1.00 1.61 191. 1.00 2.08 168. 1.00
B/M 0.85 130. 1.00 1.16 166. 1.00 1.49 161. 1.50 1.76 189. 1.00 2.27 171. 1.50
2 D&C 3.60 10.7 1.50 3.71 14.0 1.50 3.38 17.3 2.00 3.81 19.7 2.00 3.42 22.9 2.00
RFQR 2.35 62.5 2.50 3.52 157. 3.50 3.27 177. 3.00 3.01 80.5 3.00 3.71 149. 4.00
QR 6.15 19.0 2.25 8.46 25.6 3.00 10.3 28.4 4.06 12.0 37.9 3.69 13.0 47.1 6.13
Q-LAG 0.70 72.5 0.50 1.02 72.1 0.50 1.27 112. 0.50 1.50 137. 0.50 1.78 111. 0.50
B/M 1.13 77.9 1.00 1.13 76.3 1.00 1.16 77.1 1.00 1.39 109. 1.00 1.69 160. 1.00
3 D&C 3.67 12.1 1.00 3.80 14.8 1.50 3.50 18.5 1.13 3.99 20.9 2.00 3.93 24.2 1.00
RFQR 2.21 64.5 2.00 2.65 99.3 2.50 3.29 64.7 3.00 2.79 90.4 2.50 4.23 108. 2.50
QR 6.57 20.1 2.00 9.06 26.9 2.63 10.7 26.5 3.75 12.7 36.0 3.19 13.9 38.2 4.00
Q-LAG 0.71 11.3 0.16 0.83 14.4 0.08 1.81 51.2 0.85 1.44 47.7 0.09 1.13 48.0 0.51
B/M 0.92 12.1 0.64 1.60 45.0 1.28 1.81 72.3 0.85 1.44 56.3 1.28 2.14 121. 1.02
4 D&C 5.22 11.4 2.56 6.33 12.8 2.56 6.61 15.6 1.71 7.64 16.2 3.20 8.81 19.6 2.05
RFQR 11.0 13.1 10.9 18.0 21.5 17.9 26.9 63.1 27.3 20.6 49.5 20.5 23.1 60.4 13.3
QR 14.9 19.1 10.2 23.9 31.6 17.9 27.0 35.0 8.53 30.7 37.5 14.1 32.2 46.4 14.8
Q-LAG 0.48 21.1 0.10 0.53 19.1 0.05 0.66 23.9 0.52 0.68 18.8 0.52 0.97 32.5 0.52
B/M 0.89 24.5 0.82 0.89 54.2 0.82 0.84 68.3 1.05 0.91 40.4 1.05 1.08 56.5 1.05
5 D&C 4.05 12.2 1.64 3.78 16.9 1.64 2.98 13.6 3.15 3.94 17.5 3.15 3.22 16.1 3.15
RFQR 2.91 28.1 2.87 4.94 29.0 4.92 6.61 38.8 5.76 11.5 42.0 5.76 5.76 44.8 5.76
QR 9.06 24.7 5.74 11.1 25.8 6.55 10.8 29.6 6.29 14.4 35.9 6.29 14.6 47.3 6.29
Q-LAG 0.31 3.54 — 0.28 5.57 —— 0.38 4.03 — 0.28 5.55 — 1.02 4.53 —
B/M 1.06 4.05 — 0.92 5.57 — 0.85 4.53 — 1.36 5.56 — 1.02 4.55 —
6 D&C 2.29 16.7 — 3.22 19.6 — 3.06 15.9 —— 2.07 20.2 — 3.01 10.2 —
RFQR 6.56 7.01 — 13.3 6.02 — 18.7 4.54 — 30.6 5.54 — 19.9 4.05 —
QR 12.5 74.5 — 9.06 224. -— 16.3 360. — 14.3 486. — 12.1 632. —
Q-LAG 0.77 4.85 — 0.92 10.8 — 1.56 30.0 — 1.49 30.8 — 1.43 18.7 —
B/M 1.13 14.8 — 0.85 51.6 — 0.99 30.0 — 1.16 30.9 — 1.07 18.4 —
7 D&C 7.29 26.2 — 6.23 26.7 — 5.70 13.2 — 10.0 22.3 — 7.58 13.4 —
RFQR 10.6 5.33 — 12.3 23.8 — 24.3 14.4 — 31.1 66.6 — 24.3 23.4 —-
QR 12.0 84.4 — 20.1 180. — 25.3 249. — 31.0 253. — 31.7 343. —
Q-LAG 0.84 20.8 — 0.92 37.9 — 1.17 78.9 — 1.45 112. — 1.39 127. —
B/M 0.92 14.6 —— 0.96 37.1 — 1.36 88.1 — 1.61 104. — 1.19 109. —
8 D&C 8.11 13.7 — 6.60 16.2 — 7.90 21.1 —- 7.79 20.0 — 8.30 23.5 —
RFQR 4.66 13.5 — 6.56 40.4 — 7.30 57.3 — 6.20 82.2 — 28.9 55.8 —-
QR 12.0 24.1 — 16.4 27.6 —- 20.9 42.1 — 22.3 36.9 —— 37.3 71.3 —
Q-LAG 1.29 9.98 — 7.44 26.9 — 24.5 5.75 — 5.44 15.4 — 97.4 45.8 —
B/M 0.69 4.18 — 1.25 5.28 — 2.15 5.64 —- 10.7 47.6 — 20.5 15.6 —
9 D&C 4.26 18.2 — 6.76 30.1 — 3.94 32.3 — 5.61 28.1 — 4.74 18.0 —
RFQR 4.46 3.56 — 4.97 7.10 — 7.13 11.2 — 9.17 9.94 — 8.10 10.9 —
QR 4.72 42.3 —— 4.82 80.6 — 8.41 123. — 10.7 157. — 10.6 206. —
Q-LAG 1.50 213. — 0.91 3.14 — 2.11 39.2 — 2.99 9607 — 3.33 628. ——
B/M 0.58 257. — 1.40 45.4 — 1.28 1800 — 2.01 484. — 2.88 5683 _-
10 D&C 1.14 19.4 — 1.31 16.2 — 1.36 13.5 -— 1.41 18.0 — 1.49 15.3 —
RFQR 0.01 4.91 —- 0.22 9.09 — 0.08 10.7 — 0.17 11.4 — 0.02 13.4 ——
QR 0.03 42.3 —— 0.41 67.0 — 0.58 95.5 — 0.17 194. — 0.49 210. —
Q-LAG 0.45 9.57 — 0.98 28.9 — 0.60 42.5 — 0.69 83.4 — 0.83 61.8 —
B/M 0.81 9.64 — 1.02 40.6 -— 0.98 40.0 — 1.00 68.7 — 0.91 71.3 —
11 D&C 4.06 13.8 — 3.62 16.5 — 10.9 23.0 — 3.29 18.1 — 4.48 24.8 —
RFQR 8.48 11.7 — 3.09 25.8 — 10.9 23.0 — 14.1 54.9 — 16.7 35.0 —
QR 8.63 19.7 — 9.44 29.7 — 14.5 47.6 — 18.0 53.8 — 18.0 81.7 —
Q-LAG 1.74 4.73 — 1.75 2.17 — 1.11 13.7 —- 1.75 135. — 3.16 3.38 ——
B/M 1.23 102. — 6.48 3496 — 1.48 18.5 — 0.70 92.6 —— 0.80 43.0 —
12 D&C 0.37 17.2 —— 0.28 17.4 -— 0.52 14.0 — 1.01 20.0 —— 1.53 11.9 —
RFQR 0.23 5.25 —— 0.28 10.2 -- 0.48 12.3 — 0.02 16.1 — 0.06 14.8 -—
QR 0.45 38.6 —— 0.64 80.8 —- 0.49 148. — 0.04 228. — 0.64 358. —
Table 3: Accuracy on all matrices (for types 1, 2, 3, we choose a = 2,b = 1). All

errors in the table are relative to machine precision 6. Note that only the first 5 types

of matrices have valid direct error D.

4 THE COMPLEX QUASI—LAGUERRE ITERATION 54
4 The complex quasi-Laguerre iteration

In the rest of the work, we shall turn our attention to solving polynomial
with complex zeros. Although the quasi-Laguerre iteration is derived
for solving polynomials with only real zeros, it can be used to solve
those with complex zeros in the same way Laguerre’s iteration is used.

We shall only focus on approximating simple zeros.

4. 1 The formulation

Given a polynomial f (2:) of degree n with all of its zeros being real, the
quasi-Laguerre iteration is formulated in Section 2.1 as follows: If 2 is

a simple zero of f and is the closest zero to the right of two iterates

11:33-1) < (C(13), then 2:33“), obtained from (31) with m = 1, is the next

quasi-Laguerre iterate. Specifically, let a = $(k—l) and b = 2213), then
211‘ +1) is the larger solution of u of

 

 

 

a—u + a—v = q(a) : f(a) (55)
1 n—l f’ b
0:; b—v : q(b) : T((b—)1

since it is the solution to the optimization problem (P2) with m = 1
in Chapter 1. Similarly, if z is the closest simple zero to the left of two

iterates :r(_k) < xiii—1), then 21"“), obtained from (32) with m = 1, is

the next quasi-Laguerre iterate. If we let a = 29°) and b = 29.1), then

.779c +1) is in fact the smaller solution of u of equations (55) since it is
the solution to the optimization problem (P1) with m = 1 . We shall
use this essential property to establish the quasi-Laguerre iteration for
solving polynomials with complex zeros.

Let f (2:) be a complex polynomial of degree n (n 2 3), and 2 be a

simple zero, real or complex, of f (11:) Given two complex iterates yk_1

4 THE COMPLEX QUASI—LAGUERRE ITERATION 55

and yk for z, the next quasi-Laguerre iterate yk+1 is defined to be one

of the solutions of u for the following equations:

 

 

 

 

1 n—l _ _ f’gyk)

1].—u + yk—v — qr...) — 701.)
1 n—l _ _ f'(yk_1)
yk—l—U + 3114—1—11 _ q(y1_1) AA f(yk—l)

By solving the above equations,

a = 91-
2171 - (yr—1 — yr)q(yk—1)]
(ilk—1 — 9103 + 29041—1) :1: \/R[(yk_1— yk)2R + 4(n — 1)]

where
<>

There are two values for u in the above formula. We will choose yk+1

 

 

(56)

 

to be the one which is closer to yk, i.e., ka+1 -— ykl is smaller. More

precisely, let
A2 = Rl(yk—1 — 91023 + 4(n — 1)] and Q =(yk_1— yk)R + 2q(yk—1)-

Since

lo :l: M" = |0|2 + |A|2 :l: 2Re(r'iA),

to make lyk+1 — ykl smaller or IO :1: A] larger, we choose

 

yr. _ 2ln-(yk_§12;yAk)7(yk—1)l, if Re(§A) Z 0;

yk+1 = (57)

3}]: _ 2[n_(yk—§12:yAk)q(yk—l)l, if R€(§A) < O

 

To start the iteration, two initial points are required. For given yo, we
always apply one step of Newton’s iteration at yo to obtain the second

initial point; namely,
91 = 340 — f(yo)
f '(310)

 

4 THE COMPLEX QUASI—LAGUERRE ITERATION 56

4.2 Convergence properties

Theorem 4.1 The quasi-Laguerre iteration {yk} defined in the last

section converges to 2 when yo is close to z. The ultimate convergence
rate is x/2 + 1.

Proof. Let
yk — Z

yk-i - Z,
Clearly {yk} converges to 2 if for some r < 1, 16k] < r, Vk. From (56),

 

VkZl.

6k:

yk+1—Z :1_

 

 

 

 

 

 

 

 

 

 

€k+1 =
31k - Z
2 [n - (yr—1 — yr)q(yk—1)l . (58)
(ilk-Z) {lyre—1-yk)R+24(yk—1)i\/Rl(yk-1—yk)QR+4(n—1)l}
Let 2(1), - - - , 2("‘1) be the zeros of f besides 2. Then,
f’(yk)_ _ 1 . _ "—1 1
q(yk)= m”)— yk _ Z+Q(yk) With 62(9) = 2:31 W
_ f’—(yk 1)_ 21
(Am—1) — f(yk-l) ilk—1 - Z + Q0114)
and
R = n (dz/1:); E(::—l)) - (1(yk)(1(yk—1)
_ n -1 _ Qty/H) _ Q(yk) _
— (yr—Z )(yk 1_Z) yk—Z yk_1_an(yk)Q(yk—1)
— 1
+nR1(:l/k— 1, Me) With Rl(yk— 1, 31k) En i=1 l_(yk l—zli))(yk—z(ii)°
Now,

 

(a) (ilk—1 - yk)€1(yk—1) = (yk-1_ ilk) (yr—1 _ z + Q(yk—1))
=(1—6k)+Q(3/k)(1-€k)(yk—1-Z)=1-€k+0(yk—1—Z)
(b) (yk—Z)(yk—1—Z)R = "—1—(311—1- Z)Q(yk—1)

4 THE COMPLEX QUASI—LAGUERRE ITERATION 57

—(yk - Z)Q QZU k—i Q(yk)Q (yrs—1)(yk—Z)(yk—1—Z)

+7701}. — z)(y1_1— 2)R1(yk_1,y1)= Ti -1+ 0(71 1— Z)
(C) (ilk - Z)(yk—1 - 310R = (1 — €k)(yk — Z)(3/k—1 — ZlR

= (1 — €1)(n — 1) + 001-1 — 2) (because of (b) above)
((1) (yr: — Z)2R = 6143/4 - Z)(yk—1 — Z)R

= 6k(n — 1) + O(yk_1 — 2) (because of (b) above)
(8) (yr: - Z)QR [(yk—I - 31021? + 4(n — 1)]

= [(y1. - Z)(yk—1 - 341:)ng + 4(" -1)(y1 - Z)QR

= (l—ek)2(n—1)2 + 4(n—1)€k(n—l) + O(yk_1—z)

= (n — 1)2(1+ 61)2 + O(yk—r - Z)

(because of (c) and (d) above).

If yo —) z and y1 —+ 2, then

 

 

 

 

n—l 1 n—l 1
= . —) E . < ,
Q<y0) 1:1 370 _ 2(2) Q {22:1 2 _ 2(1) 00
n—l 1 R n—l 1
R , = . —) * E . < .
“yo 9‘) I; (90 - 40(41- 2(2)) .2 (z - 2(7)? 00

In particular, if go is close to z and '61] < 1, then |y1— 2] < Iyo — z| and

|Q(yo)l < 2162] and |R1(3/anl)| < 2112*]-

In addition, if [6,] < 1 for i = 2, - - - ,k, then

lyk—Z|<“'<|3/1—Z|<|yo—Z|-

Hence,
lQ(yk)l < 2IQI and 1R1(i/k—1,yk)| < 2|R*|.

It follows that the bounds of O(yk_1—z) in all the above expressions (0)
through (e) are independent of k. Substituting (0), (c), (e) into (58)

4 THE COMPLEX QUASI-LAGUERRE ITERATION 58

yields

2(n— 1) +2ek+oz1
(n — 1)(1— ck) +26;c :l: (n — 1)(1 +61) +a2’

where 011 = O(yk_1— z) and a2 = O(yk_1 — 2); namely, 3C1 > 0 and
C2 > 0, independent of k, such that

 

€k+1 =1—

IOI1| _<_ Cllyk—l — 2| and |a2| S CQka—l — 2|-
Let

A = (n — 1)(1 — 6k) +2ek:l: (n — 1)(1 +61)
2(n — 1) + 26],. if ‘+’ is chosen;
(4 — 2n)ek if ‘—’ is chosen.
In (57), since
(yr. — 201931: A) = 444—02,
thus, making |yk+1 — ykl smaller is equivalent to making |A] larger.
When lekl < 1, |2(n — l) + 2ek| > [(4 — 2n)ek|. Hence, to have |A|

larger, we always choose ‘+’ for A, making

A=2(n—1)+2€k

 

 

and
|A| = |2(n — 1) + 2ek| Z 2(n — 1) — 2|ck| 2 2n — 4.
So,
A+ozl 612—01
:1. z . 59

€k+1 A—I-CYQ A-i-ag ( )

If |61| < 1, I62] < 1,~-,|ek| < 1 and lyo - Z] < "(32" then
|072| S C2|y1_1— 2] < Cglyo — 2| < 1.

Hence,

|A+072|Z |A|—|oz2| 22n—5>0.

4 THE COMPLEX QUASI-LAGUERRE ITERATION 59
From (59),

|€k+1| =

M<(§ifi|ly _,
|A+072| _ 272—5 [CA-1 A

In summary, we have shown the following:

Lemma 4.2 If yo is close to 2, then there exists a constant C > 0 such
that V14: 2 1, if |€1| < 1,---,|€k| < 1, then

|€k+1| _<_ Clyk—l — 3| (60)
Now, given a real number r < 1, when lgo — 2| < %, by (60) above,
|€2| S Clyo — 2| < r or lyg — 2| < rlyl — 2|.

Notice that from the way yl is chosen; namely, by applying one step of

Newton’s iteration on given yo, we always have
191 — 2| < |y0 — 2|, or equivalently |€1| < 1

when yo is close to 2, say for instance, yo is in the convergence neigh-

borhood of Newton’s iteration at 2. By (60) again,

I153] S Clyl — 2| < Clyo — 2| < r, and so, |y3 — 2| < r|y2 - 2|.
It follows recursively that

|ek+1| S Clyk_1— 2| < r and |yk+1— 2| < rlyk — 2|, We 21

since C is independent of 16. Therefore {yk} converges to 2. Further-

more, from (60),
|ek+1| S C|y1_1— 2| —> 0, as k —-> 00.

Thus the convergence of {ykfiozl to 2 is super-linear.
The rest of the proof follows exactly the same argument as the proof

of part (ii) of Theorem 2.3 with m = M = 1. D

REFERENCES 60

References

[1] E. ANDERSON, Z. 3.41, C. BISCHOF, J. DEMMEL, J. DONGARRA, J. Du
CRoz, A. GREENBAUM, S. HAMMARLING, A. MCKENNEY, S. OSTROU-
CHOV, and D. SORENSON, LAPACK User’s Guide, SIAM, Philadelphia, 1992.

[2] J. J. M. CUPPEN, A divide and conquer method for the symmetric tridiagonal
eigenproblem, Numer. Math., 36 (1981), pp. 177—195.

[3] J. J. DONGARRA AND D. C. SORENSEN, A fully parallel algorithm for the
symmetric eigenvalue problem, SIAM J. Sci. Stat. Comput., 8 (1987), pp. 139—
154.

[4] G. H. GOLUB AND C. F. VAN LOAN, Matrix Computations, 2nd Ed., The
Johns Hopkins University Press, Baltimore, MD, 1989.

[5] R. T. GREGORY AND D. L. KARNEY, A Collection of Matrices for Testing

Computational Algorithms, Robert E. Krieger Publishing Company, Hunting-
ton, New York, 1978.

[6] M. GU AND S. C. EISENSTAT, A divide-and-conquer algorithm for the sym-
metric tridiagonal eigenproblem, SIAM J. Matrix Anal. Appl., Vol. 16, No. 1
(1995), pp. 172-191.

[7] W. KAHAN, Notes On Laguerre’s Iteration, preprint, University of California,
Berkeley (1992).

[8] T. Y. Li AND Z. ZENG, Laguerre’s iteration in solving the symmetric tridiago-
nal eigenproblem — revisited, SIAM J. Sci. Comput., Vol. 15, No. 5 (1994), pp.
1145-1173.

[9] B. N. PARLETT, The Symmetric Eigenvalue Problem, Prentice-Hall, Englewood
Cliffs, NJ, 1980.

[10] B. N. PARLETT, The use ofa refined error bound when updating eigenvalues of
tridiagonals, Lin. Alg. & Appls., Vol. 68 (1985), pp. 179-219.

[11] J. RUTTER,A serial implementation of Cuppen’s divide and conquer algorithm
for the symmetric eigenvalue problem, LAPACK lawn 69 (1994).

REFERENCES 61

[12] D. C. SORENSEN AND P. T. P. TANG, On the orthogonality of eigenvectors
computed by divide-and-conquer techniques, SIAM. J. Numer. Anal., 28 (1991),
pp. 1752 - 1775.

[13] C. TREFFTZ, C. C. HUANG, P. MCKINLEY, T. Y. L1, AND Z. ZENG, A
scalable eigenvalue solver for symmetric tridiagonal matrices, to appear, Parallel

Comput.

[14] J. H. WILKINSON, The Algebraic Eigenvalue Problem, Oxford University Press,
Oxford, 1965.