Rational Krylov Algorithms for Eigenvalue Computation ... - CiteSeerX
Rational Krylov Algorithms for Eigenvalue Computation and Model Reduction Axel Ruhe1 ? and Daniel Skoogh1 Department of Mathematics, Chalmers Institute of Technology and the University of Goteborg, S-41296 Goteborg, Sweden fruhe,[email protected].
Abstract. Rational Krylov is an extension of the Lanczos or Arnoldi
eigenvalue algorithm where several shifts (matrix factorizations) are performed in one run. A variant has been developed, where these factorizations are performed in parallel. It is shown how Rational Krylov can be used to nd a reduced order model of a large linear dynamical system. In Electrical Engineering, it is important that the reduced model is accurate over a wide range of frequencies, and then Rational Krylov with several shifts comes to advantage. Results for numerical examples coming from Electrical Engineering applications are demonstrated.
1 Introduction, Eigenvalues and Linear Dynamic Systems In many scienti c and engineering contexts one needs to solve an algebraic eigenvalue problem, (A ? B )x = 0 ;
(1)
for eigenvalues k and eigenvectors xk . Its solution is used as the rst step on the way from a static to a dynamic description of a system of interacting entities. Think of a network of electric components or a mechanical construction of masses and springs! Small perturbations of a static, equilibrium, position will be formed as eigenvectors and move as the eigenvalues indicate. The foremost mathematical modeling tool is the Laplace transform, where we map a description in time and space to one involving decay rates and oscillation frequences, these are the real and imaginary parts of the eigenvalues. In this realm we use eigenvalues to get the solution of a linear system of ordinary di erential equations,
x_ = Ax ; x(0) = a ; ?
Talk given at PARA 98, Umea Sweden, June 17, 1998. Partial support given by TFR, the Swedish Research Council for Engineering Sciences, Dnr 222, 96-555, and from the Royal Society of Arts and Sciences in Goteborg.
2
Axel Ruhe and Daniel Skoogh
over time as a linear combination of eigensolutions,
x(t) = k k xk ek t : Today we will concentrate on how to nd reduced order models of large linear dynamic systems. Let us specialize, and consider a model of an electric RLC network, that is a network containing only passive components, resistances (R) inductances (L) and capacitances (C). Such a network follows the equations,
C x_ = ?Gx + ru ; y = lT x
(2)
where x is the state vector, u is the vector of inputs (controls) and y is the vector of outputs. G is the matrix of DC components (resistances) and C is the matrix of components with memory (inductances and capacitances). The external connections for input and output are given by r and l, which for the single input single output case are just column vectors. For an RLC network, the state vector and matrices are partitioned into,
x = vi C = C0 L0 G = ?GB T B0 ; (3) where v is the vector of node voltages, i is the vector of currents in those branches that have inductors, C is the matrix of capacitances, L of inductances and G of conductances (reciprocal of resistance). The rectangular matrix B is just the
incidence matrix of the inductor branches. In an RLC network, the matrices C , L and G are all symmetric positive semide nite. In realistic cases we may also safely assume that the pencil (G; C ) is regular. Make a Laplace transform of the system (2) and get,
sC X = ?GX + rU ; Y = lT X
(4)
giving the input output relation,
Y = lT (G + sC )?1 rU ; = H (s)U which de nes the transfer function H (s).
(5)
We see that the transfer function has poles in the complex plane at the eigenvalues sk of the pencil (G + sC )x = 0 :
(6)
Purely imaginary eigenvalues correspond to persistent oscillations, and we are often interested in the norm of H (s) along the imaginary axis. An eigenvalue close to the imaginary axis shows up as a peak of this norm, and an excitation
Rational Krylov algorithms
3
at this frequency will give rise to a slowly decaying oscillation. This is, of course, provided the system is stable, which for our choice of sign in the pencil (6) occurs when all eigenvalues are in the right half plane. The task is now to nd a new reduced order model, with a state space of a much smaller dimension, chosen so that the transfer function H~ (s) of the reduced model keeps all the important properties of the original transfer function (5). Let us assume that we have found a subspace spanned by the rst k columns of the nonsingular matrix S , and express an arbitrary vector x in the state space by,
x = S x~ : The linear system (2) is now transformed into, S ?1 C S x~_ = ?S ?1GS x~ + S ?1 ru ; y= lT S x~
(7)
and we get the system of the reduced order model by keeping the appropriate leading k k submatrices of these. Many di erent ways to choose a basis S have been tried. An evident choice is to take the eigenvectors that correspond to eigenvalues in the frequency range of interest, this is standard practice in Mechanical Engineering. Here one is most often interested in the lowest frequencies, or in frequencies close to a slow external load, like the hum of an engine or the shaking of an earthquake. Recently Krylov spaces, computed by Lanczos or Arnoldi algorithms applied to a shifted and inverted matrix pencil (8) have caught on, especially in Electrical Engineering, (Freund and Feldmann at Bell Laboratories [4, 5], Elfadel and Ling et al at IBM [2]) but also in Mechanical Engineering (E. Grimme et al [6]) and even in Quantum Chemistry, the tridiagonal matrix from Lanczos has been used to predict photoionization cross sections ( H. Karlsson [7]). A Rational Krylov approach in the spirit we are going to discuss today has already been tried at several places, see [6] and [1]!
2 Rational Krylov A standard practice to nd eigenvalues of large matrix pencils (1) is to apply a Krylov space algorithm, Lanczos or Arnoldi, to a shift invert spectral transformation,
C = (A ? B )?1 B :
(8)
It will compute eigenvalues close to the shift point in the complex plane[3]. The Rational Krylov algorithm [9] is a further development of this approach, where an orthonormal basis Vj is built up one column at a time. In each step a shifted and inverted matrix operator (8) is applied, and the resulting vector is orthogonalized to the already computed part of the basis by means of Gram
4
Axel Ruhe and Daniel Skoogh
Schmidt. Rational Krylov di ers from shifted and inverted Arnoldi, in that it may use di erent shifts j in di erent steps j . Let us rst give the algorithm step by step, then derive the basic recursion between the computed vectors, and last show how it is equivalent to a certain shifted and inverted Arnoldi iteration. Algorithm RKS
1. Start with vector v1 , where kv1 k = 1. Set t1 = e1 . 2. For j = 1; 2; : : : until convergence (a) Choose shift j . (b) Set r = Vj tj , (continuation vector). (c) Operate r := (A ? j B )?1 Br (d) Orthogonalize r := r ? Vj hj where hj = VjH r, (Gram Schmidt). (e) Get new vector vj+1 = r=hj+1;j , where hj+1;j = krk, (normalize). (f) Compute approximate solutions, (ij) and yi(j) , and test for convergence. (g) Get tj+1 , (continuation combination). 3. Compute eigenvectors . End.
Let us follow what happens during a typical step j of Algorithm RKS. Eliminate the intermediate vector r, used in steps 2b { 2e, and get
Vj+1 hj = (A ? j B )?1 BVj tj ; the vector hj now has length j + 1. Multiply from the left by (A ? j B ), (A ? j B )Vj+1 hj = BVj tj : (9) Separate terms with A to the left and B to the right, AVj+1 hj = BVj+1 (hj j + tj ) ; with a zero added to the bottom of the vector tj to give it length j + 1. This is the relation for the j th step, now put the corresponding vectors from
the previous steps in front of this and get,
AVj+1 Hj+1;j = BVj+1 Kj+1;j ; (10) with two (j + 1) j Hessenberg matrices, Hj+1;j = [h1 h2 : : : hj ], containing the Gram Schmidt orthogonalization coecients, and
Kj+1;j = Hj+1;j diag(i ) + Tj+1;j (11) with the upper triangular matrix Tj+1;j = [t1 t2 : : : tj ], built up from the continuation combinations used in step 2b.
Rational Krylov algorithms
5
If we choose all zero shifts, j = 0, and continue with the latest vector, tj = ej , in all steps j , we get back the shifted and inverted Arnoldi, A?1 BVj = Vj+1 Hj+1;j : In the general case, it is also possible to bring the recursion (10) into an Arnoldi recursion by careful application of orthogonal transformations to the Hessenberg matrices H and K , bringing one of them to triangular form. We nd orthogonal matrices Qj+1;j+1 and Zj;j so that the transformed matrices, K~ = QHj+1;j+1 Kj+1;j Zj;j is Hessenberg, R = QHj+1;j+1 Hj+1;j Zj;j is upper triangular, Wj+1 = Vj+1 Qj+1;j+1 is a new orthogonal basis. This transformation is done one column at a time, it is easiest to understand if we follow one minor step in detail. Assume that j = 4 and that the rst 2 columns are already in the desired form,
3 3 2 2 xxxx xxxx 66x x x x77 66 x x x77 7 6 H = 6 x x7 ; K = 66 x x x77 : 4 x x5 4 y x5 x x Apply a rotation R3;4 in the (3; 4) plane from the left to zero out the element y at h4;3 . However this rotation lls k4;2 , which has to be chased up left along the main
diagonal. Let us write this in a table, where the arrows denote multiplication from left or right, R3;4 ! zeros h4;3 lls k4;2 S2;3 zeros k4;2 lls h3;2 R2;3 ! zeros h3;2 lls k3;1 S1;2 zeros k3;1 lls h2;1 R1;2 ! zeros h2;1 : After this sweep both matrices have the same nonzero pattern as before, except that the element y at h4;3 is put to zero. After we had made one such sweep for each column up to column j , we have the desired form, (12) AWj Rj;j = BWj+1 K~ j+1;j ; and can get an approximate eigensolution from approximate eigenvalue i(j) ; K~ j;j yi(j) = Rj;j yi(j) i(j) approximate eigenvector x~i : x~i = Wj Rj;j yi(j) Note that this pencil has the same nonzero pattern as the intermediate matrix in the QZ algorithm. Note also that this could have been done on a shifted pencil,
6
Axel Ruhe and Daniel Skoogh
and yield a more accurate recursion than we described in [8]. For an appropriate choice of shift, wj+1 will be a good continuation vector as described in [9]. Note, however also that even if (12) is an Arnoldi recursion, it does not start at the original starting vector v1 but at w1 , which is a linear combination of all the vectors computed during the whole RKS algorithm.
3 Reduced Order Model Now return to the linear dynamic system (2) and its response function H (s) de ned by (5). Let us specialize to the single input single output case where both matrices l and r are single vectors and we get,
H (s) = lT (G + sC )?1 r :
(13)
Our task is now to evaluate this for many values of the frequency parameter s. We see that for each s we need to solve the system (G + sC )x(s) = r
(14)
for the vector x(s) and then do a scalar multiplication with l. This can be done with few factorizations of shifted large matrices (14) in the following way. Apply the Rational Krylov algorithm to the matrix pencil (6) i. e. set A = G, B = C and shift with j = ?sj , some appropriate frequencies for which the matrices G + si C are nonsingular, and chosen so that the frequency range of interest is covered. Start at the vector v1 0 = (G + s1 C )?1 r, and note that in this case we do not use a random start. After j steps we may compute an approximate solution to x(s) (14) as,
x~(s) = Vj+1 (Kj+1;j + s1 Hj+1;j ) y~(s) ; where,
y~(s) = (Kj;j + sHj;j )?1 e1 0 ;
(15)
a solution to a j j linear system with a coecient matrix of Hessenberg form. To see how this comes about, shift the basic Rational Krylov recursion (10) and get, (G + sC )V (K + si H ) = (G + si C )V (K + sH ) ; and multiply with a shifted and inverted matrix to get (G + si C )?1 (G + sC )V (K + si H ) = V (K + sH ) :
(16)
This is useful to get an expression for the premultiplied residual,
r(s) = (G + si C )?1 f(G + sC ) x~(s) ? rg :
(17)
Rational Krylov algorithms
7
The choice to start at v1 0 = (G + s1 C )?1 r and the relation (16) will now for
i = 1 imply that,
r(s) = (G + s1 )?1 (G + sC )Vj+1 (Kj+1;j + s1 Hj+1;j ) y~(s) ? v1 0 = Vj+1 f(Kj+1;j + sHj+1;j ) y~(s) ? e1 0 g = vj+1 (kj+1;j + shj+1;j )( y~(s))j = vj+1 hj+1;j (s ? sj )( y~(s))j ;
(18)
a unit length vector vj+1 multiplied by the last Gram Schmidt orthogonalization coecient hj+1;j , the di erence s ? sj between the evaluation frequency s and the last shift sj , and nally the last element ( y~(s))j of the solution to the j j system (15). This last factor is likely to get small when we have used shifts si reasonably close to s during the Rational Krylov process. Compare to the similar expressions one gets when applying Lanczos or Arnoldi to an eigenvalue problem or GMRES to a linear system of equations. Our way (15) to get the approximate solution actually corresponds to the FOM (Full Orthogonalization Method), which is closely related to GMRES. One advantage of using FOM, instead of the more accurate GMRES, is that we can represent the reduced order model as,
H~ (s) = ~lT (G~ + sC~ )?1 r~ ;
(19)
with ~lT = lT Vj+1 (Kj+1;j + s1 Hj+1;j ) G~ = Kj;j C~ = Hj;j r~ = e1 0 : The error of the approximate solution x~(s) can be expressed as
x~(s) = x(s) ? x~(s) = (G + sC )?1 (G + si C )r(s) = (I + (si ? s)(G + sC )?1 C )r(s);
(20)
where r(s) is the premultiplied residual (17). An upper bound for the error of the approximate transfer function H~ (s) is
j(s)j k l k j (s)j(1 + jsi ? sj k (G + sC )? C vj k): 1
+1
(21)
where (s) = hj+1;j (s ? sj )( y~(s))j , see (18). In order to compute an ecient and reliable error estimate we need to estimate k (G + sC )?1 C vj+1 k. How this is done is described in [10].
8
Axel Ruhe and Daniel Skoogh
4 Two Numerical Examples Let us report tests on two models of medium size radio frequency circuits. The rst test example is provided by Zhaojun Bai and is of dimension n = 256. The rational Krylov algorithm pis applied with the shifts 1 = ?j 108, 2 = ?j 109 and 3 = ?j 1010, where j = ?1. The shifts are applied 5, 7 and 7 times respectively. In gure 2 we show the approximate transfer function H~ (s), and in gure 3 we show the error and an error estimate of H~ (s). Note that the true error is below the error estimate. The eigenvalue distribution is given in gure 1. Compare the eigenvalue distribution with the plots for the transfer function and error estimate, and note that the transfer function has a top near the imaginary part of each eigenvalue. We did another run with the same matrices but the shifts closer together at 1 = ?j 109, 2 = ?2j 109 and 3 = ?7j 109, used 4, 4 and 8 times respectively. Look at gures 4 and 5 and note that now the errors are much smaller at the higher frequencies, where there are many poles. It is evident that an adaptive strategy where the shifts are chosen where we want small errors will be of advantage. The second test example is provided by Peter Feldmann and is of dimension n = 1841. Now the rational Krylov algorithm is applied with the shifts 1 = ?j 108, 2 = ?j 109 and 3 = ?4j 109, applied 4, 8 and 13 times respectively. In gure 6 we show the approximate transfer function H~ (s), and in gure 7 the error and an error estimate of H~ (s). The eigenvalue distribution is given in gure 8. 4
0
x 10
−0.5 −1 −1.5
real
−2 −2.5 −3 −3.5 −4 −4.5 −5 8 10
9
10 imaginary
10
10
Fig. 1. The eigenvalue distribution of the rst test example given by Bai, n = 256.
Rational Krylov algorithms
9
3
10
2
|H(j 2 π f)|
10
1
10
0
10
−1
10
8
10
9
10 2π f
10
10
Fig. 2. The rst test example of dimension n = 256. The gure shows the approximate transfer function H~ (s). The rational Krylov algorithm is applied with the shifts 1 =
?j 108 , 2 = ?j 109 and 3 = ?j 1010 , where j = p?1. The shifts are applied 5, 7 and
7 times respectively. 5
10
0
error
10
−5
10
−10
10
−15
10
8
10
9
10 2π f
10
10
Fig. 3. The error and error estimate of the rst test example. The error estimate is plotted with a solid line ?, and the true error is plotted with a dashed line ? ?
10
Axel Ruhe and Daniel Skoogh 3
10
2
|H(j 2 π f)|
10
1
10
0
10
−1
10
8
10
9
10 2π f
10
10
Fig. 4. First test example is of dimension n = 256. Approximate transfer function H~ (s9)
for rational Krylov algorithm is applied with the shifts closer together at 1 = ?j 10 , 2 = ?2j 109 and 3 = ?7j 109 , applied 4, 4 and 8 times respectively. 5
10
0
error
10
−5
10
−10
10
−15
10
8
10
9
10 2π f
10
10
Fig. 5. The error and error estimate of the rst test example with shifts closer together,
corresponding to the gure 4. The error estimate is plotted with a solid line ?, and the true error is plotted with a dashed line ? ?
Rational Krylov algorithms
11
0
10
−1
|H(j f)|
10
−2
10
−3
10
−4
10
8
10
9
10 f
10
10
Fig. 6. Second test example, dimension n = 1841. Approximate transfer function H~ (s) 8 9 for the rational Krylov algorithm applied with shifts = ? j 10 , 1 2 = ?j 10 and p 9 3 = ?4j 10 , where j = ?1. The shifts are applied 4, 8 and 13 times respectively. 5
10
0
error
10
−5
10
−10
10
−15
10
8
10
9
10 f
10
10
Fig. 7. The error and error estimate of the second test example. The error estimate is plotted with a solid line ?, and the true error is plotted with a dashed line ? ?
12
Axel Ruhe and Daniel Skoogh 9
0
x 10
real
−2 −4 −6 −8 8 10
9
10 imaginary/(2*π)
10
10
Fig. 8. The eigenvalue distribution of the second test example given by Feldmann, n = 1841.
References 1. I. M. Elfadel and D. D. Ling, A block rational Arnoldi algorithm for multipoint passive model-order reduction of multiport RLC networks, in Technical Digest of the 1997 IEEE/ACM International Conference on Computer-Aided Design, IEEE Computer Society Press, 1997, pp. 66{71. 2. , Zeros and passivity of Arnoldi-reduced-order models for interconnect networks, in Proc 34nd IEEE/ACM Design Automation Conference, ACM, New York, 1997, pp. 28{33. 3. T. Ericsson and A. Ruhe, The spectral transformation Lanczos method for the numerical solution of large sparse generalized symmetric eigenvalue problems, Mathematics of Computation, 35 (1980), pp. 1251{1268. 4. R. W. Freund, Circuit simulation techniques based on Lanczos-type algorithms, in Systems and Control in the Twenty-First Century, C. I. Byrnes, B. N. Datta, D. S. Gilliam, and C. F. Martin, eds., Birkauser, 1997, pp. 171{184. 5. , Reduced-order modeling techniques based on Krylov subspaces and their use in circuit simulation, Applied and Computational Control, Signals, and Circuits, (1998, to appear). 6. K. Gallivan, E. Grimme, and P. Van Dooren, A rational Lanczos algorithm for model reduction, Numerical Algorithms, 12 (1996), pp. 33{64. 7. H. O. Karlsson, Atomic and Molecular Density-of-States by Direct Lanczos Methods, PhD thesis, Uppsala University, Department of Quantum Chemistry, 1994. 8. A. Ruhe, Eigenvalue algorithms with several factorizations { a uni ed theory yet? to appear, 1998. 9. , Rational Krylov, a practical algorithm for large sparse nonsymmetric matrix pencils, SIAM J. Sci. Comp., 19 (1998), pp. 1535{1551. 10. D. Skoogh, Model order reduction by the rational Krylov method. work in progress, 1998.
Abstract. Rational Krylov is an extension of the Lanczos or Arnoldi
eigenvalue algorithm where several shifts (matrix factorizations) are performed in one run. A variant has been developed, where these factorizations are performed in parallel. It is shown how Rational Krylov can be used to nd a reduced order model of a large linear dynamical system. In Electrical Engineering, it is important that the reduced model is accurate over a wide range of frequencies, and then Rational Krylov with several shifts comes to advantage. Results for numerical examples coming from Electrical Engineering applications are demonstrated.
1 Introduction, Eigenvalues and Linear Dynamic Systems In many scienti c and engineering contexts one needs to solve an algebraic eigenvalue problem, (A ? B )x = 0 ;
(1)
for eigenvalues k and eigenvectors xk . Its solution is used as the rst step on the way from a static to a dynamic description of a system of interacting entities. Think of a network of electric components or a mechanical construction of masses and springs! Small perturbations of a static, equilibrium, position will be formed as eigenvectors and move as the eigenvalues indicate. The foremost mathematical modeling tool is the Laplace transform, where we map a description in time and space to one involving decay rates and oscillation frequences, these are the real and imaginary parts of the eigenvalues. In this realm we use eigenvalues to get the solution of a linear system of ordinary di erential equations,
x_ = Ax ; x(0) = a ; ?
Talk given at PARA 98, Umea Sweden, June 17, 1998. Partial support given by TFR, the Swedish Research Council for Engineering Sciences, Dnr 222, 96-555, and from the Royal Society of Arts and Sciences in Goteborg.
2
Axel Ruhe and Daniel Skoogh
over time as a linear combination of eigensolutions,
x(t) = k k xk ek t : Today we will concentrate on how to nd reduced order models of large linear dynamic systems. Let us specialize, and consider a model of an electric RLC network, that is a network containing only passive components, resistances (R) inductances (L) and capacitances (C). Such a network follows the equations,
C x_ = ?Gx + ru ; y = lT x
(2)
where x is the state vector, u is the vector of inputs (controls) and y is the vector of outputs. G is the matrix of DC components (resistances) and C is the matrix of components with memory (inductances and capacitances). The external connections for input and output are given by r and l, which for the single input single output case are just column vectors. For an RLC network, the state vector and matrices are partitioned into,
x = vi C = C0 L0 G = ?GB T B0 ; (3) where v is the vector of node voltages, i is the vector of currents in those branches that have inductors, C is the matrix of capacitances, L of inductances and G of conductances (reciprocal of resistance). The rectangular matrix B is just the
incidence matrix of the inductor branches. In an RLC network, the matrices C , L and G are all symmetric positive semide nite. In realistic cases we may also safely assume that the pencil (G; C ) is regular. Make a Laplace transform of the system (2) and get,
sC X = ?GX + rU ; Y = lT X
(4)
giving the input output relation,
Y = lT (G + sC )?1 rU ; = H (s)U which de nes the transfer function H (s).
(5)
We see that the transfer function has poles in the complex plane at the eigenvalues sk of the pencil (G + sC )x = 0 :
(6)
Purely imaginary eigenvalues correspond to persistent oscillations, and we are often interested in the norm of H (s) along the imaginary axis. An eigenvalue close to the imaginary axis shows up as a peak of this norm, and an excitation
Rational Krylov algorithms
3
at this frequency will give rise to a slowly decaying oscillation. This is, of course, provided the system is stable, which for our choice of sign in the pencil (6) occurs when all eigenvalues are in the right half plane. The task is now to nd a new reduced order model, with a state space of a much smaller dimension, chosen so that the transfer function H~ (s) of the reduced model keeps all the important properties of the original transfer function (5). Let us assume that we have found a subspace spanned by the rst k columns of the nonsingular matrix S , and express an arbitrary vector x in the state space by,
x = S x~ : The linear system (2) is now transformed into, S ?1 C S x~_ = ?S ?1GS x~ + S ?1 ru ; y= lT S x~
(7)
and we get the system of the reduced order model by keeping the appropriate leading k k submatrices of these. Many di erent ways to choose a basis S have been tried. An evident choice is to take the eigenvectors that correspond to eigenvalues in the frequency range of interest, this is standard practice in Mechanical Engineering. Here one is most often interested in the lowest frequencies, or in frequencies close to a slow external load, like the hum of an engine or the shaking of an earthquake. Recently Krylov spaces, computed by Lanczos or Arnoldi algorithms applied to a shifted and inverted matrix pencil (8) have caught on, especially in Electrical Engineering, (Freund and Feldmann at Bell Laboratories [4, 5], Elfadel and Ling et al at IBM [2]) but also in Mechanical Engineering (E. Grimme et al [6]) and even in Quantum Chemistry, the tridiagonal matrix from Lanczos has been used to predict photoionization cross sections ( H. Karlsson [7]). A Rational Krylov approach in the spirit we are going to discuss today has already been tried at several places, see [6] and [1]!
2 Rational Krylov A standard practice to nd eigenvalues of large matrix pencils (1) is to apply a Krylov space algorithm, Lanczos or Arnoldi, to a shift invert spectral transformation,
C = (A ? B )?1 B :
(8)
It will compute eigenvalues close to the shift point in the complex plane[3]. The Rational Krylov algorithm [9] is a further development of this approach, where an orthonormal basis Vj is built up one column at a time. In each step a shifted and inverted matrix operator (8) is applied, and the resulting vector is orthogonalized to the already computed part of the basis by means of Gram
4
Axel Ruhe and Daniel Skoogh
Schmidt. Rational Krylov di ers from shifted and inverted Arnoldi, in that it may use di erent shifts j in di erent steps j . Let us rst give the algorithm step by step, then derive the basic recursion between the computed vectors, and last show how it is equivalent to a certain shifted and inverted Arnoldi iteration. Algorithm RKS
1. Start with vector v1 , where kv1 k = 1. Set t1 = e1 . 2. For j = 1; 2; : : : until convergence (a) Choose shift j . (b) Set r = Vj tj , (continuation vector). (c) Operate r := (A ? j B )?1 Br (d) Orthogonalize r := r ? Vj hj where hj = VjH r, (Gram Schmidt). (e) Get new vector vj+1 = r=hj+1;j , where hj+1;j = krk, (normalize). (f) Compute approximate solutions, (ij) and yi(j) , and test for convergence. (g) Get tj+1 , (continuation combination). 3. Compute eigenvectors . End.
Let us follow what happens during a typical step j of Algorithm RKS. Eliminate the intermediate vector r, used in steps 2b { 2e, and get
Vj+1 hj = (A ? j B )?1 BVj tj ; the vector hj now has length j + 1. Multiply from the left by (A ? j B ), (A ? j B )Vj+1 hj = BVj tj : (9) Separate terms with A to the left and B to the right, AVj+1 hj = BVj+1 (hj j + tj ) ; with a zero added to the bottom of the vector tj to give it length j + 1. This is the relation for the j th step, now put the corresponding vectors from
the previous steps in front of this and get,
AVj+1 Hj+1;j = BVj+1 Kj+1;j ; (10) with two (j + 1) j Hessenberg matrices, Hj+1;j = [h1 h2 : : : hj ], containing the Gram Schmidt orthogonalization coecients, and
Kj+1;j = Hj+1;j diag(i ) + Tj+1;j (11) with the upper triangular matrix Tj+1;j = [t1 t2 : : : tj ], built up from the continuation combinations used in step 2b.
Rational Krylov algorithms
5
If we choose all zero shifts, j = 0, and continue with the latest vector, tj = ej , in all steps j , we get back the shifted and inverted Arnoldi, A?1 BVj = Vj+1 Hj+1;j : In the general case, it is also possible to bring the recursion (10) into an Arnoldi recursion by careful application of orthogonal transformations to the Hessenberg matrices H and K , bringing one of them to triangular form. We nd orthogonal matrices Qj+1;j+1 and Zj;j so that the transformed matrices, K~ = QHj+1;j+1 Kj+1;j Zj;j is Hessenberg, R = QHj+1;j+1 Hj+1;j Zj;j is upper triangular, Wj+1 = Vj+1 Qj+1;j+1 is a new orthogonal basis. This transformation is done one column at a time, it is easiest to understand if we follow one minor step in detail. Assume that j = 4 and that the rst 2 columns are already in the desired form,
3 3 2 2 xxxx xxxx 66x x x x77 66 x x x77 7 6 H = 6 x x7 ; K = 66 x x x77 : 4 x x5 4 y x5 x x Apply a rotation R3;4 in the (3; 4) plane from the left to zero out the element y at h4;3 . However this rotation lls k4;2 , which has to be chased up left along the main
diagonal. Let us write this in a table, where the arrows denote multiplication from left or right, R3;4 ! zeros h4;3 lls k4;2 S2;3 zeros k4;2 lls h3;2 R2;3 ! zeros h3;2 lls k3;1 S1;2 zeros k3;1 lls h2;1 R1;2 ! zeros h2;1 : After this sweep both matrices have the same nonzero pattern as before, except that the element y at h4;3 is put to zero. After we had made one such sweep for each column up to column j , we have the desired form, (12) AWj Rj;j = BWj+1 K~ j+1;j ; and can get an approximate eigensolution from approximate eigenvalue i(j) ; K~ j;j yi(j) = Rj;j yi(j) i(j) approximate eigenvector x~i : x~i = Wj Rj;j yi(j) Note that this pencil has the same nonzero pattern as the intermediate matrix in the QZ algorithm. Note also that this could have been done on a shifted pencil,
6
Axel Ruhe and Daniel Skoogh
and yield a more accurate recursion than we described in [8]. For an appropriate choice of shift, wj+1 will be a good continuation vector as described in [9]. Note, however also that even if (12) is an Arnoldi recursion, it does not start at the original starting vector v1 but at w1 , which is a linear combination of all the vectors computed during the whole RKS algorithm.
3 Reduced Order Model Now return to the linear dynamic system (2) and its response function H (s) de ned by (5). Let us specialize to the single input single output case where both matrices l and r are single vectors and we get,
H (s) = lT (G + sC )?1 r :
(13)
Our task is now to evaluate this for many values of the frequency parameter s. We see that for each s we need to solve the system (G + sC )x(s) = r
(14)
for the vector x(s) and then do a scalar multiplication with l. This can be done with few factorizations of shifted large matrices (14) in the following way. Apply the Rational Krylov algorithm to the matrix pencil (6) i. e. set A = G, B = C and shift with j = ?sj , some appropriate frequencies for which the matrices G + si C are nonsingular, and chosen so that the frequency range of interest is covered. Start at the vector v1 0 = (G + s1 C )?1 r, and note that in this case we do not use a random start. After j steps we may compute an approximate solution to x(s) (14) as,
x~(s) = Vj+1 (Kj+1;j + s1 Hj+1;j ) y~(s) ; where,
y~(s) = (Kj;j + sHj;j )?1 e1 0 ;
(15)
a solution to a j j linear system with a coecient matrix of Hessenberg form. To see how this comes about, shift the basic Rational Krylov recursion (10) and get, (G + sC )V (K + si H ) = (G + si C )V (K + sH ) ; and multiply with a shifted and inverted matrix to get (G + si C )?1 (G + sC )V (K + si H ) = V (K + sH ) :
(16)
This is useful to get an expression for the premultiplied residual,
r(s) = (G + si C )?1 f(G + sC ) x~(s) ? rg :
(17)
Rational Krylov algorithms
7
The choice to start at v1 0 = (G + s1 C )?1 r and the relation (16) will now for
i = 1 imply that,
r(s) = (G + s1 )?1 (G + sC )Vj+1 (Kj+1;j + s1 Hj+1;j ) y~(s) ? v1 0 = Vj+1 f(Kj+1;j + sHj+1;j ) y~(s) ? e1 0 g = vj+1 (kj+1;j + shj+1;j )( y~(s))j = vj+1 hj+1;j (s ? sj )( y~(s))j ;
(18)
a unit length vector vj+1 multiplied by the last Gram Schmidt orthogonalization coecient hj+1;j , the di erence s ? sj between the evaluation frequency s and the last shift sj , and nally the last element ( y~(s))j of the solution to the j j system (15). This last factor is likely to get small when we have used shifts si reasonably close to s during the Rational Krylov process. Compare to the similar expressions one gets when applying Lanczos or Arnoldi to an eigenvalue problem or GMRES to a linear system of equations. Our way (15) to get the approximate solution actually corresponds to the FOM (Full Orthogonalization Method), which is closely related to GMRES. One advantage of using FOM, instead of the more accurate GMRES, is that we can represent the reduced order model as,
H~ (s) = ~lT (G~ + sC~ )?1 r~ ;
(19)
with ~lT = lT Vj+1 (Kj+1;j + s1 Hj+1;j ) G~ = Kj;j C~ = Hj;j r~ = e1 0 : The error of the approximate solution x~(s) can be expressed as
x~(s) = x(s) ? x~(s) = (G + sC )?1 (G + si C )r(s) = (I + (si ? s)(G + sC )?1 C )r(s);
(20)
where r(s) is the premultiplied residual (17). An upper bound for the error of the approximate transfer function H~ (s) is
j(s)j k l k j (s)j(1 + jsi ? sj k (G + sC )? C vj k): 1
+1
(21)
where (s) = hj+1;j (s ? sj )( y~(s))j , see (18). In order to compute an ecient and reliable error estimate we need to estimate k (G + sC )?1 C vj+1 k. How this is done is described in [10].
8
Axel Ruhe and Daniel Skoogh
4 Two Numerical Examples Let us report tests on two models of medium size radio frequency circuits. The rst test example is provided by Zhaojun Bai and is of dimension n = 256. The rational Krylov algorithm pis applied with the shifts 1 = ?j 108, 2 = ?j 109 and 3 = ?j 1010, where j = ?1. The shifts are applied 5, 7 and 7 times respectively. In gure 2 we show the approximate transfer function H~ (s), and in gure 3 we show the error and an error estimate of H~ (s). Note that the true error is below the error estimate. The eigenvalue distribution is given in gure 1. Compare the eigenvalue distribution with the plots for the transfer function and error estimate, and note that the transfer function has a top near the imaginary part of each eigenvalue. We did another run with the same matrices but the shifts closer together at 1 = ?j 109, 2 = ?2j 109 and 3 = ?7j 109, used 4, 4 and 8 times respectively. Look at gures 4 and 5 and note that now the errors are much smaller at the higher frequencies, where there are many poles. It is evident that an adaptive strategy where the shifts are chosen where we want small errors will be of advantage. The second test example is provided by Peter Feldmann and is of dimension n = 1841. Now the rational Krylov algorithm is applied with the shifts 1 = ?j 108, 2 = ?j 109 and 3 = ?4j 109, applied 4, 8 and 13 times respectively. In gure 6 we show the approximate transfer function H~ (s), and in gure 7 the error and an error estimate of H~ (s). The eigenvalue distribution is given in gure 8. 4
0
x 10
−0.5 −1 −1.5
real
−2 −2.5 −3 −3.5 −4 −4.5 −5 8 10
9
10 imaginary
10
10
Fig. 1. The eigenvalue distribution of the rst test example given by Bai, n = 256.
Rational Krylov algorithms
9
3
10
2
|H(j 2 π f)|
10
1
10
0
10
−1
10
8
10
9
10 2π f
10
10
Fig. 2. The rst test example of dimension n = 256. The gure shows the approximate transfer function H~ (s). The rational Krylov algorithm is applied with the shifts 1 =
?j 108 , 2 = ?j 109 and 3 = ?j 1010 , where j = p?1. The shifts are applied 5, 7 and
7 times respectively. 5
10
0
error
10
−5
10
−10
10
−15
10
8
10
9
10 2π f
10
10
Fig. 3. The error and error estimate of the rst test example. The error estimate is plotted with a solid line ?, and the true error is plotted with a dashed line ? ?
10
Axel Ruhe and Daniel Skoogh 3
10
2
|H(j 2 π f)|
10
1
10
0
10
−1
10
8
10
9
10 2π f
10
10
Fig. 4. First test example is of dimension n = 256. Approximate transfer function H~ (s9)
for rational Krylov algorithm is applied with the shifts closer together at 1 = ?j 10 , 2 = ?2j 109 and 3 = ?7j 109 , applied 4, 4 and 8 times respectively. 5
10
0
error
10
−5
10
−10
10
−15
10
8
10
9
10 2π f
10
10
Fig. 5. The error and error estimate of the rst test example with shifts closer together,
corresponding to the gure 4. The error estimate is plotted with a solid line ?, and the true error is plotted with a dashed line ? ?
Rational Krylov algorithms
11
0
10
−1
|H(j f)|
10
−2
10
−3
10
−4
10
8
10
9
10 f
10
10
Fig. 6. Second test example, dimension n = 1841. Approximate transfer function H~ (s) 8 9 for the rational Krylov algorithm applied with shifts = ? j 10 , 1 2 = ?j 10 and p 9 3 = ?4j 10 , where j = ?1. The shifts are applied 4, 8 and 13 times respectively. 5
10
0
error
10
−5
10
−10
10
−15
10
8
10
9
10 f
10
10
Fig. 7. The error and error estimate of the second test example. The error estimate is plotted with a solid line ?, and the true error is plotted with a dashed line ? ?
12
Axel Ruhe and Daniel Skoogh 9
0
x 10
real
−2 −4 −6 −8 8 10
9
10 imaginary/(2*π)
10
10
Fig. 8. The eigenvalue distribution of the second test example given by Feldmann, n = 1841.
References 1. I. M. Elfadel and D. D. Ling, A block rational Arnoldi algorithm for multipoint passive model-order reduction of multiport RLC networks, in Technical Digest of the 1997 IEEE/ACM International Conference on Computer-Aided Design, IEEE Computer Society Press, 1997, pp. 66{71. 2. , Zeros and passivity of Arnoldi-reduced-order models for interconnect networks, in Proc 34nd IEEE/ACM Design Automation Conference, ACM, New York, 1997, pp. 28{33. 3. T. Ericsson and A. Ruhe, The spectral transformation Lanczos method for the numerical solution of large sparse generalized symmetric eigenvalue problems, Mathematics of Computation, 35 (1980), pp. 1251{1268. 4. R. W. Freund, Circuit simulation techniques based on Lanczos-type algorithms, in Systems and Control in the Twenty-First Century, C. I. Byrnes, B. N. Datta, D. S. Gilliam, and C. F. Martin, eds., Birkauser, 1997, pp. 171{184. 5. , Reduced-order modeling techniques based on Krylov subspaces and their use in circuit simulation, Applied and Computational Control, Signals, and Circuits, (1998, to appear). 6. K. Gallivan, E. Grimme, and P. Van Dooren, A rational Lanczos algorithm for model reduction, Numerical Algorithms, 12 (1996), pp. 33{64. 7. H. O. Karlsson, Atomic and Molecular Density-of-States by Direct Lanczos Methods, PhD thesis, Uppsala University, Department of Quantum Chemistry, 1994. 8. A. Ruhe, Eigenvalue algorithms with several factorizations { a uni ed theory yet? to appear, 1998. 9. , Rational Krylov, a practical algorithm for large sparse nonsymmetric matrix pencils, SIAM J. Sci. Comp., 19 (1998), pp. 1535{1551. 10. D. Skoogh, Model order reduction by the rational Krylov method. work in progress, 1998.
