Sylvester Equation Solver

• Sylvester equation solvers compute the matrix X satisfying AX + XB = C under disjoint spectral conditions for uniqueness.
• They employ inverse-free iterative methods using orthogonal polynomial recurrences to avoid costly matrix inversions.
• The approach guarantees geometric convergence and efficiency in both dense and low-rank settings, with applications in control and PDE discretizations.

A Sylvester equation solver computes the matrix $X$ that satisfies $AX + XB = C$ for given matrices $A$, $B$, and $C$, possibly with structural or spectral constraints. These solvers are critical in control, signal processing, numerical linear algebra, and systems theory. Recent developments extend traditional direct and iterative approaches—such as Bartels–Stewart and ADI—with new paradigms that are inverse-free or exploit structure, spectral localization, low-rank properties, and fast transforms.

1. Problem Setting and Spectral Preconditions

The classical Sylvester equation seeks $X \in \mathbb{C}^{n \times m}$ solving $A X + X B = C$, where $A \in \mathbb{C}^{n \times n}$, $B \in \mathbb{C}^{m \times m}$, $C \in \mathbb{C}^{n \times m}$. The question of existence and uniqueness is governed by the disjointness of spectra: $\sigma(A)\cap(-\sigma(B)) = \emptyset$ is both necessary and sufficient for uniqueness when $A$ and $B$ are diagonalizable.

A key paradigm in (Ballew et al., 21 Mar 2025) restricts to the case where $A$ and $-B$ have real, diagonalizable spectra on disjoint intervals: $\sigma(A)\subset [\beta_2,\gamma_2]$, $\sigma(-B)\subset[\beta_1,\gamma_1]$, and $\gamma_1 < \beta_2$. This structural information is leveraged to recast Sylvester equation solvers as iterations involving orthogonal polynomials on two intervals and to enable inverse-free computation.

2. Akhiezer Iteration Framework and Matrix Function Expansions

The Akhiezer iteration employs orthogonal polynomials $\{p_k\}$ constructed on $$\Sigma = [\beta_1,\gamma_1] \cup [\beta_2,\gamma_2]$$ with a canonical weight:

$$w(x) \propto \begin{cases} \sqrt{(x - \beta_j)(\gamma_j - x)}, & x \in [\beta_j,\gamma_j],\ j=1,2, \ 0, & \text{otherwise} \end{cases}$$

These polynomials satisfy a three-term recurrence:

$$x p_0(x) = a_0 p_0(x) + b_0 p_1(x),\quad x p_k(x) = b_{k-1} p_{k-1}(x) + a_k p_k(x) + b_k p_{k+1}(x),\; k \geq 1,$$

with $b_k > 0$.

Given a matrix $M$ with $\sigma(M)\subset\Sigma$ and $f$ analytic on a neighborhood of $\Sigma$, polynomial approximations $f(M)\approx \sum_{k=0}^{K-1} \alpha_k p_k(M)$ (where $\alpha_k=\int_\Sigma f(x)p_k(x)w(x)dx$) enable iterative computation of matrix functions without direct inversion.

The block matrix

$$H = \begin{pmatrix} A & 0 \ C & -B \end{pmatrix}$$

is the fundamental object: extracting certain block components of functions of $H$ yields the solution of $AX + XB = C$ (via the sign-function or other analytic $f$).

3. Inverse-Free Algorithmic Approaches

Two computational approaches are defined:

3.1 Block Sign-Function Akhiezer Iteration

Based on the observation that the matrix sign-function applied to $H$ yields

$$\operatorname{sign}(H) = \begin{pmatrix} I & 0 \ 2X & -I \end{pmatrix}$$

so that $X = \frac{1}{2} [\operatorname{sign}(H)]_{21}$. The Akhiezer iteration expands

$$\operatorname{sign}(H) \approx \sum_{k=0}^{K-1} \alpha_k p_k(H)$$

with each $p_k(H)$ built recursively. The solution $X_K$ at iteration $K$ is given as $X_K = \frac{1}{2} [F_K]_{21}$, where $F_K$ is the degree-$K$ partial sum. This method only requires forming matrix-matrix products and three-term recurrences; no inverses or solutions of linear systems are required (Ballew et al., 21 Mar 2025).

3.2 Decoupled Direct Sylvester Recurrence

A block structure allows avoiding the explicit $H \in \mathbb{C}^{(n+m)\times(n+m)}$ operations. In this formulation:

• $p_k(H)$ has block form with $p_k(A)$, $p_k(-B)$, and a $(2,1)$ block $G_k$ satisfying its own three-term recurrence:

$$G_0 = -C, \quad G_1 = \frac{1}{b_0}(G_0 A + (a_0+1) C), \quad G_k = \frac{1}{b_{k-1}}\bigl( G_{k-1}A + p_{k-1}(-B) C -a_{k-1}G_{k-1} - b_{k-2}G_{k-2} \bigr)$$

for $k\geq2$. The polynomial recurrences for $p_k(A)$ and $p_k(-B)$ are standard three-term recursions.

The solution is accumulated as

$$X_{k+1} = X_k + \frac{\alpha_k}{2} ( C p_k(A) + G_k )$$

Each step requires only matrix multiplications of matching sizes and is readily adapted for low-rank $C$ (Ballew et al., 21 Mar 2025).

4. Convergence Guarantees and Rate Analysis

Geometric convergence at rate $\varrho^{-1}$ is guaranteed, where $\varrho$ is computed from potential theory based on the two intervals. Explicitly, for the Akhiezer coefficients:

$|\alpha_k| \leq C \varrho^{-k}$

The error in the sign function approximation is (Theorem 4.8 in (Ballew et al., 21 Mar 2025)):

$$\|\operatorname{sign}(H) - F_K\|_2 \leq D_H \frac{\varrho^{-K}}{1 - \varrho^{-1}},$$

implying

$$\|X_K - X\|_2 \leq \frac{1}{2} D_H \frac{\varrho^{-K}}{1 - \varrho^{-1}}.$$

Here, $D_H$ depends on the conditioning of the eigenvector matrix of $H$ and $\varrho$ is associated with a generalized Green's function determined by the spectral gap.

The convergence is thus computable and geometric; the required $K$ can be predicted for a given tolerance.

5. Computational Complexity, Structure Exploitation, and Algorithm Selection

• Block sign-function approach: Each iteration requires $O((n+m)^3)$ operations, with $K=O(\log(n+m))$ for convergence, total complexity $O((n+m)^3\log(n+m))$.
• Direct recurrence approach (dense $C$): $O(n^3 + m^3)$ per step, overall $O((n^3 + m^3)\log(n+m))$.
• Low-rank $C$: If $C=UV$ is rank $r\ll n,m$, then each step costs $O(r(n^2 + m^2))$ (matrix–low-rank-matrix products), and compression via QR+SVD of $O((r+k)^3)$ every $k$ steps, yielding $O(K r (n^2+m^2))$ overall.
• No inverses or triangular solves: Only matrix multiplications and three-term recurrences are performed.

The approaches excel when $A$ and $B$ are dense and costly to invert or when $C$ is low-rank and $A$, $B$ are accessed via fast mat-vec routines or hierarchical representations.

6. Practical Implementation and Parameter Guidance

• For spectral intervals, the $\varrho$ parameter is computed from the maximal value of the Green's function across the gap between spectral intervals.
• $K$ can be set as $$K = \lceil -\log_\varrho(\epsilon(1-\varrho^{-1})/D_H) \rceil$$ for target error tolerance $\epsilon$ and estimated $D_H$ (the paper suggests $D_H\leq 10(n+m)$ as rule of thumb).
• For low-rank right-hand side $C$, perform all steps in low-rank factorized form and compress $G_k$ at tolerance proportional to $\alpha_k$ decay.
• For dense $C$ and moderate $n+m$, use the block scheme.
• For very large $n$, $m$, and $C$ numerically low-rank, the low-rank recurrence variant is superior.

7. Applications, Benchmarks, and Limitations

• Matrix equations in operator-theoretic PDE discretizations: Collocation of certain integral equations (e.g., convolution-type separable operators) yields Sylvester equations with low-rank $C$ and dense $A,B$, where the inverse-free Akhiezer method yields significant speedup over Bartels–Stewart (Ballew et al., 21 Mar 2025).
• Fréchet derivatives of matrix functions: The recurrence and polynomial expansion formalism allows computation of $L_f(A,E)$ for matrix function $f$ by applying recurrences to augmented blocks.
• Benchmarks: For instances with $n\sim 2000$ and $\operatorname{rank} C=1$, the Akhiezer low-rank solver runs in seconds, outperforming ADI and Bartels–Stewart methods.
• Limitations: For banded $A,B$ or cases where fast $(A-\sigma I)^{-1}$ applies (rational Krylov/ADI), or for very small $n$, direct (factorization-based) methods can be superior. The spectral interval disjointness is essential for convergence and for the three-term recurrences to be well-behaved.

---

The Akhiezer iteration framework thus constitutes a rigorously convergent, explicit, inverse-free family of Sylvester equation solvers, with provable geometric rates, especially advantageous for dense or large-scale problems with spectral gaps and/or low-rank structure (Ballew et al., 21 Mar 2025).