3.4. Numerical examples 63
64 Linear Systems
10 20 30 40 50
10−5 10−4 10−3 10−2 10−1 100
Rank of Pˆn
||P−Pˆn||F/||P||F SVD IRKA
10 20 30 40 5010−4 10−3 10−2 10−1 100
Rank of Pnˆ
||P−Pˆn||L/||P||L
SVD IRKA
Figure 3.1: Steel profile with n = 1357.
we always compare the constructed low rank approximations with the singular value decomposition of the true solution, the dimensions of the considered matrix equations are only medium-sized such that a fast and reliable computation of the exact solution can actually be obtained by the lyap function from the MATLAB Control System Toolbox.
At this point, keep in mind that this comparison allows to point out the actual difference between our approximations and the best rank-ˆn approximation which is given by the SVD.
Energy norm for the Lyapunov equation
The first example is a semi-discretized heat transfer problem from the Oberwolfach benchmark collection2. Here, we use the coarsest discretization leading to symmetric matrices E,A ∈ R1357×1357 together with the input matrix B ∈ R1357×7, see [30]. In Figure3.1, we present a comparison between the SVD-based approximation of the exact solution with the approximation given by the rational Krylov subspace obtained by IRKA. As expected from Theorem3.3.1, the relative error in the Frobenius norm is better when the SVD approximation is used. However, the slope of the IRKA approximation is almost parallel to that and for an approximation of rank 50, the relative error (≈8·10−5) is almost as good as the best approximation (≈5·10−5) given by the SVD. On the other hand, we see that IRKA outperforms the SVD for every rank ˆn when the relative error is measured in terms of the L-norm. Although we have discussed a possible way of minimizing the residual, for the symmetric case we do not include a comparison here.
As it has already been obtained in [125], the residual norm is a worse estimator of the true error than the energy norm. Moreover, the corresponding minimizers may be suboptimal and, thus, for symmetric state space systems it does not seem to make sense to minimize the residual at all. Note that the phenomenon of IRKA producing nearly optimal low
2http://portal.uni-freiburg.de/imteksimulation/downloads/benchmark
3.4. Numerical examples 65
10 20 30 40 50
10−6 10−4 10−2 100
Rank ofPnˆ
||P−Pˆn||F/||P||F SVD IRKA
10 20 30 40 50
10−6 10−4 10−2 100
Rank ofPnˆ
||P−Pˆn||L/||P||L
SVD IRKA
Figure 3.2: Tunable optical filter with n = 1668.
rank approximations has already been numerically observed and investigated in [45,59].
The second example is also quite common in the context of model order reduction. The symmetric system matrices E,A ∈ R1668×1668 and B ∈ R1668×5 stem from the finite element discretization of a thermal model of a filter device and thus are sparse, see [79].
Similar to the previous example, from Figure 3.2 we can again conclude that the SVD approximation dominates the performance with respect to the Frobenius norm while the IRKA approach performs better when the L-norm is taken as a basis for judgment.
Moreover, starting from values ˆn = 10, the error resulting from IRKA approximations follows the slope of the error of the SVD-based approximation showing that the subspaces seem to be very close to the optimal ones. Interestingly enough, for dimensions ˆn <10, the energy norm does not seem to be a reasonable error estimator. In particular, although the IRKA approximations outperform the SVD with respect to the energy norm, the error almost stagnates with respect to the Frobenius norm. However, since the norms are not equivalent, we cannot expect an exact one-to-one correspondence between the errors.
Residual norm for the Lyapunov equation
Next, we consider examples that result in unsymmetric dynamical systems with complex system poles that do not allow for a definition of an energy norm. Nevertheless, as we have seen, we can still try to locally minimize the residual R = APˆn+PnˆAT +BBT by the iteration specified in Algorithm 3.3.1. The first example was introduced in [110]
and is one of the SLICOT benchmarks2. The transfer function exhibits three peaks corresponding to six complex system poles while the rest of the poles is completely real.
In Figure 3.3, we compare the results given by the SVD of the true solution with the
2http://www.slicot.org/index.php?site=benchmodred
66 Linear Systems
5 10 15 20 25 30
10−16 10−11 10−6 10−1
Rank of Pˆn
||P−Pˆn||F/||P||F SVD MinRes
5 10 15 20 25 3010−14 10−9 10−4 101
Rank of Pˆn
||R||F/||BBT ||F
SVD MinRes
Figure 3.3: Fom model with n= 1006.
low rank approximations obtained by Algorithm 3.3.1 and abbreviated with MinRes.
Similar to the symmetric case, we see that the latter approach follows the error slope resulting from the SVD. However, there seems to occur some stagnation for larger ranks ˆ
n. On the other hand, we see that our method indeed yields residuals that are smaller than those we get from the SVD. Again, for larger ranks, the approximations seem to be suboptimal and are outperformed by the SVD. Recall Remark 3.3.4 and Remark 3.3.5.
Since we are essentially working with the normal equations, the condition number of the problem is squared and we thus cannot expect the results to be as accurate as in the symmetric case.
The second example is the CD player model we have discussed in Chapter1. The system matrices A and B are part of the SLICOT benchmark collection and are of dimension n = 120 with m = 2 inputs. Again, the matrix A is unsymmetric and exhibits a complex spectrum. In Figure 3.4, we compute the results for low rank approximations varying from ˆn = 1, . . . ,20. While for ranks up to ˆn = 15, we do not observe problems with convergence of Algorithm 3.3.1, the convergence criterion is not fulfilled for larger approximations Pˆn. However, we still obtain the same results as in the previous case.
While the new method performs worse than the SVD when the Frobenius norm is used, we obtain smaller residuals, even if the algorithm does not converge at all. Nevertheless, note that for ˆn= 8 and ˆn = 10, the residuals are larger than those from the SVD. This might indicate that minimizing the residual is not as robust as a minimization of the energy norm as has already been pointed out in [125].
Energy norm for the Sylvester equation
Let us now briefly draw our attention to the more general case specified in (3.47). Analog to the Lyapunov case, we consider an example given by the process of optimal cooling
3.4. Numerical examples 67
5 10 15 20
10−6 10−4 10−2 100
Rank of Pˆn
||P−Pˆn||F/||P||F SVD MinRes
5 10 15 20
10−3 10−1 101
Rank of Pnˆ
||R||F/||BBT ||F
SVD MinRes
Figure 3.4: CD player with n= 120.
10 20 30 40 50 60
10−11 10−8 10−5 10−2
Rank of Xnˆ
||X−Xˆn||F/||X||F SVD (Sy)2IRKA
10 20 30 40 50 60
10−7 10−5 10−3 10−1
Rank ofXˆn
||X−Xˆn||LS/||X||LS
SVD (Sy)2IRKA
Figure 3.5: Steel profiles with discretizations n = 5177 andn = 1357.
of steel profiles. In order to end up with a generalized Sylvester equation including different matrix dimensions, we use a matrix set (A,E,M,H,B,C), where A,M,B are as specified for the Lyapunov case, while E,H,Care obtained by a finer resolution with mesh size m = 5177. In Figure 3.5, we see a comparison between the rational Krylov subspace approximation computed by Algorithm3.3.2 which is abbreviated with (Sy)2IRKA and the SVD-based approximation. Due to the lack of an exact solver for the original Sylvester equation, we use our new method with an approximation of rank 250 for reference values. It should be mentioned that the relative residual for this approach is smaller than 10−13 and thus should be sufficient for comparison. Again, we see that (Sy)2IRKA is dominated by the SVD approximation if the quality is measured in terms of the Frobenius norm while it performs constantly better if we use the LS-norm.
68 Linear Systems