Inexact Linear Solves In Model Reduction of Bilinear Dynamical Systems

TL;DR

This paper analyzes the stability of the TBIRKA algorithm for bilinear model order reduction when using inexact linear solves, proposing new conditions and an intuitive methodology for ensuring stability in large-scale systems.

Contribution

It derives stability conditions for TBIRKA with inexact solves and introduces a novel, more intuitive approach for achieving stability, extending previous work on BIRKA.

Findings

Established stability conditions for TBIRKA with inexact linear solves

Proposed a new, more intuitive methodology for stability assurance

The techniques can be extended to other MOR methods like balanced truncation

Abstract

Bilinear dynamical systems are commonly used in science and engineering because they form a bridge between linear and non-linear systems. However, simulating them is still a challenge because of their large size. Hence, a lot of research is currently being done for reducing such bilinear dynamical systems (termed as bilinear model order reduction or bilinear MOR). Bilinear iterative rational Krylov algorithm (BIRKA) is a very popular, standard and mathematically sound algorithm for bilinear MOR, which is based upon interpolatory projection technique. An efficient variant of BIRKA, Truncated BIRKA (or TBIRKA) has also been recently proposed. Like for any MOR algorithm, these two algorithms also require solving multiple linear systems as part of the model reduction process. For reducing very large dynamical systems, which is now-a-days becoming a norm, scaling of such linear systems with respect to input dynamical system size is a bottleneck. For efficiency, these linear systems are often solved by an iterative solver, which introduces approximation errors. Hence, stability analysis of MOR algorithms with respect to inexact linear solves is important. In our past work, we have shown that under 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 novel contribution is the more intuitive methodology for achieving this. This approach exploits the fact that in TBIRKA a bilinear dynamical system can be represented by a finite set of functions, which was not possible in BIRKA (because infinite such functions were needed there). The stability analysis techniques that we propose here can be extended to many other methods for doing MOR of bilinear dynamical systems, e.g., using balanced truncation or the ADI methods.

Equations

$${\Vert H_k\Vert }_{H_2}^2={\left(\frac{1}{2\pi }\right)}^k{\int }_{-\infty }^{\infty }\dots {\int }_{-\infty }^{\infty }{\Vert H_k\left(i{\omega }_1,\dots ,i{\omega }_k\right)\Vert }_F^{2}d{\omega }_1\dots d{\omega }_k,$$

$${\Vert H_k\Vert }_{H_{\infty }}=\underset{{\omega }_1,\dots ,{\omega }_k\in \mathbb{R}}{\max }{\Vert H_k\left(i{\omega }_1,\dots ,i{\omega }_k\right)\Vert }_2.$$

$$P\otimes Q=\left[\begin{array}{ccc}{p}_{11}Q& \cdots & p_{1n}Q\\ ⋮& \ddots & ⋮\\ p_{m1}Q& \cdots & p_{mn}Q\end{array}\right],$$

$$vec(P)={\left(p_{11},\dots ,p_{m1},p_{12},\dots ,p_{m2},\dots \dots ,p_{1n},\dots ,p_{mn}\right)}^T.$$

$$\zeta:\left\{\begin{array}[]{ll}\dot{x}(t)=Ax(t)+Nx(t)u(t)+bu(t),\\ y(t)=cx(t),\end{array}\right.$$

$$\zeta =\underset{k\to \infty }{\mathrm{Lim}}{\zeta }^k,$$

$$H_k\left(s_1,\dots ,s_k\right)=c{\left(s_{k}I-A\right)}^{-1}N{\left(s_{k-1}I-A\right)}^{-1}\dots N{\left(s_{1}I-A\right)}^{-1}b.$$

$$\zeta_{r}:\left\{\begin{array}[]{ll}\dot{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),\end{array}\right.$$

$$\begin{array}{rl}{W}_r^T& \left(V_r{\dot{x}}_r(t)-A{V}_r{x}_r(t)-N{V}_r{x}_r(t)u(t)-bu(t)\right)=0,\\ y(t)& =c{V}_r{x}_r(t).\end{array}$$

$$\begin{array}{rl}{A}_r& ={\left(W_r^T{V}_r\right)}^{-1}{W}_r^{T}A{V}_r,N_r={\left(W_r^T{V}_r\right)}^{-1}{W}_r^{T}N{V}_r,\\ b_r& ={\left(W_r^T{V}_r\right)}^{-1}{W}_r^{T}b,\text{and}{c}_r=c{V}_r.\end{array}$$

$$\sum _{k=1}^{\infty }\sum _{l_1=1}^r\dots \sum _{l_k=1}^r{\varphi }_{l_1,l_2,\dots ,l_k}{H}_k\left(-{\lambda }_{l_1},-{\lambda }_{l_2},\dots ,-{\lambda }_{l_k}\right)=$$

$$\sum _{k=1}^{\infty }\sum _{l_1=1}^r\dots \sum _{l_k=1}^r{\varphi }_{l_1,l_2,\dots ,l_k}{H}_{r_k}\left(-{\lambda }_{l_1},-{\lambda }_{l_2},\dots ,-{\lambda }_{l_k}\right),\text{and}$$

$$\sum _{k=1}^{\infty }\sum _{l_1=1}^r\dots \sum _{l_k=1}^r{\varphi }_{l_1,l_2,\dots ,l_k}\left(\sum _{j=1}^k\frac{\partial }{\partial s_j}{H}_k\left(-{\lambda }_{l_1},-{\lambda }_{l_2},\dots ,-{\lambda }_{l_k}\right)\right)=$$

$$\sum _{k=1}^{\infty }\sum _{l_1=1}^r\dots \sum _{l_k=1}^r{\varphi }_{l_1,l_2,\dots ,l_k}\left(\sum _{j=1}^k\frac{\partial }{\partial s_j}{H}_{r_k}\left(-{\lambda }_{l_1},-{\lambda }_{l_2},\dots ,-{\lambda }_{l_k}\right)\right),$$

$${\zeta }^M=\left\{H_1\left(s_1\right),H_2\left(s_1,s_2\right),H_3\left(s_1,s_2,s_3\right),\dots ,H_M\left(s_1,\dots ,s_M\right)\right\},$$

$$\sum _{k=1}^M\sum _{l_1=1}^r\dots \sum _{l_k=1}^r{\varphi }_{l_1,l_2,\dots ,l_k}{H}_k\left(-{\lambda }_{l_1},-{\lambda }_{l_2},\dots ,-{\lambda }_{l_k}\right)=$$

$$\sum _{k=1}^M\sum _{l_1=1}^r\dots \sum _{l_k=1}^r{\varphi }_{l_1,l_2,\dots ,l_k}{H}_{r_k}\left(-{\lambda }_{l_1},-{\lambda }_{l_2},\dots ,-{\lambda }_{l_k}\right),\text{and}$$

$$\sum _{k=1}^M\sum _{l_1=1}^r\dots \sum _{l_k=1}^r{\varphi }_{l_1,l_2,\dots ,l_k}\left(\sum _{j=1}^k\frac{\partial }{\partial s_j}{H}_k\left(-{\lambda }_{l_1},-{\lambda }_{l_2},\dots ,-{\lambda }_{l_k}\right)\right)=$$

$$\sum _{k=1}^M\sum _{l_1=1}^r\dots \sum _{l_k=1}^r{\varphi }_{l_1,l_2,\dots ,l_k}\left(\sum _{j=1}^k\frac{\partial }{\partial s_j}{H}_{r_k}\left(-{\lambda }_{l_1},-{\lambda }_{l_2},\dots ,-{\lambda }_{l_k}\right)\right),$$

$$vec\left({\mathbf{V}}_1\right)={\left(-\Lambda \otimes I_n-I_r\otimes A\right)}^{-1}\left({\check{\stackrel{ˇ}{b}}}^T\otimes b\right),$$

$$vec\left({\mathbf{W}}_1\right)={\left(-\Lambda \otimes I_n-I_r\otimes A^T\right)}^{-1}\left({\check{\stackrel{ˇ}{c}}}^T\otimes c^T\right).$$

$$vec\left({\mathbf{V}}_j\right)={\left(-\Lambda \otimes I_n-I_r\otimes A\right)}^{-1}\left({\check{\stackrel{ˇ}{N}}}^T\otimes N\right)vec\left({\mathbf{V}}_{j-1}\right),$$

$$vec\left({\mathbf{W}}_j\right)={\left(-\Lambda \otimes I_n-I_r\otimes A^T\right)}^{-1}\left(\check{\stackrel{ˇ}{N}}\otimes N^T\right)vec\left({\mathbf{W}}_{j-1}\right).$$

$$\sum _{j=1}^M{W}_j\perp R_{b_j}\text{and}\sum _{j=1}^M{V}_j\perp R_{c_j},$$

$${\Vert {\zeta }^M-{\tilde{\zeta }}^M\Vert }_{H_2}=\mathcal{O}\left({\Vert F\Vert }_2\right),$$

$${\zeta }^M=\left\{H_1\left(s_1\right),H_2\left(s_1,s_2\right),\dots ,H_M\left(s_1,\dots ,s_M\right)\right\},$$

$$\begin{array}{rl}& H_k\left(s_1,\dots ,s_k\right)=c{\left(s_{k}I-A\right)}^{-1}\\ & N{\left(s_{k-1}I-A\right)}^{-1}\dots N{\left(s_{1}I-A\right)}^{-1}b,\end{array}$$

$${\tilde{\zeta }}^M=\left\{{\tilde{H}}_1\left(s_1\right),{\tilde{H}}_2\left(s_1,s_2\right),\dots ,{\tilde{H}}_M\left(s_1,\dots ,s_M\right)\right\},$$

$$\begin{array}{rl}& {\tilde{H}}_k\left(s_1,\dots ,s_k\right)=c{\left(s_{k}I-\left(A+F\right)\right)}^{-1}\\ & N{\left(s_{k-1}I-\left(A+F\right)\right)}^{-1}\dots N{\left(s_{1}I-\left(A+F\right)\right)}^{-1}b,\end{array}$$

$$\displaystyle\begin{split}&\left\|\zeta^{M}-\widetilde{\zeta}^{M}\right\|_{H_{2}}^{2}=\left(\begin{bmatrix}c&-c\end{bmatrix}\otimes\begin{bmatrix}c&-c\end{bmatrix}\right)\sum_{j=0}^{M}\Biggl{[}\biggl{(}-\begin{bmatrix}A&0\\ 0&A+F\end{bmatrix}\otimes\\ &\quad\begin{bmatrix}I_{n}&0\\ 0&I_{n}\end{bmatrix}-\begin{bmatrix}I_{n}&0\\ 0&I_{n}\end{bmatrix}\otimes\begin{bmatrix}A&0\\ 0&A+F\end{bmatrix}\biggr{)}^{-1}\begin{bmatrix}N&0\\ 0&N\end{bmatrix}\otimes\begin{bmatrix}N&0\\ 0&N\end{bmatrix}\Biggr{]}^{j}\\ &\quad\biggl{(}-\begin{bmatrix}A&0\\ 0&A+F\end{bmatrix}\otimes\begin{bmatrix}I_{n}&0\\ 0&I_{n}\end{bmatrix}-\begin{bmatrix}I_{n}&0\\ 0&I_{n}\end{bmatrix}\otimes\begin{bmatrix}A&0\\ 0&A+F\end{bmatrix}\biggr{)}^{-1}\\ &\quad\left(\begin{bmatrix}b\\ b\end{bmatrix}\otimes\begin{bmatrix}b\\ b\end{bmatrix}\right).\end{split}$$

$$\displaystyle\left\|\zeta-\widetilde{\zeta}\right\|_{H_{2}}^{2}=\left(\begin{bmatrix}c&-c\end{bmatrix}\otimes\begin{bmatrix}c&-c\end{bmatrix}\right)\biggl{(}-\begin{bmatrix}A&0\\ 0&A+F\end{bmatrix}\otimes\begin{bmatrix}I_{n}&0\\ 0&I_{n}\end{bmatrix}-$$

$$\displaystyle\qquad\begin{bmatrix}I_{n}&0\\ 0&I_{n}\end{bmatrix}\otimes\begin{bmatrix}A&0\\ 0&A+F\end{bmatrix}-\begin{bmatrix}N&0\\ 0&N\end{bmatrix}\otimes\begin{bmatrix}N&0\\ 0&N\end{bmatrix}\biggr{)}^{-1}\left(\begin{bmatrix}b\\ b\end{bmatrix}\otimes\begin{bmatrix}b\\ b\end{bmatrix}\right).$$

$$\begin{array}{rl}{\Vert H_1\left(s_1\right)-{\tilde{H}}_1\left(s_1\right)\Vert }_{H_2}& \propto \mathcal{O}\left({\Vert F\Vert }_2\right),\\ {\Vert H_2\left(s_1,s_2\right)-{\tilde{H}}_2\left(s_1,s_2\right)\Vert }_{H_2}& \propto \mathcal{O}\left({\Vert F\Vert }_2\right),\\ & ⋮\\ {\Vert H_M\left(s_1,\dots ,s_M\right)-{\tilde{H}}_M\left(s_1,\dots ,s_M\right)\Vert }_{H_2}& \propto \mathcal{O}\left({\Vert F\Vert }_2\right),\end{array}$$

$$H_M\left(s_1,\dots ,s_M\right)=c{\left(s_M{I}_n-A\right)}^{-1}N{\left(s_{M-1}{I}_n-A\right)}^{-1}\dots N{\left(s_1{I}_n-A\right)}^{-1}b$$

$${\tilde{H}}_M\left(s_1,\dots ,s_M\right)=$$

Full text

Like for any MOR algorithm, these two algorithms also require solving multiple linear systems as part of the model reduction process. For reducing very large dynamical systems, which is now-a-days becoming a norm, scaling of such linear systems with respect to input dynamical system size is a bottleneck. For efficiency, these linear systems are often solved by an iterative solver, which introduces approximation errors. Hence, stability analysis of MOR algorithms with respect to inexact linear solves is important. In our past work, we have shown that under 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 novel contribution is the more intuitive methodology for achieving this. This approach exploits the fact that in TBIRKA a bilinear dynamical system can be represented by a finite set of functions, which was not possible in BIRKA (because infinite such functions were needed there). The stability analysis techniques that we propose here can be extended to many other methods for doing MOR of bilinear dynamical systems, e.g., using balanced truncation or the ADI methods.

Rajendra Choudhary and Kapil Ahuja

1 Introduction

A dynamical system, usually represented by a set of differential equations, can be linear or non-linear [1, 2, 3]. Linear dynamical systems have been studied more than the non-linear ones because of the obvious ease in working with them. Bilinear dynamical systems form a good bridge between the linear and the non-linear cases, and are usually approximated by a varying degree of bilinearity [4, 5, 6]. In this manuscript, we focus on bilinear dynamical systems.

Dynamical systems coming from different real world applications are very large in size. Thus, simulation and computation with such systems is prohibitively expensive. Model Order Reduction (MOR) techniques provide a smaller system that besides being cheaper to work with, also replicates the input-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 recently studied, 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 projection method. BIRKA (Bilinear Iterative Rational Krylov Algorithm) [15] is a very popular algorithm based upon this technique for reducing first-order bilinear dynamical systems111First-order implies that the highest derivative of the state variable in the dynamical system is one. Second-order and higher-orders are similarly defined.. BIRKA’s biggest drawback is that it does not scale well in time (with respect to increase in the size of the input dynamical system). A cheaper variant of BIRKA, called TBIRKA (Truncated Bilinear Iterative Rational Krylov Algorithm) [21, 22] has also been proposed.

Like in any other MOR algorithm, in BIRKA and TBIRKA also, people often use direct methods like LU-factorization, etc., to solve the arising linear systems, which are too expensive (i.e., computational complexity of $\mathcal{O}(n^{3})$, where $n$ is the original system size). A common solution to this scaling problem is to use iterative methods like the Krylov subspace methods, etc., which have a reduced computational complexity (i.e., $\mathcal{O}(n\times nnz)$, where $nnz$ is the number of nonzeros in the system matrix) [23, 24]. Although iterative methods are cheap, they are inexact too. Hence, studying stability of the underlying MOR algorithm (here BIRKA and TBIRKA) with respect to such approximate linear solves becomes important [25, 26].

One of the first works that performed such a stability analysis focused on popular MOR algorithms for first-order linear dynamical systems [27]. Here, the authors briefly mention that their analysis would be easily carried from the first-order to the second-order case. A detailed stability analysis focusing on reducing second-order linear dynamical systems has been done in [28]. A different kind of stability analysis for MOR of second-order linear dynamical systems has been done in [29]. In this, the authors first show that the SOAR algorithm (second order Arnoldi) is unstable with respect to the machine precision errors (and not inexact linear solves). Then, they propose a Two-level orthogonal Arnoldi (TOAR) algorithm that cures this instability of SOAR. An extended stability analysis for BIRKA (as above, a popular MOR algorithm for first-order bilinear dynamical systems) has been recently done in [30]. For rest of this manuscript, whenever stability analysis is referred, we mean it with respect to inexact linear solves.

We follow the stability analysis framework of BIRKA from [30] and propose equivalent theorems for TBIRKA. The approach here is slightly different, which forms our most novel contribution. Norm of the dynamical system plays an important role in stability analysis (the kind of norm is discussed later). In BIRKA stability analysis, a single expression for bilinear dynamical system norm is used (involving a Volterra series). In TBIRKA stability analysis, a similar single expression (involving a truncated Volterra series) leads to complications. Alternatively, in TBIRKA, because of truncation, the bilinear dynamical system can be represented by a finite set of functions. This was not possible in BIRKA where infinite such functions were needed. Thus, in TBIRKA stability analysis, we use norm of all such functions leading to easier derivations. Our stability analysis, as done for BIRKA earlier and for TBIRKA here, can be easily extended to other MOR algorithms for bilinear dynamical 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, we first give a brief overview of MOR for bilinear dynamical systems using a projection 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. We analyze 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 of this paper, we use the terms and notations as listed below.

1. a.

The $H_{2}-$norm is a functional norm defined as [27, 21, 22]

[TABLE]

where $i$ denotes $\sqrt{-1}$. Here, we assume that all $H_{2}-$norms computed further exist. In other words, the improper integrals defined by the $H_{2}-$norm give finite value. This is a reasonable assumption because this happens often in practice (see [27], where stability analysis of IRKA is done as well as [30], where stability analysis of BIRKA is done). 2. b.

The $H_{\infty}-$norm is also a functional norm, defined as [27, 21, 22]

[TABLE] 3. c.

The Kronecker product between two matrices $P$ (of size $m\times n$), and $Q$ (of size $s\times t$) is defined as

[TABLE]

where $p_{ij}$ is an element of matrix $P$ and order of $P\otimes Q$ is $ms\times nt$. 4. d.

$vec$ operator on a matrix $P$ is defined as

[TABLE] 5. e.

Also, $I_{n}$ denotes an identity matrix of size $n\times n$ and $\mathbb{R}$ denotes the set of real numbers.

2 Background

A first-order bilinear dynamical system is usually represented as [14, 15]

[TABLE]

where $A,N\in\mathbb{R}^{n\times n}$, $b\in\mathbb{R}^{n\times 1}\ \textnormal{and}\ c\in\mathbb{R}^{1\times n}$. Also, $u(t):\mathbb{R}\rightarrow\mathbb{R},\ y(t):\mathbb{R}\rightarrow\mathbb{R}\ \textnormal{and}\ x(t):\mathbb{R}\rightarrow\mathbb{R}^{n}$, represent input, output and state of the bilinear dynamical system, respectively. This is a Single Input Single Output (SISO) system, which we have chosen for ease of our analysis. We plan to look at Multiple Input Multiple Output (MIMO) systems as part of our future work.

A bilinear dynamical system can also be represented by a series of transfer functions, i.e.,

[TABLE]

where $\zeta^{k}=\left\{H_{1}\left(s_{1}\right),\ H_{2}\left(s_{1},\ s_{2}\right),\ \ldots,\ H_{k}\left(s_{1},\ s_{2},\ \ldots,\ s_{k}\right)\right\}$. Here,

$H_{k}\left(s_{1},\ s_{2},\ \ldots,\ s_{k}\right)$ is called the $k^{th}$ order transfer function of the bilinear dynamical system and is defined as

[TABLE]

After reduction, the bilinear dynamical system (2) can be represented as

[TABLE]

where $A_{r},\ N_{r}\in\mathbb{R}^{r\times r}$, b_{r}\in\mathbb{R}^{r\times 1}\ \textnormal{and}\ c_{r}\in\mathbb{R}^{1\times r}\ with $r\ll n$. The main goal of model reduction is to approximate $\zeta$ by $\zeta_{r}$ in an appropriate norm, such that for all admissible inputs, $y_{r}(t)$ is nearly same to $y(t)$.

As mentioned earlier, we use interpolatory projection for performing model reduction. The two common and standard algorithms here, BIRKA and TBIRKA, use a Petrov-Galerkin projection. Let $\mathcal{V}_{r}$ and $\mathcal{W}_{r}$ be the two $r$-dimensional subspaces, such that $\mathcal{V}_{r}=Range(V_{r})$ and $\mathcal{W}_{r}=Range(W_{r})$, where $V_{r}\in\mathbb{R}^{n\times r}$ and $W_{r}\in\mathbb{R}^{n\times r}$ are matrices. Also, let $\left(W_{r}^{T}V_{r}\right)$ be invertible222Obtaining such an invertible matrix is not difficult [15, 22].. 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), we get the reduced system as

[TABLE]

Comparing the above two equations with their respective equations in (5), we get a relation between the original system matrices and the reduced system matrices, i.e.,

[TABLE]

One way of obtaining subspaces $\mathcal{V}_{r}$ and $\mathcal{W}_{r}$ is to use Volterra series interpolation. Further, to decide where to interpolate so as to obtain an optimal reduced model, an $H_{2}-$optimization problem is commonly solved (Theorem 4.7 from [21]).

Theorem 1**.**

[21*]**

Let $\zeta$ be a bilinear system of order n. Let $\zeta_{r}$ be an $H_{2}-$ optimal approximation of order r. Then, $\zeta_{r}$ satisfies the following multi-point Volterra series interpolation conditions:*

[TABLE]

where $\phi_{l_{1},\ l_{2},\ \ldots,\ l_{k}}$ and $\lambda_{l_{1}},\ \lambda_{l_{2}},\ \ldots,\ \lambda_{l_{k}}$ are residues and poles of the transfer function $H_{r_{k}}$ associated with $\zeta_{r}$, respectively.

BIRKA is designed in such a way that at convergence, the conditions of Theorem 1 are satisfied leading to a locally $H_{2}-$optimal reduced model. Algorithm 1 lists BIRKA.

TBIRKA is similar to BIRKA in most aspects except that it performs a truncated Volterra series interpolation. Here, instead of $\zeta$ in (2)-(3), they work with $\zeta^{M}$, which is defined as

[TABLE]

with $H_{k}\left(s_{1},\ \ldots,\ s_{k}\right)$ for $k=1,\ \ldots,\ M$ is given by (4). Equivalent of Theorem 1 above is as follows (Theorem 4.8 from [21]):

Theorem 2**.**

[21*]**

Let $\zeta=(A,N,b,c)$ be an order $n$ bilinear system and $\zeta^{M}$ be the polynomial system determined by $\zeta$. Let $\zeta_{r}=(A_{r},N_{r},b_{r},c_{r})$ be a bilinear system of order $r$, and define $\zeta_{r}^{M}$ as the polynomial system determined by $\zeta_{r}$. Suppose that $\zeta_{r}^{M}$ is an $H_{2}-$optimal approximation to $\zeta^{M}$. Then $\zeta_{r}^{M}$ satisfies*

[TABLE]

[TABLE]

where $\phi_{l_{1},\ l_{2},\ \ldots,\ l_{k}}$ and $\lambda_{l_{1}},\ \lambda_{l_{2}},\ \ldots,\ \lambda_{l_{k}}$ are residues and poles of the transfer function $H_{r_{k}}$ associated with $\zeta_{r}^{M}$, respectively.

Algorithm 2 lists TBIRKA.

Both BIRKA and TBIRKA in turn require solving large sparse linear systems of equations. If we compare Algorithm 1 and 2, we realize that the cost of linear solvers at each step of the While loop in the former is $\mathcal{O}\left(nr\times nr\right)$ and in the latter is $\mathcal{O}\left(M\times nr\times nr\right)$. This makes it seem that TBIRKA is more expensive than BIRKA. However, TBIRKA is implemented in such a way that the Kronecker products are avoided leading to cost of linear solves at each step of the While loop to be $\mathcal{O}\left(M\times r\times\left(n\times n\right)\right)$. Thus, if $M<r$, which is usually the case (e.g, $M=3\ \textnormal{or}\ 4$ and $r\gg M$), then TBIRKA is more efficient than BIRKA. 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 systems introduces approximation errors. We have done a detailed stability analysis of BIRKA with respect to the inexact linear solves in [30], and we briefly revisit this next. Generally, accuracy is the metric that tells about the correctness in the output of an inexact algorithm. Due to unavailability of the exact output, it is not possible to determine accuracy [25, 30]. A more easier metric is stability. Backward stability is one such notation, which says “A backward stable algorithm gives exactly the right output to nearly the right input” [25]. In our context, theoretically we obtain two reduced systems. One by applying an inexact MOR algorithm (with iterative linear solves) on the original full model, and other by applying the same MOR algorithm but exactly (with direct linear solves) on a perturbed full model (the perturbation is introduced in the original full model as part of stability analysis, and is an unknown quantity). If these two reduced systems are equal (first condition), with the difference between the original full model and the perturbed full model equal to the order of perturbation (second condition), then the MOR algorithm under consideration is called backward stable. The two theorems summarizing this 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 BIRKA satisfies the first condition of backward stability with respect to these solves.

Theorem 4**.**

[30]** Let \widehat{Q}=\Biggl{(}-\begin{bmatrix}A&0\\ 0&A\end{bmatrix}\otimes\begin{bmatrix}I_{n}&0\\ 0&I_{n}\end{bmatrix}-\begin{bmatrix}I_{n}&0\\ 0&I_{n}\end{bmatrix}\otimes\begin{bmatrix}A&0\\ 0&A\end{bmatrix}-\begin{bmatrix}N&0\\ 0&N\end{bmatrix}\otimes\begin{bmatrix}N&0\\ 0&N\end{bmatrix}\Biggr{)}, where $I_{n}$ is an identity matrix of size $n\times n$ and $\otimes$ denotes the standard Kronecker product. Also, let $\widehat{\widehat{F}}=\left(I_{2n}\otimes\widehat{F}+\widehat{F}\otimes I_{2n}\right)$ with $\widehat{F}=\begin{bmatrix}0&0\\ 0&F\end{bmatrix}$, where $F$ is the perturbation introduced in $A$ matrix of the input dynamical system and $I_{2n}$ is an identity matrix of size $2n\times 2n$. If $\widehat{Q}$ is invertible, $\left\|\widehat{Q}^{-1}\right\|_{2}<1$, and $\left\|\widehat{\widehat{F}}\right\|_{2}<1$, then BIRKA satisfies the second condition of backward stability with respect to the inexact linear solves.

Next, we analyze the backward stability of TBIRKA. Here, the first condition is satisfied in a way similar to that of BIRKA except that some extra orthogonality conditions are imposed on the linear solver (discussed below).

Theorem 5**.**

Let the inexact linear solves in TBIRKA (lines 3b. and 3c. of Algorithm 2) are solved using a Petrov-Galerkin framework with

[TABLE]

where $R_{b_{1}}$ and $R_{b_{j}}$ are the residuals in the first equations of 3b. and 3c. of Algorithm 2, respectively; $R_{c_{1}}$ and $R_{c_{j}}$ are the residuals in the second equations of 3b. and 3c. of Algorithm 2, respectively; and $j=2,\ \ldots,\ M.$ Then, TBIRKA satisfies the first condition of backward stability with respect to these solves.

Proof.

Follows the same pattern as the proof for Theorem 2.1 in [30]. ∎

From the above theorem, we infer that the underlying iterative solver should firstly be based upon a Petrov-Galerkin framework. Since BiConjugate Gradient (i.e., BiCG) is one such algorithm [23, 24], we propose its use in TBIRKA. This is exactly same as for BIRKA. Secondly, this particular solver should also satisfy the extra orthogonalities (7).

These orthogonalities have a form similar to the orthogonalities required while reducing second order linear dynamical systems ((23) and (24) in [28]), and can be easily satisfied by using a recycling variant of the underlying iterative 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 satisfy the similar orthogonalities there, without any extra cost, the authors use Recycling Conjugate Gradient (i.e., RCG) [33]. Since here BiCG is the ideal iterative solver (as discussed above), we propose the use of Recycling BiConjugate Gradient (i.e., RBiCG) [32, 31], which would ensure that orthogonalities given by (7) are satisfied without any extra cost.

To satisfy the second condition of backward stability of TBIRKA, we need to show that

[TABLE]

where $\zeta^{M}$ is given by (6) or

[TABLE]

[TABLE]

and assuming perturbation in $A$ matrix of the input dynamical system (as for BIRKA stability; see Theorem 4 earlier).

One way to satisfy (8) is to use the definition of the $H_{2}-$norm of $\zeta^{M}-\widetilde{\zeta}^{M}$, i.e., from Lemma 5.1 of [22]

[TABLE]

This approach is followed in satisfying the second condition of backward stability of BIRKA, but for TBIRKA it turns out to be more challenging. The reason for this is that the definition of the $H_{2}-$norm of $\zeta-\widetilde{\zeta}$ used in BIRKA is different from (LABEL:eq:H_2_norm_error_system_truncated_Voterra)333Recall, in BIRKA we work with $\zeta$ rather than $\zeta^{M}$., i.e., from Corollary 4.1 of [15] or Theorem 4.5 of [21]

[TABLE]

Another way to satisfy (8) in case of TBIRKA, which turns out to be more easier, is to show that

[TABLE]

where $H_{k}\left(s_{1},\ \ldots,\ s_{k}\right)$ for $k=1,\ \ldots,\ M$ is given by (4) and (9b), and $\widetilde{H}_{k}\left(s_{1},\ \ldots,\ s_{k}\right)$ for $k=1,\ \ldots,\ M$ is given by (10b). This way was not possible in BIRKA because there $M\rightarrow\infty$ (see (2)-(4)).

3 Satisfying the Second Condition of Backward Stability

We use induction to prove (12). To prove this condition, we first abstract out the term containing the perturbation $F$ from the normed difference between the transfer function of the original system $\left(H_{M}\left(s_{1},\ \ldots,\ s_{M}\right)\right)$ and the transfer function of the perturbed system $\left(\widetilde{H}_{M}\left(s_{1},\ \ldots,\ s_{M}\right)\right)$ (in Lemma 6). Next, we use induction to prove that this abstracted term is $\mathcal{O}\left(\left\|F\right\|_{2}\right)$ (in Lemma 7). Finally, we use the result of these two lemmas to prove (12) (in Theorem 8). Note, that in all our subsequent derivations, we assume that all inverses used exist. This is an acceptable assumption because the inverse of matrices arising here are of the form as in [27] and [30] (the papers that discuss stability of IRKA and BIRKA, respectively).

Lemma 6**.**

If

[TABLE]

and

[TABLE]

then the $H_{2}-$norm of their difference

[TABLE]

where, $\mathcal{K}\left(s_{i}\right)=\left(s_{i}I_{n}-A\right)$ for $i=1,\ \ldots,\ M$ and

[TABLE]

Proof.

Using the definition of $H_{2}-$norm (1), we get

[TABLE]

[TABLE]

Using $U\left(s_{1},\ \ldots,\ s_{M}\right)$ given by (13), $\left\|XYZ\right\|_{F}\leq\left\|X\right\|_{F}\left\|YZ\right\|_{F}$, $\left\|YZ\right\|_{F}\leq\left\|Y\right\|_{F}\left\|Z\right\|_{2}$, and comparison integral inequality444This inequality says if $f\left(x\right)$ and $g\left(x\right)$ are integrable over $\left[a,b\right]$ and $f\left(x\right)\leq g\left(x\right)$, then $\int_{a}^{b}f\left(x\right)dx\leq\int_{a}^{b}g\left(x\right)dx$. Note that although we have improper integrals here, this inequality still holds because of the earlier assumption that such integrals give a finite value. [34] for any matrices $X,\ Y$ and $Z$, in the above equation, we have

[TABLE]

From the mean value theorem of integration [34] we know

[TABLE]

[TABLE]

Using this property in (14) we get555As mentioned in Footnote 4, the improper integral does not affect application of this mean value theorem because all such integrals are assumed to give a finite value.

[TABLE]

∎

Lemma 7**.**

*If $U_{M}=U(s_{1},\ \ldots,\ s_{M})$ defined in (13), $\left\|F\right\|_{2}<1$, and

$\left\|\mathcal{K}^{-1}\left(s\right)\right\|_{H_{\infty}}<1$, where $\mathcal{K}\left(s\right)=\left(sI_{n}-A\right)$. Then we have*

[TABLE]

Proof.

We prove this by mathematical induction.

Base Case

1. (i)

*First subsystem

*This is the linear system case, which has been already proved in [27] (see Theorem 4.3 of [27]). 2. (ii)

Second subsystem

Substituting $M=2$ in (13), we get

[TABLE]

If $\left\|F\mathcal{K}^{-1}\left(s_{2}\right)\right\|_{H_{\infty}}<1$ and $\left\|\mathcal{K}^{-1}\left(s_{1}\right)F\right\|_{H_{\infty}}<1$, then by the Neumann series, we get666From [35, page 527], we know $\left(I-A\right)^{-1}=\sum\limits_{k=0}^{\infty}A^{k}$ when $\left\|A\right\|<1$ for any matrix norm. Here, for the first inequality we have $\left\|F\mathcal{K}^{-1}\left(s_{2}\right)\right\|_{H_{\infty}}<1$ or $\underset{\omega_{2}\in\mathbb{R}}{max}\left\|F\mathcal{K}^{-1}\left(i\omega_{2}\right)\right\|_{2}<1$, and hence, the applicable matrix norm is $2-$norm. Similarly for the second inequality.

[TABLE]

Taking $H_{\infty}-$norm on both sides, we get

[TABLE]

Using $\left\|XY\right\|_{2}\leq\left\|X\right\|_{2}\left\|Y\right\|_{2}$ and $\left\|X+Y\right\|_{2}\leq\left\|X\right\|_{2}+\left\|Y\right\|_{2}$, for any two matrices $X$ and $Y$, in the above equation, we get

[TABLE]

Technically by definition of the $H_{\infty}-$norm and how $\mathcal{K}\left(s\right)$ is defined in our hypotheses, $\left\|\mathcal{K}\left(s_{1}\right)\right\|_{H_{\infty}}=\left\|\mathcal{K}\left(s_{2}\right)\right\|_{H_{\infty}}=\left\|\mathcal{K}\left(s\right)\right\|_{H_{\infty}}$, however, for sake of exposition, we keep them separate. Similarly for the $H_{\infty}-$norm of inverses of $\mathcal{K}\left(s_{1}\right)$ and $\mathcal{K}\left(s_{2}\right)$.

To abstract $\left\|F\right\|_{2}$ out from the above inequality, let us look at $\underset{\omega_{2}\in\mathbb{R}}{max}\left\|\left(I_{n}-F\mathcal{K}^{-1}\left(i\omega_{2}\right)\right)^{-1}\right\|_{2}$ separately. Recall, while applying Neumann series we assumed that $\left\|F\mathcal{K}^{-1}\left(s_{2}\right)\right\|_{H_{\infty}}<1$ or $\underset{\omega_{2}\in\mathbb{R}}{max}\left\|F\mathcal{K}^{-1}\left(i\omega_{2}\right)\right\|_{2}<1$. Since the maximum of such a norm is less than one, we have for all $\omega_{2}\in\mathbb{R}$, $\left\|F\mathcal{K}^{-1}\left(i\omega_{2}\right)\right\|_{2}<1$. Using this along with Lemma 2.3.3 from [36]777If $F\in\mathbb{R}^{n\times n}$ and $\left\|F\right\|_{p}<1$, then $I-F$ is nonsingular and $\left(I-F\right)^{-1}=\sum\limits_{k=0}^{\infty}F^{k}$ with $\left\|\left(I-F\right)^{-1}\right\|_{p}\leq\dfrac{1}{1-\left\|F\right\|_{p}}$. in the above expression, we get

[TABLE]

If we assume $\left\|F\right\|_{2}<1$ and $\left\|\mathcal{K}^{-1}\left(s\right)\right\|_{H_{\infty}}<1$ (as in our hypotheses), then using earlier used matrix norm properties, we get

[TABLE]

as assumed for applying Neumann series earlier as well as Lemma 2.3.3 from [36] above. Thus, no extra assumptions beyond those in hypotheses are needed. Further, we also get

[TABLE]

Similarly, we can bound the last term of (15) as follows:

[TABLE]

Substituting (16)-(17) and (18)-(19) in (15), we get

[TABLE]

From the above inequality it is clear that if $\left\|F\right\|_{2}\left\|\mathcal{K}^{-1}\left(s_{2}\right)\right\|_{H_{\infty}}<1$ and $\left\|\mathcal{K}^{-1}\left(s_{1}\right)\right\|_{H_{\infty}}\left\|F\right\|_{2}<1$, which are true from our hypotheses, then

[TABLE] 3. (iii)

*Third subsystem

*Using $M=3$ in (13), we get

[TABLE]

Following the same steps as in the case of the second subsystem we get

[TABLE]

Again, from the above inequality it is clear that if $\left\|F\right\|_{2}\left\|\mathcal{K}^{-1}\left(s_{2}\right)\right\|_{H_{\infty}}$, $\left\|F\right\|_{2}\left\|\mathcal{K}^{-1}\left(s_{3}\right)\right\|_{H_{\infty}}$ and $\left\|\mathcal{K}^{-1}\left(s_{1}\right)\right\|_{H_{\infty}}\left\|F\right\|_{2}$ are all less than $1$, which are true from our hypotheses, then

[TABLE]

Induction Hypothesis

From (13), we know

[TABLE]

Let

[TABLE]

Induction Step

Again, from (13), we know

[TABLE]

We first write $U_{M+1}$ in terms of $U_{M}$. Thus, in the above equation using Neumann series, we get

[TABLE]

In the above equation, taking $N\mathcal{K}^{-1}\left(s_{M}\right)$ common from the first two terms of the bigger bracket, we have

[TABLE]

Now we look at expression of $U_{M}$. Multiplying $\mathcal{K}^{-1}\left(s_{M-1}\right)\ldots\mathcal{K}^{-1}\left(s_{1}\right)$ on both the sides of (20) from left, we get

[TABLE]

Substituting (22) in (21), we get

[TABLE]

Taking $H_{\infty}-$norm on both sides, we get

[TABLE]

As earlier, using the norm inequality properties in the above equation, we get

[TABLE]

Similar to (16) and (17), here also, using Lemma 2.3.3 from [36] we get

[TABLE]

From induction hypothesis we know $\left\|U_{M}\right\|_{H_{\infty}}\propto\mathcal{O}\left(\left\|F\right\|_{2}\right)$. Using this we get

[TABLE]

∎

Theorem 8**.**

If

[TABLE]

$\left\|F\right\|_{2}<1$, and $\left\|\mathcal{K}^{-1}\left(s\right)\right\|_{H_{\infty}}<1$, where $\mathcal{K}\left(s\right)=\left(sI_{n}-A\right)$, then

[TABLE]

TBIRKA satisfies the second condition of backward stability with respect to inexact linear solves.

4 Conclusions & Future Work

In this paper, we apply iterative linear solvers during model order reduction (MOR) of bilinear dynamical systems. Since such solvers are inexact, the stability of the underlying MOR algorithm, with respect to these approximation errors, is important. Here, we extend the earlier stability analysis done for BIRKA in [30], which is a standard algorithm for obtaining a reduced bilinear dynamical system, to TBIRKA, a cheaper variant of BIRKA.

Proving that an algorithm is stable, typically requires satisfying two conditions. In TBIRKA, fulfilling the first condition for stability leads to constraints on the iterative linear solver, which are similar to those obtained during BIRKA’s stability analysis. The second condition for a stable TBIRKA is satisfied using an approach different than the one used in BIRKA, and is more intuitive.

Our first future direction is to extend our analysis from SISO (Single Input Single Output) to MIMO (Multiple Input Multiple Output) systems. The stability analysis as done for BIRKA earlier and TBIRKA here, all give us sufficiency conditions for a stable underlying MOR algorithm. Hence, second, we plan to derive the necessary conditions for the same. In recent years, there have been a lot of efforts in performing data-driven MOR algorithm (specially using Leowner framework [19]). Similar linear systems and challenges associated with the scaling arise here as well. Our third future direction is to apply this stability analysis to such classes of algorithms as well.