# Akiko Fukuda

Akiko Fukuda Received on October ... ilarity transformation by D yields L ... we compare our routine dhLV, which is a programming code of the dhLV alg...

Journal of Math-for-Industry, Vol. 3 (2011A-4), pp. 47–52

A note on eigenvalue computation for a tridiagonal matrix with real eigenvalues D案

Akiko Fukuda Received on October 6, 2010 / Revised on February 7, 2011 Abstract. The target matrix of the dhLV algorithm is already shown to be a class of nonsymmetric band matrix with complex eigenvalues. In the case where the band width M = 1 in the dhLV algorithm, it is applicable to a tridiagonal matrix, with real eigenvalues, whose upper and lower subdiagonal entries are restricted to be positive and 1, respectively. In this paper, we ﬁrst clarify that the dhLV algorithm is also applicable to the eigenvalue computation of nonsymmetric tridiagonal matrix with relaxing the restrictions for subdiagonal entries. We also demonstrate that the wellknown packages are not always desirable for computing nonsymmetric eigenvalues with respect to numerical accuracy. Through some numerical examples, it is shown that the tridiagonal eigenvalues computed by the dhLV algorithm are to high relative accuracy. Keywords. matrix eigenvalues, tridiagonal matrix, discrete hungry Lotka-Volterra system

1.

are real. Additionally, in , the upper subdiagonal entries are required to be positive, in order to guarantee convergence of the dhLV algorithm. It is noted here that such tridiagonal matrices always have real eigenvalues. The ﬁrst purpose of this paper is to expand the applicable range of the dhLV algorithm with M = 1. The target matrix of the dhLV algorithm with M = 1 is shown to be the nonsymmetric tridiagonal matrix whose upper and lower subdiagonal entries are not restricted to be positive and 1, respectively. The second is to demonstrate that the computed eigenvalues by the dhLV algorithm are to high relative accuracy, which is not discussed precisely in , in the case where the eigenvalues are both real and complex, through some numerical examples. In this paper, we ﬁrst review in Section 2 the dhLV algorithm with M = 1 brieﬂy, and then we clarify that it is also applicable to the eigenvalue computation of nonsymmetric tridiagonal matrix which is not a class of the target matrix in . In Section 3, by some numerical examples, we observe that MATLAB and LAPACK are not always desirable with respect to the numerical accuracy of the computed eigenvalues. In the comparison with the numerical results by MATLAB and LAPACK, the computed eigenvalues by the dhLV algorithm are shown to have almost high relative accuracy. Finally, in Section 4, we give concluding remarks.

Introduction

Several routines for nonsymmetric eigenvalues are provided in MATLAB  which is an interactive software for matrixbased computation, LAPACK  which is the famous numerical linear algebra package, and so on. The QR algorithm [3, 4] is the most standard algorithm for nonsymmetric eigenvalues, so is adopted as the LAPACK routines. However, it is not diﬃcult to ﬁnd an example such that the eigenvalues computed by LAPACK are not to high relative accuracy. The author designs in  an algorithm, named the dhLV algorithm, for complex eigenvalues of a certain nonsymmetric band matrix. The dhLV algorithm is based on the integrable discrete hungry Lotka-Volterra (dhLV) system (1)

(n+1)

uk

(n)

= uk

M ∏

1 + δ (n) uk+j

j=1

1 + δ (n+1) uk−j

(n)

(n+1)

,

k = 1, 2, . . . , Mm ,

Mk := (k − 1)M + k,

(n) u1−j

≡ 0,

≡ 0,

(n) uMm +j

j = 1, 2, . . . , M, (n)

where M , m are given positive integers, and δ (n) > 0, uk denote the values of δ and uk at the discrete step n, respectively. The dhLV system (1) is a discretized version of the continuous-time hungry Lotka-Volterra system , which is the prey-predator model in mathematical biology. The parameter M originally denotes the number of species which a species can prey, whereas, in the dhLV algorithm, M gives the location of positive entries in the target matrix. In the case where M = 1, the target matrix of the dhLV algorithm becomes a nonsymmetric tridiagonal matrix whose lower subdiagonal entries are 1 and eigenvalues

2. The dhLV algorithm for nonsymmetric tridiagonal matrix We ﬁrst explain the dhLV algorithm with M = 1.

47

E案

48

Journal of Mathematics for Industry, Vol. 3 (2011A-4)

Let us introduce two kinds of 2m × 2m matrices,   (n) 0 U1   (n)   1 0 U2     .. .. (2) L(n) :=  , . . 1     .. .. (n)  . . U2m−1  1 0  (n)  V1   (n) V2  0    ..  (n)  (n) . R :=  δ (3) , 0     .. .. .. . . .   (n) (n) δ 0 V2m where (n)

:= uk (1 + δ (n) uk−1 ),

(n)

:= (1 + δ (n) uk )(1 + δ (n) uk−1 ).

Uk Vk

(n)

(n)

(n)

(n)

Then, the matrix representation for the dhLV system with M = 1 is given by R(n) L(n+1) = L(n) R(n) .

(4)

The equality in each entry of (4) leads to the dhLV system (0) (0) (1) with M = 1. Let us assume that u1 > 0, u2 > (0) 0, . . . , u2m−1 > 0, then it is obvious from (1) that, for all n, (n) (n) (n) u1 > 0, u2 > 0, . . . , u2m−1 > 0. By taking account that (n) Vk > 1, we see that R(n) is nonsingular. So, (4) can be transformed as L(n+1) = (R(n) )−1 L(n) R(n) . This implies that the eigenvalues of L(n) are invariant under the time evolution from n to n+1 of the dhLV system (1) with M = 1. Hence, the matrices L(0) and L(1) , L(2) , . . . are similar to each other. For the unit matrix I and an arbitrary constant d, the matrices L(0) + dI and L(1) + dI, L(2) + dI, . . . are also similar. The asymptotic behavior as n → ∞ of the dhLV variables are given as

Hence, the characteristic polynomial of L∗ becomes ∗

det(L − λI) =

(5) (6)

(n)

= ck ,

lim u2k = 0,

n→∞

k = 1, 2, . . . , m,

k = 1, 2, . . . , m − 1,

where c1 , c2 , . . . , cm are positive constants such that c1 > c2 > · · · > cm . See  for the proof of (5) and (6). From (5) and (6), we see that the characteristic polynomial of L∗ := limn→∞ L(n) + dI coincides with that of the block diagonal matrix diag(L1 , L2 , . . . , Lm ), where Lk is the 2 × 2 matrix ( ) d ck Lk = . 1 d

[(d − λ)2 − ck ].

k=1

Consequently, 2m eigenvalues of L(0) + dI are given by d±

(7)

ck ,

k = 1, 2, . . . , m.

Namely, all the eigenvalues are real for nonsymmetric tridi(0) agonal matrix L(0) + dI. If each entry Uk in L(0) + dI (0) is given, the initial value uk in the dhLV system (1) with (0) (0) M = 1 is set as Uk /(1 + δ (n) uk−1 ). Since, for suﬃciently √ (N ) (n) large N , u2k−1 is an approximation of ck , d + uMk leads to the approximation of the eigenvalues of L(0) + dI. We here expand the applicable range of the dhLV algorithm. Let us introduce the diagonal matrix D := diag(1, α1 , α1 α2 , . . . , (α1 α2 · · · α2m−1 )), with arbitrary constants α1 , α2 , . . . , α2m−1 . Then the similarity transformation by D yields ˆ (n) + dI := D(L(n) + dI)D−1 L  ˆ (n) d U 1  ˆ (n)  α1 d U 2   .. = (8) . α2   ..  .

..

.

..

.

α2m−1

ˆ (n) U 2m−1 d

    ,   

ˆ (n) = U (n) /αk . Obviously, the eigenvalues of L ˆ (n) where U k k (n) + dI coincide with those of L + dI. In other words, the applicable range of the dhLV algorithm with M = 1 is the ˆ (0) + dI are matrix given as (8). Hence the eigenvalues of L (0) (0) (0) given as (7), if the initial values U1 , U2 , . . . , U2m−1 are ˆ (0) , U ˆ (0) , . . . , U ˆ (0) and α1 , α2 , set, in accordance with U 1 2 2m−1 . . . , α2m−1 , as (0)

(n) lim u n→∞ 2k−1

m ∏

Uk

(0)

ˆ αk , =U k

k = 1, 2, . . . .2m − 1.

Note that the convergence of the dhLV algorithm is guaran(0) ˆ (0) αk > 0 for k = 1, 2, . . . , 2m − teed if Uk > 0, namely, U k ˆ (0) αk > 0 implies that U ˆ (0) and αk have 1. The positivity U k k the same sign for each k. It is remarkable that the nonˆ (n) + dI with U ˆ (0) αk > 0 symmetric tridiagonal matrix L k has only real eigenvalues. In Table 1, we present the dhLV algorithm for the eigenˆ (0) + dI. The values of nonsymmetric tridiagonal matrix L lines from the 7th to the 11th are repeated until maxk u2k ≤ eps or n > nmax is satisﬁed for suﬃciently small eps > 0.

3.

Numerical experiments

In this section, we show some numerical results.

49

Akiko Fukuda

Table 1: The dhLV algorithm with M = 1. 01: 02: 03: 04: 05: 06: 07: 08: 09: 10: 11: 12: 13: 14:

for k := 1, 2, . . . , 2m − 1 do (0) ˆ (0) αk Uk = U k end for for k := 1, 2, . . . , 2m − 1 do (0) (0) (0) uk = Uk /(1 + δ (0) uk−1 ) end for for n := 1, 2, . . . , nmax do for k := 1, 2, . . . , 2m − 1 do (n+1) (n) (n) (n+1) uk := uk (1 + δ (n) uk+1 )/(1 + δ (n+1) uk−1 ) end for end for for k := 1,√ 2, . . . , m do λk = d ± end for

(n+1) u2k−1

whose eigenvalues are theoretically given by √ kπ , 2 ℓ cos 101

15 10 5 0 −5 −10 −15 −20 −40

−20

0

20

40

Figure 1: A graph of the real part (x-axis), and the imaginary part (y-axis) of the eigenvalues of T1 , given by (9) (plotted by ×), computed by eig (plotted by △) and computed by dhseqr (plotted by ¤).

Totally nonnegative (TN) matrix is a class of nonsymmetric matrices whose eigenvalues are computable to high relative accuracy. Here TN matrix is a nonsymmetric matrix with real and positive eigenvalues. Koev proposes in  an algorithm for computing eigenvalues of a TN matrix. Watkins claims in  “This is the ﬁrst example of a class of (mostly) nonsymmetric matrices whose eigenvalues can be determined to high relative accuracy”. In other words, it is not easy to compute eigenvalues of nonsymmetric matrices to high relative accuracy, except for a TN matrix. The readers should pay attention to that the following example matrices are not TN. Numerical experiments have been carried out on our computer with CPU: Intel (R) CPU L2400 @ 1.66GHz, RAM: 2GB. The dhLV algorithm is implemented by the compiler: Microsoft(R) C/C++ Optimizing Compiler Version 15.00.30729.01. We also use MATLAB R2009b (Version 7.9.0.529) and LAPACK-3.2.1 with complier: gcc-4.3.2. In this section, with respect to the numerical accuracy of computed eigenvalues, we compare our routine dhLV, which is a programming code of the dhLV algorithm, with the MATLAB routine eig and the LAPACK routine dhseqr for nonsymmetric Hessenberg matrix and dsterf for symmetric tridiagonal matrix. According to , the MATLAB routine eig is based on the LAPACK routine dgeev. In dhLV, we set δ (n) = 1.0 for n = 0, 1, . . . . Example 1. The ﬁrst example is the 100 × 100 nonsymmetric matrix   0 1  ℓ 0 1      . . .. .. , ℓ T1 =      . . .. .. 1   ℓ 0

(9)

20

k = 1, 2, . . . , 100.

5

10

0

10

−5

10

−10

10

−15

10

−20

10

0

20

40

60

80

100

Figure 2: A graph of the index k of the computed eigenˆ k of T1 (x-axis), and the relative error rk (y-axis) value λ in eig (plotted by △), dhseqr (plotted by ¤) and dhLV (plotted by ⃝).

Figure 1 shows the eigenvalues of T1 with ℓ = 100, given by (9), computed by eig, and computed by dhseqr in double precision arithmetic. Obviously, all the eigenvalues of T1 with ℓ = 100 are real. We, however, observe that some computed eigenvalues by eig and dhseqr are not real. The routines eig and dhseqr give rise to 30 and 16 complex eigenvalues, respectively. Only in dhLV, all the computed eigenvalues are real. And, if the computed eigenvalues by dhLV are plotted by ◦ in Figure 1, the mark ◦ approximately overlaps with ×, given in (9). So, in Figure 2, we clarify ˆ k − λk |/|λk | where λk with |λ1 | > the relative error rk := |λ ˆ k with |λ ˆ 1 | > |λ ˆ 2 | > · · · > |λ ˆ 100 | |λ2 | > · · · > |λ100 | and λ denote the eigenvalues given by (9) and computed by eig, dhseqr or dhLV, respectively. Not only Figure 1 but also Figure 2 imply that the relative errors r1 , r2 , . . . , r100 in eig and dhseqr are not small. Figure 2 also demonstrates that r1 , r2 , ..., r100 in dhLV are much smaller than those in

50

Journal of Mathematics for Industry, Vol. 3 (2011A-4)

Table 2: The number of irrelevant complex eigenvalues. ℓ eig dhseqr

10−10

10−5

10−1

100

101

105

1010

26 26

22 44

0 14

0 0

76 92

20 76

20 52

Table 3: The average of the relative error rave . ℓ

eig

10−10 10−5 10−1 100 101 105 1010

1.18 × 100 1.11 × 100 1.48 × 10−8 6.07 × 10−16 2.03 × 100 1.07 × 100 1.14 × 100

dhseqr 1.28 × 100 1.13 × 100 2.88 × 100 1.66 × 10−15 2.62 × 100 1.39 × 100 1.58 × 100

dhLV 8.60 × 10−10 4.94 × 10−13 1.85 × 10−15 1.40 × 10−15 2.42 × 10−15 1.31 × 10−15 2.14 × 10−15

ple, without changing the eigenvalues, T1 with ℓ = 100 is symmetrized as T˜1 := D1 T1 (D1 )−1  0 10  10 0 10   . 10 . . =   ..  .

−14

10

0

20

40

60

80

100

Figure 3: A graph of the index k (x-axis) of the computed ˆ k , and the relative error rk (y-axis) in eig eigenvalue λ (plotted by △) and dhseqr (plotted by ¤).

both eig and dhseqr. Moreover, let us consider the cases where ℓ = 10−10 , −5 10 , 10−1 , 100 , 101 , 105 and 1010 . Let us recall here that T1 with positive ℓ has only real eigenvalues. Table 2 gives the number of irrelevant complex eigenvalues by eig and dhseqr. Similarly to the case where ℓ = 100, some computed eigenvalues by eig and dhseqr are complex in every case. Of course, all the computed eigenvalues by dhLV are not complex. ∑100 Table 3 lists the average of relative errors rave := k=1 rk /100 in eig, dhseqr and dhLV. Except in the case where ℓ = 100 , rave in eig and dhseqr are not small. Though rave in dhLV tends to be larger as ℓ becomes smaller, the eigenvalues by dhLV have small relative errors in the comparison with those by eig and dhseqr. It is hence concluded that the eigenvalues of T1 are computed by dhLV with high relative accuracy. It is worth noting that T1 with positive ℓ becomes the symmetric matrix by similarity transformation. For exam-

. 10

1 0

..

.

1

..

.

..

.

10

−16

.

..

   ,   10  0

example is the 100 × 100 non-

−15

10

..

√ √ where D1 = diag(1, 10−2 , . . . , 10−198 ). All the computed eigenvalues by eig and dsterf are real in the case where T˜1 is the target matrix. As is shown in Figure 3, the computed eigenvalues by eig and dsterf are also high relative accuracy. This is an example such that the symmetrization process is useful to improve the numerical accuracy of computed eigenvalues. Example 2. The second symmetric matrix  0  8  10   T2 =    

−13

10

1 .. . 1

    ,   1  0

whose two eigenvalues are ±10000 and the others have the absolute values in [0.03, 2.0]. In the case where T2 is the target matrix, the computed eigenvalues by eig and dhseqr become real without the symmetrization process. Figure 4 shows that the relative errors in dhLV and eig are almost smaller than those in dhseqr. On the other hand, by imposing the symmetrization for T2 beforehand and using dsterf which is for symmetric tridiagonal matrices, the relative errors in dsterf are small, whereas those in eig become larger, which is shown in Figure 5. Example 3. The third example is  0 1  8 ..  10 . 0   .. ..  . . 1  8  10 0 1 T3 (k) =    1 0   . ..  |

{z

}

.. ..

.

. 1

      ,      1  0

k

whose matrix size is 100. The relative errors in eig and dsterf are evaluated only in the case where T3 (k) is symmetrized beforehand. First, let us consider the case where k = 98. Two eigenvalues of T3 (98) are ±0.14, and the others have the absolute values in [628, 19990]. Since the relative errors r97 , r98 , r99

51

Akiko Fukuda −8

10

10

−10

−10

10

10

−12

−12

10

10

−14

−14

10

−16

10

0

10 20

40

60

80

100

Figure 4: A graph of the index k (x-axis) of the computed ˆ k , and the relative error rk (y-axis) in eig eigenvalue λ (plotted by △), dhseqr (plotted by ¤) and dhLV (plotted by ⃝). −11

10

−16

0

20

40

60

80

100

Figure 6: A graph of the index k (x-axis) of the computed ˆ k , and the relative error rk (y-axis) in eig eigenvalue λ (plotted by △), dsterf (plotted by ¤) and dhLV (plotted by ⃝). 10

−10

−12

10

10

−13

−12

10

−14

10

10

−14

−15

10

−16

10

0

10 20

40

60

80

100

Figure 5: A graph of the index k (x-axis) of the computed ˆ k , and the relative error rk (y-axis) in eig eigenvalue λ (plotted by △) and dsterf (plotted by ¤).

and r100 in dhLV are smaller than 1.0 × 10−16 , so they are not plotted in Figure 6. We observe from Figure 6 that r1 , r2 , . . . , r98 in eig and dsterf are suﬃciently small, whereas r99 and r100 are not small. By comparing with Figures 2 and 4, we also see that the relative errors in dhLV are small similar to the case where T1 and T2 are the target matrices. Next, let k = 50. The matrix T3 (50) has the absolute values of 50 eigenvalues in [1207, 19964] and the others in [0.03, 2.0]. Figure 7 claims that almost half of the computed eigenvalues by eig and dsterf are not high relative accuracy, even in the case of employing the symmetrization ˆ1, λ ˆ2, . . . , λ ˆ 100 with process. The routine dhLV brings to λ r1 , r2 , . . . , r100 < 10−14 . It is numerically veriﬁed that the relative errors in dhLV are smaller than O(10−14 ) in almost every example. It is also observed that there are some cases where the relative errors in eig, dhseqr and dsterf are larger than O(10−12 ).

−16

0

20

40

60

80

100

Figure 7: A graph of the index k (x-axis) of the computed ˆ k , and the relative error rk (y-axis) in eig eigenvalue λ (plotted by △), dsterf (plotted by ¤) and dhLV (plotted by ⃝).

It is therefore concluded that the dhLV algorithm is better than MATLAB and LAPACK in order not to get the computed eigenvalues with larger relative errors.

4.

Concluding remarks

In this paper, we ﬁrst review the dhLV algorithm with M = 1 which is designed from the dhLV system, and then we expand the applicable class of nonsymmetric tridiagonal matrix with real eigenvalues. Even in the case where the MATLAB and the LAPACK routines are not desirable with respect to numerical accuracy of computed eigenvalues, the dhLV algorithm enables us to compute eigenvalues with high relative accuracy.

52

Journal of Mathematics for Industry, Vol. 3 (2011A-4)

Acknowledgments The author would like to thank Dr. E. Ishiwata, Dr. M. Iwasaki, Prof. Y. Nakamura, Prof. H. Hasegawa, Prof. G. L. G. Sleijpen and Mr. K. Yadani for many fruitful discussions and helpful advices in this work.

References  The MathWorks, Inc.: MATLAB Reference Guide, Natick, MA, 1992.  Anderson E., Bai Z., Bischof C., Demmel J., Dongarra J., Du Croz J., Greenbaum A., Hammarling S., McKenney A., Ostrouchov S., and Sorensen D.: LAPACK Users’ Guide, 2nd ed., SIAM, Philadelphia,1995. See also LAPACK: http://www.netlib.org/lapack/.  Francis J. G. F.: The QR transformation: a unitary analogue to the LR transformation I, Comput. J. 4 (1961/1962) 265–271.  Kublanovskaya V. N.: On some algorithms for the solution of the complete eigenvalue problem, U.S.S.R. Comput. Math. and Math. Phys. 3 (1961) 637–657.  Fukuda A., Ishiwata E., Iwasaki M., and Nakamura Y.: The discrete hungry Lotka-Volterra system and a new algorithm for computing matrix eigenvalues, Inverse Problems 25 (2009) 015007.  Itoh Y.: Integrals of a Lotka-Volterra system of odd number of variables, Prog. Theor. Phys. 78 (1987) 507–510.  P. Koev P.: Accurate eigenvalues and SVDs of totally nonnegative matrices, SIAM J. Matrix Anal. Appl. 27 (2005) 1–23.  Watkins D. S.: Product eigenvalue problems, SIAM Review 47 (2005) 3–40. Akiko Fukuda Graduate School of Science, Tokyo University of Science, Tokyo 162-8601, Japan E-mail: j1409704(at)ed.kagu.tus.ac.jp