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. At each iteration, IRKA does an Hermite type interpolation of the original system transfer function. Each interpolation requires solving [math]\displaystyle{ r }[/math] shifted pairs of linear systems, each of size [math]\displaystyle{ n \times n }[/math]; where [math]\displaystyle{ n }[/math] is the original system order, and [math]\displaystyle{ r }[/math] is the desired reduced model order (usually [math]\displaystyle{ r \ll n }[/math]). The algorithm was first introduced by Gugercin, Antoulas and Beattie in 2008. It is based on a first order necessary optimality condition, initially investigated by Meier and Luenberger in 1967. The first convergence proof of IRKA was given by Flagg, Beattie and Gugercin in 2012, for a particular kind of systems.

MOR as an optimization problem

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

[math]\displaystyle{ \begin{cases} \dot{x}(t) = A x(t) + b v(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. }[/math]

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

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

Assume [math]\displaystyle{ G }[/math] is stable. Given [math]\displaystyle{ r \lt n }[/math], MOR tries to approximate the transfer function [math]\displaystyle{ G }[/math], by a stable rational transfer function [math]\displaystyle{ G_r }[/math], of order [math]\displaystyle{ r }[/math]:

[math]\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. }[/math]

A possible approximation criterion is to minimize the absolute error in [math]\displaystyle{ H_{2} }[/math] norm:

[math]\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 . }[/math]

This is known as the [math]\displaystyle{ H_{2} }[/math] optimization problem. This problem has been studied extensively, and it is known to be non-convex; 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 [math]\displaystyle{ H_{2} }[/math] problem, is of great importance for the IRKA algorithm.

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

[math]\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 . }[/math]

Note that the poles [math]\displaystyle{ \lambda_i(A_r) }[/math] are the eigenvalues of the reduced [math]\displaystyle{ r \times r }[/math] matrix [math]\displaystyle{ A_r }[/math].

Hermite interpolation

An Hermite interpolant [math]\displaystyle{ G_r }[/math] of the rational function [math]\displaystyle{ G }[/math], through [math]\displaystyle{ r }[/math] distinct points [math]\displaystyle{ \sigma_1, \ldots, \sigma_r \in \mathbb{C} }[/math], has components:

[math]\displaystyle{ A_r = W_r^* A V_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; }[/math]

where the matrices [math]\displaystyle{ V_r = ( v_1 \mid \ldots \mid v_r ) \in \mathbb{C}^{n \times r} }[/math] and [math]\displaystyle{ W_r = ( w_1 \mid \ldots \mid w_r ) \in \mathbb{C}^{n \times r} }[/math] may be found by solving [math]\displaystyle{ r }[/math] dual pairs of linear systems, one for each shift [Theorem 1.1]:

[math]\displaystyle{ (\sigma_i I-A) v_i=b, \quad (\sigma_i I-A)^* w_i=c, \quad \forall \, i=1,\ldots,r . }[/math]

IRKA algorithm

As can be seen from the previous section, finding an Hermite interpolator [math]\displaystyle{ G_r }[/math] of [math]\displaystyle{ G }[/math], through [math]\displaystyle{ r }[/math] 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 [math]\displaystyle{ r }[/math] arbitrary interpolation points (closed under conjugation), and then, at each iteration [math]\displaystyle{ m }[/math], it imposes the first order necessary optimality condition of the [math]\displaystyle{ H_2 }[/math] problem:

1. find the Hermite interpolant [math]\displaystyle{ G_r }[/math] of [math]\displaystyle{ G }[/math], through the actual [math]\displaystyle{ r }[/math] shift points: [math]\displaystyle{ \sigma_1^m,\ldots,\sigma_r^m }[/math].

2. update the shifts by using the poles of the new [math]\displaystyle{ G_r }[/math]: [math]\displaystyle{ \sigma_i^{m+1} = -\lambda_i(A_r), \, \forall \, i=1,\ldots,r . }[/math]

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:

[math]\displaystyle{ \frac{ |\sigma_i^{m+1}-\sigma_i^m| }{|\sigma_i^m|} \lt \text{tol}, \, \forall \, i=1,\ldots,r . }[/math]

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

[math]\displaystyle{ (\sigma_i^m I-A) v_{i} = b, \quad (\sigma_i^m I-A)^* w_i = c, \quad \forall \, i=1,\ldots,r . }[/math]

Also, updating the shifts requires finding the [math]\displaystyle{ r }[/math] poles of the new interpolant [math]\displaystyle{ G_r }[/math]. That is, finding the [math]\displaystyle{ r }[/math] eigenvalues of the reduced [math]\displaystyle{ r \times r }[/math] matrix [math]\displaystyle{ A_r }[/math].

Pseudocode

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

algorithm IRKA
    input: [math]\displaystyle{ A,b,c }[/math], [math]\displaystyle{ \text{tol}\gt 0 }[/math], [math]\displaystyle{ \sigma_1,\ldots,\sigma_r }[/math] closed under conjugation
        [math]\displaystyle{ (\sigma_i I-A)v_i=b, \, \forall \, i=1,\ldots,r }[/math] % Solve primal systems
        [math]\displaystyle{ (\sigma_i I-A)^* w_i=c, \, \forall \, i=1,\ldots,r }[/math] % Solve dual systems

    while relative change in {[math]\displaystyle{ \sigma_{i} }[/math]} > tol
        [math]\displaystyle{ A_{r} = W_r^* AV_r }[/math] % Reduced order matrix
        [math]\displaystyle{ \sigma_i = -\lambda_i(A_r), \, \forall \, i=1,\ldots,r }[/math] % Update shifts, using poles of [math]\displaystyle{ G_{r} }[/math]
        [math]\displaystyle{ (\sigma_i I-A)v_i=b, \, \forall \, i=1,\ldots,r }[/math] % Solve primal systems
        [math]\displaystyle{ (\sigma_i I-A)^{*}w_{i}=c, \, \forall \, i=1,\ldots,r }[/math] % Solve dual systems
    end while

    return [math]\displaystyle{ A_r=W_r^* AV_r, \, b_r=W_r^{*}b, \, c_r^T=c^T V_r }[/math] % Reduced order model

Convergence

A SISO linear system is said to have symmetric state space (SSS), whenever: [math]\displaystyle{ A=A^{T}, \, b=c . }[/math] 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. For SSS systems with distinct poles, the following convergence result has been proven: "IRKA is a locally convergent fixed point iteration to a local minimizer of the [math]\displaystyle{ H_{2} }[/math] 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.

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 [Remark 4.1].

See also

Model order reduction

External links

• Model Order Reduction Wiki