A Low-Rank BUG Method for Sylvester-Type Equations

TL;DR

The paper tackles solving Sylvester-type equations $A X + X B^T = C$ by a low-rank BUG-inspired approach that decomposes the problem into smaller reduced tasks, enabling efficient computation when the solution is low-rank or when coefficients are sparse. It introduces fixed-rank and rank-adaptive BUG Sylvester solvers with $X=USV^T$, derives reduced equations for $K=US$, $L=VS^T$, and a small Galerkin system for $S$, achieving a complexity of $O( k r (n^2+m^2+mn+r^2) )$. The methodology is extended to Tucker decompositions for tensor Sylvester equations, solving per-mode reduced problems and employing a rank-adaptive truncation to update the core and factor matrices without forming dense intermediates. Numerical experiments on Poisson problems and random matrices/tensors demonstrate rapid convergence and accurate low-rank approximations, especially when coefficient matrices exhibit larger spectral separation. While empirical results are promising, a rigorous convergence theory remains an open direction for future work.

Abstract

We introduce a low-rank algorithm inspired by the Basis-Update and Galerkin (BUG) integrator to efficiently approximate solutions to Sylvester-type equations. The algorithm can exploit both the low-rank structure of the solution as well as any sparsity present to reduce computational complexity. Even when a standard dense solver, such as the Bartels-Stewart algorithm, is used for the reduced Sylvester equations generated by our approach, the overall computational complexity for constructing and solving the associated linear systems reduces to O(kr(n^2+m^2 +mn + r^2)), for X in R^{m \times n}, where k is the number of iterations and r the rank of the approximation.

Figures (5)

• Figure 1: Left figure: Convergence of the rank-adaptive algorithm for matrices. Middle figure: The singular values of the solution $X$ given by a direct numerical solution, over its largest singular value, alongside the singular values of the low-rank approximation and the threshold $\vartheta = 10^{-10}\|\Sigma_Y\|$. Right figure: The distance between the singular values of the exact solution to the low-rank approximation.
• Figure 2: Left figure: Convergence of the rank-adaptive algorithm for matrices. Middle figure: The singular values of the solution $X$ given by a direct numerical solution, over its largest singular value, alongside the singular values of the low-rank approximation and the threshold $\vartheta = 10^{-10}\|\Sigma_Y\|_F$. Right figure: The distance between the singular values of the exact solution to the low-rank approximation.
• Figure 3: Convergence of the rank-adaptive method for Poisson problem with Dirichlet boundary conditions. In both cases, the truncation tolerance was set to $\vartheta = 10^{-10}\|C_Y\|_F$.
• Figure 4: Left figure: Convergence of the rank-adaptive algorithm for matrices. Middle figure: The singular values of the solution $X$ given by a direct numerical solution, over its largest singular value, alongside the singular values of the low-rank approximation and the threshold $\vartheta = 10^{-10}\|\Sigma_X\|$. Right figure: The distance between the singular values of the exact solution to the low-rank approximation.
• Figure 5: Convergence of the rank-adaptive algorithm for Tucker tensors with random coefficient matrices of prescribed spectra. The truncation tolerance was set to $\vartheta = 10^{-10}\|C_Y\|_F$.