Numerical Solution of Large-Scale Low-Rank Coupled Generalized Sylvester Equations

1. Introduction

• Descriptor-system conversion requires solving coupled generalized Sylvester equations to obtain equivalent ODE state-space representations.The paper identifies controller design and numerical simulation of DAE-based models as motivations.Source: §1 · p. 1
• The equations have low-rank right-hand-side structure, and the paper assumes their solutions are numerically low-rank.Rapid singular-value decay is stated to permit accurate low-rank approximations.Source: §1 · p. 1
• The paper assumes nonsingular E1 and E2 and disjoint spectra for E1^-1A1 and E2^-1A2, which ensure a unique solution for every right-hand side.The stated uniqueness consequence is tied to these assumptions.Source: §1 · p. 1

2.1. Petrov–Galerkin Projection Framework

• Projected equations have a unique solution for every right-hand side under the stated condition, while disjoint reduced-pencil spectra give the necessary-and-sufficient criterion.The sufficient condition is not necessary; the exact criterion is disjoint spectra, equivalently coprimeness of the two reduced determinant polynomials.Source: §2.1 · pp. 2–3

2.2. Choice of Projection Matrices

• The solution's SVD-factor ranges are contained in enlarged bases built from resolvent directions, motivating projection spaces that approximate dominant transfer-function dynamics.When the reduced system has a unique solution and the exact factor ranges are contained in the constructed bases, the lifted approximations recover X and Y exactly.Source: §2.2 · pp. 4–5
• Nearly uncontrollable and unobservable poles have vanishingly small residues, explaining the rapidly decaying singular values that produce numerical low-rank solutions.Strictly uncontrollable or unobservable poles have zero residues and can reduce the exact rank.Source: §2.2 · p. 4

2.3. Extended Krylov-based Projection Method

• The extended Krylov method uses fixed interpolation points at s = 0 and s = ∞ and requires A1 and A2 to be invertible.Fixed points avoid shift selection; sparse factorizations are reused, and the method stops when normalized residuals fall below the prescribed tolerance.Source: §2.3 · p. 5

2.4. Rational Krylov-based Projection Method

• The rational Krylov algorithm selects shifts automatically, avoiding user intervention in a method whose performance depends strongly on shift choice.The paper contrasts its flexible shift selection with the fixed points used by extended Krylov methods.Source: §2.4 · p. 5

10 end

• The trial bases are chosen using the structural properties established in Proposition 2.2.The supplied passage does not specify the properties themselves.Source: §10 · p. 6
• The rational Krylov algorithm estimates unknown controllable poles using reduced-pencil Ritz values and selects dominant values for expansion.Dominant Ritz values are selected using the measure used in dominant pole estimation algorithms.Source: §10 · p. 6
• Deflation removes right-hand-side components already captured by current subspaces, while residue normalization retains useful directions when selected values are weakly controllable or observable.Such values often occur near convergence, when unnormalized residue norms may be negligible.Source: §10 · p. 6
• Block expansion uses the full right-hand-side matrices, whereas tangential expansion uses two columns per update.Tangential updates are particularly useful when m is not small, because block updates can rapidly increase projection dimensions and make solving equation (4) computationally expensive or infeasible.Source: §10 · p. 6

3. Numerical Results

• The evaluation uses MATLAB R2025b on a Windows 11 laptop with 32 GB of RAM and an Intel Core Ultra 9 285H processor running at 2.9 GHz.The MATLAB code for reproducing the results is publicly available.Source: §3 · p. 6
• The examples use 20 dominant poles for one pencil and 10 for the other, with the remaining poles designated nondominant.The passage says the examples are constructed so that X and Y are numerically low-rank, but the supplied text cuts off before completing that statement.Source: §3 · p. 6
• The algorithm takes the coupled-equation matrices, a tolerance, and a block-or-tangential setting, and returns low-rank approximations of X and Y.The stopping criteria compare normalized Frobenius residuals for both coupled equations against the tolerance.Source: §3 · p. 7
• In block mode, the update matrices B and C are set to the full perpendicular right-hand-side matrices.The algorithm then solves the projected coupled generalized Sylvester equations to compute the reduced approximations.Source: §3 · p. 7

19 end

• The experiments use a stopping tolerance of τ = 10^-8 and assess accuracy through the maximum normalized residuals of the two equations.The residuals correspond to the equations involving A1, A2 and E1, E2, respectively.Source: §19 · p. 7

3.1. Model Description

• The test models use user-specified real poles and complex-conjugate pole pairs; the first nd pairs are weighted more heavily and are the most controllable and observable.The constructed coefficient matrices are sparse, and the real poles are nearly uncontrollable and unobservable.Source: §3.1 · p. 8

3.2. Example 1

• In Example 1, both methods converged rapidly and produced low-rank approximations with relative differences on the order of 10^-8 from the reference solution.The direct solution served as the accuracy reference.Source: §3.2 · pp. 8–9, Figure 1 · p. 9

3.3. Example 2

• In Example 2, adaptive shifts gave RKPA-C-SYLV lower-rank approximations, but tangential updates slowed convergence.Block RKPA-C-SYLV was fastest, and all three methods converged in less than 15 seconds.Source: §3.3 · p. 9, Table 1 · p. 10

3.4. Example 3

• In Example 3, tangential RKPA-C-SYLV converged with the most compact approximation, while block RKPA-C-SYLV failed and EKPA-C-SYLV converged rapidly.EKPA-C-SYLV was fastest despite less compact approximations; block RKPA-C-SYLV terminated after updates led to linear dependence.Source: §3.4 · p. 10, Table 2 · p. 10
• In Example 3, RKPA-C-SYLV required new factorizations each iteration, whereas EKPA-C-SYLV reused its sparse LU factorizations.The block RKPA-C-SYLV variant took longest and terminated without convergence.Source: §3.4 · p. 10, Table 2 · p. 10

• The dominant subspaces of the exact solutions lie in spaces spanned by specific resolvent directions, providing the theoretical basis for the methods.Source: §4 · p. 11
• Projected coupled generalized Sylvester equations are uniquely solvable under a field-of-values separation condition.Source: §4 · p. 11
• EKPA-C-SYLV uses fixed shifts at s = 0 and s = ∞, while RKPA-C-SYLV adaptively selects shifts from reduced-pencil Ritz values.Source: §4 · p. 11
• Numerical experiments on large-scale problems show that both algorithms compute accurate low-rank solutions efficiently.Source: §4 · p. 11

• The equation (33) has at most one solution because it shares the left-hand Sylvester operator of (8), whose solution V11 is unique.Source: Appendix · p. 12
• The equation (38) has at most one solution because it shares the left-hand Sylvester operator of (10), whose solution W21 is unique.Source: Appendix · p. 14