A new restarting method in the harmonic projection algorithm for computing the eigenvalues of a nonsymmetric matrix

A new restarting method in the harmonic projection algorithm for computing the eigenvalues of a nonsymmetric matrix

Available online at www.sciencedirect.com Applied Mathematics and Computation 198 (2008) 143–149 www.elsevier.com/locate/amc A new restarting method...

347KB Sizes 0 Downloads 6 Views

Available online at www.sciencedirect.com

Applied Mathematics and Computation 198 (2008) 143–149 www.elsevier.com/locate/amc

A new restarting method in the harmonic projection algorithm for computing the eigenvalues of a nonsymmetric matrix H. Saberi Najafi *, H. Moosaei Department of Mathematics, Faculty of Sciences, Guilan University, P.O. Box 41335-1914, Rasht, Iran

Abstract The harmonic projection method can be used to find interior eigenpairs of large matrices. Given a target point or shift s to which the needed interior eigenvalues are close, the desired interior eigenpairs are the eigenvalues nearest s and the associated eigenvectors. However, it has been shown that the harmonic projection method may converge erratically and even may fail to do so. In this paper, we present a new restarting method in the harmonic projection algorithm for computing the eigenvalues of a nonsymmetric matrix. The implementation of the algorithm has been tested by numerical examples, the results show that the algorithm converges fast and works with high accuracy.  2007 Elsevier Inc. All rights reserved. Keywords: Eigenvalue; Nonsymmetric matrix; Harmonic; Arnoldi; Restarting

1. Introduction The eigenvalue problem is one of the most important subjects in the applied sciences and engineering. One of the most important and practical topics in computational mathematics is computing some of the interior eigenvalues close to target point or shift s and the associated eigenvectors. The harmonic projection method has been accepted to be one of commonly used method for computing interior eigenpairs. In this paper, we present a new restarting method in the harmonic projection algorithm, for computing some of the interior eigenvalues close to s of the eigenvalue problem, AXi = kiXi, where A is an n · n matrix, and (ki, Xi) is referred to as an eigenpair of A with kXik = 1, here the norm used is the Euclidean norm. The implementation of the algorithm has been tested by numerical examples, the results show that the algorithm converges fast and works with high accuracy.

*

Corresponding author. E-mail addresses: hnajafi@guilan.ac.ir (H. Saberi Najafi), [email protected] (H. Moosaei).

0096-3003/$ - see front matter  2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2007.08.069

144

H. Saberi Najafi, H. Moosaei / Applied Mathematics and Computation 198 (2008) 143–149

Definition 1. Let A be a n · n matrix and v 2 Rn is a vector, then the subspace generated by vectors v, Av, . . . , Am1v is called as Krylov subspace and denoted by km(A, v), i.e. k m ðA; vÞ ¼ spanfv; Av; . . . ; Am1 vg Algorithm 1 (Arnoldi MGS process) Choose a vector v1 of norm 1 For j = 1, . . . , m do w :¼ Avj For i = 1, . . . , j do hij :¼ (w, vi) w :¼ w  hi,jvi End do hj+1,j :¼ kwk2 w vjþ1 :¼ hjþ1;j End do.

Theorem 1.1. The vectors v1, v2, . . . , vm produced by the Arnoldi algorithm form an orthonormal basis of the subspace km = span{v1, Av1, . . . , Am1v1}. Proof. In [2].

h

Theorem 1.2. Denote by Vm a n · m matrix with column vectors v1, v2, . . . , vm and by Hm a m · m Hessenberg matrix whose nonzero entries are defined by the algorithm. Then the following relations hold: AV m ¼ V m H m þ hmþ1;m vmþ1 eHm ; V Hm AV m ffi H m : Proof. In [2].

h

2. The harmonic projection method If s is not an eigenvalue of A then we have from AXi = kiXi that 1 Xi ðA  sIÞ1 X i ¼ ki  s Therefore, the interior eigenvalues near s are transformed into exterior ones with largest magnitudes of (A  sI)1. e i Þ satisfying For the given s and a subspace km(A, v) the harmonic projection method seeks the pairs ð~ki ; X the harmonic projection (for more details see [3,4]), ( e i 2 k m ðA; v1 Þ X e i ? ðA  sIÞk m ðA; v1 Þ ei  ~ AX ki X and uses them to approximate some eigenvalues of A near s and the associated eigenvectors. Algorithm 2 (Harmonic projection algorithm) 1. Input, s, v1 with kv1k = 1, l: the numbers of desired eigenvalues. 2. Run Algorithm 1. (For computing Vm+1 = [v1, v2, . . . , vm+1] and Hm.)

H. Saberi Najafi, H. Moosaei / Applied Mathematics and Computation 198 (2008) 143–149

145

3. If Hm  sI is nonsingular then Solving ððH m  sIÞ þ ðH m  sIÞ

H

em hHmþ1;m hmþ1;m eHm Þgi ¼ ð~ki  sÞgi

else solving H H ððH m  sIÞ ðH m  sIÞ þ em hHmþ1;m hmþ1;m eHm Þgi ¼ ð~ki  sÞðH m  sIÞ gi

end if; (For computing ð~ ki ; gi Þ, i = 1, . . . , m.) 4. Select ~ ki ’s with respect to the smallest value of ð~ki  sÞ’s is to approximate the desired eigenvalues i = 1, . . . , l. e i ¼ V m gi Þ, i = 1, . . . , l as approximations. 5. Take the harmonic Ritz pairs ð~ ki ; X 3. Restarting method in harmonic projection algorithm In the harmonic projection algorithm the initial vector v1 is chosen as an arbitrary vector. Since this vector can increase the accuracy of the result, it would be more suitable if it was chosen in a better way, which is using a process called restarting. This is the main idea in the restarting which can be seen in the following algorithm: Algorithm 3 (Harmonic projection restarting process) Choose initial vector v(0) such that kv(0)k = 1 For S = 0, 1, . . . do v1 = v(s) Run harmonic projection algorithm with initial v1 ðsÞ ðsÞ e ðsÞ Compute the Ritz vectors X i ¼ V m g i i ¼ 1; l . . . ; l Pharmonic ðsÞ l ðsþ1Þ Set v ¼ i¼1 X i and normalized. End do. 3.1. Numerical test 1 Matrix A is a banded matrix select from [6], which use in Quantum Chemistry, i.e.: 2

1 6 0:11 6 6 6 6 0:12 6 6 6 0 6 6 6 0:34 6 A¼6 6 0:45 6 6 6 6 6 6 6 6 6 4

0:21 1:2 2 0:21 0:11

3

0:21

0:12 0:11 0 0:34 0:45

0:12 0 .. ..

0 1:2

4 ..

.

..

.

.

..

.

.

..

.

..

.

0:13 0 1:2 .. . .. . .. . .. . .. . 0:34 0:45

1:42 0:13

3

..

.

1:42 .. . .. . .. .

..

.

0:21

1:2

0

0:11

197

0:21

1:2

0:12

0:11

198

0:21

0 0:34

0:12 0

0:11 0:12

199 0:11

0 ..

.

..

.

..

.

..

..

.

0:13

.

7 7 7 7 7 7 7 7 7 7 1:42 7 7 7 0:13 7 7 7 0 7 7 7 7 1:2 7 7 7 0:21 5 200

200200

We apply Algorithm 3 for m = 11, l = 3, s = 0. The results are shown in Fig. 1. The Fig. 1 shown that Algorithm 3 does not have a good convergence.

146

H. Saberi Najafi, H. Moosaei / Applied Mathematics and Computation 198 (2008) 143–149

Fig. 1. shows the behavior of the average values of ri for three eigenvalues. Horizontal axis: the number restarts. Vertical axis: average values of ri at each restart.

4. New restarting method for harmonic projection method All restarting methods described in [1,5,7,8] can also be used for computing the harmonic Ritz pair, in this section a new restart method for computing the harmonic Ritz pair is described. This method has a very simple algorithm and can be used on personal computers. First, we prove that following theorem and then H-R algorithm will be described. ðmÞ e ðmÞ Theorem 4.1. Let X ¼ V m gi be the harmonic Ritz vector, ~ki be the harmonic Ritz value then: i " # ðmÞ ðmÞ ~ ðH  k IÞg m ðmÞ ðmÞ i i ðA  ~ ki IÞX i ¼ V mþ1 ðmÞ hmþ1;m eHm gi

Proof. Since we have AV m ¼ V m H m þ hmþ1;m vmþ1 eHm then ðmÞ

AX i

ðmÞ

¼ AV m gi

ðmÞ

¼ V m H m gi

ðmÞ

þ hmþ1;m vmþ1 eHm gi

So ðmÞ

ðmÞ

ðA  ~ ki IÞX i

ðmÞ

ðmÞ ðmÞ ðmÞ þ hmþ1;m vmþ1 eHm gi  k~i V m gi " # ðmÞ ðmÞ ðH m  ~ki IÞgi ðmÞ ðmÞ ðmÞ H ¼ V m ðH m  ~ ki IÞgi þ hmþ1;m vmþ1 em gi ¼ V mþ1 : ðmÞ hmþ1;m eHm gi

¼ AV m gi

ðmÞ

ðmÞ

~ ki V m g i

ðmÞ

¼ V m H m gi

Algorithm 4 (New restarting for Harmonic projection process (H-R)). Choose initial vector v(0) such that kv(0)k = 1 For k = 1,2, . . . , m do For s = 1,2,. . . do



H. Saberi Najafi, H. Moosaei / Applied Mathematics and Computation 198 (2008) 143–149

147

Run harmonic projection algorithm with initial v1 ðsÞ Compute harmonic Ritz value ~ ki and harmonic Ritz vector e ðsÞ X i i ¼ 1; . . . ; l ðsÞ

ðsÞ

ðsÞ

m Set c ¼ minfri ¼ kAX i  kðsÞ i X i kgi¼k for P = k, k + 1, . . . , m do ðsÞ e ðsÞ If c ¼ X P then, set v1 ¼ X P and normalized. end if; end do; end do; ðsÞ ðsÞ XP $ Xk ðsÞ U k ¼P Xk P ðsÞ v1 ¼ kj¼1 U j þ lj¼kþ1 X j and normalized End do

4.1. Numerical test 2 Let A be a 200 · 200 matrix used in numerical test 1. We obtain three eigenvalue close to s of matrix A. Step 1. Let m = 11, l = 3, s = 0 in the H-R algorithm, then we obtain three eigenvalues namely k1, k2, k3. These eigenvalues are approximations of three eigenvalues of A, and residual values are r1, r2, r3: ð1Þ

k1 ¼ 6:035;

ð1Þ

r1 ¼ 9213:495;

ð1Þ

ð1Þ

ð1Þ

ð1Þ

k2 ¼ 18:072; r2 ¼ 26493:277; k3 ¼ 35:302; r3 ¼ 52175:339: ð1Þ

ð1Þ

Since the residual value of k1 is minimum, therefore, we select vector V 1 ¼ X 1 for the next iteration of harmonic projection algorithm.After 26 iterations the values are ð26Þ

k1

ð26Þ k2 ð26Þ k3

ð26Þ

¼ 0:842; r1 ¼ 3:517; ¼ 7:686;

ð26Þ r2 ð26Þ r3

¼ 0:00000274; ¼ 2:948; ¼ 6:703:

According to the results shown, the value of k1 has improved very well (see Fig. 2). ð26Þ ð26Þ ð26Þ Step 2. To improve the other eigenvalues, set v1 = U1 + X2 + X3 where U 1 ¼ X 1 ; X 2 ¼ X 2 ; X 3 ¼ X 3 after 75 iterations, the eigenvalue k2 improved very well (see Fig. 3.) ð101Þ

k2

ð101Þ

¼ 5:946; r2

¼ 0:00527:

Step 3. As the results show, the eigenvalue k3 improved, so for the next iteration set v1 = U1 + U2 + X2 where ð26Þ ð101Þ ð101Þ U 1 ¼ X 1 ; U 2 ¼ X 2 ; X 2 ¼ X 3 and doing H-R algorithm again.

Fig. 2. Showing the improvement of the first eigenvalue (k1) in the first step of H-R restarting after a number of restarts.

148

H. Saberi Najafi, H. Moosaei / Applied Mathematics and Computation 198 (2008) 143–149

Fig. 3. Showing the improvement of the second eigenvalue (k2) in the second step of H-R restarting after a number of restarts.

Fig. 4. Showing the improvement of the third eigenvalue (k3) in the third step of H-R restarting after a number of restarts.

Table 1 Show the outputs of H-M algorithm for m = 6, l = 2, s = 0 Step 1 k1

ð76Þ

¼ 0:842

r1

ð76Þ k2

ð76Þ

¼ 6:735

r2

¼ 0:0000954

ð76Þ

¼ 3:177

Step 2 ð213Þ

ð213Þ

¼ 0:0000215

ð213Þ

¼ 0:0000257

k1

¼ 8:259

r1

ð213Þ k2

¼ 20:161

r2

Method

Iteration

Error

Algorithm 1 H-R algorithm

600 213

No convergence 0.000074

After 76 iterations, the eigenvalue k3 improved very well (see Fig. 4.) ð177Þ

k3

ð177Þ

¼ 24:003 r3

¼ 0:00000707

As the results show the new restarting of harmonic projection algorithm works better than Algorithm 3 and it gives the results with high accuracy. 4.2. Numerical test 3 Let A be a 200 · 200 matrix used in numerical test 1. We obtain three eigenvalue close to s of matrix A. See Table 1 for results. As the results show the new restarting of harmonic projection (Algorithm 4) works better than Algorithm 3 and it gives the results with high accuracy. 5. Comments and conclusions 1. The first method (Algorithm 3) has a very easy and simple algorithm but does not have a good convergence, as was shown before; this method gives a little improvement to the approximated eigenvalues.

H. Saberi Najafi, H. Moosaei / Applied Mathematics and Computation 198 (2008) 143–149

149

2. The H-R algorithm (Algorithm 4) has the following properties: (a) This algorithm is very simple to run. (b) It can be run on any PC. (c) The eigenvalue problem can be computed with a high accuracy and less consuming time. (d) The approximation obtained by this method can be improved more easily compared to the other method. References [1] Y. Sadd, Numerical method for large eigenvalue problem, in: Algorithm and Architectures for Advanced Scientific Computing, Manchester University Press, 1992. [2] Y. Sadd, Iterative Method for Sparse Linear System, PWS Publishing Company, A Division of International Thomson Publishing Inc., USA, 1996. [3] G. Wu, A modified harmonic block Arnoldi Algorithm with adaptive shifts for large interior eigenproblems, J. Comput. Appl. Math. 205 (1) (2007) 343–365. [4] R.B. Morgan, M. Zeng, Harmonic projection method for large non-symmetric eigenvalue problem, Numer. Linear Algebra Appl. 5 (1998) 33–55. [5] H. Saberi Najafi, E. Khaleghi, A new restarting method in the Arnoldi algorithm for computing the eigenvalues of a nonsymmetric matrix, Appl. Math. 156 (2004) 59–71. [6] S.H. Lam, The role of CSP in reduce chemistry modeling, in: IMA Symposium, Minnesota University, 13–19 October 1999. [7] H. Saberi Najafi, H. Ghazvini, Weighted restarting method in the weighted Arnoldi algorithm for computing the eigenvalues of a nonsymmetric matrix, Appl. Math. 175 (2006) 1279–1287. [8] H. Saberi Najafi, A. Refahi, A new restarting method in Lanczos algorithm for generalized eigenvalue problem, Appl. Math. Comput. 184 (2007) 421–428.