Accurate eigenvalue decomposition of arrowhead
matrices and applications
N. Jakovčević Stora,1,∗, I. Slapničara,1 , J. Barlowb,2a Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture,

University of Split, Rudjera Boškovića 32, 21000 Split, Croatiab 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 symmetric arrowhead matrix. The algorithm computes all eigenvalues and all components of the corresponding eigenvectors with high relative accuracy in O(n2 )operations. The algorithm is based on a shift-and-invert approach. Double
precision is eventually needed to compute only one element of the inverse of
the shifted matrix. Each eigenvalue and the corresponding eigenvector can becomputed separately, which makes the algorithm adaptable for parallel computing. Our results extend to Hermitian arrowhead matrices, real symmetricdiagonal-plus-rank-one matrices and singular value decomposition of real triangular arrowhead matrices.
Keywords: eigenvalue decomposition, arrowhead matrix, high relative
accuracy, singular value decomposition
2000 MSC: 65F15

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. Since
eigenvalues are invariant under similarity transformations, we can symmetricallypermute the rows and the columns of the given matrix. Therefore, we assume
without loss of generality that the matrix A is a n×n real symmetric arrowhead
matrix of the form
∗ Corresponding author

Email addresses: nevena@fesb.hr (N. Jakovčević Stor), ivan.slapnicar@fesb.hr (I.
Slapničar), barlow@cse.psu.edu (J. Barlow)
1 The research of Ivan Slapničar and Nevena Jakovčević Stor was supported by the Ministryof Science, Education and Sports of the Republic of Croatia under grant 023-0372783-1289.
2 The research of Jesse L. Barlow was supported by the National Science Foundation under
grant CCF-1115704.

Preprint submitted to Elsevier

March 1, 2013


A=

D
zT

z
α


,

(1)

where
D = diag(d1 , d2 , . . . , dn−1 )is diagonal matrix of order n − 1,
z=



ζ1

···

ζ2

ζn−1

T

(2)

is a vector and α is a scalar.
Such matrices arise in the description of radiationless transitions in isolated
molecules [3], oscillators vibrationally coupled with a Fermi liquid [8], quantum
optics [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 inO(n2 ) operations.
Without loss of generality we may assume that A is irreducible, that is,ζi 6= 0, for all ianddi 6= dj , for all i 6= j, i, j = 1, . . . , n − 1.
If A has a zero in the last column, say ζi = 0, then the diagonal element di
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 di = dj , then di is eigenvalue of matrix A (this follows from the interlacing
property (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 thatd1 > d2 > · · · > dn−1 .
(3)Hence, we will consider only ordered and irreducible arrowhead matrices. Without loss of generality we can also assume that ζi > 0 for all i, which can be attained 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(λ1 , λ2 , . . . , λn )is a diagonal matrix whose diagonal elements are the eigenvalues of A, and
V = v1 · · · vn
2

is 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−1
X

ζi2
= α − λ − z T (D − λI)−1 z,
d
−
λ
i
i=1

and the corresponding eigenvectors are given by


−1
xi
(D − λi I) z,
xi =, i = 1, . . . , n.
vi =−1
kxi k2(5)(6)

Diagonal elements of the matrix D, di , 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
λ1 > d1 > λ2 > d2 > · · · > dn−2 > λn−1 > dn−1 > λn .

(7)Since A is symmetric, its eigenvalues may be computed by invoking anyof a number of standard programs (LAPACK [1]). However, these programs
usually 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(n2 ) computations
and O(n) storage. Although the idea is conceptually simple and in fact has
been 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 vi may not be
sufficiently orthogonal (see Example 3). The existing algorithms for arrowhead
matrices [11, 16] obtain orthogonal eigenvectors with the following procedure:- compute the eigenvalues λ̃i of A by solving (5);- construct a new matrix


Ã =

D
z̃ T

z̃
α̃



by solving inverse problem with the prescribed eigenvalues λ̃, and diagonal
matrix D, that is, compute new z̃ and α̃ as
v




u
u
i
n−1

Y
λ̃j − diY λ̃j − diuζ̃i = t di − λ̃n λ̃1 − di,(dj−1 − di ) j=i+1 (dj − di )j=2α̃ = λ̃n +

n−1
X
j=1

- compute eigenvectors of Ã by (6).
3λ̃j − dj .

Since the formulas for ζ˜i involve only multiplications, division and subtractions of exact quantities, each ζi is computed with relative error of O(εM ),where εM denotes the machine precision.3 Therefore, Ã = A + δA, wherekδAk2 = O(M ). Here k · k2 denotes the spectral matrix norm. We conclude
that 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 eigenvalues and there is no need for follow-up orthogonalization. The algorithm isbased on shift-and-invert technique. Basically, the eigenvalue λ is computed as
the largest or the smallest eigenvalue of the inverse of the matrix shifted to thepole di which is nearest to λ, that is,
λ=1+ di ,
ν(8)

where ν is either smallest or largest eigenvalue of the matrix
A−1≡ (A − di I)−1 .
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,

−1 
×
×
×



×
×





0
×

 = ×
× ×
× × × × ×



×
× ×
× × ×
× ×
×



×
,

0

and the inverse of the full arrowhead matrix is a diagonal-plus-rank-one (DPR1)
matrix,

−1 
×
×
×



×
×
×




 =
×
×





× ×
× × × × × ± uuT .



×
×
0

3 The machine precision ε
M is defined as a smallest positive number such that in the
floating-point arithmetic 1 + εM 6= 1. In Matlab or FORTRAN REAL(8) arithmetic εM =2.2204 · 10−16 , thus the floating-point numbers have approximately 16 significant decimal
digits. The term “double of the working precision” means that the computations are performed
with numbers having approximately 32 significant decimal digits, or with the machine precisionequal to ε2M .

4Our algorithm is completely parallel, since the computation of one eigenvalue and its eigenvector is completely independent of the computation of othereigenvalues and eigenvectors.
In Section 2 we describe the basic idea of our algorithm named aheig (ArrowHead EIGenvalues). In Section 3 we discuss the accuracy of the algorithm.
In Section 4 we present the complete algorithm which uses double of the working
precision, if necessary. In Section 5 we illustrate algorithm with few examples
and in Section 6 we apply our results to eigenvalue decomposition of Hermitian
arrowhead matrix, singular value decomposition of real triangular arrowhead
matrix and eigenvalue decomposition of real symmetric diagonal-plus-rank-one
matrix. 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 di be the pole which is closest to λ.
Clearly, from (7) it follows that either λ = λi or λ = λi+1 . Let Ai be the shiftedmatrix


D1 0
0 z1
 0
0
0 ζi 
,Ai = A − d i I = (9) 0
0 D2 z2 TT
z1 ζi z2awhere
D1 = diag(d1 − di , . . . , di−1 − di ),D2 = diag(di+1 − di , . . . , dn−1 − di ),T
z1 = ζ1 ζ2 · · · ζi−1,Tz2 = ζi+1 ζi+2 · · · ζn−1,
a = α − di .
Notice that D1 (D2 ) is positive (negative) definite.
Obviously, if λ is an eigenvalue of A, thenµ = λ − diis an eigenvalue of Ai , and vice versa, and they both have the same eigenvector.
The inverse of Ai is
 −1

D1
w1
0
0
 w1T
b
w2T 1/ζi 
,
A−1
=
(10)
i
 0
w2 D2−1
0 
0
1/ζi
0
0

5

where1,ζi1w2 = −D2−1 z2 ,ζi1b = 2 −a + z1T D1−1 z1 + z2T D2−1 z2 .
ζiw1 = −D1−1 z1(11)Notice thatb = f¯ (di ) /ζi2
where

−1
z̄
f¯ (di ) = α − di − z̄ T D̄ − di I

where D̄ is the diagonal matrix D without di and z̄ is z without ζi .
The eigenvector x from (6) is given by−1(D1 − µI) z1
 
x1


ζi

−
 ..  .

x= . =µ (D2 − µI)−1 z2 
xn
−1

(12)

If λ is an eigenvalue of A which is closest to the pole di , then µ is the
eigenvalue of matrix Ai which is closest to zero andν=1= ± A−1.
i2µIn this case, if all entries of A−1
are computed with high relative accuracy,ithen, according to standard perturbation theory, ν is computed to high relative
accuracy (by any reasonable algorithm). In Section 3 we show that all entriesof A−1are indeed computed to high relative accuracy, except possibly b (seei,(11)). If b is not computed to high relative accuracy and it influences A−1i2
it is sufficient to compute it in double of the working precision (see Section 4).
Further, if µ is not the eigenvalue of Ai which is closest to zero, then |ν| <A−1, and the quantity
i
2
A−1
i
2
Kν =
(13)
|ν|
tells us how far is ν from the absolutely largest eigenvalue of A−1
i . If Kν  1,then the standard perturbation theory does not guarantee that the eigenvalueµ will be computed with high relative accuracy. Remedies of this situation are
described in Remark 2.
Finally, the exterior eigenvalues λ1 and λn of A are excluded from the above
reasoning. The absolutely larger of those is equal to kAk2 and is, according
to standard perturbation theory, computed with high relative accuracy. The
6other exterior eigenvalue can also be computed with high relative accuracy if it
is of the same order of magnitude as the larger one. If this is not the case, theabsolutely smaller eigenvalue is computed as interior eigenvalues of A.
With this approach the componentwise high relative accuracy of the eigenvectors computed by (12) follows from high relative accuracy of the computed
eigenvalues (see Theorem 3). Componentwise high relative accuracy of the computed eigenvectors implies, in turn, their orthogonality.
The described procedure is implemented in algorithm aheig basic (Algorithm 1). The computation of the inverse of the shifted matrix, A−1
i , according
to formulas (10) and (11), is implemented in Algorithm 2. Algorithm 3 computes 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. We
assume tha standard model of floating point arithmetic where subtraction is
preformed 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 their
proofs we shall use the standard first order approximations, that is, we neglect
the terms of order O(ε2M ) or higher. Moreover, we assume that neither overflow
or underflow occurs during the computation.
We shall use the following notation:
MatrixA
Ai
ei = f l(Ai )
A
A−1
i
−1
^
(A
) = f l(A−1 )iiExact eigenvalue
λ
µ
µ
b
ν

Computed eigenvalueeλνbνe = f l(bν)−µ
e = f l(b
µ)
−

(14)

Here


D1 (I + E1 ) 000ei = f l (Ai ) = 
A

0
0
z1T
ζi


0
z1

0
ζi
,

D2 (I + E2 )z2Ta (1 + εa )z2where E1 and E2 are diagonal matrices whose elements are bounded by εM inabsolute values and |εa | ≤ εM .

7

Algorithm 1[λ, v] = aheig basic (D, z, α, k)% Computes the k -th eigenpair of an irreducible arrowhead matrix% A = [diag (D) z; z 0 α]
n = max(size(D)) + 1% Set the shift index i initially to zero
i=0
% Compute exterior eigenvalue (k = 1 or k = n)if k == 1
λ = bisect(D, z, α, 0 R0 )λn = bisect(D, z, α, 0 L0 )if |λn /λ| < 10
% the eigenvalue is accuratev = vect(D, z, λ)else% determine the shift σ , the shift index i, and whether λ is on the left
% or the right side of the nearest poleσ = d1i=1side = 0 R0end
elseif k == nλ = bisect(D, z, α, 0 L0 )λ1 = bisect(D, z, α, 0 R0 )if |λ1 /λ| < 10
% the eigenvalue is accuratev = vect(D, z, λ)else% determine the shift, shift index, and sideσ = dn−1i=n−1side = 0 L0end

8Algorithm 1 (Continued)else% Compute interior eigenvalue (k ∈ {2, . . . , n − 1})% Determine the shift (the nearest pole) and the shift index

Dtemp = D − dk
atemp = α − dk
middle = Dtempk−1 /2P
F middle = atemp − middle − (z 2 ./(Dtemp − middle))
if F middle < 0
σ = dki=kside = 0 R0else
σ = dk−1i=k−1
side = 0 L0end
end
if i 6= 0% Compute the inverse of the shifted matrix, A−1
i[invD1 , invD2 , w1 , w2 , wζ , b] = invA(D, z, α, i)% Compute the leftmost or the rightmost eigenvalue of A−1iν = bisect([invD1 ; 0; invD2 ], [w1 ; wζ ; w2 ], b, side)% Compute the corresponding eigenvectorµ = 1/ν
v = vect(D − σ, z, µ)% Shift the eigenvalue back

λ=µ+σ
end

Algorithm 2[invD1 , invD2 , w1 , w2 , wζ , b] = invA (D, z, α, i)% Computes the inverse of an arrowhead matrix A = [diag(D − di ) z; z 0 α − di ]% according to (10) and (11).

n = max(size(D)) + 1
D = D − dia = α − diw1 = −z1:i−1 ./D1:i−1 /zi
w2 = −zi+1:n−1 ./Di+1:n−1 /ziwζ = 1/zi
invD1 = 1./D1:i−1
invD2 = 1./Di+1:n−1
b = (−a + sum(z1:i−1 .^2./D1:i−1 ) + sum(zi+1:n−1 .^2./Di+1:n−1 ))/zi ^2

9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 0 α] by bisection.

n = max(size(D)) + 1% Determine the starting interval for bisection, [lef t, right]

if side == 0 L0
lef t = min{D − |z|, α − kzk1 }right = min dielseright = max{D + |z|, α + kzk1 }lef t = max diend% Bisection

middle = (lef t + right)/2
while (right − lef t)/abs(middle) > 2 ∗ eps
F middle = α − middle − sum(z.^2./ (D − middle))
if F middle > 0
lef t = middle
elseright = middle
end
middle = (lef t + right)/2end% Eigenvalue

λ = right

Algorithm 4v = vect (D, z, λ)% Computes the eigenvector of an arrowhead matrix A = [diag(D) z; z 0 α]
% which corresponds to the eigenvalue λ by using (6).

v = [z./(D − λ); −1]v = v/kvk2

10

Further we define the quantities κλ , κµ and κb as follows:e = f l (λ) = λ (1 + κλ εM ) ,
λ

(15)

µ
e = f l (µ) = µ (1 + κµ εM ) ,eb = f l (b) = b (1 + κb εM ) .

(16)
(17)

We also define the quantity
Kb =|a| + z1T D1−1 z1 + z2T D2−1 z2−a + z1T D1−1 z1 + z2T D2−1 z2

.

(18)

3.1. Connection between accuracy of λ and µ
Letλ = µ + di
be an eigenvalue of the matrix A, where µ is the corresponding eigenvalue ofthe shifted matrix Ai = A − di from which λ is computed. Lete = f l(eλµ + di )ebe the computed eigenvalue. Theorem 1 gives us dependency of accuracy of λin (15) upon accuracy of µe in (16).
e from (15) and µ and µ
Theorem 1. For λ and λe from (16) we have|κλ | ≤

|di | + |µ|(|κµ | + 1) .
|λ|

(19)

Proofs of this theorem and subsequent theorems are given in Appendix A.
e depends on κµ and the size
From Theorem 1 we see that the accuracy of λof the quotient|di | + |µ|
.
(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 (di ) = sign (µ) ,then|di | + |µ|
= 1.
|λ|

(ii) If λ is between two poles of the same sign and sign (di ) 6= sign (µ) (seeFigure 1 (ii)), then|di | + |µ|
≤ 3.
|λ|
11

di+1

µ

µ

z}|{
s

z}|{
s

λ

di

0

di

λ

(i)

di−1

0

(ii)
Figure 1: Typical situations from Theorem 2

Theorem 2 does not cover the following cases:(a) If d1 < 0, then µ > 0. If, further, |d1 | ≈ |µ|, then λ1 is near zero, and(|d1 | + |µ|)/ |λ1 |  1 (see Figure 2 (a)).4(b) If dn > 0, then µ < 0. If, further, |dn | ≈ |µ|, then λn is near zero, andagain (|dn | + |µ|)/ |λn |  1.5(c) If λ is between two poles of the different signs and sign (di ) 6= sign (µ), theneither di+1 < 0 < di and µ < 0, or di < 0 < di−1 and µ > 0. In both cases,
if, additionally, |di | ≈ |µ|, then λ is near zero, and (|di | + |µ|)/ |λ|  1 (seeFigure 2 (c)).

µ

z

}|

µ

{

d1

z }| {
s

s

λ1 0

di+1

(a)

0

λ

di

(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|di | + |µ|≤ 3.
|λ|If one of the above cases does occur, remedies are given in the following
remark.
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.
4 In this case λ is computed as a difference of two close quantities and cancellation can
1occur.
5 In this case λ is computed as a difference of two close quantities and cancellation can
noccur.

12

Notice that the inverse is of an unreduced arrowhead matrix with non-zero shaft
is a diagonal-plus-rank-one (DPR1) matrix of the form
 −1

D
A−1 =
+ ρuuT ,0where
u = z T D−1

−1

T

,

ρ=

1
.
a − z T D−1 z

Eigenvalues of A−1 are zeros of (see [2, 14])ϕ (λ) = 1 + ρ

n
X

u2j
.
d −λ
j=1 j

Since the absolutely largest eigenvalue of A−1 is computed accurately according
to standard perturbation theory, and 1/|λ| = kA−1 k2 , λ is also computed withhigh relative accuracy. In computing matrix A−1 , eventually ρ needs to becomputed in higher precision. For more details see Remark 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.
3.2. Accuracy of the eigenvectors
Since the eigenvector is computed by (12), its accuracy depends on the accuracy of µe as described by the following theorem:
Theorem 3. Let (16) hold and let




−1
(D1 (I + E1 ) − µeI) z1xe1ζi .. 


−


)
x
e =  .  = f l(
µ
e




−1
 (D2 (I + E2 ) − µeI) z2 
x
en
−1

(21)

be the computed un-normalized eigenvector corresponding to µ and λ. Then
xej = xj 1 + εxj , εxj ≤ 3 (|κµ | + 3) εM , j = 1, . . . , n.
In other words, if κµ is small, then all components of the eigenvector aree and xcomputed to high relative accuracy. Since the accuracy of λe depends
on the accuracy of µe (on the size of κµ ) in the next three subsections tells we
discuss the accuracy of µe. Since µ
e is computed as an inverse of the eigenvalueof the matrix f l(A−1
),
we
first
discuss
the accuracy of that matrix.
i13

3.3. Accuracy of the matrix A−1i
We have the following theorem:
Theorem 4. For the computed elements of the matrix A−1
from (10) and (11)ifor all (j, k) 6= (i, i) we have



^
A−1
= f l A−1
= A−1
(1 + εjk ), |εjk | ≤ 3εM .
i
i
i
jk
jk
jk

For the computed element b ≡ A−1
from (17) we haveiii|κb | ≤ 3Kb ,where Kb is defined by (18).
The above theorem states that all elements of the matrix A−1
are computediwith high relative accuracy except possibly b. Therefore, we have to monitor
whether b is computed accurately, and, if not, it needs to be computed in doubleof the working precision (see Section 4 for details).
3.4. Accuracy of bisection
Let λmax be the absolutely largest eigenvalue of a symmetric arrowhead
emax be the eigenvalue computed by bisection as implemented
matrix A, an let λin Algorithm 3. The error bound from [16, Section 3.1] immediately implies
thatemax − λmax
λ√
= κbis εM , κbis ≤ 1.06n n + 1 .
(22)|λmax |
Notice that the similar error bound holds for all eigenvalues which are of thesame order of magnitude as |λmax |.
3.5. Accuracy of exterior eigenvalues of A−1i
The desired interior eigenvalue and, in some cases, also absolutely smaller
exterior eigenvalue λ of A is in Algorithm 1 computed by (8), where ν is one of
the exterior eigenvalues of the matrix A−1
i .
The following theorem covers the case when ν is the absolutely largest eigenvalue of A−1, and gives two different bounds.
i
2
Theorem 5. Let A−1
be defined by (10) and let ν be its eigenvalue such thati.
|ν| = A−1
i
2

(23)



. Let
Let νb be the exact eigenvalue of the computed matrix ^A−1= f l A−1iiνb = ν (1 + κν εM ) .
Then

(24)


n−1

√
√
1 X
|κν | ≤ min 3 n max{1, Kb }, 6 n +|ζk | ,
|ζi |


k=1
k6=i

where Kb is defined by (18).
14

(25)

3.6. Final error bounds
All previous error bounds are summarized as follows.
e be the computed eigenvalue of an unreduced arrowhead maTheorem 6. Let λfi from (9), and let νe be thetrix A, let µe be computed eigenvalue of the matrix Acorresponding computed eigenvalue of the matrix ^A−1 from (10). If µ is theieigenvalue of Ai closest to zero (or, equivalently, if (23) holds), then the error
e is given by (15) within the computed eigenvalue λ|κλ | ≤ 3(|κν | + κbis ) + 4,

(26)

and the error in the computed un-normalized eigenvector x
e is given by Theorem
3 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 guarantee high relative accuracy of the computed eigenvalue and high componentwise
relative accuracy of the computed eigenvector if ν is such that |ν| = O(kA−1
i k2 )and it is computed accurately. This is certainly fulfilled if the following conditions are met:
C1. The quantity Kν from (13) is moderate, andC2. (i) either the quantity Kb from (18) is small, or
n−1
P|ζk | from (25) is of order O(n).
(ii) the quantity |ζ1i |
k=1
k6=i

The condition C1 implies that ν will be computed accurately according to
the standard perturbation theory. The conditions C2 (i) or C2 (ii) imply thatκν from (25) is small, which, together with C1, implies that ν is computed
accurately.
If the condition C1 does not hold, that is, if Kν  1, remedies are given
in 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 Section
4.
Remark 2. If λ = λ1 or λ = λn , it can be computed from the starting matrix
A as in Algorithm 1. We have two more possibilities:(a) we can compute λ by shifting to another neighboring pole provided thatKν is in this case small (shifting to the pole di−1 instead of di in Figure
3 (a)),(b) if shifting to another neighboring pole is not possible (Kν  1, see Figure
3 (b)), we can invert A − σI, where shift σ is chosen near λ, and σ ∈/{λ, di , di−1 }. This results in a DPR1 matrix(D − σI)−1−1(A − σI) =+ ρuuT ,015whereu = z T (D − σI)−1

T
−1 ,

ρ=

1
a − z T (D − σI)−1 z.

Eigenvalues of this matrix are zeros ofϕ (λ) = 1 + ρ

n
X

u2j
,
(dj − σ) − λj=1

and the absolutely largest eigenvalue is computed accurately. Eventually,ρ needs to be computed in higher precision.6

di
s

s

λi+1

λ(λi )

di−1

s

di
s

s

λi−1

λi+1

λ(λi )

(a)

di−1
s
λi−1

(b)
Figure 3: Typical situations from Remark 2

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−1
needs to be computed in higher precision. The following theoremi
implies that if 1  Kb ≤ O(1/εM ), it is sufficient to evaluate (11) in double of
the working precision.7 8Theorem 7. If −a > 0 in (11), set
P = −a + z1T D1−1 z1 ,Q = −z2T D2−1 z2 ,and if −a < 0 in (11) set
P = z1T D1−1 z1 ,Q = a − z2T D2−1 z2 .

6 Determining whether ρ needs to be computed in higher precision is done similarly as
determining whether element b of A−1
needs to be computed in higher precision, which isidescribed in Section 4. Further, Theorem 7 implies that it suffices to compute ρ in double of
the working precision.
7 If K ≥ O(1/ε ), that is, if K = 1/εM
E for some εE < εM , then, in view of Theorem 7,bbb needs to be computed with extended precision εE .
8 Usage of higher precision in conjunction with the eigenvalue computation for DPR1 matrices is analysed 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.

16

Notice that in both cases P, Q ≥ 0 and b = (P − Q)/ζi2 . Let Pe = f l(P ) ande = f l(Q) be evaluated in standard precision, εM . Assume that Pe 6= Q
e andQKb ≤ O(1/εM ). If P , Q and b are all evaluated in double of the working
precision, ε2M , then (17) holds with |κb | ≤ O(n).
We summarize the above results in one, complete algorithm, aheig. The
algorithm 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 is
not the case, the quantity Kb is computed, and if Kb  1, the eigenpair (λ, v)
is computed by Algorithm 1 but with evaluation of b in double of the working
precision. At the end, the quantity Kν is computed, and if Kν  1, one of theremedies from Remark 2 is applied.
Algorithm 5[λ, v] = aheig (D, z, α, k)% Computes the k -th eigenpair of an irreducible arrowhead matrix% A = [diag (D)z; z 0 α]
 n−1P
if the quantity|ζk | /|ζi | from (25) is of O(n)
k=1
k6=i

% standard precision is enough[λ, v] = aheig basic(D, z, α, k)elsecompute the quantity Kb from (18)
if Kb  1% double precision is necessary[λ, v] = aheig basic(D, z, α, k) with evaluation of b in double precisionelse% standard precision is enough[λ, v] = aheig basic(D, z, α, k)
end
endcompute the quantity Kν from (13)if Kν  1apply one of the remedies from Remark 2
end

4.1. On implementing double precision
Implementation of the double of the working precision depends upon whether
the 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- digits(d) specifies the number of significant decimal digits d used to do
variable precision arithmetic vpa,
17

- 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 have
non-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 decimal
number, which is, for example, assumed by extended precision computation in
Mathematica [20]. In this case, the options in Matlab are to either use symbolic
computation, or to cast the input to a string, and then convert it to extended
precision:
a1=vpa(num2str(a,16),32)
In this case, the the decimal interpretation of the variable a1 has all zero
entries after 16-th significant decimal digit, but the binary representation of
the 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) a1
character(25) string...
write(string,*) aread(string,*) a1
If the input consists of numbers for which decimal and binary representation
are equal (for example, integers, as in Example 3 below), then the two above
approaches give the same results.
5. Numerical Examples
We illustrate out algorithm with four numerically demanding examples. Examples 1 and 2 illustrate Algorithm 1, Example 3 illustrates the use of double
precision arithmetic, and Example 4 illustrates and application of higher dimension.

18

Example 1. In this example both quantities Kν from (13) and Kb from (18) arefor all eigenvalues approximately equal to 1, so we guarantee that all eigenvalues
and all components of their corresponding eigenvectors are computed with highrelative accuracy by Algorithm 5 (aheig) using only standard machine precision.
Let
2 · 10−3
0
0
0
0
107

0
10−7 0
0
0
107 



0
0
0
0
0
1 

.
A=
0
0
0 −10−7
0
107 



0
0
0
0
−2 · 10−3 107 
107
107 1
107
107
1020
The eigenvalues computed by Matlab [13] routine eig, Algorithm 5 and Mathematica [20] with 100 digits precision, are, respectively:λ(eig)
1.000000000000000 · 1020
1.999001249000113 · 10−3
4.987562099695390 · 10−9
−1.000644853973479 · 10−20−2.004985562101759 · 10−6
−2.001001251000111 · 10−3λ(aheig)
1.000000000000000 · 1020
1.999001249000113 · 10−3
4.987562099722817 · 10−9−9.999999999980001 · 10−20−2.004985562101717 · 10−6
−2.001001251000111 · 10−3λ(M ath)
1.000000000000000 · 1020
1.999001249000113 · 10−3
4.987562099722817 · 10−9−9.999999999980001 · 10−20−2.004985562101717 · 10−6
−2.001001251000111 · 10−3

We see that even the tiniest eigenvalues λ3 and λ4 , 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 eigenvectors
computed by Algorithm 5 are componentwise accurate up to machine precision,and therefore, orthogonal up to machine precision. For example:
(eig)

v4
4.999993626151683 · 10−11
9.999999962328609 · 10−7
9.999999999990000 · 10−1
−9.999999964673912 · 10−7
−5.000006338012225 · 10−11−9.999999963825105 · 10−21(aheig)v4−4.999999999985000 · 10−11−9.999999999969000 · 10−7
−9.999999999989999 · 10−1
9.999999999970999 · 10−7
4.999999999985000 · 10−11
9.999999999970000 · 10−21(M ath)v4−4.999999999985000 · 10−11−9.999999999969000 · 10−7
−9.999999999989999 · 10−1
9.999999999970999 · 10−7
4.999999999985000 · 10−11
9.999999999969999 · 10−21 .

Example 2. In this example, despite very close diagonal elements, we again
guarantee that all eigenvalues and all components of their corresponding eigenvectors are computed with high relative accuracy, without deflation. Let
1 + 4εM00010
1 + 3εM
0
0
2 



0
0
1 + 2εM
0
3 
A=


0
0
0
1 + 1εM 4 
1
2
3
4
0
where εM = 2 · 2−53 = 2.2204 · 10−16 . For this matrix the quantities Kν
and Kb are again of order one for all eigenvalues, so Algorithm 5 uses only

19

standard working precision. The eigenvalues computed by Matlab, Algorithm
5 and Mathematica with 100 digits precision, are, respectively:λ(eig)
6.000000000000000
1.000000000000001
1.000000000000001
1.000000000000000−5.000000000000000

λ(aheig)
6.000000000000001
1.000000000000001
1.000000000000001
1.000000000000000−4.999999999999999

λ(M ath)
6.000000000000000
1.000000000000001
1.000000000000001
1.000000000000000−5.000000000000000

Although the eigenvalue computed by Matlab appear to be accurate, they arenot. Namely,(aheig)λ2(eig)(aheig)− λ3= 2.220446049250313 · 10−16 ,(eig)while λ2
− λ3
= 0, so the eigenvalues computed by Matlab do not satisfy
the interlacing property. Notice that despite of very close eigenvalues, Algorithm
5 works without deflation. Due to the accuracy of the computed eigenvalues,
the eigenvectors computed by Algorithm 5 are componentwise accurate up to
the machine precision, and are therefore orthogonal. For example:9
(eig)

v3
−4.999999999985000 · 10−1
−9.999999999969000 · 10−2
−9.999999999989999 · 10−2
9.999999999970999 · 10−10(aheig)v3−7.758019882294899 · 10−2
−9.085893014435115 · 10−1
3.53462724159426 · 10−1
2.085926573079233 · 10−1
2.077382436267882 · 10−17(M ath)v3−7.758019882294900 · 10−2
−9.085893014435116 · 10−1
3.53462724159426 · 10−1
2.085926573079233 · 10−1
2.077382436267882 · 10−17

Example 3. In this example we can guarantee all eigenvalues and eigenvectors, componentwise will be computed with high relative accuracy only if we
use double of the working precision when computing b from (11) in matrices
−1
−1
−1
A−1
2 , A3 , A4 and A5 . Let
1010 0 0 0 0 1010
 0
4 0 0 01  0
0 3 0 0
1 
.
A=
 0
0 0 2 01  0
0 0 0 11 1010 1 1 1 1 1010

9 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.2204460492503130808472633361816 · 10−16 .

20The quantities Kν and Kb are:10Kν−
9.999999090793056 · 10−1
9.999996083428923 · 10−1
1.000000117045544 · 100
9.999998561319470 · 10−1
7.941165469988994 · 100

Kb−3.243243243540540 · 109
3.636363636818182 · 109
4.444444445000000 · 109
5.217390439488477 · 109
5.217390439488477 · 109

It is clear, from the condition numbers, that the element b in each of the ma−1−1−1trices A−1
2 , A3 , A4 and A5 needs to be computed in double of the working
precision. For example, 1010 − 4 0 0
0
0
1010


0
0 0
0
0
1




0
0 −1 001.
A2 = A − d 2 I = 0
0 0 −2 0
1




0
0 0
0 −3
1
10
10
10
1 1
1
1 10 − 4The element b = A−1
computed by Algorithm 2 gives b = 6.16666666667,
2
22
Matlab’s routine inv yields b = 6.166665889418350, while computing b in doubleof the working precision gives the correct value b = 6.166666668266667.
Eigenvalues computed by Algorithm 1 (aheig basic, using only standard
working precision), Algorithm 5 (aheig, using double of the working precision tocompute respective b’s) and Mathematica with 100 digits precision, respectively,are:
λaheig basic
2.000000000000000 · 1010
4.150396802313551 · 100
3.161498641452035 · 100
2.188045596352105 · 100
1.216093560005649 · 100−7.160348702977373 · 10−1

λaheig2.000000000000000 · 1010
4.150396802279712 · 100
3.161498641430967 · 100
2.188045596339914 · 100
1.216093584948579 · 100−7.160346250991725 · 10−1

λM ath2.000000000000000 · 1010
4.150396802279713 · 100
3.161498641430967 · 100
2.188045596339914 · 100
1.216093584948579 · 100−7.160346250991725 · 10−1

The eigenvectors computed by Algorithm 5 are componentwise accurate to machine precision and therefore orthogonal.
Example 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,- di is a frequency of the i-th optical mode, and
10 Algorithm 5 does not compute K and K for the first eigenvalue, since it is an absolutelyνblargest one.

21

- ζ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 103 and104 . We ran a test example for n = 2501 where, typically,di ∈ [5.8769928900036225 · 1014 , 1.3709849013800450 · 1015 ],ζi ∈ [3.2087698694339995 · 106 , 4.9584253488898976 · 106 ],α = 9.7949881500060375 · 1014 .
For this matrix the condition number Kν ∼ 1 for all eigenvalues, and the components of the vector z are all of the same order of magnitude. Therefore,
conditions C1 and C2 (ii) from Section 3 are fulfilled, so all eigenvalues andall components of all eigenvectors are computed with high relative accuracy by
Algorithm 5 using only standard working precision. On the other hand about
half of the eigenvalues computed by the Matlab routine eig do not satisfy theinterlacing property.
6. Applications
In this section we extend our results to eigenvalue decompositions of Hermitian arrowhead matrices, singular value decompositions of real triangular arrowhead matrices and eigenvalue decompositions of real symmetric diagonalplus-rank-one matrices.
6.1. Hermitian arrowhead matrices
Let

C=

D
z∗

z
α


,

where
D = diag(d1 , d2 , . . . , dn−1 ),is a real diagonal matrix of order n − 1,z = ζ1 ζ2 · · ·

ζn−1

∗

,

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 eigenvalue
decomposition of C is given by
C = U ΛU ∗where Λ = diag(λ1 , . . ., λn ) ∈ Rn×n is a diagonal matrix of eigenvalues, andU = u1 u2 · · · un 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:22A = Φ∗ CΦ =where
Φ = diag



D
|z|T

|z|
α


,

ζ1 ζ2
ζn−1
,
,...,
,1
|ζ1 | |ζ2 ||ζn−1 |(28)


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 3. Similarly, for irreducible non-symmetric arrowhead matrix


D z̆
G=
,
z̊ T αwhere sign(ζ̊i ) = sign(ζ̆i ), i = 1, . . . , n − 1, we define the diagonal matrix


s
s
ζ̊
ζ̊
1
n−1
, . . . , sign(ζ̊n−1 ), 1 .
Ψ = diag sign(ζ̊1 )
ζ̆1
ζ̆n−1
The matrix
−1

A=Ψ


GΨ =

D
zT

z
α


,

q
where ζi = ζ̆1 ζ̊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 needed
to compute b in Algorithm 5, the elements ζi need to be computed in double of
the working precision, as well.
6.2. Singular value decomposition of a triangular arrowhead matrix
Let


D z
B=
,
0 α
be an irreducible upper triangular arrowhead matrix, that is, di 6= dj for i 6= jand ζi 6= 0 for all i. The matrix


D2
Dz
A = BT B =,zT D α + zT z
is an irreducible symmetric arrowhead matrix.
23

When applying Algorithm 5 to the matrix A, we must ensure that all components of A−1
in (10) are computed to high relative accuracy. This is obviouslyitrue for elements of the vectors wi and w2 . Diagonal elements, except b, arecomputed with high relative accuracy as differences of squares of original quantities,1, j 6= i.
[A−1i ]jj =
(dj − di )(dj + di )11
The element b = [A−1i ]ii from (11) is computed as


2 2
X
ζ
d
1 
j j
.
b = 2 2 −α − z T z + di +di ζi(dj − di )(dj + di )j6=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(σ1 , . . . , σn ) are the singular values, the columns of V are the corresponding
right singular vectors and the columns of U are the corresponding left singular
vectors.
√ 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 corresponding
right 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


 T
 D z
 T

u1:n−1 unvn ,= σ v1:n−10 αwhich implies
u1:n−1 = σv1:n−1 D−1 .
From the relation BV = U Σ for the k-th column we have





D z v1:n−1u= σ 1:n−1 ,
0 α
vn
un
which implies

αvn.
σComponents of u are computed by multiplication and division of quantities
which are accurate to almost full machine precision, so the are accurate toalmost full machine precision, as well.
un =6.3. Diagonal-plus-rank-one matricesLet
M = D + uuT ,11 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.

24

whereD = diag(d1 , . . . , dn ), d1 > d2 > · · · > dn ,Tu = u1 · · · un , ui =6 0, i = 1, . . . , n,be a n × n irreducible ordered real symmetric diagonal-plus-rank-one (DPR1)
matrix. Let
D̄ = diag(d1 , . . . , dn−1 ),∆ = (D̄ − dn )1/2 ,Tū = u1 · · · un−1 ,un ∆−1 0.
L=−ūT ∆−1 1
Then
−1

A=L



D̄
ML = T
z


z
,
α

where
z = ū∆,α = dn + uT u,

is an irreducible real symmetric arrowhead matrix.
When applying Algorithm 5 to the matrix A, we must ensure that all components of A−1
in (10) are computed to high relative accuracy. This is obviouslyitrue for elements of the vectors wi and w2 . Diagonal elements, except b, arecomputed with high relative accuracy as differences of original quantities, andthe element b = [A−1i ]ii from (11) is computed as


X ζj
1 
.
b=
−dn − uT u + di +ζidj − dij6=i

If double precision is needed in Algorithm 5, all entries of A need to be computedin double precision.
Let M = QΛQT and A = V ΛV T be the eigenvalue decompositions of M
and A, respectively. Since M is by assumption irreducible, its eigenvalues satisfy
interlacing property
λ1 > d1 > λ2 > d2 > · · · > λn > dn .

(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−1 M LV = Λ implies that the columns of the matrix X = LV
are the unnormalized eigenvectors of the matrix M . Further, since, by (29), alleigenvalues are simple, we conclude that X = QΣ, where Σ = diag(σ1 , . . . , σn )is a positive definite matrix. Notice that QΣV T = L is, in fact, singular value
decomposition of L.
25

Equating k-th columns of the equation X = LV gives
 

 
un ∆−1 0 v̄x̄x=
= Lv =,
xn−ūT ∆−1 1 vn
where x and v are partitioned according to L. This immediately implies thatx̄ = un ∆−1 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 xn accurately.
Let
 
q̄
q=
qn
be the k-th column of Q and let σ = Σkk . Equating k-th rows of the equation
X −1 = Σ−1 QT = V T L−1
gives for the n-th elementqn

11= x n 2 = vn .
σ
σ

Thus,
x n = σ 2 vn
and, in order to compute xn , it is necessary to compute σ 2 . From X = U Σ = LV
it follows that V T LT LV = Σ2 , or, equivalently, LV = L−T V Σ2 . Equating k-th
columns of this equation gives
1∆−1 v̄ un = ∆v̄
+ ∆−1 ū vn σ 2 .
unThis gives n − 1 equations for σ 2 , and we can choose the numerically mostaccurate one.
Therefore, xn will be computed to almost full machine precision, as are theentries of x̄, and it remains to normalize x and obtain q = x/σ.
Remark 4. Notice that DPR1 matrices of the form D −uuT cannot be reduced
to symmetric arrowhead matrix by the procedure described in this section. By
using ideas from this paper, it is possible to derive highly accurate algorithm
for DPR1 matrices without prior transformation to arrowhead form. This algorithm, which is a topic of our forthcoming paper, covers more general DPR1
matrices of the form
D + ρuuT , ρ ∈ R.

26

Appendix A. Proofs
Proof of Theorem 1.
e be defined by (14). ThenLet µe and λe ≡ f l (di + µλe) = (di + µe) (1 + ε1 ) .
By simplifying the equality(di + µ (1 + κµ εM )) (1 + ε1 ) = λ (1 + κλ εM )and using λ = µ + di , we havedi ε1 + µ (κµ εM + ε1 ) = λκλ εM .
Taking absolute value gives|κλ | ≤

|di | + |µ|(|κµ | + 1) . |λ|Proof of Theorem 2.
(i) The assumption sign (di ) = sign (µ) immediately implies|di + µ||di | + |µ|
=
= 1.
|λ|
|di + µ|(ii) The assumptions imply that either
0 < di+1 < λ < di ,µ < 0,

di < λ < di−1 < 0,

µ > 0.

orIn the first case λ is closest to the pole di and|di | + 12 |di − di+1 |di + 12 di − 12 di+1
|di | + |µ|
≤
≤
1
1
1
|λ|
2 |di + di+1 |2 di + 2 di+1
3

di − 21 di+1
3di
≤
= 3.
1
di
d
+
d
i
i+1
2
2

≤ 21

Here we used the inequalities |µ| ≤ 12 |di − di+1 | and |λ| ≥ 12 |di + di+1 | for the
first inequality, di − di+1 > 0 and di + di+1 > 0 for the second inequality anddi+1 > 0 for the fourth inequality, respectively.
The proof for the second case is analogous. 27

Proof of Theorem 3.
Let x and x
e be defined by (12) and (21), respectively. The theorem obviously
holds for xn = x
en = −1. For xei we have


ζi
ζi
x
ei = f l −=−(1 + ε1 ) = xi (1 + εxi ) .
µeµ (1 + κµ εM )By using (16) and (21), the first order approximation gives|εxi | ≤ (|κµ | + 1) εM .
For j ∈
/ {i, n}, by solving the equality
x
ej =

ζj
ζj
(1 + ε3 ) =(1 + εx )((dj − di ) (1 + ε1 ) − µ (1 + κµ εM )) (1 + ε2 )dj − λfor εx , using (16) and λ = µ + di , and ignoring higher order terms, we haveεx =(dj − di ) (ε1 + ε2 + ε3 ) − µ (κµ εM + ε2 + ε3 ).
dj − λTherefore,|εx | ≤

|dj − di | + |µ|(|κµ | + 3) εM .
|dj − λ|(A.1)

To complete the proof we need to analyze two cases. If
sign (dj − di ) = − sign µ,
then|dj − di − µ||dj − λ||dj − di | + |µ|=== 1.
|dj − λ||dj − λ||dj − λ|Ifsign (dj − di ) = sign µ,then, since di is pole closest to λ, we have |µ| ≤ 0.5 |dj − di | and3|dj − di ||dj − di | + |µ||dj − di | + |µ|
≤
≤ 21
= 3.
|dj − λ||dj − di | − |µ|
2 |dj − di |Finally, the theorem follows by inserting this into (A.1).



Proof of Theorem 4.
For the non-zero computed elements of the matrix A−1
from (10) and (11)i
we have:
1)=(1 + ε2 ) , j ∈/ {i, n},f l( A−1ijj(dj − di ) (1 + ε1 )




−ζj
f l( A−1) = f l( A−1)=(1 + ε5 ) , j ∈/ {1, n},iiji
ij
(dj − di ) (1 + ε3 ) ζi (1 + ε4 )1−1
f l([A−1(1 + ε6 ),i ]ni ) = f l([Ai ]in ) =ζi28

where |εk | ≤ εM for all indices k. The first statement of the theorem now follows
by using standard first order approximations.
Similar analysis of the formula (11) yieldsef l([A−1
i ]ii ) = b = b + δb,
where|δb| ≤

3|−a| + z1T D1−1 z1 + z2T D2−1 z2 εM .
2
ζi

(A.2)

This, in turn, implies (17) with|κb | ≤

|a| + z1T D1−1 z1 + z2T D2−1 z2|δb| 1=3= 3Kb ,
|b| |εM |−a + z1T D1−1 z1 + z2T D2−1 z2where Kb is defined by (18).



Proof of Theorem 5.
Let

−1
^
A−1
= A−1
i
i + δAi ,or, equivalently,
νb − ν = δA−1
i .
From this, (24), and Theorem 4, we have−1|νκν εM | ≤ kδA−1
i k2 ≤ 3k|Ai |k2 max{1, Kb }εM .
√nkA−1
Since k|A−1
i k2 , we havei |k2 ≤√
|κν | ≤ 3 n max{1, Kb }.

(A.3)

For the second part of the proof, we may assume that the element a =α − di is not perturbed. More precisely, we represent the perturbation in aas perturbations in D1 and D2 by pre- and post-multiplication of the matrix
ei = f l(Ai ) with the diagonal matrix
A
"
#
1/2
(1 + εa ) In−1
0
.
−1/2
0
(1 + εa )This introduces additional errors in D1 and D2 bounded by εM . Therefore, we
can write−1|νκν εM | ≤ kδA−1(A.4)
i k2 ≤ 6k|Ai |k2 + |δb|.
By modifying (A.2) to accommodate lack of perturbation in a, we have|δb| ≤

6
ζi2

z1T D1−1 z1 + z2T D2−1 z2

29



εM .

By inserting this bound into (A.4) and rearranging, we have√
1 z1T D1−1 z1 + z2T D2−1 z2|κν | ≤ 6
n+
.
|ν|
ζi2

(A.5)

Since
−1
A−1
= |ν| = max A−1
i x 2 ≥ Ai e k 2i2
kxk2 =1
s
ζk2
1
|ζk |
=
2 ≥ |ζ | |d − d | ,
2 + 2
i
k
i
(dk − di )ζi (dk − di )by simply dividing each term
ζk22
ζi |dk − di |in (A.5) with the corresponding quotient|ζk |,|ζi | |dk − di |we obtain|κν | ≤ 6


√


n−1
1 X
|ζk | .
n+
ζi

(A.6)

k=1
k6=i

The bound (25) now follows from (A.3) and (A.6).



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. The
bound (26) now follows by inserting (27) into Theorems 1 and 2. 
Proof of Theorem 7.
Let the assumptions of the theorem hold. Let b be computed in double of theworking precision, ε2M , and then stored in the standard precision. The standard
floating-point error analysis with neglecting higher order terms gives P −Q
P 1 + κP ε2M − Q 1 + κQ ε2M
1 + κ1 ε2M =(1 + κb εM )2ζiζi2
≡ b (1 + κb εM ) ,30

where |κP |, |κQ | ≤ (n + 1) and |κ1 | ≤ 3. Solving the above equality for κb ,neglecting higher order terms, and taking absolute values gives|κb | ≤

|P | + |Q|(n + 4) εM ≡ Kb (n + 4)εM .
|P − Q|

Since, by assumption, Kb ≤ O(1/εM ), this implies|κb | ≤ O(n),as desired.



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 eigenvalue
problem, 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 Method
for the Generalized Real Symmetric Definite Tridiagonal Eigenproblem, in
Numerical Linear Algebra and Scientific Computing, L. Reichel, A. Ruttan
and 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 tridiagonal 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 Theory, Phys. Rev. B, 24 (1981) 1651-1663.
[9] D. Goldberg, What Every Computer Scientist Should Know About
Floating-Point Arithmetic, ACM Computing Surveys, 23:1 (1991) 5-48.
[10] G. H. Golub and C. F. Van Loan, Matrix Computations, The John Hopkins
University Press, Baltimore, 3rd ed. (1996).
[11] M. Gu and S. C. Eisenstat, A divide-and-conquer algorithm for the symmetric tridiagonal eigenproblem, SIAM J. Matrix Anal. Appl., 16 (1995)
79-92.

31

[12] Intel Fortran Compiler, http://software.intel.com/en-us/fortran-compilers[13] MATLAB. The MathWorks,
http://www.mathworks.com.

Inc.,

Natick,Massachusetts,USA,[14] A. Melman, Numerical solution of a secular equation, Numer. Math., 69(1995) 483-493.
[15] D. Mogilevtsev, A. Maloshtan, S. Kilin, L. E. Oliveira and S. B. Cavalcanti,
Spontaneous emission and qubit transfer in spin-1/2 chains, J. Phys. B: At.
Mol. Opt. Phys., 43 095506 (2010).
[16] D. P. O’Leary and G.W. Stewart, Computing the eigenvalues and eigenvectors of symmetric arrowhead matrices, J. Comput. Phys. 90, 2 (1990)
497-505.
[17] S. Oliveira, A new parallel chasing algorithm for transforming arrowhead
matrices to tridiagonal form, Math. Comp., 221 (1998) 221–235.
[18] B. N. Parlett, The Symmetric Eigenvalue Problem, Prentice-Hall, Englewood Cliffs, (1980).
[19] J. H. Wilkinson, The Algebraic Eigenvalue Problem, Clarendon Press, Oxford, (1965).
[20] Wolfram Mathematica, Documentation Center,
http://reference.wolfram.com/mathematica/guide/Mathematica.html

32