Inexact Linear Solves In Model Reduction of Bilinear Dynamical Systems
aa r X i v : . [ m a t h . NA ] M a r An. S¸t. Univ. Ovidius Constant¸a Vol. xx ( x ), , INEXACT LINEAR SOLVES IN MODELREDUCTION OF BILINEAR DYNAMICALSYSTEMS
Rajendra Choudhary and Kapil Ahuja
Abstract
Bilinear dynamical systems are commonly used in science and engi-neering because they form a bridge between linear and non-linear sys-tems. However, simulating them is still a challenge because of their largesize. Hence, a lot of research is currently being done for reducing suchbilinear dynamical systems (termed as bilinear model order reduction orbilinear MOR). Bilinear iterative rational Krylov algorithm (BIRKA) isa very popular, standard and mathematically sound algorithm for bilin-ear MOR, which is based upon interpolatory projection technique. Anefficient variant of BIRKA, Truncated BIRKA (or TBIRKA) has alsobeen recently proposed.Like for any MOR algorithm, these two algorithms also require solv-ing multiple linear systems as part of the model reduction process. Forreducing very large dynamical systems, which is now-a-days becominga norm, scaling of such linear systems with respect to input dynamicalsystem size is a bottleneck. For efficiency, these linear systems are oftensolved by an iterative solver, which introduces approximation errors.Hence, stability analysis of MOR algorithms with respect to inexactlinear solves is important. In our past work, we have shown that un-der mild conditions, BIRKA is stable (in the sense as discussed above).Here, we look at stability of TBIRKA in the same context.Besides deriving the conditions for a stable TBIRKA, our other novelcontribution is the more intuitive methodology for achieving this. Thisapproach exploits the fact that in TBIRKA a bilinear dynamical systemcan be represented by a finite set of functions, which was not possible
Key Words: Model order reduction, Bilinear dynamical systems, Interpolatory projec-tion, Stability analysis, Backward stability, Perturbation analysis, Volterra series interpolation.2010 Mathematics Subject Classification: Primary 34D10, 34D20, 41A05, 65F10. NEXACT LINEAR SOLVES IN MODEL REDUCTION OF BILINEARDYNAMICAL SYSTEMS in BIRKA (because infinite such functions were needed there). Thestability analysis techniques that we propose here can be extended tomany other methods for doing MOR of bilinear dynamical systems, e.g.,using balanced truncation or the ADI methods. A dynamical system, usually represented by a set of differential equations,can be linear or non-linear [1, 2, 3]. Linear dynamical systems have beenstudied more than the non-linear ones because of the obvious ease in workingwith them. Bilinear dynamical systems form a good bridge between the linearand the non-linear cases, and are usually approximated by a varying degreeof bilinearity [4, 5, 6]. In this manuscript, we focus on bilinear dynamicalsystems.Dynamical systems coming from different real world applications are verylarge in size. Thus, simulation and computation with such systems is pro-hibitively expensive. Model Order Reduction (MOR) techniques provide asmaller system that besides being cheaper to work with, also replicates theinput-output behaviour of the original system to a great extent [7, 8, 9, 10,11, 12, 13, 14, 15, 16]. Since, bilinear dynamical systems have been recentlystudied, the techniques for reducing them are also recent.Out of the many methods available for performing bilinear MOR [14, 15,17, 18, 19, 20, 21, 22], we focus on a commonly used interpolatory projectionmethod. BIRKA (Bilinear Iterative Rational Krylov Algorithm) [15] is a verypopular algorithm based upon this technique for reducing first-order bilineardynamical systems ∗ . BIRKA’s biggest drawback is that it does not scale wellin time (with respect to increase in the size of the input dynamical system).A cheaper variant of BIRKA, called TBIRKA (Truncated Bilinear IterativeRational Krylov Algorithm) [21, 22] has also been proposed.Like in any other MOR algorithm, in BIRKA and TBIRKA also, peopleoften use direct methods like LU-factorization, etc., to solve the arising linearsystems, which are too expensive (i.e., computational complexity of O ( n ),where n is the original system size). A common solution to this scaling problemis to use iterative methods like the Krylov subspace methods, etc., whichhave a reduced computational complexity (i.e., O ( n × nnz ), where nnz isthe number of nonzeros in the system matrix) [23, 24]. Although iterativemethods are cheap, they are inexact too. Hence, studying stability of theunderlying MOR algorithm (here BIRKA and TBIRKA) with respect to suchapproximate linear solves becomes important [25, 26]. ∗ First-order implies that the highest derivative of the state variable in the dynamicalsystem is one.
Second-order and higher-orders are similarly defined.
Rajendra Choudhary and Kapil Ahuja
One of the first works that performed such a stability analysis focused onpopular MOR algorithms for first-order linear dynamical systems [27]. Here,the authors briefly mention that their analysis would be easily carried fromthe first-order to the second-order case. A detailed stability analysis focusingon reducing second-order linear dynamical systems has been done in [28]. Adifferent kind of stability analysis for MOR of second-order linear dynamicalsystems has been done in [29]. In this, the authors first show that the SOARalgorithm (second order Arnoldi) is unstable with respect to the machine pre-cision errors (and not inexact linear solves). Then, they propose a Two-levelorthogonal Arnoldi (TOAR) algorithm that cures this instability of SOAR. Anextended stability analysis for BIRKA (as above, a popular MOR algorithmfor first-order bilinear dynamical systems) has been recently done in [30]. Forrest of this manuscript, whenever stability analysis is referred, we mean it withrespect to inexact linear solves.We follow the stability analysis framework of BIRKA from [30] and pro-pose equivalent theorems for TBIRKA. The approach here is slightly differ-ent, which forms our most novel contribution. Norm of the dynamical systemplays an important role in stability analysis (the kind of norm is discussedlater). In BIRKA stability analysis, a single expression for bilinear dynamicalsystem norm is used (involving a Volterra series). In TBIRKA stability anal-ysis, a similar single expression (involving a truncated Volterra series) leadsto complications. Alternatively, in TBIRKA, because of truncation, the bi-linear dynamical system can be represented by a finite set of functions. Thiswas not possible in BIRKA where infinite such functions were needed. Thus,in TBIRKA stability analysis, we use norm of all such functions leading toeasier derivations. Our stability analysis, as done for BIRKA earlier and forTBIRKA here, can be easily extended to other MOR algorithms for bilineardynamical systems, e.g., projection based [14], implicit Volterra series [17],balanced truncation [18], gramian based [20], etc.The rest of the paper is divided into three more parts. In Section 2, wefirst give a brief overview of MOR for bilinear dynamical systems using aprojection method. Next, we review the stability analysis of BIRKA from[30]. A stability analysis typically involves satisfying two conditions. Hence,finally in this section, we study the first condition for a stable TBIRKA. Weanalyze the second condition of stability in TBIRKA’s context in Section 3.In Section 4, we give conclusions and discuss the future work. For the rest ofthis paper, we use the terms and notations as listed below.a. The H − norm is a functional norm defined as [27, 21, 22] k H k k H = (cid:18) π (cid:19) k Z ∞−∞ . . . Z ∞−∞ k H k ( iω , . . . , iω k ) k F dω . . . dω k , (1) NEXACT LINEAR SOLVES IN MODEL REDUCTION OF BILINEARDYNAMICAL SYSTEMS i denotes √−
1. Here, we assume that all H − norms computedfurther exist. In other words, the improper integrals defined by the H − norm give finite value. This is a reasonable assumption becausethis happens often in practice (see [27], where stability analysis of IRKAis done as well as [30], where stability analysis of BIRKA is done).b. The H ∞ − norm is also a functional norm, defined as [27, 21, 22] k H k k H ∞ = max ω , ..., ω k ∈ R k H k ( iω , . . . , iω k ) k . c. The Kronecker product between two matrices P (of size m × n ), and Q (of size s × t ) is defined as P ⊗ Q = p Q · · · p n Q ... . . . ... p m Q · · · p mn Q , where p ij is an element of matrix P and order of P ⊗ Q is ms × nt .d. vec operator on a matrix P is defined as vec ( P ) = ( p , . . . , p m , p , . . . , p m , . . . . . . , p n , . . . , p mn ) T . e. Also, I n denotes an identity matrix of size n × n and R denotes the setof real numbers. A first-order bilinear dynamical system is usually represented as [14, 15] ζ : (cid:26) ˙ x ( t ) = Ax ( t ) + N x ( t ) u ( t ) + bu ( t ) ,y ( t ) = cx ( t ) , (2)where A, N ∈ R n × n , b ∈ R n × and c ∈ R × n . Also, u ( t ) : R → R , y ( t ) : R → R and x ( t ) : R → R n , represent input, output and state of the bilineardynamical system, respectively. This is a Single Input Single Output (SISO)system, which we have chosen for ease of our analysis. We plan to look atMultiple Input Multiple Output (MIMO) systems as part of our future work.A bilinear dynamical system can also be represented by a series of transferfunctions, i.e., ζ = Lim k →∞ ζ k , (3) Rajendra Choudhary and Kapil Ahuja where ζ k = { H ( s ) , H ( s , s ) , . . . , H k ( s , s , . . . , s k ) } . Here, H k ( s , s , . . . , s k ) is called the k th order transfer function of the bilineardynamical system and is defined as H k ( s , . . . , s k ) = c ( s k I − A ) − N ( s k − I − A ) − . . . N ( s I − A ) − b. (4)After reduction, the bilinear dynamical system (2) can be represented as ζ r : (cid:26) ˙ x r ( t ) = A r x r ( t ) + N r x r ( t ) u ( t ) + b r u ( t ) ,y r ( t ) = c r x r ( t ) , (5)where A r , N r ∈ R r × r , b r ∈ R r × and c r ∈ R × r with r ≪ n . The main goalof model reduction is to approximate ζ by ζ r in an appropriate norm, suchthat for all admissible inputs, y r ( t ) is nearly same to y ( t ).As mentioned earlier, we use interpolatory projection for performing modelreduction. The two common and standard algorithms here, BIRKA andTBIRKA, use a Petrov-Galerkin projection. Let V r and W r be the two r -dimensional subspaces, such that V r = Range ( V r ) and W r = Range ( W r ),where V r ∈ R n × r and W r ∈ R n × r are matrices. Also, let (cid:0) W Tr V r (cid:1) be in-vertible † . Applying the projection x ( t ) = V r x r ( t ), and enforcing the Petrov-Galerkin conditions [15, 22] on the original bilinear dynamical system (2), weget the reduced system as W Tr ( V r ˙ x r ( t ) − AV r x r ( t ) − N V r x r ( t ) u ( t ) − bu ( t )) = 0 ,y ( t ) = cV r x r ( t ) . Comparing the above two equations with their respective equations in (5), weget a relation between the original system matrices and the reduced systemmatrices, i.e., A r = (cid:0) W Tr V r (cid:1) − W Tr AV r , N r = (cid:0) W Tr V r (cid:1) − W Tr N V r ,b r = (cid:0) W Tr V r (cid:1) − W Tr b, and c r = cV r . One way of obtaining subspaces V r and W r is to use Volterra series inter-polation. Further, to decide where to interpolate so as to obtain an optimalreduced model, an H − optimization problem is commonly solved (Theorem4.7 from [21]). Theorem 1. [21] Let ζ be a bilinear system of order n. Let ζ r be an H − optimal approximation of order r. Then, ζ r satisfies the following multi-point † Obtaining such an invertible matrix is not difficult [15, 22].
NEXACT LINEAR SOLVES IN MODEL REDUCTION OF BILINEARDYNAMICAL SYSTEMS Volterra series interpolation conditions: ∞ X k =1 r X l =1 . . . r X l k =1 φ l , l , ..., l k H k ( − λ l , − λ l , . . . , − λ l k ) = ∞ X k =1 r X l =1 . . . r X l k =1 φ l , l , ..., l k H r k ( − λ l , − λ l , . . . , − λ l k ) , and ∞ X k =1 r X l =1 . . . r X l k =1 φ l , l , ..., l k k X j =1 ∂∂s j H k ( − λ l , − λ l , . . . , − λ l k ) = ∞ X k =1 r X l =1 . . . r X l k =1 φ l , l , ..., l k k X j =1 ∂∂s j H r k ( − λ l , − λ l , . . . , − λ l k ) , where φ l , l , ..., l k and λ l , λ l , . . . , λ l k are residues and poles of the transferfunction H r k associated with ζ r , respectively. BIRKA is designed in such a way that at convergence, the conditions ofTheorem 1 are satisfied leading to a locally H − optimal reduced model. Al-gorithm 1 lists BIRKA. Algorithm 1
BIRKA [15, 21, 22] Given an input bilinear dynamical system
A, N, b, c . Select an initial guess for the reduced system as ˇ A, ˇ N , ˇ b, ˇ c . Also selectstopping tolerance btol . While (cid:0) relative change in eigenvalues of ˇ A ≥ btol (cid:1) a. R Λ R − = ˇ A, ˇˇ b = ˇ b T R − T , ˇˇ c = ˇ cR, ˇˇ N = R T ˇ N R − T .b. vec ( V ) = (cid:16) − Λ ⊗ I n − I r ⊗ A − ˇˇ N T ⊗ N (cid:17) − (cid:16) ˇˇ b T ⊗ b (cid:17) .c. vec ( W ) = (cid:16) − Λ ⊗ I n − I r ⊗ A T − ˇˇ N ⊗ N T (cid:17) − (cid:0) ˇˇ c T ⊗ c T (cid:1) .d. V r = orth ( V ) , W r = orth ( W ).e. ˇ A = ( W Tr V r ) − W Tr A V r , ˇ N = (cid:0) W Tr V r (cid:1) − W Tr N V r , ˇ b = (cid:0) W Tr V r (cid:1) − W Tr b, ˇ c = c V r . A r = ˇ A, N r = ˇ N , b r = ˇ b, c r = ˇ c .TBIRKA is similar to BIRKA in most aspects except that it performs atruncated Volterra series interpolation. Here, instead of ζ in (2)-(3), they work Rajendra Choudhary and Kapil Ahuja with ζ M , which is defined as ζ M = { H ( s ) , H ( s , s ) , H ( s , s , s ) , . . . , H M ( s , . . . , s M ) } , (6)with H k ( s , . . . , s k ) for k = 1 , . . . , M is given by (4). Equivalent of Theorem1 above is as follows (Theorem 4.8 from [21]): Theorem 2. [21] Let ζ = ( A, N, b, c ) be an order n bilinear system and ζ M bethe polynomial system determined by ζ . Let ζ r = ( A r , N r , b r , c r ) be a bilinearsystem of order r , and define ζ Mr as the polynomial system determined by ζ r .Suppose that ζ Mr is an H − optimal approximation to ζ M . Then ζ Mr satisfies M X k =1 r X l =1 . . . r X l k =1 φ l , l , ..., l k H k ( − λ l , − λ l , . . . , − λ l k ) = M X k =1 r X l =1 . . . r X l k =1 φ l , l , ..., l k H r k ( − λ l , − λ l , . . . , − λ l k ) , and M X k =1 r X l =1 . . . r X l k =1 φ l , l , ..., l k k X j =1 ∂∂s j H k ( − λ l , − λ l , . . . , − λ l k ) = M X k =1 r X l =1 . . . r X l k =1 φ l , l , ..., l k k X j =1 ∂∂s j H r k ( − λ l , − λ l , . . . , − λ l k ) , where φ l , l , ..., l k and λ l , λ l , . . . , λ l k are residues and poles of the transferfunction H r k associated with ζ Mr , respectively. Algorithm 2 lists TBIRKA.Both BIRKA and TBIRKA in turn require solving large sparse linear sys-tems of equations. If we compare Algorithm 1 and 2, we realize that the costof linear solvers at each step of the
While loop in the former is O ( nr × nr ) andin the latter is O ( M × nr × nr ). This makes it seem that TBIRKA is moreexpensive than BIRKA. However, TBIRKA is implemented in such a way thatthe Kronecker products are avoided leading to cost of linear solves at each stepof the While loop to be O ( M × r × ( n × n )). Thus, if M < r , which is usuallythe case (e.g, M = 3 or 4 and r ≫ M ), then TBIRKA is more efficient thanBIRKA. For further details on this see chapter 4 in [21] and Section 5.3 in [22].These implementation details do not affect our stability analysis, and hence,we use Algorithm 2 in the current form as our base.As mentioned earlier, using iterative methods for solving such linear sys-tems introduces approximation errors. We have done a detailed stability anal-ysis of BIRKA with respect to the inexact linear solves in [30], and we briefly NEXACT LINEAR SOLVES IN MODEL REDUCTION OF BILINEARDYNAMICAL SYSTEMS Algorithm 2
TBIRKA [21, 22] Given an input bilinear dynamical system
A, N, b, c . Select an initial guess for the reduced system as ˇ A, ˇ N , ˇ b, ˇ c . Also selectthe truncation index M and stopping tolerance tbtol . While (cid:0) relative change in eigenvalues of ˇ A ≥ tbtol (cid:1) a. R Λ R − = ˇ A, ˇˇ b = ˇ b T R − T , ˇˇ c = ˇ cR, ˇˇ N = R T ˇ N R − T .b. Compute vec ( V ) = ( − Λ ⊗ I n − I r ⊗ A ) − (cid:16) ˇˇ b T ⊗ b (cid:17) ,vec ( W ) = (cid:0) − Λ ⊗ I n − I r ⊗ A T (cid:1) − (cid:0) ˇˇ c T ⊗ c T (cid:1) . c. For j = 2 , . . . , M , solve vec ( V j ) = ( − Λ ⊗ I n − I r ⊗ A ) − (cid:16) ˇˇ N T ⊗ N (cid:17) vec ( V j − ) ,vec ( W j ) = (cid:0) − Λ ⊗ I n − I r ⊗ A T (cid:1) − (cid:16) ˇˇ N ⊗ N T (cid:17) vec ( W j − ) . d. V = M P j =1 V j , W = M P j =1 W j .e. V r = orth ( V ) , W r = orth ( W ).f. ˇ A = ( W Tr V r ) − W Tr A V r , ˇ N = (cid:0) W Tr V r (cid:1) − W Tr N V r , ˇ b = (cid:0) W Tr V r (cid:1) − W Tr b, ˇ c = c V r . A r = ˇ A, N r = ˇ N , b r = ˇ b, c r = ˇ c .revisit this next. Generally, accuracy is the metric that tells about the correct-ness in the output of an inexact algorithm. Due to unavailability of the exactoutput, it is not possible to determine accuracy [25, 30]. A more easier metricis stability. Backward stability is one such notation, which says “A backwardstable algorithm gives exactly the right output to nearly the right input” [25].In our context, theoretically we obtain two reduced systems. One by applyingan inexact MOR algorithm (with iterative linear solves) on the original fullmodel, and other by applying the same MOR algorithm but exactly (withdirect linear solves) on a perturbed full model (the perturbation is introducedin the original full model as part of stability analysis, and is an unknownquantity). If these two reduced systems are equal ( first condition ), with thedifference between the original full model and the perturbed full model equal Rajendra Choudhary and Kapil Ahuja to the order of perturbation ( second condition ), then the MOR algorithm un-der consideration is called backward stable. The two theorems summarizingthis stability analysis for BIRKA are listed below.
Theorem 3. [30] If the inexact linear solves in BIRKA (lines 3b. and 3c.of Algorithm 1) are solved using a Petrov-Galerkin framework, then BIRKAsatisfies the first condition of backward stability with respect to these solves.
Theorem 4. [30] Let b Q = − (cid:20) A A (cid:21) ⊗ (cid:20) I n I n (cid:21) − (cid:20) I n I n (cid:21) ⊗ (cid:20) A A (cid:21) − (cid:20) N N (cid:21) ⊗ (cid:20) N N (cid:21)! , where I n is an identity matrix of size n × n and ⊗ denotes the standard Kronecker product. Also, let bb F = (cid:16) I n ⊗ b F + b F ⊗ I n (cid:17) with b F = (cid:20) F (cid:21) , where F is the perturbation introduced in A matrix of theinput dynamical system and I n is an identity matrix of size n × n . If b Q is invertible, (cid:13)(cid:13)(cid:13) b Q − (cid:13)(cid:13)(cid:13) < , and (cid:13)(cid:13)(cid:13)(cid:13) bb F (cid:13)(cid:13)(cid:13)(cid:13) < , then BIRKA satisfies the secondcondition of backward stability with respect to the inexact linear solves. Next, we analyze the backward stability of TBIRKA. Here, the first con-dition is satisfied in a way similar to that of BIRKA except that some extraorthogonality conditions are imposed on the linear solver (discussed below).
Theorem 5.
Let the inexact linear solves in TBIRKA (lines 3b. and 3c. ofAlgorithm 2) are solved using a Petrov-Galerkin framework with M X j =1 W j ⊥ R b j and M X j =1 V j ⊥ R c j , (7) where R b and R b j are the residuals in the first equations of 3b. and 3c. ofAlgorithm 2, respectively; R c and R c j are the residuals in the second equa-tions of 3b. and 3c. of Algorithm 2, respectively; and j = 2 , . . . , M. Then,TBIRKA satisfies the first condition of backward stability with respect to thesesolves.Proof.
Follows the same pattern as the proof for Theorem 2.1 in [30].From the above theorem, we infer that the underlying iterative solvershould firstly be based upon a Petrov-Galerkin framework. Since BiConju-gate Gradient (i.e., BiCG) is one such algorithm [23, 24], we propose its use in
NEXACT LINEAR SOLVES IN MODEL REDUCTION OF BILINEARDYNAMICAL SYSTEMS
Secondly, this particular solvershould also satisfy the extra orthogonalities (7).These orthogonalities have a form similar to the orthogonalities requiredwhile reducing second order linear dynamical systems ((23) and (24) in [28]),and can be easily satisfied by using a recycling variant of the underlying itera-tive solver. In [28], the ideal iterative solver to be used is Conjugate Gradient(i.e., CG) [23, 24] (due to the use of Galerkin projection). Hence, to sat-isfy the similar orthogonalities there, without any extra cost, the authors useRecycling Conjugate Gradient (i.e., RCG) [33]. Since here BiCG is the idealiterative solver (as discussed above), we propose the use of Recycling BiConju-gate Gradient (i.e., RBiCG) [32, 31], which would ensure that orthogonalitiesgiven by (7) are satisfied without any extra cost.To satisfy the second condition of backward stability of TBIRKA, we needto show that (cid:13)(cid:13)(cid:13) ζ M − e ζ M (cid:13)(cid:13)(cid:13) H = O ( k F k ) , (8)where ζ M is given by (6) or ζ M = { H ( s ) , H ( s , s ) , . . . , H M ( s , . . . , s M ) } , (9a)with H k ( s , . . . , s k ) for k = 1 , . . . , M given by (4) or H k ( s , . . . , s k ) = c ( s k I − A ) − N ( s k − I − A ) − . . . N ( s I − A ) − b, (9b) e ζ M = n e H ( s ) , e H ( s , s ) , . . . , e H M ( s , . . . , s M ) o , (10a)with for k = 1 , . . . , M e H k ( s , . . . , s k ) = c ( s k I − ( A + F )) − N ( s k − I − ( A + F )) − . . . N ( s I − ( A + F )) − b, (10b)and assuming perturbation in A matrix of the input dynamical system (as forBIRKA stability; see Theorem 4 earlier).One way to satisfy (8) is to use the definition of the H − norm of ζ M − e ζ M ,0 Rajendra Choudhary and Kapil Ahuja i.e., from Lemma 5.1 of [22] (cid:13)(cid:13)(cid:13) ζ M − e ζ M (cid:13)(cid:13)(cid:13) H = (cid:0)(cid:2) c − c (cid:3) ⊗ (cid:2) c − c (cid:3)(cid:1) M X j =0 "(cid:18) − (cid:20) A A + F (cid:21) ⊗ (cid:20) I n I n (cid:21) − (cid:20) I n I n (cid:21) ⊗ (cid:20) A A + F (cid:21)(cid:19) − (cid:20) N N (cid:21) ⊗ (cid:20) N N (cid:21) j (cid:18) − (cid:20) A A + F (cid:21) ⊗ (cid:20) I n I n (cid:21) − (cid:20) I n I n (cid:21) ⊗ (cid:20) A A + F (cid:21)(cid:19) − (cid:18)(cid:20) bb (cid:21) ⊗ (cid:20) bb (cid:21)(cid:19) . (11)This approach is followed in satisfying the second condition of backward sta-bility of BIRKA, but for TBIRKA it turns out to be more challenging. Thereason for this is that the definition of the H − norm of ζ − e ζ used in BIRKAis different from (11) ‡ , i.e., from Corollary 4.1 of [15] or Theorem 4.5 of [21] (cid:13)(cid:13)(cid:13) ζ − e ζ (cid:13)(cid:13)(cid:13) H = (cid:0)(cid:2) c − c (cid:3) ⊗ (cid:2) c − c (cid:3)(cid:1) (cid:18) − (cid:20) A A + F (cid:21) ⊗ (cid:20) I n I n (cid:21) − (cid:20) I n I n (cid:21) ⊗ (cid:20) A A + F (cid:21) − (cid:20) N N (cid:21) ⊗ (cid:20) N N (cid:21)(cid:19) − (cid:18)(cid:20) bb (cid:21) ⊗ (cid:20) bb (cid:21)(cid:19) . Another way to satisfy (8) in case of TBIRKA, which turns out to be moreeasier, is to show that (cid:13)(cid:13)(cid:13) H ( s ) − e H ( s ) (cid:13)(cid:13)(cid:13) H ∝ O ( k F k ) , (cid:13)(cid:13)(cid:13) H ( s , s ) − e H ( s , s ) (cid:13)(cid:13)(cid:13) H ∝ O ( k F k ) , ... (cid:13)(cid:13)(cid:13) H M ( s , . . . , s M ) − e H M ( s , . . . , s M ) (cid:13)(cid:13)(cid:13) H ∝ O ( k F k ) , (12)where H k ( s , . . . , s k ) for k = 1 , . . . , M is given by (4) and (9b), and e H k ( s , . . . , s k ) for k = 1 , . . . , M is given by (10b). This way was notpossible in BIRKA because there M → ∞ (see (2)-(4)). We use induction to prove (12). To prove this condition, we first abstractout the term containing the perturbation F from the normed difference be- ‡ Recall, in BIRKA we work with ζ rather than ζ M . NEXACT LINEAR SOLVES IN MODEL REDUCTION OF BILINEARDYNAMICAL SYSTEMS H M ( s , . . . , s M )) and thetransfer function of the perturbed system (cid:16) e H M ( s , . . . , s M ) (cid:17) (in Lemma 6).Next, we use induction to prove that this abstracted term is O ( k F k ) (inLemma 7). Finally, we use the result of these two lemmas to prove (12) (inTheorem 8). Note, that in all our subsequent derivations, we assume that allinverses used exist. This is an acceptable assumption because the inverse ofmatrices arising here are of the form as in [27] and [30] (the papers that discussstability of IRKA and BIRKA, respectively). Lemma 6. If H M ( s , . . . , s M ) = c ( s M I n − A ) − N ( s M − I n − A ) − . . . N ( s I n − A ) − b and e H M ( s , . . . , s M ) = c ( s M I n − ( A + F )) − N ( s M − I n − ( A + F )) − . . . N ( s I n − ( A + F )) − b, then the H − norm of their difference (cid:13)(cid:13)(cid:13) H M ( s , . . . , s M ) − e H M ( s , . . . , s M ) (cid:13)(cid:13)(cid:13) H ≤ (cid:13)(cid:13) c K − ( s M ) (cid:13)(cid:13) H (cid:13)(cid:13) K − ( s M − ) (cid:13)(cid:13) H . . . (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H k U ( s , . . . , s M ) k H ∞ (cid:13)(cid:13) K − ( s ) b (cid:13)(cid:13) H ∞ where, K ( s i ) = ( s i I n − A ) for i = 1 , . . . , M and U ( s , . . . , s M ) = K ( s ) . . . K ( s M − ) N K − ( s M − ) . . . N K − ( s ) N − (cid:0) I n − F K − ( s M ) (cid:1) − N K − ( s M − ) (cid:0) I n − F K − ( s M − ) (cid:1) − . . .N K − ( s ) (cid:0) I n − F K − ( s ) (cid:1) − N (cid:0) I n − K − ( s ) F (cid:1) − ! . (13) Proof.
Using the definition of H − norm (1), we get (cid:13)(cid:13)(cid:13) H M ( s , . . . , s M ) − e H M ( s , . . . , s M ) (cid:13)(cid:13)(cid:13) H = (cid:18) π (cid:19) M Lim m →∞ Z m − m . . . Z m − m (cid:13)(cid:13) c K − ( iω M ) N K − ( iω M − ) . . . N K − ( iω ) b − c ( K ( iω M ) − F ) − N ( K ( iω M − ) − F ) − . . . N ( K ( iω ) − F ) − N ( K ( iω ) − F ) − b (cid:13)(cid:13)(cid:13) F dω . . . dω M Rajendra Choudhary and Kapil Ahuja = (cid:18) π (cid:19) M Lim m →∞ Z m − m . . . Z m − m (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) c K − ( iω M ) N K − ( iω M − ) . . . N K − ( iω ) N − (cid:0) I n − F K − ( iω M ) (cid:1) − N K − ( iω M − ) (cid:0) I n − F K − ( iω M − ) (cid:1) − . . .N K − ( iω ) (cid:0) I n − F K − ( iω ) (cid:1) − N (cid:0) I n − K − ( iω ) F (cid:1) − ! K − ( iω ) b (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) F dω . . . dω M = (cid:18) π (cid:19) M Lim m →∞ Z m − m . . . Z m − m (cid:13)(cid:13) c K − ( iω M ) K − ( iω M − ) . . . K − ( iω ) K ( iω ) . . . K ( iω M − ) N K − ( iω M − ) . . . N K − ( iω ) N − (cid:0) I n − F K − ( iω M ) (cid:1) − N K − ( iω M − ) (cid:0) I n − F K − ( iω M − ) (cid:1) − . . .N K − ( iω ) (cid:0) I n − F K − ( iω ) (cid:1) − N (cid:0) I n − K − ( iω ) F (cid:1) − ! K − ( iω ) b (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) F dω . . . dω M . Using U ( s , . . . , s M ) given by (13), k XY Z k F ≤ k X k F k Y Z k F , k Y Z k F ≤k Y k F k Z k , and comparison integral inequality § [34] for any matrices X, Y and Z , in the above equation, we have (cid:13)(cid:13)(cid:13) H M ( s , . . . , s M ) − e H M ( s , . . . , s M ) (cid:13)(cid:13)(cid:13) H ≤ (cid:18) π (cid:19) M Lim m →∞ Z m − m . . . Z m − m (cid:13)(cid:13) c K − ( iω M ) (cid:13)(cid:13) F (cid:13)(cid:13) K − ( iω M − ) (cid:13)(cid:13) F . . . (cid:13)(cid:13) K − ( iω ) (cid:13)(cid:13) F k U ( iω , . . . , iω M ) k (cid:13)(cid:13) K − ( iω ) b (cid:13)(cid:13) dω . . . dω M . (14)From the mean value theorem of integration [34] we know Z m − m Z m − m f ( iω ) g ( iω , iω ) h ( iω ) dω dω = Z m − m f ( iω ) (cid:18)Z m − m g ( iω , iω ) h ( iω ) dω (cid:19) dω § This inequality says if f ( x ) and g ( x ) are integrable over [ a, b ] and f ( x ) ≤ g ( x ),then R ba f ( x ) dx ≤ R ba g ( x ) dx . Note that although we have improper integrals here, thisinequality still holds because of the earlier assumption that such integrals give a finitevalue. NEXACT LINEAR SOLVES IN MODEL REDUCTION OF BILINEARDYNAMICAL SYSTEMS ≤ Z m − m f ( iω ) (cid:18) max c ∈ R ( g ( ic, iω )) Z m − m h ( iω ) dω (cid:19) dω ≤ (cid:18)Z mm f ( iω ) max c ∈ R ( g ( ic, iω )) dω (cid:19) Z m − m h ( iω ) dω ≤ max c,d ∈ R ( g ( ic, id )) Z m − m f ( iω ) dω Z m − m h ( iω ) dω . Using this property in (14) we get ¶ (cid:13)(cid:13)(cid:13) H M ( s , . . . , s M ) − e H M ( s , . . . , s M ) (cid:13)(cid:13)(cid:13) H ≤ (cid:18) π (cid:19) M Lim m →∞ Z m − m . . . Z m − m (cid:13)(cid:13) c K − ( iω M ) (cid:13)(cid:13) F (cid:13)(cid:13) K − ( iω M − ) (cid:13)(cid:13) F . . . (cid:13)(cid:13) K − ( iω ) (cid:13)(cid:13) F dω . . . dω M max ω , ..., ω M ∈ R k U ( iω , . . . , iω M ) k max ω ∈ R (cid:13)(cid:13) K − ( iω ) b (cid:13)(cid:13) ≤ (cid:13)(cid:13) c K − ( s M ) (cid:13)(cid:13) H (cid:13)(cid:13) K − ( s M − ) (cid:13)(cid:13) H . . . (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H k U ( s , . . . , s M ) k H ∞ (cid:13)(cid:13) K − ( s ) b (cid:13)(cid:13) H ∞ . Lemma 7. If U M = U ( s , . . . , s M ) defined in (13) , k F k < , and (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ < , where K ( s ) = ( sI n − A ) . Then we have k U M k H ∞ ∝ O ( k F k ) . Proof.
We prove this by mathematical induction.
Base Case (i)
First subsystem
This is the linear system case, which has been already proved in [27] (seeTheorem 4.3 of [27]).(ii)
Second subsystem
Substituting M = 2 in (13), we get U = K ( s ) (cid:16) N − (cid:0) I n − F K − ( s ) (cid:1) − N (cid:0) I n − K − ( s ) F (cid:1) − (cid:17) . ¶ As mentioned in Footnote § , the improper integral does not affect application of thismean value theorem because all such integrals are assumed to give a finite value. Rajendra Choudhary and Kapil Ahuja If (cid:13)(cid:13) F K − ( s ) (cid:13)(cid:13) H ∞ < (cid:13)(cid:13) K − ( s ) F (cid:13)(cid:13) H ∞ <
1, then by the Neumannseries, we get k U = K ( s ) (cid:18) N − (cid:16) I n + F K − ( s ) + (cid:0) F K − ( s ) (cid:1) + · · · (cid:17) N (cid:16) I n + K − ( s ) F + (cid:0) K − ( s ) F (cid:1) + · · · (cid:17) (cid:19) = K ( s ) (cid:18) N − N − N K − ( s ) F (cid:0) I n + K − ( s ) F + · · · (cid:1) − F K − ( s ) (cid:0) I n + F K − ( s ) + · · · (cid:1) N (cid:16) I n + K − ( s ) F + (cid:0) K − ( s ) F (cid:1) + · · · (cid:17) (cid:19) = K ( s ) (cid:18) − N K − ( s ) F (cid:0) I n − K − ( s ) F (cid:1) − − F K − ( s ) (cid:0) I n − F K − ( s ) (cid:1) − N (cid:0) I n − K − ( s ) F (cid:1) − (cid:19) = K ( s ) (cid:16) − N K − ( s ) F − F K − ( s ) (cid:0) I n − F K − ( s ) (cid:1) − N (cid:17)(cid:0) I n − K − ( s ) F (cid:1) − . Taking H ∞ − norm on both sides, we get k U k H ∞ = (cid:13)(cid:13)(cid:13) K ( s ) (cid:16) − N K − ( s ) F − F K − ( s ) (cid:0) I n − F K − ( s ) (cid:1) − N (cid:17)(cid:0) I n − K − ( s ) F (cid:1) − (cid:13)(cid:13)(cid:13) H ∞ = max ω ,ω ∈ R (cid:13)(cid:13)(cid:13) K ( iω ) (cid:16) − N K − ( iω ) F − F K − ( iω ) (cid:0) I n − F K − ( iω ) (cid:1) − N (cid:17) (cid:0) I n − K − ( iω ) F (cid:1) − (cid:13)(cid:13)(cid:13) . Using k XY k ≤ k X k k Y k and k X + Y k ≤ k X k + k Y k , for any two k From [35, page 527], we know ( I − A ) − = ∞ P k =0 A k when k A k < (cid:13)(cid:13) F K − ( s ) (cid:13)(cid:13) H ∞ < max ω ∈ R (cid:13)(cid:13) F K − ( iω ) (cid:13)(cid:13) <
1, andhence, the applicable matrix norm is 2 − norm. Similarly for the second inequality. NEXACT LINEAR SOLVES IN MODEL REDUCTION OF BILINEARDYNAMICAL SYSTEMS X and Y , in the above equation, we get k U k H ∞ ≤ max ω ,ω ∈ R k K ( iω ) k (cid:18) k N k (cid:13)(cid:13) K − ( iω ) (cid:13)(cid:13) k F k + k F k (cid:13)(cid:13) K − ( iω ) (cid:13)(cid:13) (cid:13)(cid:13)(cid:13)(cid:0) I n − F K − ( iω ) (cid:1) − (cid:13)(cid:13)(cid:13) k N k (cid:19)(cid:13)(cid:13)(cid:13)(cid:0) I n − K − ( iω ) F (cid:1) − (cid:13)(cid:13)(cid:13) ! ≤ k K ( s ) k H ∞ k N k k F k (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ + (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ max ω ∈ R (cid:13)(cid:13)(cid:13)(cid:0) I n − F K − ( iω ) (cid:1) − (cid:13)(cid:13)(cid:13) ! max ω ∈ R (cid:13)(cid:13)(cid:13)(cid:0) I n − K − ( iω ) F (cid:1) − (cid:13)(cid:13)(cid:13) . (15)Technically by definition of the H ∞ − norm and how K ( s ) is defined inour hypotheses, k K ( s ) k H ∞ = k K ( s ) k H ∞ = k K ( s ) k H ∞ , however, forsake of exposition, we keep them separate. Similarly for the H ∞ − normof inverses of K ( s ) and K ( s ).To abstract k F k out from the above inequality, let us lookat max ω ∈ R (cid:13)(cid:13)(cid:13)(cid:0) I n − F K − ( iω ) (cid:1) − (cid:13)(cid:13)(cid:13) separately. Recall, while apply-ing Neumann series we assumed that (cid:13)(cid:13) F K − ( s ) (cid:13)(cid:13) H ∞ < max ω ∈ R (cid:13)(cid:13) F K − ( iω ) (cid:13)(cid:13) <
1. Since the maximum of such a norm is lessthan one, we have for all ω ∈ R , (cid:13)(cid:13) F K − ( iω ) (cid:13)(cid:13) <
1. Using this alongwith Lemma 2.3.3 from [36] ∗∗ in the above expression, we get max ω ∈ R (cid:13)(cid:13)(cid:13)(cid:0) I n − F K − ( iω ) (cid:1) − (cid:13)(cid:13)(cid:13) ≤ max ω ∈ R − k F K − ( iω ) k ≤ − max ω ∈ R k F K − ( iω ) k ≤ − k F K − ( s ) k H ∞ . (16) ∗∗ If F ∈ R n × n and k F k p <
1, then I − F is nonsingular and ( I − F ) − = ∞ P k =0 F k with (cid:13)(cid:13)(cid:13) ( I − F ) − (cid:13)(cid:13)(cid:13) p ≤ − k F k p . Rajendra Choudhary and Kapil Ahuja
If we assume k F k < (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ < (cid:13)(cid:13) F K − ( s ) (cid:13)(cid:13) H ∞ = max ω ∈ R (cid:13)(cid:13) F K − ( iω ) (cid:13)(cid:13) ≤ k F k max ω ∈ R (cid:13)(cid:13) K − ( iω ) (cid:13)(cid:13) ≤ k F k (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ ≤ , as assumed for applying Neumann series earlier as well as Lemma 2.3.3from [36] above. Thus, no extra assumptions beyond those in hypothesesare needed. Further, we also get1 − k F k (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ ≤ − (cid:13)(cid:13) F K − ( s ) (cid:13)(cid:13) H ∞ or11 − k F K − ( s ) k H ∞ ≤ − k F k k K − ( s ) k H ∞ . (17)Similarly, we can bound the last term of (15) as follows: max ω ∈ R (cid:13)(cid:13)(cid:13)(cid:0) I n − K − ( iω ) F (cid:1) − (cid:13)(cid:13)(cid:13) ≤ − k K − ( s ) F k H ∞ and (18)11 − k K − ( s ) F k H ∞ ≤ − k K − ( s ) k H ∞ k F k . (19)Substituting (16)-(17) and (18)-(19) in (15), we get k U k H ∞ ≤ k K ( s ) k H ∞ k N k k F k " (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ + (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ − k F k k K − ( s ) k H ∞ − k K − ( s ) k H ∞ k F k ! . From the above inequality it is clear that if k F k (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ < (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ k F k <
1, which are true from our hypotheses, then k U k H ∞ = O ( k F k ) . (iii) Third subsystem
Using M = 3 in (13), we get U = K ( s ) K ( s ) (cid:16) N K − ( s ) N − (cid:0) I n − F K − ( s ) (cid:1) − N K − ( s ) (cid:0) I n − F K − ( s ) (cid:1) − N (cid:0) I n − K − ( s ) F (cid:1) − (cid:17) . NEXACT LINEAR SOLVES IN MODEL REDUCTION OF BILINEARDYNAMICAL SYSTEMS k U k H ∞ ≤ k K ( s ) k H ∞ k K ( s ) k H ∞ k N k k F k (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ " (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ + (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ − k F k k K − ( s ) k H ∞ + (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ (cid:0) − k F k k K − ( s ) k H ∞ (cid:1) (cid:0) − k F k k K − ( s ) k H ∞ (cid:1) − k K − ( s ) k H ∞ k F k ! . Again, from the above inequality it is clear that if k F k (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ , k F k (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ and (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ k F k are all less than 1, whichare true from our hypotheses, then k U k H ∞ = O ( k F k ) . Induction Hypothesis
From (13), we know U M = K ( s ) . . . K ( s M − ) N K − ( s M − ) . . . N K − ( s ) N − (cid:0) I n − F K − ( s M ) (cid:1) − N K − ( s M − ) (cid:0) I n − F K − ( s M − ) (cid:1) − . . .N K − ( s ) (cid:0) I n − F K − ( s ) (cid:1) − N (cid:0) I n − K − ( s ) F (cid:1) − ! . (20)Let k U M k H ∞ ≤ k K ( s ) k H ∞ . . . k K ( s M − ) k H ∞ k N k M − k F k (cid:13)(cid:13) K − ( s M − ) (cid:13)(cid:13) H ∞ (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ " (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ + (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ − k F k k K − ( s ) k H ∞ ...+ (cid:13)(cid:13) K − ( s M ) (cid:13)(cid:13) H ∞ (cid:0) − k F k k K − ( s M ) k H ∞ (cid:1) . . . (cid:0) − k F k k K − ( s ) k H ∞ (cid:1) − k K − ( s ) k H ∞ k F k ! or k U M k H ∞ = O ( k F k ) . Rajendra Choudhary and Kapil Ahuja
Induction Step
Again, from (13), we know U M +1 = K ( s ) . . . K ( s M ) N K − ( s M ) . . . N K − ( s ) N − (cid:0) I n − F K − ( s M +1 ) (cid:1) − N K − ( s M ) (cid:0) I n − F K − ( s M ) (cid:1) − . . .N K − ( s ) (cid:0) I n − F K − ( s ) (cid:1) − N (cid:0) I n − K − ( s ) F (cid:1) − ! . We first write U M +1 in terms of U M . Thus, in the above equation usingNeumann series, we get U M +1 = K ( s ) . . . K ( s M ) N K − ( s M ) . . . N K − ( s ) N − (cid:16) I n + F K − ( s M +1 ) + (cid:0) F K − ( s M +1 ) (cid:1) + · · · (cid:17) N K − ( s M ) (cid:0) I n − F K − ( s M ) (cid:1) − . . . N K − ( s ) (cid:0) I n − F K − ( s ) (cid:1) − N (cid:0) I n − K − ( s ) F (cid:1) − ! = K ( s ) . . . K ( s M ) N K − ( s M ) . . . N K − ( s ) N − N K − ( s M ) (cid:0) I n − F K − ( s M ) (cid:1) − . . . N K − ( s ) (cid:0) I n − F K − ( s ) (cid:1) − N (cid:0) I n − K − ( s ) F (cid:1) − − F K − ( s M +1 ) (cid:0) I n − F K − ( s M +1 ) (cid:1) − N K − ( s M ) (cid:0) I n − F K − ( s M ) (cid:1) − . . . N K − ( s ) (cid:0) I n − F K − ( s ) (cid:1) − N (cid:0) I n − K − ( s ) F (cid:1) − ! . In the above equation, taking N K − ( s M ) common from the first two terms NEXACT LINEAR SOLVES IN MODEL REDUCTION OF BILINEARDYNAMICAL SYSTEMS K ( s ) . . . K ( s M ) N K − ( s M ) (cid:18) N K − ( s M − ) . . . N K − ( s ) N − (cid:0) I n − F K − ( s M ) (cid:1) − N K − ( s M − ) (cid:0) I n − F K − ( s M − ) (cid:1) − . . .N K − ( s ) (cid:0) I n − F K − ( s ) (cid:1) − N (cid:0) I n − K − ( s ) F (cid:1) − (cid:19) − F K − ( s M +1 ) (cid:0) I n − F K − ( s M +1 ) (cid:1) − N K − ( s M ) (cid:0) I n − F K − ( s M ) (cid:1) − . . . N K − ( s ) (cid:0) I n − F K − ( s ) (cid:1) − N (cid:0) I n − K − ( s ) F (cid:1) − ! . (21)Now we look at expression of U M . Multiplying K − ( s M − ) . . . K − ( s ) onboth the sides of (20) from left, we get K − ( s M − ) . . . K − ( s ) U M = N K − ( s M − ) . . . N K − ( s ) N − (cid:0) I n − F K − ( s M ) (cid:1) − N K − ( s M − ) (cid:0) I n − F K − ( s M − ) (cid:1) − . . .N K − ( s ) (cid:0) I n − F K − ( s ) (cid:1) − N (cid:0) I n − K − ( s ) F (cid:1) − ! . (22)Substituting (22) in (21), we get U M +1 = K ( s ) . . . K ( s M ) N K − ( s M ) (cid:18) K − ( s M − ) . . . K − ( s ) U M (cid:19) − F K − ( s M +1 ) (cid:0) I n − F K − ( s M +1 ) (cid:1) − N K − ( s M ) (cid:0) I n − F K − ( s M ) (cid:1) − . . . N K − ( s ) (cid:0) I n − F K − ( s ) (cid:1) − N (cid:0) I n − K − ( s ) F (cid:1) − ! . Rajendra Choudhary and Kapil Ahuja
Taking H ∞ − norm on both sides, we get k U M +1 k H ∞ = (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) K ( s ) . . . K ( s M ) N K − ( s M ) K − ( s M − ) . . . K − ( s ) U M − F K − ( s M +1 ) (cid:0) I n − F K − ( s M +1 ) (cid:1) − N K − ( s M ) (cid:0) I n − F K − ( s M ) (cid:1) − . . . N K − ( s ) (cid:0) I n − F K − ( s ) (cid:1) − N (cid:0) I n − K − ( s ) F (cid:1) − !(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) H ∞ . As earlier, using the norm inequality properties in the above equation, we get k U M +1 k H ∞ ≤ max ω , ..., ω M +1 ∈ R " k K ( iω ) k . . . k K ( iω M ) k (cid:18) k N k (cid:13)(cid:13) K − ( iω M ) (cid:13)(cid:13) . . . (cid:13)(cid:13) K − ( iω ) (cid:13)(cid:13) k U ( iω , . . . , iω M ) k + k F k (cid:13)(cid:13) K − ( iω M +1 ) (cid:13)(cid:13) (cid:13)(cid:13)(cid:13)(cid:0) I n − F K − ( iω M +1 ) (cid:1) − (cid:13)(cid:13)(cid:13) k N k (cid:13)(cid:13) K − ( iω M ) (cid:13)(cid:13) (cid:13)(cid:13)(cid:13)(cid:0) I n − F K − ( iω M ) (cid:1) − (cid:13)(cid:13)(cid:13) . . . k N k (cid:13)(cid:13) K − ( iω ) (cid:13)(cid:13) (cid:13)(cid:13)(cid:13)(cid:0) I n − F K − ( iω ) (cid:1) − (cid:13)(cid:13)(cid:13) k N k (cid:13)(cid:13)(cid:13)(cid:0) I n − K − ( iω ) F (cid:1) − (cid:13)(cid:13)(cid:13) (cid:19) . Similar to (16) and (17), here also, using Lemma 2.3.3 from [36] we get k U M +1 k H ∞ ≤ k K ( s ) k H ∞ . . . k K ( s M ) k H ∞ k N k (cid:13)(cid:13) K − ( s M ) (cid:13)(cid:13) H ∞ . . . (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ " (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ k U M k H ∞ + k N k M − (cid:13)(cid:13) K − ( s M +1 ) (cid:13)(cid:13) H ∞ (cid:0) − k F k k K − ( s M +1 ) k H ∞ (cid:1) . . . (cid:0) − k F k k K − ( s ) k H ∞ (cid:1) ·k F k − k K − ( s ) k H ∞ k F k . From induction hypothesis we know k U M k H ∞ ∝ O ( k F k ). Using this we get k U M +1 k H ∞ ∝ O ( k F k ) . NEXACT LINEAR SOLVES IN MODEL REDUCTION OF BILINEARDYNAMICAL SYSTEMS Theorem 8. If H M ( s , . . . , s M ) = c ( s M I n − A ) − N ( s M − I n − A ) − . . . N ( s I n − A ) − b, e H M ( s , . . . , s M ) = c ( s M I n − ( A + F )) − N ( s M − I n − ( A + F )) − . . . N ( s I n − ( A + F )) − b, k F k < , and (cid:13)(cid:13) K − ( s ) (cid:13)(cid:13) H ∞ < , where K ( s ) = ( sI n − A ) , then (cid:13)(cid:13)(cid:13) H M ( s , . . . , s M ) − e H M ( s , . . . , s M ) (cid:13)(cid:13)(cid:13) H = O (cid:16) k F k (cid:17) or TBIRKA satisfies the second condition of backward stability with respect toinexact linear solves.
In this paper, we apply iterative linear solvers during model order reduction(MOR) of bilinear dynamical systems. Since such solvers are inexact, the sta-bility of the underlying MOR algorithm, with respect to these approximationerrors, is important. Here, we extend the earlier stability analysis done forBIRKA in [30], which is a standard algorithm for obtaining a reduced bilineardynamical system, to TBIRKA, a cheaper variant of BIRKA.Proving that an algorithm is stable, typically requires satisfying two con-ditions. In TBIRKA, fulfilling the first condition for stability leads to con-straints on the iterative linear solver, which are similar to those obtained dur-ing BIRKA’s stability analysis. The second condition for a stable TBIRKAis satisfied using an approach different than the one used in BIRKA, and ismore intuitive.Our first future direction is to extend our analysis from SISO (Single In-put Single Output) to MIMO (Multiple Input Multiple Output) systems. Thestability analysis as done for BIRKA earlier and TBIRKA here, all give ussufficiency conditions for a stable underlying MOR algorithm. Hence, second,we plan to derive the necessary conditions for the same. In recent years, therehave been a lot of efforts in performing data-driven MOR algorithm (speciallyusing Leowner framework [19]). Similar linear systems and challenges associ-ated with the scaling arise here as well. Our third future direction is to applythis stability analysis to such classes of algorithms as well.
References [1] H. Kaper, H. Engler.
Mathematics And Climate . SIAM, Philadelphia,PA, USA, 2013.2
Rajendra Choudhary and Kapil Ahuja [2] J. D. Meiss.
Differential Dynamical Systems, Revised Edition . SIAM,Philadelphia, PA, USA, 2017.[3] J. T. Stuart. Taylor-Vortex flow: A dynamical system.
SIAM Review ,28(3):315–342, 1986.[4] P. D’Alessandro, A. Isidori, and A. Ruberti. Realization and structuretheory of bilinear dynamical systems.
SIAM J. Control Opt. , 12(3):517–535, 1974.[5] R. J. Wilson.
Nonlinear System Theory: The Volterra/ Wiener Approach .Johns Hopkins Series in Information Sciences and Systems, Johns HopkinsUniversity Press, Baltimore, 1981.[6] F. Carravetta. Global exact quadratization of continuous–time nonlinearcontrol systems.
SIAM J. Control Opt. , 53(1):235–261, 2015.[7] A. C. Antoulas.
Approximation of Large-Scale Dynamical Systems . Ad-vances in Design and Control, SIAM, Philadelphia, PA, USA, 2005.[8] P. Benner, D. C. Sorensen, and V. Mehrmann.
Dimension Reductionof Large-Scale Systems . Lecture Notes in Computational Science andEngineering, Springer, Berlin, Germany, 2005.[9] E. J. Grimme. Krylov Projection Methods for Model Reduction.
Ph.D.thesis , University of Illinois at Urbana-Champaign, Urbana, IL, USA,1997.[10] Z. Bai. Krylov subspace techniques for reduced-order modeling of large-scale dynamical systems.
Appl. Num. Math. , 43(1):9–44, 2002.[11] K. Willcox and J. Peraire. Balanced model reduction via the properorthogonal decomposition.
AIAA Journal , 40(11):2323–2330, 2002.[12] S. Gugercin, A. C. Antoulas, and C. Beattie. H model reductionfor large-scale linear dynamical systems. SIAM J. Matrix Anal. Appl. ,30(2):609–638, 2008.[13] A. Bunse-Gerstner, D. Kubali´nska, G. Vossen, and D. Wilczek. h -normoptimal model reduction for large scale discrete dynamical MIMO sys-tems. J. Comput. Appl. Math. , 233(5):1202–1216, 2010.[14] Z. Bai and D. Skoogh. A projection method for model reduction of bilineardynamical systems.
Linear Algebra Appl. , 415(2-3):406–425, 2006.
NEXACT LINEAR SOLVES IN MODEL REDUCTION OF BILINEARDYNAMICAL SYSTEMS H -model reduction ofbilinear control systems. SIAM J. Matrix Anal. Appl. , 33(3):859–885,2012.[16] T. Bonin, H. Faßbender, A. Soppa, and M. Zaeh. A fully adaptive ratio-nal global Arnoldi method for the model-order reduction of second-orderMIMO systems with proportional damping.
Math. Comp. Simulation. ,122(C):1–19, 2016.[17] M. I. Ahmad, U. Baur, and P. Benner. Implicit Volterra series interpo-lation for model reduction of bilinear systems.
J. Comput. Appl. Math. ,316:15–28, 2017.[18] P. Benner and T. Damm. Lyapunov equations, energy functionals, andmodel order reduction of bilinear and stochastic systems.
SIAM J. ControlOpt. , 49(2):686–711, 2011.[19] A. C. Antoulas, I. V. Gosea, and A. C. Ionita. Model reduction of bilinearsystems in the Loewner framework.
SIAM J. Control Opt. , 38(5):B889–B916, 2016.[20] K. L. Xu, Y. L. Jiang, and Z. X. Yang. H optimal model order reductionby two-sided technique on Grassmann manifold via the cross-Gramian ofbilinear systems. Inter. J. Control , 90(3):616–626, 2017.[21] G. M. Flagg. Interpolation Methods for the Model Reduction of BilinearSystems.
Ph.D. thesis , Virginia Polytechnic Institute and State Univer-sity, Blacksburg, VA, USA, 2012.[22] G. M. Flagg and S. Gugercin. Multipoint Volterra series interpolationand H optimal model reduction of bilinear systems. SIAM J. MatrixAnal. Appl. , 36(2):549–579, 2015.[23] A. Greenbaum.
Iterative Methods for Solving Linear Systems, Frontiersin Applied Mathematics . SIAM, Philadelphia, PA, USA, 1997.[24] Y. Saad.
Iterative Methods for Sparse Linear Systems . SIAM, Philadel-phia, PA, USA, 2003.[25] L. N. Trefethen and D. Bau.
Numerical Linear Algebra . SIAM, Philadel-phia, PA, USA, 1997.[26] J. W. Demmel.
Applied Numerical Linear Algebra . SIAM, Philadelphia,PA, USA, 1997.4
Rajendra Choudhary and Kapil Ahuja [27] C. Beattie, S. Gugercin, and S. Wyatt. Inexact solves in interpolatorymodel reduction.
Linear Algebra Appl. , 436(8):2916–2943, 2012.[28] N. P. Singh and K. Ahuja. Stability analysis of inexact solves in momentmatching based model reduction.
Preprint arXiv:1803.09283 , 2018.[29] D. Lu, Y. Su, and Z. Bai. Stability analysis of the two-level orthogonalArnoldi procedure.
SIAM J. Matrix Anal. Appl. , 37(1):195–214, 2016.[30] R. Choudhary and K. Ahuja. Stability analysis of bilinear iterative ratio-nal Krylov algorithm.
Linear Algebra Appl. , 538:56–88, 2018.[31] K. Ahuja. Recycling Krylov subspaces and preconditioners.
Ph.D. the-sis , Virginia Polytechnic Institute and State University, Blacksburg, VA,USA, 2011.[32] K. Ahuja, E. D. Sturler, S. Gugercin, and E. R. Chang. RecyclingBiCG with an application to model reduction.
SIAM J. Sci. Comput. ,34(4):A1925–A1949, 2012.[33] M. L. Parks, E. D. Sturler, G. Mackey, D. D. Johnson, and S. Maiti.Recycling Krylov subspaces for sequences of linear systems.
SIAM J. Sci.Comput. , 28(5):1651-1674, 2006.[34] A. Jeffrey.
Handbook of Mathematical Formulas and Integrals, Ed. 3 .Academic Press, 2003.[35] C. D. Meyer.
Matrix Analysis and Applied Linear Algebra . SIAM,Philadelphia, PA, USA, 2000.[36] G. H. Golub and C. F. Van Loan.
Matrix Computations, Vol. 3 . JohnsHopkins University Press, 2012.. JohnsHopkins University Press, 2012.