Iterative rational Krylov algorithm

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.

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

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

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.}$

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}}$.

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 while

    return ${\displaystyle A_{r}=W_{r}^{*}AV_{r},\,b_{r}=W_{r}^{*}b,\,c_{r}^{T}=c^{T}V_{r}}$ % Reduced order model

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]

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].

Model order reduction

• Model Order Reduction Wiki