Iterative rational Krylov algorithm

The iterative rational Krylov algorithm (IRKA), is an iterative algorithm, useful for model order reduction (MOR) of single-input single-output (SISO) linear time-invariant dynamical systems. ^[1] At each iteration, IRKA does an Hermite type interpolation of the original system transfer function. Each interpolation requires solving ${\displaystyle r}$ shifted pairs of linear systems, each of size ${\displaystyle n\times n}$; where ${\displaystyle n}$ is the original system order, and ${\displaystyle r}$ is the desired reduced model order (usually ${\displaystyle r\ll n}$).

The algorithm was first introduced by Gugercin, Antoulas and Beattie in 2008. ^[2] It is based on a first order necessary optimality condition, initially investigated by Meier and Luenberger in 1967. ^[3] The first convergence proof of IRKA was given by Flagg, Beattie and Gugercin in 2012, ^[4] for a particular kind of systems.

MOR as an optimization problem

Consider a SISO linear time-invariant dynamical system, with input ${\displaystyle v(t)}$, and output ${\displaystyle y(t)}$:

${\displaystyle {\begin{cases}{\dot {x}}(t)=Ax(t)+bv(t)\\y(t)=c^{T}x(t)\end{cases}}\qquad A\in \mathbb {R} ^{n\times n},\,b,c\in \mathbb {R} ^{n},\,v(t),y(t)\in \mathbb {R} ,\,x(t)\in \mathbb {R} ^{n}.}$

Applying the Laplace transform, with zero initial conditions, we obtain the transfer function ${\displaystyle G}$, which is a fraction of polynomials:

${\displaystyle G(s)=c^{T}(sI-A)^{-1}b,\quad A\in \mathbb {R} ^{n\times n},\,b,c\in \mathbb {R} ^{n}.}$

Assume ${\displaystyle G}$ is stable. Given ${\displaystyle r<n}$, MOR tries to approximate the transfer function ${\displaystyle G}$, by a stable rational transfer function ${\displaystyle G_{r}}$, of order ${\displaystyle r}$:

${\displaystyle G_{r}(s)=c_{r}^{T}(sI_{r}-A_{r})^{-1}b_{r},\quad A_{r}\in \mathbb {R} ^{r\times r},\,b_{r},c_{r}\in \mathbb {R} ^{r}.}$

A possible approximation criterion is to minimize the absolute error in ${\displaystyle H_{2}}$ norm:

${\displaystyle G_{r}\in {\underset {\dim({\hat {G}})=r,\,{\hat {G}}{\text{ stable}}}{\operatorname {arg\min } }}\|G-{\hat {G}}\|_{H_{2}},\quad \|G\|_{H_{2}}^{2}:={\frac {1}{2\pi }}\int \limits _{-\infty }^{\infty }|G(ja)|^{2}\,da.}$

This is known as the ${\displaystyle H_{2}}$ optimization problem. This problem has been studied extensively, and it is known to be non-convex; ^[4] which implies that usually it will be difficult to find a global minimizer.

Meier–Luenberger conditions

The following first order necessary optimality condition for the ${\displaystyle H_{2}}$ problem, is of great importance for the IRKA algorithm.

Theorem ( ^[2] [Theorem 3.4] ^[4] [Theorem 1.2]) —  Assume that the ${\displaystyle H_{2}}$ optimization problem admits a solution ${\displaystyle G_{r}}$ with simple poles. Denote these poles by: ${\displaystyle \lambda _{1}(A_{r}),\ldots ,\lambda _{r}(A_{r})}$. Then, ${\displaystyle G_{r}}$ must be an Hermite interpolator of ${\displaystyle G}$, through the reflected poles of ${\displaystyle G_{r}}$:

${\displaystyle G_{r}(\sigma _{i})=G(\sigma _{i}),\quad G_{r}^{\prime }(\sigma _{i})=G^{\prime }(\sigma _{i}),\quad \sigma _{i}=-\lambda _{i}(A_{r}),\quad \forall \,i=1,\ldots ,r.}$

Note that the poles ${\displaystyle \lambda _{i}(A_{r})}$ are the eigenvalues of the reduced ${\displaystyle r\times r}$ matrix ${\displaystyle A_{r}}$.

Hermite interpolation

An Hermite interpolant ${\displaystyle G_{r}}$ of the rational function ${\displaystyle G}$, through ${\displaystyle r}$ distinct points ${\displaystyle \sigma _{1},\ldots ,\sigma _{r}\in \mathbb {C} }$, has components:

${\displaystyle A_{r}=W_{r}^{*}AV_{r},\quad b_{r}=W_{r}^{*}b,\quad c_{r}=V_{r}^{*}c,\quad A_{r}\in \mathbb {R} ^{r\times r},\,b_{r}\in \mathbb {R} ^{r},\,c_{r}\in \mathbb {R} ^{r};}$

where the matrices ${\displaystyle V_{r}=(v_{1}\mid \ldots \mid v_{r})\in \mathbb {C} ^{n\times r}}$ and ${\displaystyle W_{r}=(w_{1}\mid \ldots \mid w_{r})\in \mathbb {C} ^{n\times r}}$ may be found by solving ${\displaystyle r}$ dual pairs of linear systems, one for each shift ^[4] [Theorem 1.1]:

${\displaystyle (\sigma _{i}I-A)v_{i}=b,\quad (\sigma _{i}I-A)^{*}w_{i}=c,\quad \forall \,i=1,\ldots ,r.}$

IRKA algorithm

As can be seen from the previous section, finding an Hermite interpolator ${\displaystyle G_{r}}$ of ${\displaystyle G}$, through ${\displaystyle r}$ given points, is relatively easy. The difficult part is to find the correct interpolation points. IRKA tries to iteratively approximate these "optimal" interpolation points.

For this, it starts with ${\displaystyle r}$ arbitrary interpolation points (closed under conjugation), and then, at each iteration ${\displaystyle m}$, it imposes the first order necessary optimality condition of the ${\displaystyle H_{2}}$ problem:

1. find the Hermite interpolant ${\displaystyle G_{r}}$ of ${\displaystyle G}$, through the actual ${\displaystyle r}$ shift points: ${\displaystyle \sigma _{1}^{m},\ldots ,\sigma _{r}^{m}}$.

2. update the shifts by using the poles of the new ${\displaystyle G_{r}}$: ${\displaystyle \sigma _{i}^{m+1}=-\lambda _{i}(A_{r}),\,\forall \,i=1,\ldots ,r.}$

The iteration is stopped when the relative change in the set of shifts of two successive iterations is less than a given tolerance. This condition may be stated as:

${\displaystyle {\frac {|\sigma _{i}^{m+1}-\sigma _{i}^{m}|}{|\sigma _{i}^{m}|}}<{\text{tol}},\,\forall \,i=1,\ldots ,r.}$

As already mentioned, each Hermite interpolation requires solving ${\displaystyle r}$ shifted pairs of linear systems, each of size ${\displaystyle n\times n}$:

${\displaystyle (\sigma _{i}^{m}I-A)v_{i}=b,\quad (\sigma _{i}^{m}I-A)^{*}w_{i}=c,\quad \forall \,i=1,\ldots ,r.}$

Also, updating the shifts requires finding the ${\displaystyle r}$ poles of the new interpolant ${\displaystyle G_{r}}$. That is, finding the ${\displaystyle r}$ eigenvalues of the reduced ${\displaystyle r\times r}$ matrix ${\displaystyle A_{r}}$.

Pseudocode

The following is a pseudocode for the IRKA algorithm ^[2] [Algorithm 4.1].

algorithm IRKA     input:${\displaystyle A,b,c}$, ${\displaystyle {\text{tol}}>0}$, ${\displaystyle \sigma _{1},\ldots ,\sigma _{r}}$ closed under conjugation         ${\displaystyle (\sigma _{i}I-A)v_{i}=b,\,\forall \,i=1,\ldots ,r}$ % Solve primal systems         ${\displaystyle (\sigma _{i}I-A)^{*}w_{i}=c,\,\forall \,i=1,\ldots ,r}$ % Solve dual systems      while relative change in {${\displaystyle \sigma _{i}}$} > tol         ${\displaystyle A_{r}=W_{r}^{*}AV_{r}}$ % Reduced order matrix         ${\displaystyle \sigma _{i}=-\lambda _{i}(A_{r}),\,\forall \,i=1,\ldots ,r}$ % Update shifts, using poles of ${\displaystyle G_{r}}$${\displaystyle (\sigma _{i}I-A)v_{i}=b,\,\forall \,i=1,\ldots ,r}$ % Solve primal systems         ${\displaystyle (\sigma _{i}I-A)^{*}w_{i}=c,\,\forall \,i=1,\ldots ,r}$ % Solve dual systems     end whilereturn${\displaystyle A_{r}=W_{r}^{*}AV_{r},\,b_{r}=W_{r}^{*}b,\,c_{r}^{T}=c^{T}V_{r}}$ % Reduced order model

Convergence

A SISO linear system is said to have symmetric state space (SSS), whenever: ${\displaystyle A=A^{T},\,b=c.}$ This type of systems appear in many important applications, such as in the analysis of RC circuits and in inverse problems involving 3D Maxwell's equations. ^[4] For SSS systems with distinct poles, the following convergence result has been proven: ^[4] "IRKA is a locally convergent fixed point iteration to a local minimizer of the ${\displaystyle H_{2}}$ optimization problem."

Although there is no convergence proof for the general case, numerous experiments have shown that IRKA often converges rapidly for different kind of linear dynamical systems. ^{[1] [4]}

Extensions

IRKA algorithm has been extended by the original authors to multiple-input multiple-output (MIMO) systems, and also to discrete time and differential algebraic systems ^{[1] [2]} [Remark 4.1].

See also

Model order reduction

References

1. "Iterative Rational Krylov Algorithm". MOR Wiki. Retrieved 3 June 2021.
2. Gugercin, S.; Antoulas, A.C.; Beattie, C. (2008), ${\displaystyle H_{2}}$ Model Reduction for Large-Scale Linear Dynamical Systems, Journal on Matrix Analysis and Applications, 30, SIAM, pp. 609–638
3. L. Meier; D.G. Luenberger (1967), Approximation of linear constant systems, IEEE Transactions on Automatic Control, 12, pp. 585–588
4. G. Flagg; C. Beattie; S. Gugercin (2012), Convergence of the Iterative Rational Krylov Algorithm, Systems & Control Letters, 61, pp. 688–691

External links

• Model Order Reduction Wiki