Compensated Splitting For Generalized Lyapunov Equations

Technical Summary: Compensated Splitting for Generalized Lyapunov Equations

Problem Statement
The paper addresses the numerical solution of the generalized Lyapunov equation:
$$AX + XA^H + \sum_{j=1}^m N_j X N_j^H + BB^H = 0,$$
where $A, N_j \in \mathbb{C}^{n \times n}$ and $B \in \mathbb{C}^{n \times \ell}$ (with $\ell \ll n$). This equation arises in model order reduction for bilinear and linear parameter-varying systems, stability analysis of linear stochastic differential equations, and the linearization of stochastic continuous-time algebraic Riccati equations.

The standard approach involves a natural splitting of the linear operator into $M(X) = AX + XA^H$ and $N(X) = -\sum_{j=1}^m N_j X N_j^H$, leading to the standard fixed-point iteration (sFPI). The convergence of sFPI is guaranteed if and only if the spectral radius $\rho(M^{-1}N) < 1$. When $\rho(M^{-1}N) \ge 1$, the standard iteration diverges. Furthermore, even when convergent, the rate may be prohibitively slow if the spectral radius is close to 1. Direct solvers are infeasible for large $n$ due to $O(n^6)$ complexity, necessitating iterative methods.

Methodology
The authors propose a compensated splitting scheme to modify the iteration operator, aiming to reduce the spectral radius of the resulting iteration matrix. The core idea is to decompose the "heavy" $N$-part into a compensating term and a residual term:
$$-N(X) = (EX + XE^H) - \tilde{N}(X),$$
where $E \in \mathbb{C}^{n \times n}$ is a matrix to be determined. This leads to a new splitting:
$$\tilde{M}(X) = M(X) + EX + XE^H = (A+E)X + X(A+E)^H,$$
$$\tilde{N}(X) = \tilde{N}(X) - BB^H.$$
The resulting compensated fixed-point iteration (cFPI) solves $\tilde{M}(X_k) = \tilde{N}(X_{k-1}) - BB^H$ at each step.

To determine the optimal $E$, the authors formulate an optimization problem that minimizes the Frobenius norm of the difference between the original $N$-operator and the compensating term, subject to a regularization term to prevent $E$ from becoming too large:
$$\min_{E} \left( \left\| \sum_{j=1}^m N_j X N_j^H - (EX + XE^H) \right\|_F^2 + \mu \|EX + XE^H\|_F^2 \right),$$
where $\mu \ge 0$ is a regularization parameter.

Key Contributions and Theoretical Results

1. Derivation of Optimal $E$: The paper derives a closed-form solution for the optimal $E$ that minimizes the relaxed objective function. The solution is given by:
   $$E = \frac{1}{n(1+\mu)} \sum_{j=1}^m \text{tr}(\bar{N}_j) \left[ N_j - \frac{1}{2n}\text{tr}(N_j) \cdot I_n \right].$$
   The authors prove that this minimizer is unique in the real case and is the unique minimizer with the smallest Frobenius norm in the complex case.

2. Algorithmic Implementation: The proposed cFPI is implemented in Algorithm 2.1. A key computational advantage is that all inner Lyapunov equations share the same coefficient matrix $(A+E)$, allowing the Schur decomposition of $(A+E)$ to be computed once outside the iteration loop.

3. Theoretical Analysis of Convergence: The authors analyze the impact of the eigenvalue distribution of $A$ on the effectiveness of the compensation. They note that if $A$ has eigenvalues in the right-half plane, the compensation $E$ (which generally shifts eigenvalues rightward) moves them further away from the imaginary axis, potentially improving the invertibility of $\tilde{M}$. Conversely, if $A$ is in the left-half plane, the shift may move eigenvalues closer to zero, potentially worsening the condition of $\tilde{M}$, though the reduction in the spectral radius of the iteration operator may still yield net benefits.

Numerical Results
The paper presents three numerical experiments (two real, one complex) comparing cFPI against sFPI across varying parameters $\xi$ that control the magnitude of $N_j$.

• Example 1 & 2 (Real): In cases where sFPI diverges (spectral radius $\ge 1$), cFPI successfully converges. When sFPI converges, cFPI consistently requires fewer iterations and less CPU time.
• Eigenvalue Sensitivity: The experiments confirm the theoretical observation regarding the location of eigenvalues. When $A$ is replaced by $-A$ (shifting eigenvalues to the right-half plane), the compensated splitting yields a significantly smaller spectral radius compared to the standard splitting, leading to much faster convergence.
• Complex Case (Example 3): The method remains effective for complex matrices, demonstrating similar superiority to sFPI, including the ability to converge where sFPI diverges.
• Regularization Parameter: The authors note that in their experiments, setting $\mu = 0$ yielded the best results, though they acknowledge that a systematic method for choosing $\mu$ remains an open question.

Significance and Claims
The paper claims that the compensated splitting technique offers a practical remedy for the divergence of standard fixed-point iterations in generalized Lyapunov equations. Its primary significance lies in two areas:

1. Rescuing Divergence: It allows the fixed-point iteration to converge in scenarios where the natural splitting fails ($\rho(M^{-1}N) \ge 1$).
2. Acceleration: It accelerates convergence when the standard method is already convergent but slow.

The authors explicitly state that they do not combine this idea with advanced subspace projection techniques (like Krylov methods) in this work, focusing solely on the fixed-point iteration framework. They suggest that combining their compensation strategy with state-of-the-art subspace methods is a direction for future research. The paper concludes modestly, noting that while the current derivation of $E$ does not explicitly account for the $M$-part, numerical evidence suggests the method is robust and effective, particularly when the eigenvalues of $A$ lie in the right-half plane.