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 [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.[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 [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;[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 [math]\displaystyle{ H_{2} }[/math] problem, is of great importance for the IRKA algorithm.
Theorem ([2][Theorem 3.4] [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 [4][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 [2][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.[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 [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.[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
References
- ↑ 1.0 1.1 1.2 "Iterative Rational Krylov Algorithm". https://morwiki.mpi-magdeburg.mpg.de/morwiki/index.php/Iterative_Rational_Krylov_Algorithm.
- ↑ 2.0 2.1 2.2 2.3 Gugercin, S.; Antoulas, A.C.; Beattie, C. (2008), [math]\displaystyle{ H_{2} }[/math] Model Reduction for Large-Scale Linear Dynamical Systems, Journal on Matrix Analysis and Applications, 30, SIAM, pp. 609–638
- ↑ L. Meier; D.G. Luenberger (1967), Approximation of linear constant systems, IEEE Transactions on Automatic Control, 12, pp. 585–588
- ↑ 4.0 4.1 4.2 4.3 4.4 4.5 4.6 G. Flagg; C. Beattie; S. Gugercin (2012), Convergence of the Iterative Rational Krylov Algorithm, Systems & Control Letters, 61, pp. 688–691
External links
![]() | Original source: https://en.wikipedia.org/wiki/Iterative rational Krylov algorithm.
Read more |