A Non-Recursive, Dimension-Independent Schur-Decomposition Algorithm for $N$-Dimensional Sylvester Tensor Equations

TL;DR

This paper introduces a simple, dimension-independent direct solver for N-dimensional Sylvester tensor equations based on Schur decompositions, capable of handling very high-dimensional problems efficiently.

Contribution

The paper presents a non-recursive, dimension-independent algorithm for Sylvester tensor equations, with implementation details and applications to high-dimensional PDEs.

Findings

Successfully solved up to N=29 dimensions on a standard laptop.

Achieves comparable accuracy to state-of-the-art recursive methods.

More memory-efficient and simpler to implement, especially for large N.

Abstract

In this paper we present a non-recursive direct solver, based on the Bartels-Stewart algorithm, for $N$-dimensional Sylvester tensor equations. The method relies only on Schur decompositions of the coefficient matrices and reduces the computation to a single sequential sweep over tensor entries, making it entirely independent of the dimension $N$. Its main advantages are simplicity, a dimension-independent formulation, and the ability to solve very high-dimensional problems limited only by available memory, which is used efficiently. We successfully solve cases up to $N=29$ on a standard laptop with $32$ GB RAM. Compared with the recursive blocked method of Chen and Kressner (state of the art), both approaches achieve identical accuracy. The recursive method is faster for large coefficient matrices, whereas our solver is competitive or superior when matrices are small, especially for large $N$, where recursive methods cannot effectively exploit BLAS-3 kernels. It also uses memory more efficiently: for near-capacity problems (e.g., matrices of order $19$ with $N=7$, where solution and right-hand side occupy $\approx 28$ GB), the Chen-Kressner method exceeds available memory, while ours succeeds. The method is also significantly simpler to implement, fully independent of $N$, and correctly handles singleton dimensions. We detail the algorithm and derive accurate cost estimates. For reproducibility, we provide pseudocode and complete MATLAB implementations. As an application, we compute solutions of linear $N$-dimensional ODE systems with constant coefficients at arbitrary times, and thus of evolutionary PDEs after spatial discretization, including highly accurate solutions of an advection-diffusion equation on $\mathbb{R}^N$.