Sylvester Equation Reformulation

• Sylvester-Equation Reformulation is a method that transforms the standard Sylvester equation into a matrix function framework using block matrices and sign functions, enabling solution without direct inversion.
• It employs advanced iterative schemes such as Akhiezer polynomial and Zolotarev rational approximations to achieve rapid convergence and precise error control.
• The approach is well-suited for high-dimensional problems, leveraging low-rank structures and parallel computations to handle large-scale and sparse systems efficiently.

A Sylvester-equation reformulation is any transformation or reinterpretation of a linear, matrix-valued operator equation such as $X\,A + B\,X = C$ (or, equivalently, $A\,X + X\,B = C$, depending on conventions) that expresses the existence, uniqueness, and computation of $X$ in terms of alternative matrix-analytic or operator-theoretic objects. These reformulations underpin much of the modern numerical and analytic machinery for fast and robust solution of linear matrix equations, particularly when the problem size or structure precludes standard direct methods. Recent developments center around sign-function block formulations, optimal polynomial/rational iterative schemes, and low-rank exploitation. The following sections provide a comprehensive account of these strategies, focusing on the case where $A$ and $-B$ have spectra in disjoint real intervals and the solution can be constructed without inverting $A$ or $B$ directly (Ballew et al., 21 Mar 2025).

1. Block-Matrix and Matrix-Sign Function Formulation

Central to recent advances in Sylvester-equation reformulation is the block-matrix embedding technique. Given the Sylvester equation $X\,A+B\,X=C$, introduce the $(m+n)\times(m+n)$ block matrix

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

If $\sigma(A) \subset U_A$ and $\sigma(-B) \subset U_B$ for disjoint compact real intervals $U_A, U_B \subset \mathbb{R}$, the scalar sign function is defined as

$$\operatorname{sign}(z) = \begin{cases} +1,& z\in U_A,\ -1,& z\in U_B. \end{cases}$$

Under these conditions, $\operatorname{sign}(M)$ can be block-factorized as

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

so that $X = \frac12\, [\operatorname{sign}(M)]_{2,1}$ is the Sylvester solution. This equivalence (cf. Roberts, Higham) recasts the problem into the context of matrix functions, opening the way to polynomial and rational approximation methods acting on $M$ rather than inverting $A$ or $B$ (Ballew et al., 21 Mar 2025).

2. Akhiezer Polynomial Iteration

To approximate the matrix sign function, one employs Akhiezer-type orthogonal polynomials $p_0,p_1,\ldots$ on $\Sigma=U_B\cup U_A,$ with weight $w(x)$. These satisfy a three-term recurrence: $$x p_k(x) = b_{k-1} p_{k-1}(x) + a_k p_k(x) + b_k p_{k+1}(x),$$ with initial $p_0(x)=1$. The discontinuous sign function is expanded in this basis as

$$\operatorname{sign}(z) \approx \sum_{j=0}^n \alpha_j p_j(z), \quad \alpha_j = \int_\Sigma p_j(x)\, \operatorname{sign}(x)\, w(x)\, dx.$$

For each $n$, the matrix expansion $\sum \alpha_j p_j(M)$ approximates $\operatorname{sign}(M)$. Importantly, only selected blocks of $p_j(M)$ are needed: by compact coupled recursions using only $A,B,C$, one efficiently constructs the sequence: $$X_{k+1} = X_k + \frac{\alpha_k}{2}\left(C\, p_k(A) + G_k\right),$$ where $G_0=-C$, $G_1 = \frac{G_0 A + (a_0+1) C}{b_0}$, and for $j\ge2$,

$$G_j = \frac{1}{b_{j-1}}\big(G_{j-1}A + p_{j-1}(B)C - a_{j-1} G_{j-1} - b_{j-2} G_{j-2}\big).$$

This construct is inverse-free and avoids forming full block matrices (Ballew et al., 21 Mar 2025).

3. Rational-Approximation and Direct-Inverse Reformulation

The alternative to polynomial expansion is direct rational approximation for $L^{-1}$, where $L(X)=X A + B X$. The optimal min–max rational approximant $r_n(z)\approx 1/z$ on two intervals $\Sigma_L \supset \sigma(L)$ solves a Zolotarev problem, and is of the form: $$r_n(z) = z \prod_{j=1}^n \frac{z^2+\mu_j^2}{z^2+\nu_j^2},$$ with the coefficients $\mu_j, \nu_j$ determined from elliptic integrals. In partial-fraction form,

$r_n(z) = \sum_{j=1}^n \frac{\alpha_j}{z+\beta_j},$

so applying $r_n(L)$ to a matrix $W$ reduces to solving $n$ shifted Sylvester equations: $$X^{(j)} A + B X^{(j)} + \beta_j X^{(j)} = W,\quad j=1,\ldots,n,$$ either in parallel or sequentially. Iterative refinement

$X_{k+1} = X_k + r_n(L)\left( C - L(X_k) \right)$

reduces the Sylvester-residual by a computable geometric rate per cycle (Ballew et al., 21 Mar 2025).

4. Convergence Theory and Error Bounds

Both the polynomial (Akhiezer) and rational (Zolotarev-based) iterations exhibit geometric convergence, with explicit rates governed by the spectral gap between $\sigma(A)$ and $-\sigma(B)$. For the sign-based Akhiezer method, the expansion coefficients $\alpha_j$ admit the sharp Bernstein-type bound: $$|\alpha_j| \leq C\, \varrho^{-j}, \qquad \varrho = \exp(c^*)>1,$$ where $c^*=\max_{z\in$(gap)$}\mathrm{Re}\,g(z)$, and $g(z)$ is the Green's function for $\mathbb{C} \setminus \Sigma$. The sign-series error is thus

$$\|\operatorname{sign}(M) - F_k\|_2 \leq C' \varrho^{-k}/(1-\varrho^{-1}),$$

and for the Sylvester solution

$$\|X - X_k\|_2 \leq \frac{D}{2} \frac{\varrho^{-k}}{1-\varrho^{-1}},$$

with $D$ a condition-number dependent prefactor. For the rational method, the per-iteration contraction is $\epsilon_n = \max_{z\in \sigma(L)} |1-zr_n(z)|$, with the optimal Zolotarev value decaying like $\epsilon_n \approx 4\exp(-\pi^2 n/\log(4b/a))$, $[a,b]$ denoting interval endpoints.

5. Algorithmic Implementation and Computational Complexity

Akhiezer Polynomial Solver

For full-rank data, each iteration costs $O(n^3 + m^3)$ if $A,B$ are dense, with matrix-matrix multiplies dominating; for sparse or banded matrices, this reduces accordingly. With low-rank right-hand side $C=U V$ ($U\in \mathbb{R}^{m\times r}$, $V\in \mathbb{R}^{r\times n}$), all iterates can be maintained in compressed factored ($W Q$) form, with each step primarily but not exclusively $O(r \cdot (n^2+m^2))$, and low-rank truncations reduce storage and arithmetic costs to $O((r+k)^3)$ per iteration (for $k$ iterations).

Rational/Zolotarev Method

Each cycle requires $n$ shifted Sylvester solves, each direct or using a precomputed factorization of $A$,$B$. The overall complexity is $O(n \cdot \text{cost per Sylvester solve})$ per outer update (Ballew et al., 21 Mar 2025).

6. Practical Guidelines on Interval Selection and Degree Choice

Robust application depends on precise enclosures for $\sigma(A)$ and $-\sigma(B)$. One selects

$$[\beta_1, \gamma_1] \supset \sigma(-B), \quad [\beta_2, \gamma_2] \supset \sigma(A),$$

expanding as needed for rounding/uncertainty. The required polynomial or rational degree for error $\epsilon$ is

$$n \geq \left\lceil -\frac{\log\,[\,\epsilon(1-\varrho^{-1})/D\,]}{\log\,\varrho} \right\rceil, \quad \text{for Akhiezer;} \quad n \geq \frac{\log(4\gamma/\beta)}{\pi^2} \log\frac{4}{\epsilon}, \quad \text{for Zolotarev}.$$

Choice of method depends on the matrix structure: for structured or low-rank-compatible $A,B$, the Akhiezer approach is efficient; for settings with factorized or diagonalizable $A,B$, Zolotarev/ADI or direct inversion can be competitive.

7. Applicability and Limitations

These reformulations are effective when $\sigma(A)$ and $\sigma(-B)$ are real and separated, as in many discretized PDEs and control systems. When spectra are complex or coalescent, alternative spectral splitting or Hamiltonian/reduction strategies are needed. Low-rank structure in $C$ directly accelerates the Akhiezer and rational schemes. These iterations yield explicit geometric rates, are inverse-free, and are well suited for parallelization and structure exploitation in large or sparse systems. Implementation requires only three-term recurrence arithmetic and standard polynomial/rational weight computation; further details and generalizations can be found in (Ballew et al., 21 Mar 2025).