Download Rational Approximation and Interpolation of Functions and more Exams Advanced Calculus in PDF only on Docsity!
Solving Secular Equations
Stably and Eciently
Ren-Cang Li
Department of Mathematics
University of California at Berkeley
Berkeley, California 94720
April, 1993
Dedicated to B. N. Parlett and W. Kahan on the occasion of their 60th birthdays
Abstract A divide-and-conquer metho d for solving symmetric tridiagonal eigenproblems has evolved from work by Cupp en, Dongarra, Sorensen, Tang, and most recently Gu and Eisenstadt. At the heart of their metho ds is the solution of a so-called Secular Equation. Prop osed here is a more ecient organization of the equation{solving pro cess, including some crucial implementation details.
1 Intro duction
A continuing element in divide-and-conquer algorithms descended from Cup- p en's [4, 5 ] is the solution of an eigensystem that di ers from diagonal by
a rank{1 p erturbation. Let D = diag ( 1 ; ; n ) b e a real diagonal matrix,
a nonzero real numb er, and z = ( 1 ; ; n )T^ a real vector; the p erturb ed
system's matrix is
D +
z z T^ : (1)
Its eigenvalues are the zeros x = 1 ; ; n of the secular function (see [7])
f (x) = +
X^ n
j =
(^2) j
j x
and the eigenvector corresp onding to each k is parallel to
(D k I ) ^1 z : (3)
(Here I is the n n identity matrix.) Numerical diculties arose when these
were not computed accurately enough to b e orthogonal, but recent work [11, 9] has overcome these diculties. Ours is a small contribution to solving the secular equation f (x) = 0 as accurately as necessary but faster than b efore. Certain details necessary for robustness in the co de will b e discussed to o. Inattention to such details can cause accidents like Division by Zero , which the author has encountered while testing others' programs. The pap er is organized as follows: In Section 2, we study various his- toric ways to rationally interp olate the secular function (2) and then develop three di erent schemes for solving the equation f (x) = 0, two of which are essentially due to [3] and the other one is new and fastest among the three. The way to interp olate f (x) in Section 2 is not always ecient b ecause, in which, attention is largely paid to the p ositions of two nearby p oles. Sec- tion 3 gives a closer lo ok into cases when attention has to b e paid not only to the p ositions of most relevant p oles but also to weights over particular p oles. Imp ortant implementation issues like securing initial guesses and selecting the b est scheme are discussed in detail in Sections 4 and 6. Numerical exam- ples with detailed explanation are given in Section 7. Discussions of various stopping criteria and justi cations of our prop osed stopping criteria, are pre- sented in Section 5. Numerical examples with detailed explanation are given in Section 7. Section 8 presents our conclusions.
They will b e used in conjunction with a partition of the secular function
f (x) = + k (x) + k (x) where, for k = 1 ; 2 ; ; or n,
k (x)^ =
X^ k
j =
(^) j^2
j x
; k (x) =
X^ n
j =k +
(^) j^2
j x
(Convention sets n (x) 0.) The choice of k will dep end up on which eigen-
value k is b eing computed. It is easy to see that for k < x < k +
1 < k (x) < 0 < k (x) < + 1 :
For any approximation y to k , we shall approximate each of (^) k (x) and k (x) by a simpler form F or G chosen to match (^) k and k in value and derivative at x = y ; in other words, we p erform osculatory interp olation. So far, nothing we have said di ers from past practice [3, 5, 9 ]. What will b e novel will b e the way we cho ose which of F and G should b e used with each of (^) k and k ; in fact, we shall see that F is b est not used at all.
2.1.1 Two Ways to Rationally Interp olate (^) k (x)
Let y b e a xed approximation to k somewhere b etween k and k +1. Now we have a choice. We can cho ose parameters p and q in F (x; p; q ) so that
F (y ; p; q ) = (^) k (y ); F 0 (y ; p; q ) = (^0) k (y ):
Or we can cho ose = k to match the p ole of (^) k (x) next to y (and k ) and then cho ose parameters r and s in G(y ; k ; r; s) so that
G(y ; k ; r; s) = (^) k (y ); G^0 (y ; k ; r; s) = (^) k^0 (y ):
Here are the formulas for the parameters:
In F (x; p; q ) = q =(p x),
p = y + (^) k (y )= (^0) k (y ) (6)
=
0 k (y^ )
X^ k
j =
(^) j^2
(j y )^2
j ;
q = (^) k (y )^2 = (^0) k (y ) > 0 : (7)
In G(x; k ; r; s) = r + s=(k x),
r = k (y ) (k y ) 0 k (y ) (8)
kX 1
j =
j k
(j y )^2
j^2 0 ;
s = (k y )^2 0 k (y ) > 0 : (9)
In fact, k^2 s
Pk
j =
(^) j^2 ; if k > 1,
kP 1
j =
j^2 =(j y ) = k 1 (y ) < r < 0; and if
k = 1, r = 0.
Let these assignments to p; q ; r; s hold throughout the rest of this x 2.1.1.
Here are some of the prop erties of the two osculatory interp olating functions that our subsequent analysis will exploit.
Prop osition 1 1 < p < k < y :
This prop osition indicates that the p ole p of F (x; p; q ) lies farther from y than k do es. Therefore F (x; p; q ) may approximate (^) k (x) p o orly b etween its p ole k and y. This happ ens when k is relatively small, in which case k lies very close to k , and therefore y will do likewise. Since the interp olating rationals are unique, as we can see from the ab ove equations, they are identically equal to (^) k (x) if k = 1.
Theorem 1 If k 2 , then for k < x < + 1 and x 6 = y
G(x; k ; r; s) < (^) k (x) < F (x; p; q ):
Proof: The inequality F (x; p; q ) > k (x) for k < x < + 1 was proved by
[3]. Now we prove G(x; k ; r; s) < k (x) for k < x < + 1 in a similar way.
Let g (x) = G(x; k ; r; s) k (x). It suces to prove g (x) < 0 for k <
x < + 1. To this end, set
h(x) = g (x)
Y^ k
j =
(j x): (10)
X^ n
j =k +
j k +
(j y )^2
(^) j^2 > 0 ; (17)
s = (k +1 y )^2 ^0 k (y ) > 0 : (18)
Moreover 2 k +1 s
Pn
j =k +1 ^ 2
j ;^ if^ k^ <^ n^ ^ 1,^0 <^ r^ <^ k^ +1^ ;^ and^ if^ k^ =^ n^ ^ 1,
r = 0.
F and G are identically equal to k (x) if k = n 1.
Prop osition 2 y < k +1 < p < n :
Prop osition 2 indicates the same phenomenon as we observed in Prop o- sition 1.
Theorem 2 If k n 2 , then for y 6 = x < k +
F (x; p; q ) < k (x) < G(x; k +1 ; r; s):
2.2 Solving f (x) = 0
Three ways come to mind to nd the k th zero k of the secular function (2) by combining di erent rational osculatory interp olating functions to (^) k (x) and k (x) studied ab ove. The rst one we are going to study is the one in [3], and we describ e it as Approaching from the Left b ecause the algorithm will pro duce a sequence of monotonically increasing approximations to k provided the initial guess is b etween k and k. We repro duce it here for completeness. We describ e the second one as Approaching from the Right for the same reason. The third iteration scheme is a new and fastest one which we shall call the Midd le Way. Assume, we have a guess y to k (with k < y < k +1 ); Later, in Section 4, we shall discuss how to cho ose such a y.
2.2.1 Approaching from the Left
Assume for the moment 1 k < n.
To approach from the left, we interp olate (^) k (x) by F (x; p; q ) and k (x) by G(x; k +1 ; r; s), where p, q , r and s are determined by Equations (6), (7), (16) and (18). Instead of solving + (^) k (x) + k (x) = 0, we solve
+
q
p x
s
k +1 x
for x. Equation (19) has two ro ots x of which just one lies b etween p and
k +1 ; that one is our new approximation y + . Set k +1 = k +1 y. Then
the correction to y is
a
p
a^2 4 bc
2 c
if a 0 ; (20)
2 b a +
p
a^2 4 bc
if a > 0 ; (21)
where
a = f (y ) k +1 + k^
(y ) 0 k (y^ )
k +1 k (y ) 1 +
^0 k (y ) 0 k (y^ )
;
b = k +1 f (y ) k (y^ ) 0 k (y^ )^
;
c = + r = + k (y ) k +1 ^0 k (y ):
Eventually, as f (y )! 0, we nd a > 0, so (21) will b e used more often than
The case for k = n is sp ecial since n 0. Solving + p qx = 0 gives
= x y =
+ (^) n (y ) (^) n^0 (y ) n (y^ )^ =^
f (y ) (^) n (y ) f 0 (y )
Regarding the foregoing scheme, we have
Theorem 3 If k < y < k , then > 0 and y < y + < k ; if k < y < k +1 , then < 0 and y + < k < y.
Proof: It is clear from Theorems 1 and 2. Theorem 3 says, starting with an initial guess which is less than k , the foregoing scheme wil l yield a sequence of approximations converging mono- tonical ly upwards to k ; on the other hand if the initial guess exceeds k the
b = n 1 n f (y );
c = + r = + n 1 (y ) n 1 n^0 1 (y )
= f (y ) n 1 n^0 1 (y )
(^) n^2 n
:
(23) is usable provided f (y ) n 1 0 n 1 (y ) ^
n^2 n >^0 since^ otherwise^ y^ +^ ^ <^ n
or is in nite. Eventually as f (y )! 0, so do es f (y ) n 1 n^0 1 (y ) ^
(^2) n n b ecome p ositive. The foregoing scheme is called Approaching from the Right b ecause of this:
Theorem 4 If k < y < k +1 , then < 0 and k < y + < y. If k < y < k and k < n, then > 0 and y < k < y + . If n < y < n and y is not too close to n , then y + > n.
Proof: It is clear from Theorems 1 and 2. Remark: As we remarked after Theorem 3, in Theorem 4 it is p ossible for y + to jump ab ove k +1 when k < y < k. This must b e prevented in practice.
2.2.3 The Middle Way: a New Metho d
A new way (we call it the Midd le Way ) to solve the secular equation f (x) = 0 will b e develop ed here. Unlike the previous two ways, it need not yield a sequence of approximations that converges monotonically to k , but it converges faster as we will see. We interp olate b oth (^) k (x) and k (x) by rationals of typ e G(x; ; r; s), taking b oth nearby p oles into consideration. To b e sp eci c, when 0 < k < n we let
r + (^) s
k ^ x^
approximate to (^) k (x); and
R + (^) S
k +1 ^ x^
approximate to k (x);
where
s = ^2 k^0 k (y ) > 0 ;
r = k (y ) k 0 k (y ) 0 ;
and
S = ^2 k +1 ^0 k (y ) > 0 ;
R = k (y ) k +1 ^0 k (y ) 0 ;
and k = k y < 0 < k +1 = k +1 y as b efore. We compute our new \b etter" approximation y + to k by solving the equation
+ r +
s k x
+ R +
S k +1 x
Equation (25) has two ro ots x, one of them in nite if + r + R = 0. The ro ot we need is the one b etween k and k +1. Set = x y. Then the correction to y has forms (20) and (21) with
a = (k + k +1 )f (y ) k k +1 f 0 (y ); b = k k +1 f (y ); c = f (y ) k (^0) k (y ) k +1 ^0 k (y ):
Eventually a > 0. The case k = n is exactly the same as in x 2.2.2. Remark: Quadratic convergence of Approaching from the left was proved in [3]. A similar pro of do es likewise for Approaching from the Right , and proves that the Midd le Way converges at least quadratically.
3 Secular Equation Solvers (I I): A Close Lo ok
The Midd le Way outp erforms Approaching from the Left and Approaching from the Right b ecause it takes b oth nearby p oles into considerations while each of latter two uses one of the p oles. However, as we shall see, the Midd le Way still b ehaves badly in a few circumstances. The following is a simple
intuitive explanation for why. In
(^2) j j x ,^ ^
2 j is^ the^ weight^ over^ the^ p^ ole^ j^ and controls how much role the p ole j plays in the secular function f (x). Most of the time the simple interp olating rational r + (^) k s x approximates (^) k (x) very well, but there are situations when s overestimates (^) k^2 so much that iterations
are forced to move slowly towards the desired ro ots. Recall s (^) k^2 k x functions for the rest of the p oles in (^) k (x). Similar things happ en to k , to o. In this section, we shall provide ways to conquer this as well as diculties when just employing two nearby p oles is not enough. The case k = n will not b e handled here since the Midd le Way (Ap- proaching from the Right also) has used the two nearby p oles on the left of n already while there is no p ole on its right. We b egin with a new lo ok at our iteration formulas whose exibility will b e of big help.
Proof: It follows from (27) and (28) that
s =
^2 k k k +
(f (y ) c k +1 f 0 (y ));
S =
^2 k + k +1 k
(f (y ) c k f 0 (y )):
Substituting them into c(k + k +1 ) + s + S and ck k +1 + sk +1 + S k leads to the desired results. Prop osition 3 illustrates a surprising fact that the iteration formula by solving (29) dep ends only up on c, alone. As surely, the equations (27) and (28) are solvable for any c. Without worrying ab out global convergence, any c will give rise an ultimately at least quadratically convergent iteration scheme! In the case s > 0 and S > 0 in which we are interested, it follows from (29) that
a
p a^2 4 bc 2 c
if a 0 ; (30)
2 b a +
p a^2 4 bc
if a > 0 ; (31)
where
a = (k + k +1 )f (y ) k k +1 f 0 (y ); b = k k +1 f (y ):
In what follows, for a particular iteration scheme, we shall provide the nec- essary numb er c, as well as a pro of for s > 0 and S > 0. Di erent choices of c give rise to di erent iteration schemes. By cho osing c prop erly, a very ecient iteration scheme may b e obtained. The Midd le Way falls into the category, i.e., s > 0 and S > 0, and
c = f (y ) k (^) k^0 (y ) k +1 ^0 k (y ) = f (y ) k +1 f 0 (y ) (^0) k (y )(k k +1 ) = f (y ) k f 0 (y ) ^0 k (y )(k +1 k ):
Reorganization of c is for the reader to see its di erences from those in some other schemes to which we are ab out to get.
3.2 The Fixed Weight Metho d
We have b een aware that the weight (^) k^2 or (^) k^2 +1 may b e overestimated in the Midd le Way. Although not all overestimations are harmful, there are cases where iterations move slowly b ecause of overestimations. In order to handle these cases eciently, we prop ose the following scheme so-called the Fixed Weight Method b ecause it xes one of the weights (^) k^2 and (^) k^2 +1 while satisfying (27) and (28).
- The Case k closer to k : We set s = (^) k^2 , then
s = (^) k^2 ; (32)
S = ^2 k +1 f 0 (y )
(^) k^2 ^2 k
= (^) k^2 +1 +
X
j 6 =k ;k +
^2 k + ^2 j
(^) j^2 > (^) k^2 +1 ;
c = f (y )
(^) k^2 k
k +1 f 0 (y )
(^) k^2 ^2 k
= f (y ) k +1 f 0 (y )
(^) k^2 ^2 k
(k k +1 ): (34)
- The Case k closer to k +1 : We set S = (^) k^2 +1 , then
s = ^2 k f 0 (y )
(^) k^2 + ^2 k +
= (^) k^2 +
X
j 6 =k ;k +
^2 k ^2 j
(^) j^2 > (^) k^2 ;
S = (^) k^2 +1 ; (36)
c = f (y ) k f 0 (y )
(^) k^2 + ^2 k +
(k +1 k ): (37)
The following theorem gives a basic prop erty of the interp olating function Q(x; c; s; S ) and the next approximation y + to k.
Theorem 5 If c, s and S are de ned by (32), (33) and (34), then f (x)
Q(x; c; s; S ) for k < x < k +1 and either k < y < y + < k < k +1 or
3.4 Using Three Poles When Necessary: a Hybrid
Scheme
The Midd le Way and the Fixed Weight Method can b e combined to design more p owerful secular equation solvers by prop erly switching b etween the two. But there are cases where three p oles have to b e used in order to make
iterations go faster. For m = 2 ; ; n 1, set
fm (x) = +
X^ n
j =1;j 6 =m
(^) j^2
j x
;
which is the secular function f (x) with the mth term in the summation removed. It is easy to see that fm (x) has a zero b etween m 1 and m+. Numerically, we have discovered the following cases are p ossible dicult ones:
- k < k < k^ + 2 k +1 and fk (x) has a zero b etween k and k ;
- k^ + 2 k^ +1< k < k +1 and fk +1 (x) has a zero b etween k and k +.
To see why the rst case may b e dicult, we let k^2! 0, then, as functions
of (^) k^2 , k 1 go es monotonically upward until it hits k while k go es mono- tonically downward until it hits the zero of fk (x) b etween k and k. Now, on the contrary, let (^) k^2 go back from zero to its original value, what we will see is the exact contrary phenomena. Based on this simple observation, we can think, roughly, that k 1 dep ends largely on the p ole k and its weight (^) k^2 and, therefore, the Fixed Weight Method should b e go o d at nding it. On the other hand, k dep ends on the p ole k and its weight (^) k^2 and the zero of fk (x) where it starts from and which is controlled roughly by the the p oles k 1 and k +1 with appropriate weights. Similar intuitive arguments apply to the second case ab ove. The treatment of the two cases is almost identical. In what follows, we discuss the rst case only. A natural way to handle the rst case is to interp olate f (x) with the following simple rational:
Q^ e(x; c; s; S ) def = c + s
k 1 x
(^) k^2
k x
S
k +1 x
The parameters c, s and S can b e determined by either interp olating fk (x) = + (^) k 1 (x) + k (x) in the way that the Midd le Way do es or in the way that
the Fixed Weight Method do es dep ending on situations. Once we have (39), its zero b etween k and k +1 can b e computed in various iterative ways with
negligible cost. However, care has to b e taken in evaluating Qe(x; c; s; S ) at
a given p oint. As a matter of fact, b ecause of roundo , sometimes (though
rarely) computed Qe(y ; c; s; S ) di ers signi cantly from computed f (y ) though
the two should b e the same by interp olation in theory. It turns out that we
can evaluate Qe(x; c; s; S ) in an indirectly way to avoid this from happ ening.
Since cimputed f (y ) is available at the time we interp olate the secular func-
tion f (x), we do not compute Qe(y ; c; s; S ) at all while simply setting it to
b e f (y ). At the very next time when we need to compute Qe(x; c; s; S ) we
up date f (y ) by adding a correction to it b ecause
Q^ e(x; c; s; S ) = Qe(y ; c; s; S ) + ( Qe(x; c; s; S ) Qe(y ; c; s; S ))
= f (y ) + (x y )
s
(k 1 y ) (x y )
(^) k^2
(k y ) (x y )
S
(k +1 y ) (x y )
:
Subsequent evaluation of Qe(x; c; s; S ) is done in the same way by adding
correction to its value at the previous x. According to our exp erience, prop er treatment of p oles and their weights is crucial and delicate, esp ecially in dicult cases as we discussed ab ove. It accelerates convergence signi cantly in the sense that it reduces the over- all average numb er of iterations p er eigenvalue nding and keeps the p eak numb er of iterations for nding an eigenvalue reasonably small. which is essential for avoiding load-balancing in parallel computations. In view of this, we prop ose the following scheme{the Hybrid Scheme which combines the Midd le Way and the Fixed Weight Method in a clever way and which of course employs three p oles when necessary: Supp ose we are computing k
and supp ose k < k < (k + k +1 )=2. In practice, we compute k k instead
of k itself. We have an initial guess k + for k for which k < k + < k is guaranteed (ref. Section 4). The value of the secular function is evaluated in such a way
f (k + ) = fk (k + ) +
(^) k^2
other j 's, in which case k is very close to k or k +1. In extreme cases, as we have already observed, the iterations Approaching from the left or Right can even jump out from b etween k and k +. In what follows we will present an inexp ensive way to obtain relatively accurate initial guesses. At the same time we shall show how to decide which of k and k +1 is closes to k. This decision is imp ortant b ecause, if decided wrongly, roundo could cause an iterate to collide with an endp oint k or k +1 , or even oversho ot it. A collision would so on lead to Division by Zero; oversho oting could jeopardize subsequent convergence. To avoid these problems, we shall translate the origin temp orarily, while k is b eing computed, to whichever of k or k +1 is closer. A natural and e ective way to make the correct decision is to lo ok at the sign of f
k +k + 2
. If this value is p ositive, we know k is closer to k than to k +1 and thus the origin should b e translated to k ; otherwise k is closer to k +1 and the origin should b e translated to there. If, however, roundo
obscures the sign of f
k +k + 2
, in which case k is almost half way b etween the two p oles, the origin could b e translated to either one of them without making much di erence. The computation of f
k +k + 2
can b e done in such a way that an initial guess y is obtained as a by-pro duct.
First we consider the case 1 k < n.
Rewrite the secular function (2) as f (x) = g (x) + h(x), where
g (x) = +
X^ n
j =1;j 6 =k ;k +
(^) j^2
j x
; and h(x) =
(^) k^2
k x
(^) k^2 +
k +1 x
We cho ose our initial guess y to b e that one of the two ro ots of the equation
g
k + k + 2
lying b etween k and k +1 , where g
k +k + 2
was retained from the compu-
tation of f
k +k + 2
. (Thus it is a gift!) By sketching the graph of the sum
of the last two terms on the left hand-side of (41), we can tell without any
diculty which ro ot is needed. In case f
k +k + 2
0, Equation (41) should
b e solved for = y k , while in case f
k +k + 2
< 0, it should b e solved
for = y k +1. De ne = k +1 k and c = g
k +k + 2
. Here are the
formulas for :
= y K =
a
p
a^2 4 bc
2 c
if a 0 ; (42)
2 b a +
p
a^2 4 bc
if a > 0 ;
where if f
k +k + 2
K = k ; a = c + ( (^) k^2 + (^) k^2 +1 ); b = (^) k^2 ; (43)
and if f
k +k + 2
K = k + 1 ; a = c + ( k^2 + k^2 +1 ); b = k^2 +1 : (44)
Theorem 6 If f
k +k + 2
0 , then given by (42) and (43) satis es
k < k + < k <
k + k + 2
If f
k +k + 2
< 0 , then given by (42) and (44) satis es
k + k + 2
< k < k +1 + < k +1 :
Proof: One way to see these is to sketch graphs. It can also b e seen by subtracting one from the other of the following two equations
(^) k^2
k k
(^) k^2 +
k +1 k
= g (k );
(^) k^2
k y
(^) k^2 +
k +1 y
= g
k + k + 2
;
which pro duces
(k y )
(^) k^2
(k k )(k y )
(^) k^2 +
(k +1 k )(k +1 y )
k + k + 2
k
X
j 6 =k ;k +
(^2) j
j k^ + 2 k^ +
(j k )
