1Introduction ArevisedmomenterrorexpressionfortheAIRGAalgorithm

A revised moment error expression for the AIRGA algorithm

Heike Faßbender and Julius Mayer

Abstract

The fully adaptive rational global Arnoldi method (AIRGA) for the model- order reduction of second-order multi-input multi-output systems with propor- tional damping is revisited. The method automatically generates a reduced sys- tem approximating the transfer function. It is based on a moment-matching ap- proach. The expansion points are determined iteratively. The reduced order and the number of moments matched per expansion point are determined adaptively using a heuristic based on an error estimation. A revised moment error expres- sion is presented as well as some related findings.

1 Introduction

A continuous time-invariant second-order multi-input multi-output linear dynamical system is of the form

Mx(t) =¨ −D˙x(t)−Kx(t) +Fu(t),

y(t) =C_px(t) +C_vx(t),˙ (1) whereM, D, K∈R^n×n, F ∈R^n×m,C_p,C_v∈R^q×n are constant matrices. In (1), x(t)∈Rⁿis the state,u(t)∈R^mis the input, andy(t)∈R^qis the output. The mass matrixMand the stiffness matrixKneed not have any specific property (e.g., sym- metry, positive definiteness etc.), but only the special case of proportional damping is

Key Words: Model Order Reduction, Krylov Subspace, Global Arnoldi Algorithm, Moment Match- ing, Second Order, Proportional Damping.

2010 Mathematics Subject Classification: Primary 65F30; Secondary 70J50.

Received: March, 2017.

Revised: June, 2017.

Accepted: August, 2017.

87

considered. That is, the damping matrix is chosen asD=αM+βKfor some choice of realαandβ.

In many cases, the original system dimensionnis too large to allow for an effi- cient simulation of (1). Therefore, the goal of model reduction is to generate a low dimensional system that has, as best as possible, the same characteristics as the orig- inal system, but whose simulation requires significantly less computational effort.

The reduced system of (1) is described by

Mˆx(t) =¨ˆ −D˙ˆˆx(t)−Kˆx(t) +ˆ Fu(t),ˆ ˆ

y(t) =Cˆ_px(t) +ˆ Cˆ_vx(t),˙ˆ (2) where ˆM,K,ˆ Dˆ ∈R^r×r, ˆF∈R^r×m, ˆC_p,Cˆ_v∈R^q×randrn. In order to capture the relevant features of the original model, the damping matrix ˆDof the reduced order model is required to be ˆD=αMˆ+βK.ˆ

We will revisit the fully adaptive rational global Arnoldi method (AIRGA) for the model-order reduction of second-order multi-input multi-output systems with pro- portional damping [3]. This method uses a projection based on a moment-matching approach in order to compute the reduced order system. The AIRGA algorithm is recalled in Section 2. It makes use of a heuristic based on an error estimation of the moment error in order to adaptively determine the number of moments to be matched per expansion point. It turned out that the moment error given in [3] is not correct. In order to present our key findings, a revised moment error, the AIRGA algorithm has to be discussed in some detail. In particular, some technicalities from well-known facts are needed. For that matter, we also present some known facts whose proofs in the existing literature seem to be gappy. Our main result, a revised moment error, is given in Section 3. Some concluding remarks are given in Section 4. As all proofs are very technical, they have been moved to Section A.

2 AIRGA revisited

In this section we briefly review the AIRGA method [3]. This section is longer than usual because we have to present some well-known facts as we need some technical details of their proofs for further discussion. Moreover, we state some known facts whose proofs in the existing literature seem to be gappy.

Moment-Matching and Projection based Model Order Reduction

The objective is to generate a reduced order system (2) for which the first moments of the transfer function match those of the original system. The transfer functionH(s) of (1) is the linear mapping of the Laplace transformU(s)of the input u(t)to the

Laplace transformY(s)of the outputy(t),Y(s) =H(s)U(s).It is given by

H(s) = (C_p+sC_v)(s²M+sD+K)⁻¹F=:(C_p+sC_v)T(s). (3) Here and throughout the paper,s∈Chas to be chosen such thats²M+sD+K is nonsingular. The power series expansion ofT(s)around an expansion pointsi∈Cis given by (see, e.g., [10])

T(s) =

∞

∑

k=0

T^(k)(s_i)(s−s_i)^k, (4) where thek-th system momentsT^(k)(s_i)∈C^n×mats_iare given by

T⁽⁰⁾(s_i) =L⁻¹_i F,

T⁽¹⁾(s_i) =L⁻¹_i B_iT⁽⁰⁾(s_i), and fork=2,3, . . . T^(k)(s_i) =L⁻¹_i [B_iT^(k−1)(s_i)−MT^(k−2)(s_i)]

(5)

with

Li=s²_iM+siD+K and Bi=−(2s_iM+D). (6) From (4) we obtain

H(s) =

∞

∑

k=0

(C_p+sC_v)T^(k)(s_i)(s−s_i)^k

=

∞

∑

k=0

(C_p+s_iC_v)T^(k)(s_i)(s−s_i)^k+C_vT^(k)(s_i)(s−s_i)^k+1

=:

∞ k=0

∑

h_k(s_i)(s−s_i)^k

with the momentsh0(s_i) = (C_p+s_iCv)T⁽⁰⁾(s_i),and fork=1,2, . . . h_k(s_i) =C_vT^(k−1)(s_i) + (C_p+s_iC_v)T^(k)(s_i)∈C^q×m. Similarly, the transfer function of the reduced system (2) is given by

H(s) = (ˆ Cˆ_p+sCˆ_v)Tˆ(s), (7) with ˆT(s) = (s²Mˆ +sDˆ+K)ˆ ⁻¹Fˆ. Clearly, heres∈Chas to be chosen such that not onlyL=s²M+sD+Kis nonsingular, but also such that ˆL=s²Mˆ+sDˆ+Kˆ is nonsingular as well. In a projection based framework as considered below this will be satisfied automatically, as ˆL=V^HLV is nonsingular ifLis nonsingular andV is a n×rmatrix with linearly independent columns.

The power series expansion of ˆT(s)around an expansion points_i∈Cis given by Tˆ(s) =

∞

∑

k=0

Tˆ^(k)(s_i)(s−s_i)^k, (8)

where ˆT^(k)(s_i)∈C^r×mis defined analogously toT^(k)(s_i). The moments ˆh_k(s_i)of the reduced system are thus defined analogously toh_k(s_i)fork∈N0.

The goal of the moment-matching approach is to find a reduced order model such that the first few moments of (7) match those of (3), that is,

h_k(s_i) =hˆ_k(s_i) for k=0,1, . . . ,ki−1 for somek_i∈N.

A projection based method to generate a reduced oder model of orderrconstructs a projectionΠ=VV^†with a full rank matrixV ∈C^n×rand the pseudoinverseV^†= (V^HV)⁻¹V^H. SinceΠ=Π^Hholds,Πis an orthogonal projection. The reduced order model is given by

V^†(MVx(t) +¨ˆ DVx(t˙ˆ ) +KVx(t)ˆ −Fu(t)) =0, ˆ

y(t) = C_pVx(t) +ˆ C_vVx(t).˙ˆ (9) Thus, we have

Mˆ=V^†MV, Dˆ=V^†DV, Kˆ=V^†KV, Fˆ=V^†F,Cˆ_p=C_pV and ˆC_v=C_vV. (10) The following well-known theorem [5, 1, 9] states how to chooseV in order to achieve the desired moment-matching property. We restate the theorem as we will need the relation (12) later on.

Theorem 2.1. Assume siis chosen such that Li is nonsingular. Let V ∈C^n×r have linearly independent columns such that

colspan(V) ⊃ colspan([T⁽⁰⁾(s_i),T⁽¹⁾(s_i), . . . ,T^(ki−1)(s_i)]). (11) Then for the reduced order system (9) it holds that

T^(k)(s_i) =VTˆ^(k)(s_i) (12) and thus the moment-matching property h_k(s_i) =hˆ_k(s_i)holds for k=0,1, . . . ,ki−1.

First and second-order Krylov Subspace

Theorem 2.1 tells us how to chooseV. A numerically efficient and stable way to obtain suchV makes use of Krylov subspace methods.

A first-order Krylov subspaceKk(P,Q) of orderk∈Ngenerated by ann×n matrixPand ann×mmatrixQis the linear subspace spanned by the columns of the images ofQunder powers ofP

Kk(P,Q) =colspan([Q,PQ,P²Q, . . . ,P^k−1Q]).

A second-order Krylov subspaceGk(P₁,P₂,Q)of orderkforn×nmatricesP₁,P₂and ann×mmatrixQis defined as follows:

Gk(P₁,P₂,Q) =colspan([G₀,G1, . . . ,Gk−1])

withG₀=Q,G₁=P₁G₀andG_j=P₁Gj−1+P₂Gj−2, j=2,3, . . . ,k−1.

As already observed in [1], the system momentsT^(k)(s_i)are just the blocks of the second-order Krylov subspaceGki(L⁻¹_i B_i,L⁻¹_i M,L⁻¹_i F); that is

colspan([T⁽⁰⁾(s_i),T⁽¹⁾(s_i), . . . ,T^(ki−1)(s_i)]) =Gki(L⁻¹_i B_i,−L⁻¹_i M,L⁻¹_i F).

This also follows directly from (5).

For the special case of proportionally damped second-order systems, the second- order Krylov subspace is essentially identical to a first order Krylov subspace. This has already been observed, e.g., in [2, 5], but no discussion given so far seems to include all special cases. Let us first consider the following lemma (a similar relation has already been noted in [5, Section 3]); for a proof see the Appendix.

Lemma 2.2. Assume s_iis chosen such that L_iis nonsingular and s_iβ6=−1.Then L⁻¹_i B_i=−(γ_i,1I_n+γ_i,2L⁻¹_i M)

with

γi,1= β

s_iβ+1 and γi,2=s_i+ s_i+α s_iβ+1.

With the help of Lemma 2.2 the following theorem can be proven (see the Ap- pendix).

Theorem 2.3. Assume s_i is chosen such that L_i is nonsingular and s_iβ 6=−1.Let γi,1,γi,2be defined as in Lemma 2.2. Then, for any k_i∈N, it holds

Gki(L⁻¹_i B_i,−L_i⁻¹M,L⁻¹_i F) =Kki(−L⁻¹_i M,L⁻¹_i F), if γ_i,26=0, Gki(L⁻¹_i B_i,−L_i⁻¹M,L⁻¹_i F) =Kdk_i/2e(−L⁻¹_i M,L⁻¹_i F),if γi,2=0.

Remark 2.4. Assume s_iβ6=−1. Note thatγi,2=0holds iff eitherβ=0and s_i=−^α₂ orβ6=0and s_i=−β⁻¹±β⁻¹p

1−α β.

In general, the choices_iβ =−1 is not feasible. Assume for a moment thats_iβ=

−1.This impliess_i6=0 andβ 6=0. AsL_i=s²_iM+s_iD+K=s_i(s_i+α)M must be nonsingular, it further implies thatMhas to be nonsingular ands_i6=−α. Moreover,

Kki(−L⁻¹_i M,L⁻¹_i F) =Kki((s²_i +αs_i)⁻¹M⁻¹M,L⁻¹_i F) =Kki(I_n,L⁻¹_i F)

=colspan(M⁻¹F) and by an easy manipulation

Gki(L⁻¹_i B_i,−L⁻¹_i M,L⁻¹_i F) =Kki(M⁻¹K,M⁻¹F).

Therefore, unlessM=µK,µ∈RorK=0n×n, it follows fors_iβ =−1 Kki(−L⁻¹_i M,L⁻¹_i F)6⊃Gki(L⁻¹_i B_i,−L⁻¹_i M,L⁻¹_i F).

So, when an expansion points_iis chosen, it always has to be checked thats_iβ6=−1 as we will make use of Theorem 2.3 when constructing the matrixV.

The Global Arnoldi Method

Theorem 2.3 suggests to generate the desired matrixV fromKki (−L⁻¹_i M, L⁻¹_i F);

in caseγi,2=0,onlyKdk_i/2e(−L⁻¹_i M,L⁻¹_i F)has to be considered. Standard efficient and numerically sound methods to compute a basis (and thusV) of a Krylov subspace are, e.g., the block or the global Arnoldi algorithm [7, 6, 8].

The AIRGA method uses the global Arnoldi method. It constructs a basisW_i,1, W_i,2, . . . ,W_i,ki∈C^n×mof the Krylov subspaceKki(P_i,Q_i)withP_i=−L⁻¹_i M∈C^n×n andQ_i=L⁻¹_i F∈C^n×mwhich is block-orthonormal in the following sense

hW_i,j,W_i,pi=0 j6=p,

hW_i,j,W_i,pi=1 j=p for j,p=1, . . . ,k_i. (13) Here,<Y,Z>=trace(Y^HZ)whereY,Z∈C^n×s.The associated norm is the Frobe- nius normk · k_F.

In order to simplify the discussion, we assume that k_i is chosen such that the global Arnoldi algorithm does not break down; that is, for eachs_iit produces a ma- trixW_i= [W_i,1 · · · W_i,ki]∈C^n×ki·m,representing a block-orthonormal basis of the block Krylov subspaceKk_i(−L⁻¹_i M,L⁻¹_i F).Then the following relation holds for the block-orthonormal matrixW_i

P_iW_i=W_i(H_k⁽ⁱ⁾

i ⊗I_m) +h⁽ⁱ⁾_k

i+1,k_i[0, . . . ,0,W_i,ki+1], (14) with

W_i,1=Q_i/h⁽ⁱ⁾_1,0, h⁽ⁱ⁾_1,0=kQ_ik_F. (15) HereH_k⁽ⁱ⁾

i is an unreducedk_i×k_iupper Hessenberg matrix and⊗denotes the usual Kronecker product. Ifm=1, the global Arnoldi algorithm reduces to the standard Arnoldi algorithm.

Multiple expansion points

In order to ensure a good reduced model in the entire problem dependent frequency domain of interest, one usually employs not just one expansion point, but a set of` expansion points. That is, one considers a setS={s<sub>1</sub>, . . . ,s<sub>`</sub>}of`expansion points and the corresponding block Krylov subspaces

Kki(P_i,Q_i) =Kki(−L⁻¹_i M,L⁻¹_i F) for i=1, . . . , `

together with the associated block-orthonormal basisW_i∈C^n×ki·m(computed by the global Arnoldi algorithm such that eachW_isatisfies (14).). Recall that the expansion pointss_i,i=1, . . . , `,have to be chosen such thatL_i=s²_iM+s_iD+Kis nonsingular ands_iβ6=−1.

Generating the projectionΠ Let

W= [W₁ W2 ... W`]∈C^n×rmax, rmax=m

` i=1

∑

ki. (16) Clearly,

colspan(W)⊃colspan(W_i)⊃colspan(T^(k)(s_i))

fori=1,2, . . . , `andk=0,1, . . . ,k<sub>i</sub>−1. NowΠ=VV<sup>†</sup>can be set up using any full rank matrixV∈C<sup>n×r</sup>which has the same column space asW. Then, due to Theorem 2.1, the firstki moments at the expansion pointsi,i=1, . . . , `of the reduced order system (9) generated withV are matching those of the original system (1), that is,

h_j(s_i) =hˆ_j(s_i) holds for j=0, . . . ,k_i−1 and i=1,2, . . . , `.

Choosing the expansion points iteratively

Given the number`of expansion points, the setS={s<sub>1</sub>, . . .s<sub>`</sub>}of expansion points and the numberk_iof moments to be matched at eachs_i, the algorithm sketched so far will compute the desired reduced order model. As it is a priori not obvious how to choosek_i,the AIRGA algorithm [3] chooses thek_iadaptively given a fixed setSand the total number of number of moments to be matched,rmax/m. Thus, unlike as de- scribed so far, the algorithm does not generateWicorresponding toKki(L⁻¹_i M,L_i⁻¹F) at once. Instead, the following approach is used: The expansion points are picked iteratively. The first times_i is picked, just K1(L⁻¹_i M,L⁻¹_i F) is used to generate W_i∈C^n×m and just one moment is matched ats_i. The next times_i is picked, this is expanded toK2(L⁻¹_i M,L⁻¹_i F)andW_i∈C^n×2mmatching two moments ats_i, and so forth. Assume that the algorithm has picked the expansions points such that the firstk_imoments are matched at expansion points_i, that is,h_k(s_i) =hˆ_k(s_i)holds for

k=0,1, . . . ,k_i−1. The choice of the next expansion point to be considered is based on thek_i-th moment error at expansion points_i

||h_ki(s_i)−hˆ_ki(s_i)||_F=ε_ki(s_i). (17) The expansion points_pchosen next is the one corresponding to the maximum mo- ment error bys_p=argmax_s

iεk_i(s_i).

3 AIRGA revised

The idea for the adaptive choice of the expansion points is based on an expression which describes thek_i-th moment errorε_ki(s_i).The expression given in [3] is not correct. Here is the revised version of the result.

Theorem 3.1. Assume that si is chosen for all i=1, . . . , `, such that Li=s<sup>2</sup><sub>i</sub>M+ s<sub>i</sub>D+K is nonsingular and s<sub>i</sub>β6=−1.Let W<sub>i</sub>,i=1, . . . , `,be computed by the global Arnoldi method such that (14) and (15) hold. Let W= [W₁ W₂ · · · W_`]∈C^n×rmax be as in (16). Let V∈C^n×rbe a full rank matrix which has the same column space as W . Let the reduced order system (9) be generated via (10). Then the error of the k_i+1-th moment at s_ican be expressed as

ε_ki(s_i):=||h_ki(s_i)−hˆ_ki(s_i)||_F

=|γ_i,2^ki| ·

k_i

∏

h⁽ⁱ⁾_k+1,k

!

· ||(C_p+s_iC_v)

I_n−V V^†L_iV⁻¹ V^†L_i

W_i,ki+1||_F.

In order to be able to prove Theorem 3.1 the following observation which is in- spired by [4, Theorem 2] will be useful.

Lemma 3.2. Let P_i∈C^n×n, Q_i∈C^n×m, He_i=H_k⁽ⁱ⁾

i ⊗I_m and E_i=e_i⊗I_m, where ei ∈R^ki denotes the ki-th unit vector. Let Wi be computed by the global Arnoldi method such that (14) and (15) hold. Then it holds

P_i^kiQi=h⁽ⁱ⁾_1,0WiHe_i^kiE1+

k_i k=0

∏

h⁽ⁱ⁾_k+1,k

! Wi,k_i+1.

Theorem 3.1 gives rise to Algorithm 1. It starts with an initial set of expan- sion points and automatically and adaptively chooses the number of moments to be matched at each expansion points_ibased on Theorem 3.1. This is controlled by the inner while loop starting at line 10 whereV is computed.

One can use different methods to obtain a full rank matrixV∈C^n×rwhich has the same column space asW ∈C^n×rmax. A numerically safe way to generate the matrix V fromW is the use of the rank-revealing QR-decomposition ofW.The relevant part

Algorithm 1Revised Adaptive Iterative Rational Global Arnoldi Algorithm

1: Input:M, D,K, F,Cp,Cv, α,βsuch thatD=αM+βK, r_max, initial set of exp. pointsS={s₁, . . . ,s_`}

2: Output:r,M,ˆ D,ˆ K,ˆ F,ˆ Cˆ_p,Cˆ_v,such that ˆD=αMˆ+βKˆ

3: whileno convergencedo

4: W= [ ], seq= [ ], r_W =0, [n,m] =size(F)

5: fori=1 :`do

6: L_i=s²_iM+s_iD+K

7: R_i=L⁻¹_i F

8: h_i=||R_i||_F, R_i=R_i/h_i, h_i,Π=h_i, γ_i,Π=1, γ_i,2=s_i+_s^si+α

iβ+1

9: end for

10: whilerW≤rmax−mand no convergencedo

11: ifrW =0then

12: p=argmax_i||γ_i,Πh_i,Π(C_p+s_iC_v)R_i||_F

13: else

14: p=argmax_i||γ_i,Πh_i,Π(C_p+s_iC_v)

I_n−V V^HL_iV−1

V^HL_i R_i||_F

15: end if

16: seq= [seq,p], W= [W R_p], r_W=r_W+m

17: R_p=−L⁻¹_p M R_p

18: fork=1 : length(seq)do

19: ifseq(k) =pthen

20: h=trace(W(:,(k−1)m+1 :km)^HR_p)

21: R_p=R_p−hW(:,(k−1)m+1 :km)

22: end if

23: end for

24: h_p=||R_p||_F, h_p,Π=h_p,Π·h_p, γp,Π=γp,Π·γp,2 25: ifh_p6=0then

26: R_p=R_p/h_p

27: end if

28: DetermineV from [Q,R,E] = qr(W,0) (to deflate all linear

29: dependent columns)

30: end while

31: Choose new set of expansion pointsS={s1, . . . ,s_`}

32: end while

33: DetermineV from [Q,R,E] = qr([Re(V), Im(V)],0)

34: r = size(V,2)

35: Compute the reduced order system as in (18)

of its unitary factor is then used asV such thatV has orthonormal columns. Thus it holdsV^†=V^H and the projectionΠ=VV^†=VV^H becomes a Galerkin projection.

The system matrices of the reduced order system are given by

Mˆ =V^HMV, Dˆ =V^HDV, Kˆ=V^HKV,Fˆ=V^HF,Cˆ_p=C_pV and ˆC_v=C_vV. (18) The size of the resulting reduced system can not be predetermined. At the end,Vwill havermaxor less columns.

Please note thath_pin line 24 corresponds to a lower subdiagonal element of the associated Hessenberg matrix. In caseh_p=0, we haveR_p=0n×m.The algorithm does not break down, as this implies that the corresponding moment error is equal to zero. Thus, the corresponding expansion point will not be chosen again.

The quality of the reduced order system heavily depends on the choice of the expansion points. As a good set of expansion points is usually not available, typically such a set is determined iteratively. This is controlled by the outer while loop starting at line 3 of Algorithm 1. One starts with an initial set of expansion points, computes the corresponding reduced order model and selects a new set of expansion points, e.g., based on the eigenvalues ofλ²Mˆ+λD+ˆ K. The actual selection criterion has toˆ be based on the problem considered, see, e.g., the discussion in [3] and the references therein. This process is repeated till convergence, measured, e.g., in terms of the H2-error between the previously computed reduced order system and the current one

err=kHˆ_previous−Hˆ_currentk_H2.

HerekH−Gk_H2 = _2π^{1 R}_−∞^∞ kH(ıω)−G(ıω)k_Fdω, where the transfer functions H andGbelong to systems with the same input and output dimension. For a thorough discussion on how to determine convergence as well as a new set of expansion points in the outer loop see [3].

Finally note that allowing complex-valued expansion pointss_ileads toW∈C^n×r. ThusV and the reduced order system (9) is complex-valued, even though usually a real-valued one is desired. Using complex-conjugate pairs of expansion points, at least theoretically, the entire computations can be done in real arithmetic. A different option is to split the complex-valued matrixV into its real and imaginary part and to use a rank-revealing QR decomposition of [Re(V), Im(V)] to obtain a real matrix with orthonormal columns and the same column space. This real-valued matrix may have twice the number of columns as desired. Hence, the dimension of the reduced order system will be doubled. The number of moments matched does not change.

4 Concluding Remarks

In [3] the AIRGA algorithm for model-order reduction of second-order multi-input multi-output systems with proportional damping has been proposed. The method

relies on the moment errorεk_i(s_i)as in (17). Unfortunately the expression for the mo- ment error given in [3] is not correct. Section 3 presents our main contribution: the revised moment error given in Theorem 3.1 as well as Lemma 3.2 which is needed to proof Theorem 3.1. The idea of the AIRGA algorithm and all further details neces- sary to proof Theorem 3.1 have been summarized in Section 2. In doing so, we have rounded off some of the results employed (by explicitily stating all assumptions and findings which may have been obscured in previous publications). In [3] numerical examples have been considered. Repeating those with the revised moment error do not reveal any major differences in the results. Thus, these experiments have not been included here.

A Proofs

In this section we provide the details of the proofs for the theorems and lemmas of Sections 2 and 3. For ease of reference, the theorems and lemmas are restated except for Theorem 2.1 which is well-known. Lemma 2.2 and Theorem 2.3 appeared in slightly different form in, e.g., [2, 5], while Lemma 3.2 and Theorem 3.1 are new.

Lemma 2.2 Assume s_iis chosen such that L_iis nonsingular and s_iβ 6=−1.Then L⁻¹_i B_i=−(γi,1I_n+γ_i,2L⁻¹_i M)

with

γi,1= β

s_iβ+1 and γi,2=s_i+ s_i+α s_iβ+1. Proof.

−L⁻¹_i Bi=L⁻¹_i (2s_iM+D) =L⁻¹_i (2s_iM+αM+βK)

= (2s_i+α)L⁻¹_i M+βL⁻¹_i K

= (2s_i+α)L⁻¹_i M+ β

s_iβ+1(s_iβ+1)L⁻¹_i K

= (2s_i+α)L⁻¹_i M+ β

s_iβ+1L⁻¹_i −(s²_i +s_iα)M+ (s²_i+s_iα)M+ (s_iβ+1)K

=

2s_i+α−β(s²_i +s_iα) s_iβ+1

L⁻¹_i M+ β

s_iβ+1L⁻¹_i s²_iM+siD+K

=

s_i+ s_i+α siβ+1

L⁻¹_i M+ β

siβ+1L⁻¹_i L_i

=γi,2L⁻¹_i M+γ_i,1I_n.

Theorem 2.3Assume s_i is chosen such that L_i is nonsingular and s_iβ 6=−1. Let γi,1,γi,2be defined as in Lemma 2.2. Then for any k_i∈Nit holds

Gki(L⁻¹_i B_i,−L_i⁻¹M,L⁻¹_i F) =Kki(−L⁻¹_i M,L⁻¹_i F), if γi,26=0, Gk_i(L⁻¹_i B_i,−L_i⁻¹M,L⁻¹_i F) =Kdk_i/2e(−L⁻¹_i M,L⁻¹_i F),if γi,2=0.

Proof. Fromγi,2=s_i+_s^si+α

iβ+1it follows that forβ6=0 γi,2=s²_iβ+2s_i+α

s_iβ+1 =s²_iβ²+2s_iβ+α β β(s_iβ+1) =0 iff the numerators²_iβ²+2s_iβ+α β is zero. This yieldss_iβ=−1±p

1−α β. Assumeγi,26=0. SetP:=−L⁻¹_i M,Q:=L⁻¹_i Fsuch that the blocks of the Krylov subspaceKk(−L⁻¹_i M,L⁻¹_i F)are given asQ, PQ,P²Q, . . . ,P^k−1Q. Since

T⁽⁰⁾(s_i) =L⁻¹_i F=Q,

we haveG1(L⁻¹_i B_i,−L⁻¹_i M,L⁻¹_i F) =K1(−L⁻¹_i M,L⁻¹_i F). Next with Lemma 2.2, it holds

T⁽¹⁾(s_i) =L_i⁻¹BiL⁻¹_i F=−(γ_i,1In+γi,2L⁻¹_i M)L⁻¹_i F=−γ_i,1Q−γi,2PQ.

Thus, asγi,26=0 we haveG2(L⁻¹_i B_i,−L⁻¹_i M,L⁻¹_i F) =K2(−L⁻¹_i M,L⁻¹_i F).

Now, assumeGj(L⁻¹_i B_i,−L⁻¹_i M,L⁻¹_i F) =Kj(−L⁻¹_i M,L⁻¹_i F)for j=1,2, . . . ,p.

Then we can find µ_i^(j−1,k) ∈Cfor k=0,1, . . . ,j−1 and j=1,2, . . . ,p such that T^(j−1)(s_i) =∑_k=0^j−1µ_i^(j−1,k)P^kQ.With Lemma 2.2 it follows

T^(p)(s_i) =L⁻¹_i B_iT^(p−1)(s_i)−L⁻¹_i MT^(p−2)(s_i)

=−(γ_i,1In+γi,2L⁻¹_i M)T^(p−1)(s_i)−L⁻¹_i MT^(p−2)(s_i)

=−γ_i,1T^(p−1)(s_i)−γi,2PT^(p−1)(s_i)−PT^(p−2)(s_i)

=−γ_i,1

p−1 k=0

∑

µ_i^(p−1,k)P^kQ−γi,2P

p−1 k=0

∑

µ_i^(p−1,k)P^kQ−P

p−2 k=0

∑

µ_i^(p−2,k)P^kQ

=−γ_i,1

p−1

∑

k=0

µ_i^(p−1,k)P^kQ−γi,2 p

∑

k=1

µ_i^{(p−1,k−1)}P^kQ−

p−1

∑

k=1

µ_i^{(p−2,k−1)}P^kQ

=−γ_i,1µ_i^(p−1,0)Q−

p−1

∑

k=1

(γ_i,1µ_i^(p−1,k)+γi,2µ_i^{(p−1,k−1)}+µ_i^{(p−2,k−1)})P^kQ

−γi,2µ_i^{(p−1,p−1)}P^pQ

=:

p

∑

k=0

µ_i^(p,k)P^kQ. (19)

The above directly reveals the recursion formula µ_i^(p,0)=−γi,1µ_i^(p−1,0)

µ_i^(p,k)=−γ_i,1µ_i^(p−1,k)−γ_i,2µ_i^{(p−1,k−1)}−µ_i^{(z−2,k−1)}, fork=1,2, . . . ,p−1 µ_i^(p,p)=−γi,2µ_i^{(p−1,p−1)}

for anyp≥2 withµ_i^(0,0)=1,µ_i^(1,0)=−γi,1andµ_i^(1,1)=−γi,2. Particularly, it holds

µ_i^(p,p)= (−γ_i,2)^p. (20)

Asγi,26=0, we immediately haveGk(L⁻¹_i B_i,−L⁻¹_i M,L⁻¹_i F) =Kk(−L⁻¹_i M,L⁻¹_i F), so that the first equation of the theorem is proven by induction.

In order to prove the second statement of the theorem, assume γi,2=0. With Lemma 2.2, it follows

Gki(L⁻¹_i B_i,−L⁻¹_i M,L⁻¹_i F) =Gki(−γ_i,1I_n,−L⁻¹_i M,L⁻¹_i F) =Kdk_i/2e(−L⁻¹_i M,L⁻¹_i F).

Lemma 3.2Let P_i∈C^n×n,Q_i∈C^n×m,He_i=H_k⁽ⁱ⁾

i ⊗I_mand E_i=e_i⊗I_m, where e_i∈R^ki denotes the k_i-th unit vector. Let W_ibe computed by the global Arnoldi method such that (14) and (15) hold. Then it holds

P_i^kiQ_i=h⁽ⁱ⁾_1,0W_iHe_i^kiE₁+

ki

∏

k=0

h⁽ⁱ⁾_k+1,k

! W_i,ki+1.

Proof. Observe PiWi=Wi(H_k⁽ⁱ⁾

i ⊗Im) +h⁽ⁱ⁾_k

i+1,k_i[0, . . . ,0,W_i,ki+1] =WiH˜i+h⁽ⁱ⁾_k

i+1,k_iWi,k_i+1E_k^T_i. (21)

Multiplication from the left byP_i^ki−1and repeated use of (21) yields P_i^kiW_i=P_i^ki−2(P_iW_i)He_i+P_i^ki−1h⁽ⁱ⁾_k

i+1,k_iW_i,ki+1E_k^T

i

=P_i^ki−2

W_iHe_i+h⁽ⁱ⁾_k

i+1,k_iW_i,ki+1E_k^T

i

He_i+h⁽ⁱ⁾_k

i+1,k_iP_i^ki−1W_i,ki+1E_k^T

i

=P_i^ki−3(P_iW_i)He_i²+h⁽ⁱ⁾_k

i+1,k_i 1

∑

k=0

P^ki−1−kW_i,ki+1E_k^T

iHe_i^k

=P_i^ki−3

WiHei+h⁽ⁱ⁾_k

i+1,kiWi,k_i+1E_k^T_i

He_i² +h⁽ⁱ⁾_k

i+1,k_i 1

∑

k=0

P^ki−1−kW_i,ki+1E_k^T

iHe_i^k

=P_i^ki−4(P_iW_i)He_i³+h⁽ⁱ⁾_k

i+1,k_i 2 k=0

∑

P^ki−1−kW_i,ki+1E_k^T

iHe_i^k

=. . .

=W_iHe_i^ki+h⁽ⁱ⁾_k

i+1,k_i ki−1 k=0

∑

P^ki−1−kW_i,ki+1E_k^T

iHe_i^k.

AsQ_i=h⁽ⁱ⁾_1,0W_i,1=h⁽ⁱ⁾_1,0W_iE₁, we have

P_i^kiQ_i=h⁽ⁱ⁾_1,0P_i^kiW_iE₁=h⁽ⁱ⁾_1,0W_iHe_i^kiE₁+h⁽ⁱ⁾_1,0h⁽ⁱ⁾_k

i+1,k_i ki−1

k=0

∑

P^ki−1−kW_i,ki+1E_k^T

iHe_i^kE₁. (22)

SinceH_iis an upper Hessenberg matrix,He_i is a block upper Hessenberg matrix. It follows that

E_k^T

iHe_i^pE₁=0 for p=0,1, . . . ,k_i−2, and E_k^T

iHe_i^ki−1E₁=

k_i−1

∏

h⁽ⁱ⁾_k+1,kI_m.

Substituting this into (22) gives

P_i^kiQ_i=h⁽ⁱ⁾_1,0W_iHe_i^kiE1+h⁽ⁱ⁾_1,0h⁽ⁱ⁾_k

i+1,k_iW_i,ki+1

ki−1

∏

k=1

h⁽ⁱ⁾_k+1,kI_m

=h⁽ⁱ⁾_1,0WiHe_i^kiE1+

ki

k=0

∏

h⁽ⁱ⁾_k+1,kWi,k_i+1.

Theorem 3.1Assume that s_iis chosen for all i=1, . . . , `,such that L<sub>i</sub>=s<sup>2</sup><sub>i</sub>M+s<sub>i</sub>D+K is nonsingular and s<sub>i</sub>β 6=−1.Let W<sub>i</sub>,i=1, . . . , `,be computed by the global Arnoldi method such that (14) and (15) hold. Let W= [W₁ W₂ · · · W_`]∈C^n×rmax be as in (16). Let V∈C^n×rbe a full rank matrix which has the same column space as W . Let

the reduced order system (9) be generated via (10). Then the error of the k_i+1-th moment at s_ican be expressed as

ε_ki(s_i):=||h_ki(s_i)−hˆ_ki(s_i)||_F

=|γ_i,2^ki| ·

k_i

∏

h⁽ⁱ⁾_k+1,k

!

· ||(C_p+s_iC_v)

I_n−V V^†L_iV⁻¹ V^†L_i

W_i,ki+1||_F.

Proof. First considerk_i=0. Recall ˆT⁽⁰⁾(s_i) =Lˆ⁻¹_i F.ˆ

ε₀(s_i) =||(C_p+s_iC_v)T⁽⁰⁾(s_i)−(Cˆ_p+s_iCˆ_v)Tˆ⁽⁰⁾(s_i)||_F

=||(C_p+s_iC_v)T⁽⁰⁾(s_i)−(C_p+s_iC_v)VTˆ⁽⁰⁾||_F

=||(C_p+s_iC_v)

T⁽⁰⁾(s_i)−VLˆ_i⁻¹Fˆ

||_F

=||(C_p+s_iC_v)

T⁽⁰⁾(s_i)−V V^†L_iV⁻¹ (V^†F)

||_F

=||(C_p+s_iC_v)

T⁽⁰⁾(s_i)−V V^†L_iV⁻¹

V^†(L_iL⁻¹_i )F

||_F

=||(C_p+s_iC_v)

I_n−V V^†L_iV⁻¹ V^†L_i

T⁽⁰⁾(s_i)||_F. (23) Next considerk_i=1. Recall ˆT⁽¹⁾(s_i) =Lˆ⁻¹_i Bˆ_iTˆ⁽⁰⁾(s_i),andT⁽⁰⁾(s_i) =VTˆ⁽⁰⁾(s_i),but in generalT⁽¹⁾(s_i)6=VTˆ⁽¹⁾(s_i).

ε1(s_i) =||(C_p+s_iC_v)T⁽¹⁾(s_i)−(Cˆ_p+s_iCˆ_v)Tˆ⁽¹⁾(s_i)||_F

=||(C_p+s_iC_v)T⁽¹⁾(s_i)−(C_p+s_iC_v)VLˆ⁻¹_i Bˆ_iTˆ⁽⁰⁾(s_i)||_F

=||(C_p+s_iC_v)

T⁽¹⁾(s_i)−V V^†L_iV⁻¹

V^†B_iVTˆ⁽⁰⁾(s_i)

||_F

=||(C_p+s_iC_v)

T⁽¹⁾(s_i)−V V^†L_iV⁻¹

V^†B_i

VTˆ⁽⁰⁾(s_i)

||_F

=||(C_p+s_iC_v)

T⁽¹⁾(s_i)−V V^†L_iV−1

V^†B_i

T⁽⁰⁾(s_i)

||_F

=||(C_p+s_iC_v)

T⁽¹⁾(s_i)−V V^†L_iV⁻¹

V^†(L_iL⁻¹_i )B_iT⁽⁰⁾(s_i)

||_F

=||(C_p+s_iC_v)

I_n−V V^†L_iV⁻¹ V^†L_i

T⁽¹⁾(s_i)||_F. (24) Considerk_i>1.Recall ˆT^(ki)(s_i) =Lˆ⁻¹_i

Bˆ_iTˆ^(ki−1)(s_i)−MˆTˆ^(ki−2)(s_i)

andT^(j)(s_i) = VTˆ^(j)(s_i),for j=1, . . . ,k_i−1,but in generalT^(ki)(s_i)6=VTˆ^(ki)(s_i).

ε_ki(s_i) =||(C_p+s_iC_v)T^(ki)(s_i)−(Cˆ_p+s_iCˆ_v)Tˆ^(ki)(s_i)||_F