Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

Mixed-precision algorithms for solving the Sylvester matrix equation

Andrii Dmytryshyn, Massimiliano Fasi, Nicholas J. Higham, Xiaobo Liu

TL;DR

This work develops two mixed-precision strategies for solving the Sylvester equation $A X + X B = C$ by computing Schur factors in a low-precision format and refining to high precision. The methods rely on an iterative refinement core for perturbed quasi-triangular Sylvester equations, and stabilize the solution by either orthonormalizing unitary factors in high precision or by explicitly inverting near-unitary factors. The authors provide rigorous error and cost analyses, including a flop-based model to compare against the traditional Bartels–Stewart algorithm, and demonstrate through extensive numerical experiments that the new approaches achieve accuracy comparable to high-precision methods and can offer speedups on hardware with native low-precision capabilities. The results highlight practical guidelines for selecting precision combinations and refinement counts, and suggest promising directions for extending the framework to additional precision tiers and related matrix equations.

Abstract

We consider the solution of the Sylvester equation $AX+XB=C$ in mixed precision. We derive a new iterative refinement scheme to solve perturbed quasi-triangular Sylvester equations; our rounding error analysis provides sufficient conditions for convergence and a bound on the attainable relative residual. We leverage this iterative scheme to solve the general Sylvester equation. The new algorithms compute the Schur decomposition of the coefficient matrices $A$ and $B$ in lower than working precision, use the low-precision Schur factors to obtain an approximate solution to the perturbed quasi-triangular equation, and iteratively refine it to obtain a working-precision solution. In order to solve the original equation to working precision, the unitary Schur factors of the coefficient matrices must be unitary to working precision, but this is not the case if the Schur decomposition is computed in low precision. We propose two effective approaches to address this: one is based on re-orthonormalization in working precision, and the other on explicit inversion of the almost-unitary factors. The two mixed-precision algorithms thus obtained are tested on various Sylvester and Lyapunov equations from the literature. Our numerical experiments show that, for both types of equations, the new algorithms are at least as accurate as existing ones. Our cost analysis, on the other hand, suggests that they would typically be faster than mono-precision alternatives if implemented on hardware that natively supports low precision.

Mixed-precision algorithms for solving the Sylvester matrix equation

TL;DR

This work develops two mixed-precision strategies for solving the Sylvester equation by computing Schur factors in a low-precision format and refining to high precision. The methods rely on an iterative refinement core for perturbed quasi-triangular Sylvester equations, and stabilize the solution by either orthonormalizing unitary factors in high precision or by explicitly inverting near-unitary factors. The authors provide rigorous error and cost analyses, including a flop-based model to compare against the traditional Bartels–Stewart algorithm, and demonstrate through extensive numerical experiments that the new approaches achieve accuracy comparable to high-precision methods and can offer speedups on hardware with native low-precision capabilities. The results highlight practical guidelines for selecting precision combinations and refinement counts, and suggest promising directions for extending the framework to additional precision tiers and related matrix equations.

Abstract

We consider the solution of the Sylvester equation in mixed precision. We derive a new iterative refinement scheme to solve perturbed quasi-triangular Sylvester equations; our rounding error analysis provides sufficient conditions for convergence and a bound on the attainable relative residual. We leverage this iterative scheme to solve the general Sylvester equation. The new algorithms compute the Schur decomposition of the coefficient matrices and in lower than working precision, use the low-precision Schur factors to obtain an approximate solution to the perturbed quasi-triangular equation, and iteratively refine it to obtain a working-precision solution. In order to solve the original equation to working precision, the unitary Schur factors of the coefficient matrices must be unitary to working precision, but this is not the case if the Schur decomposition is computed in low precision. We propose two effective approaches to address this: one is based on re-orthonormalization in working precision, and the other on explicit inversion of the almost-unitary factors. The two mixed-precision algorithms thus obtained are tested on various Sylvester and Lyapunov equations from the literature. Our numerical experiments show that, for both types of equations, the new algorithms are at least as accurate as existing ones. Our cost analysis, on the other hand, suggests that they would typically be faster than mono-precision alternatives if implemented on hardware that natively supports low precision.

Paper Structure

This paper contains 23 sections, 2 theorems, 83 equations, 4 figures, 4 tables.

Table of Contents

  1. Introduction
  2. Background
  3. Error analysis
  4. Numerical algorithms
  5. The Bartels--Stewart algorithm
  6. Hermitian matrix coefficients
  7. Iterative refinement for perturbed quasi-triangular equations
  8. Computational cost and storage requirement
  9. Error analysis and attainable residual
  10. Normwise forward error analysis
  11. Normwise backward error analysis
  12. Mixed-precision algorithms for the Sylvester equation
  13. Orthonormalization of the unitary factors in high precision
  14. Computational cost and storage requirement
  15. A flop-based computational model
  16. ...and 8 more sections

Key Result

Theorem 3.1

Let alg:ir-linear-system be applied to a linear system $Mx = b$ in precision $u_h$, where $M\in\mathbb{C}^{s\times s}$ is nonsingular, and assume that the solver used on ln:linear-solve satisfies eq:assum-di and eq:assum-pertb. If is sufficiently smaller than $2$, then the normwise forward error is reduced on the $i$th iteration by a factor approximately $\phi_i$ until an iterate $\widehat{x}$ is

Figures (4)

  • Figure 4.1: Maximum number of iterations for which the iterative refinement algorithms in Section \ref{['sec:pert-sylv-algs']} will be asymptotically faster than the Bartels--Stewart algorithm run in high precision. We are plotting the quantities in \ref{['eq:funk-orth']}, \ref{['eq:funk-orth-lyap']}, \ref{['eq:funk-inv']}, and \ref{['eq:funk-inv-lyap']} against the ratio $\rho$ of the computational cost of low-precision to a high-precision operation.
  • Figure 6.1: Comparison of lyap and our MATLAB implementation of \ref{['alg:sylv-mp', 'alg:sylv-mp-inversion']} on matrix equations from the literature. The low-precision arithmetic is TensorFloat-32, for which $u_\ell = 2^{-11}$. The top and bottom row refer to Sylvester and Lyapunov equations, respectively. Left: relative residual of the computed solution. Right: number of iterative refinement steps.
  • Figure 6.2: Comparison of lyap and our MATLAB implementation of \ref{['alg:sylv-mp', 'alg:sylv-mp-inversion']} on the same test set. The low-precision arithmetic is a custom 3-byte format, for which $u_\ell = 2^{-16}$.
  • Figure A.1: Residual of solutions computed using the MATLAB function sylvester, our MATLAB implementation of GMRES-IR2, and our MATLAB implementations of \ref{['alg:sylv-mp', 'alg:sylv-mp-inversion']}. Here, $\kappa_{2}(M_f)\approx 10^{t}$, and we set $u_\ell = 2^{-24}$ (binary32) and $u_h = 2^{-53}$ (binary64).

Theorems & Definitions (2)

  • Theorem 3.1
  • Theorem 3.2