Accurate eigenvalue decomposition of arrowhead matrices and applications
aa r X i v : . [ m a t h . NA ] A p r Accurate eigenvalue decomposition of arrowheadmatrices and applications
N. Jakovˇcevi´c Stor a,1, ∗ , I. Slapniˇcar a,1 , J. Barlow b,2 a Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture,University of Split, Rudjera Boˇskovi´ca 32, 21000 Split, Croatia b Department of Computer Science and Engineering, The Pennsylvania State University,University Park, PA 16802-6822, USA
Abstract
We present a new algorithm for solving an eigenvalue problem for a real sym-metric arrowhead matrix. The algorithm computes all eigenvalues and all com-ponents of the corresponding eigenvectors with high relative accuracy in O ( n )operations. The algorithm is based on a shift-and-invert approach. Doubleprecision is eventually needed to compute only one element of the inverse ofthe shifted matrix. Each eigenvalue and the corresponding eigenvector can becomputed separately, which makes the algorithm adaptable for parallel com-puting. Our results extend to Hermitian arrowhead matrices, real symmetricdiagonal-plus-rank-one matrices and singular value decomposition of real trian-gular arrowhead matrices. Keywords: eigenvalue decomposition, arrowhead matrix, high relativeaccuracy, singular value decomposition
1. Introduction and Preliminaries
In this paper we consider eigenvalue problem for a real symmetric matrix A which is zero except for its main diagonal and one row and column. Sinceeigenvalues are invariant under similarity transformations, we can symmetricallypermute the rows and the columns of the given matrix. Therefore, we assumewithout loss of generality that the matrix A is a n × n real symmetric arrowheadmatrix of the form ∗ Corresponding author
Email addresses: [email protected] (N. Jakovˇcevi´c Stor), [email protected] (I.Slapniˇcar), [email protected] (J. Barlow) The research of Ivan Slapniˇcar and Nevena Jakovˇcevi´c Stor was supported by the Ministryof Science, Education and Sports of the Republic of Croatia under grant 023-0372783-1289. The research of Jesse L. Barlow was supported by the National Science Foundation undergrant CCF-1115704.
Preprint submitted to Elsevier April 12, 2013 = (cid:20) D zz T α (cid:21) , (1)where D = diag( d , d , . . . , d n − )is diagonal matrix of order n − z = (cid:2) ζ ζ · · · ζ n − (cid:3) T (2)is a vector and α is a scalar.Such matrices arise in the description of radiationless transitions in isolatedmolecules [3], oscillators vibrationally coupled with a Fermi liquid [8], quantumoptics [15] (see also Example 4). Such matrices also arise in solving symmetricreal tridiagonal eigenvalue problems with the divide-and-conquer method [11].In this paper we present an algorithm which computes all eigenvalues andall components of the corresponding eigenvectors with high relative accuracy in O ( n ) operations.Without loss of generality we may assume that A is irreducible, that is, ζ i = 0 , for all i and d i = d j , for all i = j, i, j = 1 , . . . , n − . If A has a zero in the last column, say ζ i = 0, then the diagonal element d i is an eigenvalue whose corresponding eigenvector is the i -th unit vector, andwe can reduce the size of the problem by deleting the i -th row and column ofthe matrix, eventually obtaining a matrix for which all elements ζ j are nonzero.If d i = d j , then d i is eigenvalue of matrix A (this follows from the interlacingproperty (7)), and we can reduce the size of the problem by annihilating ζ j witha Givens rotation in the ( i, j )-plane and proceeding as in the previous case.Further, by symmetric row and column pivoting, we can order elements of D such that d > d > · · · > d n − . (3)Hence, we will consider only ordered and irreducible arrowhead matrices. With-out loss of generality we can also assume that ζ i > i , which can be at-tained by pre- and post-multiplication of the matrix A with D = diag(sign( ζ i ))).Let A = V Λ V T (4)be the eigenvalue decomposition of A . HereΛ = diag( λ , λ , . . . , λ n )is a diagonal matrix whose diagonal elements are the eigenvalues of A , and V = (cid:2) v · · · v n (cid:3)
2s an orthonormal matrix whose columns are the corresponding eigenvectors.The eigenvalues of A are the zeros of the Pick function (see [4, 16]) f ( λ ) = α − λ − n − X i =1 ζ i d i − λ = α − λ − z T ( D − λI ) − z, (5)and the corresponding eigenvectors are given by v i = x i k x i k , x i = (cid:20) ( D − λ i I ) − z − (cid:21) , i = 1 , . . . , n. (6)Diagonal elements of the matrix D , d i , are called poles of the function f .Notice that (3) and the Cauchy interlacing theorem [10, Theorem 8.1.7]applied to matrices D and A imply the interlacing property λ > d > λ > d > · · · > d n − > λ n − > d n − > λ n . (7)Since A is symmetric, its eigenvalues may be computed by invoking anyof a number of standard programs (LAPACK [1]). However, these programsusually begin with an initial reduction of the matrix to tridiagonal form [17], oras proposed in [16], with an alternative which takes advantage of the structureof A by finding the zeros of the Pick function given in (5), for the eigenvaluesof A . This results in an algorithm which requires only O ( n ) computationsand O ( n ) storage. Although the idea is conceptually simple and in fact hasbeen used to solve other eigenvalue problems of special structure [2, 5, 6, 7],the computation is not always stable [11]. Namely, if the computed eigenvalues λ i are not accurate enough, then the computed eigenvectors v i may not besufficiently orthogonal (see Example 3). The existing algorithms for arrowheadmatrices [11, 16] obtain orthogonal eigenvectors with the following procedure:- compute the eigenvalues ˜ λ i of A by solving (5);- construct a new matrix ˜ A = (cid:20) D ˜ z ˜ z T ˜ α (cid:21) by solving inverse problem with the prescribed eigenvalues ˜ λ , and diagonalmatrix D , that is, compute new ˜ z and ˜ α as˜ ζ i = vuuut(cid:16) d i − ˜ λ n (cid:17) (cid:16) ˜ λ − d i (cid:17) i Y j =2 (cid:16) ˜ λ j − d i (cid:17) ( d j − − d i ) n − Y j = i +1 (cid:16) ˜ λ j − d i (cid:17) ( d j − d i ) ,˜ α = ˜ λ n + n − X j =1 (cid:16) ˜ λ j − d j (cid:17) .- compute eigenvectors of ˜ A by (6).3ince the formulas for ˜ ζ i involve only multiplications, division and subtrac-tions of exact quantities, each ζ i is computed with relative error of O ( ε M ),where ε M denotes the machine precision. Therefore, ˜ A = A + δA , where k δA k = O ( ǫ M ). Here k · k denotes the spectral matrix norm. We concludethat the computed eigenvalues ˜ λ i satisfy standard perturbation bounds likethose from [10, Corollary 8.1.6]. Further, since ˜ λ i are the eigenvalues of thematrix ˜ A computed to higher relative accuracy, the eigenvectors computed by(6) are orthogonal to machine precision. For details see [11, 16].Our algorithm uses a different approach. Accuracy of the eigenvectors andtheir orthogonality follows from high relative accuracy of the computed eigen-values and there is no need for follow-up orthogonalization. The algorithm isbased on shift-and-invert technique. Basically, the eigenvalue λ is computed asthe largest or the smallest eigenvalue of the inverse of the matrix shifted to thepole d i which is nearest to λ , that is, λ = 1 ν + d i , (8)where ν is either smallest or largest eigenvalue of the matrix A − i ≡ ( A − d i I ) − . Inverses of arrowhead matrices are structured in the following manner (here × stands for non-zero element): the inverse of an arrowhead matrix with zeroon the shaft is a permuted arrowhead matrix with zero on the shaft, × ×× × ×× ×× × × × × − = × ×× ×× × × × ×× ×× , and the inverse of the full arrowhead matrix is a diagonal-plus-rank-one (DPR1)matrix, × ×× ×× ×× ×× × × × × − = × × × × ± uu T . The machine precision ε M is defined as a smallest positive number such that in thefloating-point arithmetic 1 + ε M = 1. In Matlab or FORTRAN REAL(8) arithmetic ε M =2 . · − , thus the floating-point numbers have approximately 16 significant decimaldigits. The term “double of the working precision” means that the computations are performedwith numbers having approximately 32 significant decimal digits, or with the machine precisionequal to ε M . aheig (Ar-rowHead EIGenvalues). In Section 3 we discuss the accuracy of the algorithm.In Section 4 we present the complete algorithm which uses double of the workingprecision, if necessary. In Section 5 we illustrate algorithm with few examplesand in Section 6 we apply our results to eigenvalue decomposition of Hermitianarrowhead matrix, singular value decomposition of real triangular arrowheadmatrix and eigenvalue decomposition of real symmetric diagonal-plus-rank-onematrix. The proofs are given in Appendix A.
2. Basic shift-and-invert algorithm
Let λ be an eigenvalue of A , let v be its eigenvector, and let x be theunnormalized version of v from (6). Let d i be the pole which is closest to λ .Clearly, from (7) it follows that either λ = λ i or λ = λ i +1 . Let A i be the shiftedmatrix A i = A − d i I = D z ζ i D z z T ζ i z T a , (9)where D = diag( d − d i , . . . , d i − − d i ) ,D = diag( d i +1 − d i , . . . , d n − − d i ) ,z = (cid:2) ζ ζ · · · ζ i − (cid:3) T ,z = (cid:2) ζ i +1 ζ i +2 · · · ζ n − (cid:3) T ,a = α − d i . Notice that D ( D ) is positive (negative) definite.Obviously, if λ is an eigenvalue of A , then µ = λ − d i is an eigenvalue of A i , and vice versa, and they both have the same eigenvector.The inverse of A i is A − i = D − w w T b w T /ζ i w D −
00 1 /ζ i , (10)5here w = − D − z ζ i ,w = − D − z ζ i ,b = 1 ζ i (cid:0) − a + z T D − z + z T D − z (cid:1) . (11)Notice that b = ¯ f ( d i ) /ζ i where ¯ f ( d i ) = α − d i − ¯ z T (cid:0) ¯ D − d i I (cid:1) − ¯ z where ¯ D is the diagonal matrix D without d i and ¯ z is z without ζ i .The eigenvector x from (6) is given by x = x ... x n = ( D − µI ) − z − ζ i µ ( D − µI ) − z − . (12)If λ is an eigenvalue of A which is closest to the pole d i , then µ is theeigenvalue of matrix A i which is closest to zero and ν = 1 µ = ± (cid:13)(cid:13) A − i (cid:13)(cid:13) . In this case, if all entries of A − i are computed with high relative accuracy,then, according to standard perturbation theory, ν is computed to high relativeaccuracy (by any reasonable algorithm). In Section 3 we show that all entriesof A − i are indeed computed to high relative accuracy, except possibly b (see(11)). If b is not computed to high relative accuracy and it influences (cid:13)(cid:13) A − i (cid:13)(cid:13) ,it is sufficient to compute it in double of the working precision (see Section 4).Further, if µ is not the eigenvalue of A i which is closest to zero, then | ν | < (cid:13)(cid:13) A − i (cid:13)(cid:13) , and the quantity K ν = (cid:13)(cid:13) A − i (cid:13)(cid:13) | ν | (13)tells us how far is ν from the absolutely largest eigenvalue of A − i . If K ν ≫ µ will be computed with high relative accuracy. Remedies of this situation aredescribed in Remark 3.With this approach the componentwise high relative accuracy of the eigen-vectors computed by (12) follows from high relative accuracy of the computedeigenvalues (see Theorem 3). Componentwise high relative accuracy of the com-puted eigenvectors implies, in turn, their orthogonality.6he described procedure is implemented in algorithm aheig basic (Algo-rithm 1). The computation of the inverse of the shifted matrix, A − i , accordingto formulas (10) and (11), is implemented in Algorithm 2. Algorithm 3 com-putes the largest or the smallest zero of the Pick function (5) by bisection.Given eigenvalue λ , Algorithm 4 computes the corresponding eigenvector by (6)or (12), respectively.
3. Accuracy of the algorithm
We now consider numerical properties of Algorithms 1, 2, 3, and 4. Weassume tha standard model of floating point arithmetic where subtraction ispreformed with guard digit, such that [9, 18, 10, 19] f l ( a ◦ b ) = ( a ◦ b )(1 + ε ◦ ) , | ε ◦ | ≤ ε M , ◦ ∈ { + , − , ∗ , / } , where ε M is machine precision. In the statements of the theorems and theirproofs we shall use the standard first order approximations, that is, we neglectthe terms of order O ( ε M ) or higher. Moreover, we assume that neither overflowor underflow occurs during the computation.We shall use the following notation:Matrix Exact eigenvalue Computed eigenvalue A λ e λA i µ − e A i = f l ( A i ) b µ e µ = f l ( b µ ) A − i ν − ^ ( A − i ) = f l ( A − i ) b ν e ν = f l ( b ν ) (14)Here e A i = f l ( A i ) = D ( I + E ) 0 0 z ζ i D ( I + E ) z z T ζ i z T a (1 + ε a ) , where E and E are diagonal matrices whose elements are bounded by ε M inabsolute values and | ε a | ≤ ε M .Further we define the quantities κ λ , κ µ and κ b as follows: e λ = f l ( λ ) = λ (1 + κ λ ε M ) , (15) e µ = f l ( µ ) = µ (1 + κ µ ε M ) , (16) e b = f l ( b ) = b (1 + κ b ε M ) . (17)We also define the quantity K b = | a | + (cid:12)(cid:12) z T D − z (cid:12)(cid:12) + (cid:12)(cid:12) z T D − z (cid:12)(cid:12)(cid:12)(cid:12) − a + z T D − z + z T D − z (cid:12)(cid:12) . (18)7 lgorithm 1 [ λ, v ] = aheig basic ( D, z, α, k ) % Computes the k -th eigenpair of an irreducible arrowhead matrix% A = [diag ( D ) z ; z ′ α ] n = max ( size ( D )) + 1 % Determine the shift σ , the shift index i , and whether λ is on the left% or the right side of the nearest pole.% Exterior eigenvalues ( k = 1 or k = n ): if k == 1 σ = d i = 1 side = ′ R ′ elseif k == nσ = d n − i = n − side = ′ L ′ else % Interior eigenvalues ( k ∈ { , . . . , n − } ): Dtemp = D − d k atemp = α − d k middle = Dtemp k − / F middle = atemp − middle − P ( z ./ ( Dtemp − middle )) if F middle < σ = d k i = kside = ′ R ′ else σ = d k − i = k − side = ′ L ′ endend % Compute the inverse of the shifted matrix, A − i [ invD , invD , w , w , w ζ , b ] = invA ( D, z, α, i ) % Compute the leftmost or the rightmost eigenvalue of A − i ν = bisect ([ invD ; 0; invD ] , [ w ; w ζ ; w ] , b, side ) % Compute the corresponding eigenvector µ = 1 /νv = vect ( D − σ, z, µ ) % Shift the eigenvalue back λ = µ + σ lgorithm 2 [ invD , invD , w , w , w ζ , b ] = invA ( D, z, α, i ) % Computes the inverse of an arrowhead matrix A = [diag( D − d i ) z ; z ′ α − d i ] % according to (10) and (11). n = max ( size ( D )) + 1 D = D − d i a = α − d i w = − z i − ./D i − /z i w = − z i +1: n − ./D i +1: n − /z i w ζ = 1 /z i invD = 1 ./D i − invD = 1 ./D i +1: n − b = ( − a + sum ( z i − . ^ ./D i − ) + sum ( z i +1: n − . ^ ./D i +1: n − )) /z i ^ Algorithm 3 λ = bisect ( D, z, α, side ) % Computes the leftmost (for side =’L’) or the rightmost (for side =’R’) eigenvalue% of an arrowhead matrix A = [diag ( D ) z ; z ′ α ] by bisection. n = max ( size ( D )) + 1 % Determine the starting interval for bisection, [ lef t, right ] if side == ′ L ′ lef t = min { D − | z | , α − k z k } right = min d i else right = max { D + | z | , α + k z k } lef t = max d i end % Bisection middle = ( lef t + right ) / while ( right − lef t ) /abs ( middle ) > ∗ epsF middle = α − middle − sum ( z. ^ ./ ( D − middle )) if F middle > lef t = middle else right = middle end middle = ( lef t + right ) / end % Eigenvalue λ = right lgorithm 4 v = vect ( D, z, λ ) % Computes the eigenvector of an arrowhead matrix A = [diag( D ) z ; z ′ α ] % which corresponds to the eigenvalue λ by using (6). v = [ z./ ( D − λ ); − v = v/ k v k λ and µ Let λ = µ + d i be an eigenvalue of the matrix A , where µ is the corresponding eigenvalue ofthe shifted matrix A i = A − d i from which λ is computed. Let e λ = f l ( e µ + d i )be the computed eigenvalue. Theorem 1 gives us dependency of accuracy of e λ in (15) upon accuracy of e µ in (16). Theorem 1.
For λ and e λ from (15) and µ and e µ from (16) we have | κ λ | ≤ | d i | + | µ || λ | ( | κ µ | + 1) . (19)Proofs of this theorem and subsequent theorems are given in Appendix A.From Theorem 1 we see that the accuracy of e λ depends on κ µ and the sizeof the quotient | d i | + | µ || λ | . (20)Theorem 2 analyzes the quotient (20) with respect to the position of λ and signsof µ and the neighboring poles. Theorem 2.
Let the assumptions of Theorem 1 hold.(i) If (see Figure 1 (i)) sign ( d i ) = sign ( µ ) , then | d i | + | µ || λ | = 1 . (ii) If λ is between two poles of the same sign and sign ( d i ) = sign ( µ ) (seeFigure 1 (ii)), then | d i | + | µ || λ | ≤ . Theorem 2 does not cover the following cases:10 d i d i +1 λ µ z}|{ s d i − d i λ µ z}|{ (i) (ii) Figure 1: Typical situations from Theorem 2 (a) If d <
0, then µ >
0. If, further, | d | ≈ | µ | , then λ is near zero, and( | d | + | µ | ) / | λ | ≫ (b) If d n > , then µ <
0. If, further, | d n | ≈ | µ | , then λ n is near zero, andagain ( | d n | + | µ | ) / | λ n | ≫ (c) If λ is between two poles of the different signs and sign ( d i ) = sign ( µ ), theneither d i +1 < < d i and µ <
0, or d i < < d i − and µ >
0. In both cases,if, additionally, | d i | ≈ | µ | , then λ is near zero, and ( | d i | + | µ | ) / | λ | ≫ s d λ µ z }| { s λ d i d i +1 µ z }| { (a) (c) Figure 2: Typical situations for special cases
Since only one of these three cases can occur, Theorems 1 and 2 imply thatfor all eigenvalues λ ∈ σ ( A ), but eventually one, it holds | d i | + | µ || λ | ≤ . If one of the above cases does occur, remedies are given in the followingremark.
Remark 1.
If one of the cases (a), (b) or (c) occurs, then λ is an eigenvalueof A nearest to zero, and we can accurately compute it from the inverse of A .Notice that the inverse is of an unreduced arrowhead matrix with non-zero shaft In this case λ is computed as a difference of two close quantities and cancellation canoccur. In this case λ n is computed as a difference of two close quantities and cancellation canoccur.
11s a diagonal-plus-rank-one (DPR1) matrix of the form A − = (cid:20) D − (cid:21) + ρuu T , where u = (cid:2) z T D − − (cid:3) T , ρ = 1 a − z T D − z . Eigenvalues of A − are zeros of (see [2, 14]) ϕ ( λ ) = 1 + ρ n X j =1 u j d j − λ . Since the absolutely largest eigenvalue of A − is computed accurately accordingto standard perturbation theory, and 1 / | λ | = k A − k , λ is also computed withhigh relative accuracy. In computing matrix A − , eventually ρ needs to be com-puted in higher precision. For more details see Remark 3. If the denominatorin ρ is computed as zero, the matrix A is numerically singular and we can set λ = 0. Notice that all components of the the corresponding eigenvector are stillcomputed accurately. Remark 2.
Notice that Algorithm 1 (and, consequently, Algorithm 5 below)can be easily modified to return both quantities, d i and µ such that λ = d i + µ .If none of the remedies from Remark 1 were needed, these two quantities giveadditional information about λ (that is, they give a more accurate representationof λ ). An example is given in Example 2.We still need to bound the quantity κ µ from (19). This quantity essentiallydepends on the accuracy of f l ( b ). The bound for κ µ is given in Theorem 6. Since the eigenvector is computed by (12), its accuracy depends on the ac-curacy of e µ as described by the following theorem: Theorem 3.
Let (16) hold and let e x = e x ... e x n = f l ( ( D ( I + E ) − e µI ) − z − ζ i e µ ( D ( I + E ) − e µI ) − z − ) (21) be the computed un-normalized eigenvector corresponding to µ and λ . Then e x j = x j (cid:0) ε x j (cid:1) , (cid:12)(cid:12) ε x j (cid:12)(cid:12) ≤ | κ µ | + 3) ε M , j = 1 , . . . , n. In other words, if κ µ is small, then all components of the eigenvector arecomputed to high relative accuracy. Since the accuracy of e λ and e x dependson the accuracy of e µ (on the size of κ µ ) in the next three subsections tells wediscuss the accuracy of e µ . Since e µ is computed as an inverse of the eigenvalueof the matrix f l ( A − i ), we first discuss the accuracy of that matrix.12 .3. Accuracy of the matrix A − i We have the following theorem:
Theorem 4.
For the computed elements of the matrix A − i from (10) and (11)for all ( j, k ) = ( i, i ) we have ^ (cid:0) A − i (cid:1) jk = f l (cid:0) A − i (cid:1) jk = (cid:0) A − i (cid:1) jk (1 + ε jk ) , | ε jk | ≤ ε M . For the computed element b ≡ (cid:0) A − i (cid:1) ii from (17) we have | κ b | ≤ ( n + 3) K b , where K b is defined by (18). The above theorem states that all elements of the matrix A − i are computedwith high relative accuracy except possibly b . Therefore, we have to monitorwhether b is computed accurately, and, if not, it needs to be computed in doubleof the working precision (see Section 4 for details). Let λ max be the absolutely largest eigenvalue of a symmetric arrowheadmatrix A , an let e λ max be the eigenvalue computed by bisection as implementedin Algorithm 3. The error bound from [16, Section 3.1] immediately impliesthat (cid:12)(cid:12)(cid:12)e λ max − λ max (cid:12)(cid:12)(cid:12) | λ max | = κ bis ε M , κ bis ≤ . n (cid:0) √ n + 1 (cid:1) . (22)Notice that the similar error bound holds for all eigenvalues which are of thesame order of magnitude as | λ max | . A − i The desired interior eigenvalue and, in some cases, also absolutely smallerexterior eigenvalue λ of A is in Algorithm 1 computed by (8), where ν is one ofthe exterior eigenvalues of the matrix A − i .The following theorem covers the case when ν is the absolutely largest eigen-value of (cid:13)(cid:13) A − i (cid:13)(cid:13) , and gives two different bounds. Theorem 5.
Let A − i be defined by (10) and let ν be its eigenvalue such that | ν | = (cid:13)(cid:13) A − i (cid:13)(cid:13) . (23) Let b ν be the exact eigenvalue of the computed matrix ^ (cid:0) A − i (cid:1) = f l (cid:0) A − i (cid:1) . Let b ν = ν (1 + κ ν ε M ) . (24) Then | κ ν | ≤ min (cid:26) ( n + 3) √ nK b , √ n + ( n + 3) (cid:0) | ζ i | n − X k =1 k = i | ζ k | (cid:1)(cid:27) , (25) where K b is defined by (18). .6. Final error bounds All previous error bounds are summarized as follows.
Theorem 6.
Let e λ be the computed eigenvalue of an unreduced arrowhead ma-trix A , let e µ be computed eigenvalue of the matrix f A i from (9), and let e ν be thecorresponding computed eigenvalue of the matrix ^ (cid:0) A − i (cid:1) from (10). If µ is theeigenvalue of A i closest to zero (or, equivalently, if (23) holds), then the errorin the computed eigenvalue e λ is given by (15) with | κ λ | ≤ | κ ν | + κ bis ) + 4 , (26) and the error in the computed un-normalized eigenvector e x is given by Theorem3 with | κ µ | ≤ | κ ν | + κ bis + 1 , (27) where | κ ν | is bounded by (25) and κ bis is defined by (22). Since we are essentially using the shift-and-invert technique, we can guaran-tee high relative accuracy of the computed eigenvalue and high componentwiserelative accuracy of the computed eigenvector if ν is such that | ν | = O ( k A − i k )and it is computed accurately. This is certainly fulfilled if the following condi-tions are met: C1. The quantity K ν from (13) is moderate, andC2. (i) either the quantity K b from (18) is small, or(ii) the quantity | ζ i | n − P k =1 k = i | ζ k | from (25) is of order O ( n ) . The condition C1 implies that ν will be computed accurately according tothe standard perturbation theory. The conditions C2 (i) or C2 (ii) imply that κ ν from (25) is small, which, together with C1 , implies that ν is computedaccurately.If the condition C1 does not hold, that is, if K ν ≫
1, remedies are givenin Remark 2 below. If neither of the conditions
C2 (i) and
C2 (ii) holds, theremedy is to compute b in double of the working precision as described in Section4. Remark 3.
We have two possibilities:(a) we can compute λ by shifting to another neighboring pole provided that K ν is in this case small (shifting to the pole d i − instead of d i in Figure3 (a)),(b) if shifting to another neighboring pole is not possible ( K ν ≫
1, see Figure3 (b)), we can invert A − σI , where shift σ is chosen near λ , and σ / ∈{ λ, d i , d i − } . This results in a DPR1 matrix( A − σI ) − = (cid:20) ( D − σI ) − (cid:21) + ρuu T , u = (cid:2) z T ( D − σI ) − − (cid:3) T , ρ = 1 a − z T ( D − σI ) − z . Eigenvalues of this matrix are zeros of ϕ ( λ ) = 1 + ρ n X j =1 u j ( d j − σ ) − λ , and the absolutely largest eigenvalue is computed accurately. Eventually, ρ needs to be computed in higher precision. s ss d i − λ i − λ i +1 d i λ ( λ i ) s ss d i − λ i − λ i +1 d i λ ( λ i )(a) (b) Figure 3: Typical situations from Remark 3
4. Final algorithm
If neither of the conditions
C2 (i) and
C2 (ii) hold, in order to guaranteethat λ will be computed with high relative accuracy, the element b from thematrix A − i needs to be computed in higher precision. The following theoremimplies that if 1 ≪ K b ≤ O (1 /ε M ), it is sufficient to evaluate (11) in double ofthe working precision. Theorem 7. If − a > in (11), set P = − a + z T D − z , Q = − z T D − z , and if − a < in (11) set P = z T D − z , Q = a − z T D − z . Determining whether ρ needs to be computed in higher precision is done similarly asdetermining whether element b of A − i needs to be computed in higher precision, which isdescribed in Section 4. Further, Theorem 7 implies that it suffices to compute ρ in double ofthe working precision. If K b ≥ O (1 /ε M ), that is, if K b = 1 /ε E for some ε E < ε M , then, in view of Theorem 7, b needs to be computed with extended precision ε E . Usage of higher precision in conjunction with the eigenvalue computation for DPR1 ma-trices is analyzed in [2], but there the higher precision computation is potentially neededin the iterative part. This is less convenient than our approach where the higher precisioncomputation is used only to compute one element. otice that in both cases P, Q ≥ and b = ( P − Q ) /ζ i . Let e P = f l ( P ) and e Q = f l ( Q ) be evaluated in standard precision, ε M . Assume that e P = e Q and K b ≤ O (1 /ε M ) . If P , Q and b are all evaluated in double of the workingprecision, ε M , then (17) holds with | κ b | ≤ O ( n ) . We summarize the above results in one, complete algorithm, aheig . Thealgorithm first checks the components of the vector z . If they are of the sameorder of magnitude, the eigenpair ( λ, v ) is computed by Algorithm 1. If that isnot the case, the quantity K b is computed, and if K b ≫
1, the eigenpair ( λ, v )is computed by Algorithm 1 but with evaluation of b in double of the workingprecision. At the end, the quantity K ν is computed, and if K ν ≫
1, one of theremedies from Remark 3 is applied.
Algorithm 5 [ λ, v ] = aheig ( D, z, α, k ) % Computes the k -th eigenpair of an ordered irreducible arrowhead matrix% A = [diag ( D ) z ; z ′ α ]compute the shift i as in the first part of Algorithm 1 if the quantity (cid:18) n − P j =1 j = i | ζ j | (cid:19) / | ζ i | from (25) is of O ( n ) % standard precision is enough [ λ, v ] = aheig basic ( D, z, α, k ) else compute the quantity K b from (18) if K b ≫ % double precision is necessary [ λ, v ] = aheig basic ( D, z, α, k ) with evaluation of b in double precision else % standard precision is enough [ λ, v ] = aheig basic ( D, z, α, k ) endend compute the quantity K ν from (13) if K ν ≫ end Implementation of the double of the working precision depends upon whetherthe input is considered to be binary or decimal.Double standard precision in Matlab, which assumes that input is binary,is obtained by using a combination of commands vpa , digits and double [13],where 16 digits(d) specifies the number of significant decimal digits d used to dovariable precision arithmetic vpa ,- vpa(x) uses variable-precision arithmetic to compute x to d decimal digitsof accuracy,- double(x) converts x to standard precision.The assignment a1=vpa(a,32) pads the binary representation of a withzeros, which means that the decimal interpretation of the variable a1 may havenon-zero entries after 16-th significant decimal digit. The same effect is obtainedin Intel FORTRAN compiler ifort [12] by the following program segment real(8) areal(16) a1...a1=a However, the user can assume that the true input is given as a decimalnumber, which is, for example, assumed by extended precision computation inMathematica [20]. In this case, the options in Matlab are to either use symboliccomputation, or to cast the input to a string, and then convert it to extendedprecision: a1=vpa(num2str(a,16),32)
In this case, the the decimal interpretation of the variable a1 has all zeroentries after 16-th significant decimal digit, but the binary representation ofthe variable a is, in general, padded with non-zero entries. The same effectis obtained in ifort writing to and reading from a string variable as in thefollowing program segment: real(8) areal(16) a1character(25) string...write(string,*) aread(string,*) a1 If the input consists of numbers for which decimal and binary representationare equal (for example, integers, as in Example 3 below), then the two aboveapproaches give the same results.
5. Numerical Examples
We illustrate out algorithm with four numerically demanding examples. Ex-amples 1 and 2 illustrate Algorithm 1, Example 3 illustrates the use of doubleprecision arithmetic, and Example 4 illustrates an application of higher dimen-sion. 17 xample 1.
In this example both quantities K ν from (13) and K b from (18) arefor all eigenvalues approximately equal to 1, so we guarantee that all eigenvaluesand all components of their corresponding eigenvectors are computed with highrelative accuracy by Algorithm 5 ( aheig ) using only standard machine precision.Let A = · − − − − − · − . The eigenvalues computed by Matlab [13] routine eig , Algorithm 5 and Math-ematica [20] with 100 digits precision, are, respectively: λ ( eig ) λ ( aheig ) λ ( Math ) . · . · . · . · − . · − . · − . · − . · − . · − − . · − − . · − − . · − − . · − − . · − − . · − − . · − − . · − − . · − We see that even the tiniest eigenvalues λ and λ , computed by Algorithm 5, areexact to the machine precision, which is not true for the eigenvalues computedby eig . Because of the accuracy of the computed eigenvalues, the eigenvectorscomputed by Algorithm 5 are componentwise accurate up to machine precision,and therefore, orthogonal up to machine precision. For example: v ( eig )4 v ( aheig )4 v ( Math )4 . · − − . · − − . · − . · − − . · − − . · − . · − − . · − − . · − − . · − . · − . · − − . · − . · − . · − − . · − . · − . · − . Example 2.
In this example, despite very close diagonal elements, we againguarantee that all eigenvalues and all components of their corresponding eigen-vectors are computed with high relative accuracy, without deflation. Let A = ε M ε M ε M ε M
41 2 3 4 0 where ε M = 2 · − = 2 . · − . For this matrix the quantities K ν and K b are again of order one for all eigenvalues, so Algorithm 5 uses only standard18orking precision. The eigenvalues computed by Matlab and Algorithm 5 are: λ ( eig ) λ ( aheig ) . . ε M ε M ε M ε M ε M ε M − . − . The eigenvalues computed by Mathematica with 100 digits precision, properlyrounded to 32 decimal digits are : λ ( Math ) . . . . − . The eigenvalues computed by Matlab are accurate according to standard per-turbation theory, but they do not satisfy the interlacing property. Furthermore,the Matlab’s eigenvectors corresponding to λ , λ and λ only span an accurateeigenspace, and are not individually accurate. On the other hand, the eigen-values computed by Algorithm 5 are exact (they coincide with the eigenvaluescomputed by Mathematica properly rounded to 16 decimal digits). Notice thatdespite of very close eigenvalues, Algorithm 5 works without deflation. Due tothe accuracy of the computed eigenvalues, the eigenvectors computed by Al-gorithm 5 are componentwise accurate up to the machine precision, and aretherefore orthogonal.If, as suggested in Remark 2, the algorithms are modified to return d i and µ (both in standard precision), then for the eigenvalues λ , λ and λ the cor-responding pairs ( d i , µ ) give representations of those eigenvalues to 32 decimaldigits. That is, exact values d i + µ properly rounded to 32 decimal digits areequal to the corresponding eigenvalues computed by Mathematica as displayedabove. Example 3.
In this example we can guarantee all eigenvalues and eigenvec-tors, componentwise will be computed with high relative accuracy only if weuse double of the working precision when computing b from (11) in matrices Since, as described in Section 4.1, Mathematica uses decimal representation of the input,in order to obtain accurate eigenvectors we need to define ε M in Mathematica with the outputof Matlab’s command vpa(eps), ε M = 2 . · − . − , A − , A − and A − . Let A = . The quantities K ν and K b are: K ν K b − − . · − . · . · − . · . · . · . · − . · . · . · It is clear, from the condition numbers, that the element b in each of the ma-trices A − , A − , A − and A − needs to be computed in double of the workingprecision. For example, A = A − d I = − − − − − . The element b = (cid:2) A − (cid:3) computed by Algorithm 2 gives b = 6 . inv yields b = 6 . b in doubleof the working precision gives the correct value b = 6 . aheig basic , using only standardworking precision), Algorithm 5 ( aheig , using double of the working precision tocompute respective b ’s) and Mathematica with 100 digits precision, respectively,are: λ aheig basic λ aheig λ Math . · . · . · . · . · . · . · . · . · . · . · . · . · . · . · − . · − − . · − − . · − The eigenvectors computed by Algorithm 5 are componentwise accurate to ma-chine precision and therefore orthogonal. Algorithm 5 does not compute K ν and K b for the first eigenvalue, since it is an absolutelylargest one. xample 4. This example comes from the research related to decay of excitedstates of quantum dots in in real photon crystals [15]. In this case- α is quantum dot transition frequency,- d i is a frequency of the i -th optical mode, and- ζ i is an interaction constant of the quantum dot with the i -th opticalmode.The size of the matrix is changeable but, in realistic cases, it is between 10 and10 . We ran a test example for n = 2501 where, typically, d i ∈ [5 . · , . · ] ,ζ i ∈ [1 . · , . · ] ,α = 9 . · . For this matrix the condition number K ν ∼ z do not differ by much in size, thus the conditions C1 and C2 (ii) from Section 3 are fulfilled. Therefore, all eigenvalues and all compo-nents of all eigenvectors are computed with high relative accuracy by Algorithm5 using only standard working precision. On the other hand about half of theeigenvalues computed by the Matlab routine eig do not satisfy the interlacingproperty.
6. Applications
In this section we extend our results to eigenvalue decompositions of Her-mitian arrowhead matrices, singular value decompositions of real triangular ar-rowhead matrices and eigenvalue decompositions of real symmetric diagonal-plus-rank-one matrices.
Let C = (cid:20) D zz ∗ α (cid:21) , where D = diag( d , d , . . . , d n − ) , is a real diagonal matrix of order n − z = (cid:2) ζ ζ · · · ζ n − (cid:3) ∗ , is a complex valued vector and α is a real scalar. Here z ∗ denotes the conjugatetranspose of z . As in Section 1, we assume that C is irreducible. The eigenvaluedecomposition of C is given by C = U Λ U ∗ λ , . . . , λ n ) ∈ R n × n is a diagonal matrix of eigenvalues, and U = (cid:2) u u · · · u n (cid:3) is an unitary matrix of the corresponding eigenvectors.To apply Algorithm 5 to Hermitian arrowhead matrix we first transform C to real symmetric arrowhead matrix A by diagonal unitary similarity: A = Φ ∗ C Φ = (cid:20) D | z || z | T α (cid:21) , (28)where Φ = diag (cid:18) ζ | ζ | , ζ | ζ | , . . . , ζ n − | ζ n − | , (cid:19) We now compute the k -th eigenpair ( λ, v ) of A by Algorithm 5, and set u =Φ v . Since we guarantee high relative accuracy of the eigenvalue decompositionof A computed by Algorithm 5, we also guarantee high relative accuracy of theeigenvalue decomposition of C . Notice that, if double precision is needed tocompute b in Algorithm 5, the modules | ζ i | in (28) need to be computed indouble of the working precision, as well. Remark 4.
Similarly, for irreducible non-symmetric arrowhead matrix G = (cid:20) D ˘ z ˚ z T α (cid:21) , where sign(˚ ζ i ) = sign(˘ ζ i ), i = 1 , . . . , n −
1, we define the diagonal matrixΨ = diag sign(˚ ζ ) s ˚ ζ ˘ ζ , . . . , sign(˚ ζ n − ) s ˚ ζ n − ˘ ζ n − , . The matrix A = Ψ − G Ψ = (cid:20)
D zz T α (cid:21) , where ζ i = q ˘ ζ ˚ ζ i is an irreducible symmetric arrowhead matrix.We now compute the k -th eigenpair ( λ, v ) of A by Algorithm 5. The eigenpairof G is then ( λ, Ψ v ). set u = Φ v . Since we guarantee high relative accuracy ofthe eigenvalue decomposition of A , we also guarantee high relative accuracy ofthe eigenvalue decomposition of G . Notice that, if double precision is neededto compute b in Algorithm 5, the elements ζ i need to be computed in double ofthe working precision, as well. Let B = (cid:20) D z α (cid:21) ,
22e an irreducible upper triangular arrowhead matrix, that is, d i = d j for i = j and ζ i = 0 for all i . The matrix A = B T B = (cid:20) D Dzz T D α + z T z (cid:21) , is an irreducible symmetric arrowhead matrix.When applying Algorithm 5 to the matrix A , we must ensure that all com-ponents of A − i in (10) are computed to high relative accuracy. This is obviouslytrue for elements of the vectors w i and w . Diagonal elements, except b , arecomputed with high relative accuracy as differences of squares of original quan-tities, [ A − i ] jj = 1( d j − d i )( d j + d i ) , j = i. The element b = [ A − i ] ii from (11) is computed as b = 1 d i ζ i − α − z T z + d i + X j = i d j ζ j ( d j − d i )( d j + d i ) . If double precision is needed in Algorithm 5, all entries of A need to be computedin double precision.Let B = U Σ V T be the singular value decomposition of B , where Σ =diag( σ , . . . , σ n ) are the singular values, the columns of V are the correspondingright singular vectors and the columns of U are the corresponding left singularvectors. We first compute the k -th eigenpair ( λ, v ) of A by Algorithm 5. Then σ = √ λ is the corresponding singular value of B and v is the correspondingright singular vector. The value σ and all components of v are computed toalmost full accuracy. From the relation U T B = Σ V T for the k -th row we have (cid:2) u T n − u n (cid:3) (cid:20) D z α (cid:21) = σ (cid:2) v T n − v n (cid:3) , which implies u n − = σv n − D − . From the relation BV = U Σ for the k -th column we have (cid:20) D z α (cid:21) (cid:20) v n − v n (cid:21) = σ (cid:20) u n − u n (cid:21) , which implies u n = αv n σ . Components of u are computed by multiplication and division of quantitieswhich are accurate to almost full machine precision, so the are accurate toalmost full machine precision, as well. In view of Theorem 7, if double precision computation is necessary, the positive andnegative parts of this formula should be computed separately, and then added. .3. Diagonal-plus-rank-one matrices Let M = D + uu T , where D = diag( d , . . . , d n ) , d > d > · · · > d n ,u = (cid:2) u · · · u n (cid:3) T , u i = 0 , i = 1 , . . . , n, be a n × n irreducible ordered real symmetric diagonal-plus-rank-one (DPR1)matrix. Let ¯ D = diag( d , . . . , d n − ) , ∆ = ( ¯ D − d n ) / , ¯ u = (cid:2) u · · · u n − (cid:3) T ,L = (cid:20) u n ∆ − − ¯ u T ∆ − (cid:21) . Then A = L − M L = (cid:20) ¯ D zz T α (cid:21) , where z = ¯ u ∆ , α = d n + u T u, is an irreducible real symmetric arrowhead matrix.When applying Algorithm 5 to the matrix A , we must ensure that all com-ponents of A − i in (10) are computed to high relative accuracy. This is obviouslytrue for elements of the vectors w i and w . Diagonal elements, except b , arecomputed with high relative accuracy as differences of original quantities, andthe element b = [ A − i ] ii from (11) is computed as b = 1 ζ i − d n − u T u + d i + X j = i ζ j d j − d i . If double precision is needed in Algorithm 5, all entries of A need to be computedin double precision.Let M = Q Λ Q T and A = V Λ V T be the eigenvalue decompositions of M and A , respectively. Since M is by assumption irreducible, its eigenvalues satisfyinterlacing property λ > d > λ > d > · · · > λ n > d n . (29)We first compute the k -th eigenpair ( λ, v ) of A by Algorithm 5. The value λ and all components of v are computed to almost full accuracy. The relation V T AV = V T L − M LV = Λ implies that the columns of the matrix X = LV M . Further, since, by (29), alleigenvalues are simple, we conclude that X = Q Σ, where Σ = diag( σ , . . . , σ n )is a positive definite matrix. Notice that Q Σ V T = L is, in fact, singular valuedecomposition of L .Equating k -th columns of the equation X = LV gives x = (cid:20) ¯ xx n (cid:21) = Lv = (cid:20) u n ∆ − − ¯ u T ∆ − (cid:21) (cid:20) ¯ vv n (cid:21) , where x and v are partitioned according to L . This immediately implies that¯ x = u n ∆ − ¯ v. Notice that, since all components of ¯ v are computed to almost full, accuracy, thesame holds for the components of ¯ x , and it remains to compute x n accurately.Let q = (cid:20) ¯ qq n (cid:21) be the k -th column of Q and let σ = Σ kk . Equating k -th rows of the equation X − = Σ − Q T = V T L − gives for the n -th element q n σ = x n σ = v n . Thus, x n = σ v n and, in order to compute x n , it is necessary to compute σ . From X = U Σ = LV it follows that V T L T LV = Σ , or, equivalently, LV = L − T V Σ . Equating k -thcolumns of this equation gives∆ − ¯ v u n = (cid:20) ∆¯ v u n + ∆ − ¯ u v n (cid:21) σ . This gives n − σ , and we can choose the numerically mostaccurate one.Therefore, x n will be computed to almost full machine precision, as are theentries of ¯ x , and it remains to normalize x and obtain q = x/σ . Remark 5.
Notice that DPR1 matrices of the form D − uu T cannot be reducedto symmetric arrowhead matrix by the procedure described in this section. Byusing ideas from this paper, it is possible to derive highly accurate algorithmfor DPR1 matrices without prior transformation to arrowhead form. This al-gorithm, which is a topic of our forthcoming paper, covers more general DPR1matrices of the form D + ρuu T , ρ ∈ R . ppendix A. Proofs Proof of Theorem 1.
Let e µ and e λ be defined by (14). Then e λ ≡ f l ( d i + e µ ) = ( d i + e µ ) (1 + ε ) . By simplifying the equality( d i + µ (1 + κ µ ε M )) (1 + ε ) = λ (1 + κ λ ε M )and using λ = µ + d i , we have d i ε + µ ( κ µ ε M + ε ) = λκ λ ε M . Taking absolute value gives | κ λ | ≤ | d i | + | µ || λ | ( | κ µ | + 1) . (cid:3) Proof of Theorem 2.(i)
The assumption sign ( d i ) = sign ( µ ) immediately implies | d i | + | µ || λ | = | d i + µ || d i + µ | = 1 . (ii) The assumptions imply that either0 < d i +1 < λ < d i , µ < , or d i < λ < d i − < , µ > . In the first case λ is closest to the pole d i and | d i | + | µ || λ | ≤ | d i | + | d i − d i +1 | | d i + d i +1 | ≤ d i + d i − d i +112 d i + d i +1 ≤ d i − d i +112 d i + d i +1 ≤ d i d i = 3 . Here we used the inequalities | µ | ≤ | d i − d i +1 | and | λ | ≥ | d i + d i +1 | for thefirst inequality, d i − d i +1 > d i + d i +1 > d i +1 > (cid:3) roof of Theorem 3. Let x and e x be defined by (12) and (21), respectively. The theorem obviouslyholds for x n = e x n = −
1. For e x i we have e x i = f l (cid:18) − ζ i e µ (cid:19) = − ζ i µ (1 + κ µ ε M ) (1 + ε ) = x i (1 + ε x i ) . By using (16) and (21), the first order approximation gives | ε x i | ≤ ( | κ µ | + 1) ε M . For j / ∈ { i, n } , by solving the equality e x j = ζ j (( d j − d i ) (1 + ε ) − µ (1 + κ µ ε M )) (1 + ε ) (1 + ε ) = ζ j d j − λ (1 + ε x )for ε x , using (16) and λ = µ + d i , and ignoring higher order terms, we have ε x = ( d j − d i ) ( ε + ε + ε ) − µ ( κ µ ε M + ε + ε ) d j − λ . Therefore, | ε x | ≤ | d j − d i | + | µ || d j − λ | ( | κ µ | + 3) ε M . (A.1)To complete the proof we need to analyze two cases. Ifsign ( d j − d i ) = − sign µ, then | d j − d i | + | µ || d j − λ | = | d j − d i − µ || d j − λ | = | d j − λ || d j − λ | = 1 . If sign ( d j − d i ) = sign µ, then, since d i is pole closest to λ , we have | µ | ≤ . | d j − d i | and | d j − d i | + | µ || d j − λ | ≤ | d j − d i | + | µ || d j − d i | − | µ | ≤ | d j − d i | | d j − d i | = 3 . Finally, the theorem follows by inserting this into (A.1). (cid:3)
Proof of Theorem 4.
For the non-zero computed elements of the matrix A − i from (10) and (11),except the element b = [ A − i ] ii , we have: f l ( (cid:2) A − i (cid:3) jj ) = 1( d j − d i ) (1 + ε ) (1 + ε ) , j / ∈ { i, n } ,f l ( (cid:2) A − i (cid:3) ji ) = f l ( (cid:2) A − i (cid:3) ij ) = − ζ j ( d j − d i ) (1 + ε ) ζ i (1 + ε ) (1 + ε ) , j / ∈ { i, n } ,f l ([ A − i ] ni ) = f l ([ A − i ] in ) = 1 ζ i (1 + ε ) , | ε k | ≤ ε M for all indices k . The first statement of the theorem now followsby using standard first order approximations.Similar analysis of the formula (11) yields f l ([ A − i ] ii ) = e b = b + δb, where | δb | ≤ ζ i (cid:0) |− a | + (cid:12)(cid:12) z T D − z (cid:12)(cid:12) + (cid:12)(cid:12) z T D − z (cid:12)(cid:12)(cid:1) ( n + 3) ε M . (A.2)This, in turn, implies (17) with | κ b | ≤ | δb || b | | ε M | = ( n + 3) | a | + (cid:12)(cid:12) z T D − z (cid:12)(cid:12) + (cid:12)(cid:12) z T D − z (cid:12)(cid:12)(cid:12)(cid:12) − a + z T D − z + z T D − z (cid:12)(cid:12) = ( n + 3) K b , where K b is defined by (18). (cid:3) Proof of Theorem 5.
Let ^ (cid:0) A − i (cid:1) = A − i + δA − i . Therefore, | b ν − ν | = k δA − i k , which, together with (24), implies | νκ ν ε M | ≤ k δA − i k . (A.3)Theorem 4 implies that k δA − i k ≤ ( n + 3) k| A − i |k K b ε M . Since k| A − i |k ≤ √ n k A − i k and | ν | = k A − i k , from (A.3) we have | κ ν | ≤ ( n + 3) √ nK b , (A.4)which proves the first part of the bound (25).For the second part of the proof, notice that Theorem 4 also implies k δA − i k ≤ k|A|k ε M + | δb | , (A.5)where A is equal to the matrix A − i without b (that is, with A ii = 0).By bounding (A.3) with (A.5) and (A.2), and dividing the resulting inequal-ity by | νε m | , we have | κ ν | ≤ √ n + ( n + 3) (cid:18) | ν | | − a | + (cid:12)(cid:12) z T D − z (cid:12)(cid:12) + (cid:12)(cid:12) z T D − z (cid:12)(cid:12) ζ i (cid:19) . (A.6)28ince | − a | ζ i = 1 ζ i | − a + z T D − z + z T D − z − z T D − z − z T D − z |≤ | b | + 1 ζ i ( | z T D − z | + | z T D − z | ) , from (A.6) it follows | κ ν | ≤ √ n + ( n + 3) (cid:18) | b || ν | + 2 | ν | (cid:12)(cid:12) z T D − z (cid:12)(cid:12) + (cid:12)(cid:12) z T D − z (cid:12)(cid:12) ) ζ i (cid:19) . (A.7)Since | b | ≤ | ν | and (cid:13)(cid:13) A − i (cid:13)(cid:13) = | ν | = max k x k =1 (cid:13)(cid:13) A − i x (cid:13)(cid:13) ≥ (cid:13)(cid:13) A − i e k (cid:13)(cid:13) = s d k − d i ) + ζ k ζ i ( d k − d i ) ≥ | ζ k || ζ i | | d k − d i | , by simply dividing each term ζ k ζ i | d k − d i | in (A.7) with the corresponding quotient | ζ k || ζ i | | d k − d i | , we obtain | κ ν | ≤ √ n + ( n + 3) (cid:18) ζ i n − X k =1 k = i | ζ k | (cid:19) . (A.8)The bound (25) now follows from (A.4) and (A.8). (cid:3) Proof of Theorem 6.
We first prove the bound (27). Since e ν = f l ( b ν ) is computed by bisection,from (22) we have e ν = b ν (1 + κ bis ε M ) . This and (24) imply e ν = ν (1 + κ ν ε M )(1 + κ bis ε M ) . Since b µ = f l (1 / b ν ), the bound (27) follows by ignoring higher order terms. Thebound (26) now follows by inserting (27) into Theorems 1 and 2. (cid:3) roof of Theorem 7. Let the assumptions of the theorem hold. Let b be computed in double of theworking precision, ε M , and then stored in the standard precision. The standardfloating-point error analysis with neglecting higher order terms gives P (cid:0) κ P ε M (cid:1) − Q (cid:0) κ Q ε M (cid:1) ζ i (cid:0) κ ε M (cid:1) = P − Qζ i (1 + κ b ε M ) ≡ b (1 + κ b ε M ) , where | κ P | , | κ Q | ≤ ( n + 1) and | κ | ≤
3. Solving the above equality for κ b ,neglecting higher order terms, and taking absolute values gives | κ b | ≤ | P | + | Q || P − Q | ( n + 4) ε M ≡ K b ( n + 4) ε M . Since, by assumption, K b ≤ O (1 /ε M ), this implies | κ b | ≤ O ( n ) , as desired. (cid:3) References [1] E. Anderson et al.,
LAPACK Users’ Guide , SIAM 3rd ed., Philadelphia,(1999).[2] J. L. Barlow, Error analysis of update methods for the symmetric eigenvalueproblem, SIAM J. Matrix Anal. Appl., 14 (1993) 598-618.[3] M. Bixon and J. Jortner, Intramolecular radiationless transitions, J. Chem.Physics, 48 (1968) 715-726.[4] C. F. Borges, W. B. Gragg, A parallel Divide - and - Conquer Methodfor the Generalized Real Symmetric Definite Tridiagonal Eigenproblem, in
Numerical Linear Algebra and Scientific Computing , L. Reichel, A. Ruttanand R. S. Varga, eds., de Gruyter, Berlin (1993) 11-29.[5] J. R. Bunch and C. P. Nielsen, Rank-one modification of the symmetriceigenproblem, Numer. Math., 31 (1978) 31-48.[6] J. J. M. Cuppen, A divide and conquer method for the symmetric tridiag-onal eigenproblem, Numer. Math., 36 (1981) 177-195.[7] J. Dongarra and D. Sorensen, A fully parallel algorithm for the symmetriceigenvalue problem, SIAM J. Sci. Statist. Comput., 8 (1987) 139-154.[8] J. W. Gadzuk, Localized vibrational modes in Fermi liquids, general The-ory, Phys. Rev. B, 24 (1981) 1651-1663.309] D. Goldberg, What Every Computer Scientist Should Know AboutFloating-Point Arithmetic, ACM Computing Surveys, 23:1 (1991) 5-48.[10] G. H. Golub and C. F. Van Loan,
Matrix Computations , The John HopkinsUniversity Press, Baltimore, 3rd ed. (1996).[11] M. Gu and S. C. Eisenstat, A divide-and-conquer algorithm for the sym-metric tridiagonal eigenproblem, SIAM J. Matrix Anal. Appl., 16 (1995)79-92.[12]
Intel Fortran Compiler , http://software.intel.com/en-us/fortran-compilers[13]
MATLAB. The MathWorks
The Symmetric Eigenvalue Problem , Prentice-Hall, Engle-wood Cliffs, (1980).[19] J. H. Wilkinson,
The Algebraic Eigenvalue Problem , Clarendon Press, Ox-ford, (1965).[20]