Mixed-precision algorithms for solving the Sylvester matrix equation

1. Problem Statement

The paper addresses the numerical solution of the Sylvester matrix equation:
$$AX + XB = C$$
where $A \in \mathbb{C}^{m \times m}$, $B \in \mathbb{C}^{n \times n}$, and $C \in \mathbb{C}^{m \times n}$ are given, and $X$ is the unknown. A special case is the continuous Lyapunov equation ($AX + XA^* = C$).

Context:

• These equations are fundamental in control theory, model reduction, and signal processing.
• The standard solution method is the Bartels–Stewart algorithm, which involves computing Schur decompositions of $A$ and $B$, transforming the equation into a triangular system, and solving it via back-substitution.
• Motivation: With the advent of hardware supporting mixed-precision arithmetic (e.g., NVIDIA Tensor Cores), there is a need to leverage lower precision (e.g., TensorFloat-32, bfloat16) to accelerate computations while maintaining the accuracy of high precision (e.g., binary64).
• Challenge: Directly applying mixed-precision to the Bartels–Stewart algorithm fails because the unitary factors ($U_A, U_B$) computed in low precision are not unitary to high precision. Using them to recover the solution introduces errors that prevent the final result from reaching working precision.

2. Methodology

The authors propose a framework combining low-precision Schur decomposition with iterative refinement in high precision. The core strategy involves three stages:

A. Low-Precision Preprocessing

1. Compute the Schur decompositions of $A$ and $B$ in low precision ($u_\ell$):
   $$A \approx \hat{U}_A \hat{T}_A \hat{U}_A^*, \quad B \approx \hat{U}_B \hat{T}_B \hat{U}_B^*$$
   Here, $\hat{T}_A, \hat{T}_B$ are upper quasi-triangular, and $\hat{U}_A, \hat{U}_B$ are "almost-unitary" (unitary only to precision $u_\ell$).

B. Iterative Refinement for Perturbed Triangular Systems

The authors derive a new stationary iteration scheme (Algorithm 3.1) to solve the perturbed quasi-triangular Sylvester equation:
$$(\hat{T}_A + \Delta T_A)Y + Y(\hat{T}_B + \Delta T_B) = \tilde{C}$$

• This scheme treats the low-precision errors ($\Delta T$) as perturbations.
• Convergence Analysis: They prove that the iteration converges if the perturbation magnitude is smaller than the separation of the spectra ($sep(\hat{T}_A, -\hat{T}_B)$).
• Error Bounds: They establish that the attainable relative residual is bounded by the unit roundoff of the high precision ($u_h$), provided the condition number of the problem is not excessively large relative to $u_\ell$.

C. High-Precision Recovery Strategies

To recover the solution $X$ in high precision ($u_h$) from the low-precision factors, the authors propose two distinct algorithms to handle the non-unitary nature of $\hat{U}_A$ and $\hat{U}_B$:

1. Algorithm 4.1: Re-orthonormalization

   • Compute QR factorizations of the low-precision unitary factors in high precision: $\hat{U}_A = Q_A R_A$, $\hat{U}_B = Q_B R_B$.
   • Use the high-precision orthogonal matrices $Q_A, Q_B$ to transform the original equation.
   • Solve the resulting perturbed triangular system using the iterative refinement scheme.
   • Cost: Requires explicit QR factorization (approx. $2(m^3+n^3)$ flops in high precision).

2. Algorithm 4.2: Explicit Inversion

   • Instead of re-orthonormalizing, explicitly invert the "almost-unitary" factors $\hat{U}_A$ and $\hat{U}_B$ in high precision (using LU decomposition).
   • Reformulate the equation to incorporate the inversion errors as perturbations.
   • Solve the perturbed system via iterative refinement.
   • Cost: Requires LU decomposition (approx. $\frac{2}{3}(m^3+n^3)$ flops in high precision), making it computationally cheaper than Algorithm 4.1.

3. Key Contributions

• New Iterative Refinement Scheme: A rigorous derivation of a stationary iteration for solving perturbed quasi-triangular Sylvester equations, including convergence conditions and bounds on the attainable residual.
• Two Mixed-Precision Algorithms:
  • MP-Orth (Alg 4.1): Uses re-orthonormalization to ensure high-precision unitarity.
  • MP-Inv (Alg 4.2): Uses explicit inversion of low-precision factors, offering lower computational overhead.
• Theoretical Error Analysis: Detailed forward and backward error analysis showing that the algorithms can achieve backward stability in high precision ($u_h$) despite using low precision ($u_\ell$) for the Schur decomposition, provided the condition number $\kappa(M)$ satisfies $\kappa(M) \lesssim 1/u_\ell$.
• Cost Modeling: A flop-based model comparing the mixed-precision approaches against the standard high-precision Bartels–Stewart algorithm. The model introduces a ratio $\rho$ (cost of low-precision flop / cost of high-precision flop) to determine break-even points.

4. Results

The authors tested their algorithms on 19 Sylvester and 12 Lyapunov equations from literature using MATLAB and simulated hardware.

• Accuracy:

  • Both mixed-precision algorithms achieved relative residuals comparable to the standard high-precision lyap/sylvester functions (often reaching $10^{-14}$ to $10^{-16}$).
  • Failure Cases: The algorithms failed on singular problems (where spectra overlap) or ill-conditioned problems where the low-precision condition number exceeded the threshold ($10^4$ for TensorFloat-32).
  • Robustness: Using a custom 16-bit significand format (more accurate than TensorFloat-32) resolved instability issues seen with TensorFloat-32 on ill-conditioned problems, reducing the required iterations to 2–3.

• Performance & Efficiency:

  • Sylvester Equations: The mixed-precision algorithms are theoretically faster if the low-precision hardware is significantly faster ($\rho \ll 1$). For TensorFloat-32 on modern GPUs (speedup $\sim 30\times$ to $900\times$), the mixed-precision approach can be 20% to 60% faster than the standard algorithm, even with 1–5 refinement iterations.
  • Lyapunov Equations: Due to the structure of the equation (requiring fewer operations in high precision for the standard method), the mixed-precision advantage is smaller. Gains are only realized for well-conditioned problems requiring very few iterations (1–2).
  • Algorithm Comparison: Algorithm 4.2 (Inversion) consistently required fewer high-precision flops than Algorithm 4.1 (Orthonormalization), making it the preferred choice for performance.

5. Significance

• Bridging Hardware and Theory: The paper provides a concrete, mathematically sound method to exploit emerging mixed-precision hardware for dense linear algebra problems beyond simple linear systems.
• Practical Applicability: It demonstrates that one does not need to perform the entire Schur decomposition in high precision to get a high-precision solution, provided a robust iterative refinement step is used.
• Scalability: The approach is particularly relevant for large-scale dense problems where the cubic cost of the Schur decomposition dominates the total runtime.
• Future Directions: The authors highlight that the approach is most beneficial for well-conditioned problems and suggests future work on unbalanced matrices ($m \gg n$) and extending the method to T- and *-Sylvester equations.

In summary, this work establishes that mixed-precision Schur decomposition followed by high-precision iterative refinement is a viable, accurate, and potentially faster alternative to traditional high-precision solvers for Sylvester and Lyapunov equations, provided the problem is not excessively ill-conditioned relative to the low precision used.