pith. sign in

arxiv: 2412.15840 · v2 · submitted 2024-12-20 · 🧮 math.NA · cs.NA

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

Carlota M. Cuesta , Francisco de la Hoz This is my paper

Pith reviewed 2026-05-23 07:17 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords Sylvester tensor equationsSchur decompositionBartels-Stewart algorithmnon-recursive solverhigh-dimensional tensorstensor equationsdimension-independent algorithms
0
0 comments X

The pith

A non-recursive Schur-based method solves N-dimensional Sylvester tensor equations with a single sequential sweep independent of dimension.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper introduces a direct solver for N-dimensional Sylvester tensor equations that avoids recursion by relying on Schur decompositions of the coefficient matrices. The computation reduces to one sequential sweep through the tensor entries, allowing the algorithm to work for any number of dimensions without modification. This approach solves instances with up to 29 dimensions on ordinary hardware and consumes memory more efficiently than existing recursive techniques while delivering matching accuracy. The formulation is straightforward to code and correctly manages cases where some dimensions are size one. Applications include solving systems of linear ordinary differential equations in many variables and discretized evolutionary partial differential equations.

Core claim

The algorithm applies the classical Bartels-Stewart procedure by first computing Schur decompositions of each coefficient matrix and then constructing the solution tensor through a single pass that visits each entry exactly once, with no recursive calls and no dependence on the value of N.

What carries the argument

the single sequential sweep over tensor entries after computing the Schur factors of the coefficient matrices

If this is right

  • The solver handles problems with N as large as 29 using only 32 GB of RAM.
  • It matches the accuracy of the recursive Chen-Kressner method for the same problems.
  • It performs competitively or better when the coefficient matrices are small and N is large.
  • It requires less memory than recursive methods near the limit of available RAM.
  • It is simpler to implement and fully independent of N.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The method could be adapted to other multilinear tensor equations that admit Schur-based reductions.
  • High-dimensional PDE discretizations may become feasible on modest hardware with this approach.
  • Memory usage patterns suggest the technique scales better to extreme dimensions than recursion-based alternatives.

Load-bearing premise

The sequential sweep assembles the correct solution from the Schur factors without any hidden dimension-dependent corrections or stability adjustments.

What would settle it

Running both this algorithm and the recursive method on the same N=4 problem with known exact solution and checking whether the maximum entrywise difference exceeds floating-point roundoff.

Figures

Figures reproduced from arXiv: 2412.15840 by Carlota M. Cuesta, Francisco de la Hoz.

Figure 1
Figure 1. Figure 1: Elapsed time need to solve an N-dimensional problem, where X ∈ C 2×2×...×2 . 5.1 N-dimensional systems of ODEs In order to test Theorem 3.1, we have considered a 7-dimensional example. More precisely, we have chosen N = 7, X ∈ C 2×3×4×5×6×7×8 , generated randomly A1, A2, . . . , A7, B and X0, and approximated the so￾lution of Xt = P7 j=1 AjjX+B, at t = 0.1 by means of the Matlab program evolND.m, obtaining… view at source ↗
read the original 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$.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

0 major / 2 minor

Summary. The paper presents a non-recursive direct solver for N-dimensional Sylvester tensor equations based on the Bartels-Stewart algorithm. It uses Schur decompositions of the coefficient matrices and reduces the solution process to a single sequential sweep over tensor entries, yielding a formulation that is entirely independent of the dimension N. The method is claimed to solve problems up to N=29 on standard hardware with 32 GB RAM, achieve identical accuracy to the recursive Chen-Kressner method, use memory more efficiently in near-capacity cases, and include explicit cost estimates, pseudocode, and MATLAB implementations, with an application to high-dimensional linear ODE systems.

Significance. If the central claims hold, the work provides a practical advance for high-dimensional tensor linear algebra by delivering a simpler, dimension-independent solver whose memory usage scales only with the tensor size rather than recursion depth. The explicit provision of reproducible code, cost formulas, and direct comparisons to the state-of-the-art recursive method strengthens its utility for applications in evolutionary PDEs after spatial discretization.

minor comments (2)
  1. The abstract states that 'exact flop counts' are derived, yet the provided text does not display the final simplified expressions; placing these in a dedicated table or subsection would improve clarity.
  2. Section describing the sequential sweep would benefit from an explicit statement confirming that no additional N-dependent stability corrections are introduced beyond the classical Bartels-Stewart analysis.

Simulated Author's Rebuttal

0 responses · 0 unresolved

We thank the referee for the positive report and the recommendation to accept the manuscript. No major comments were raised.

Circularity Check

0 steps flagged

No significant circularity; derivation is self-contained extension of Bartels-Stewart

full rationale

The paper presents an algorithmic construction that applies the standard (non-self-cited) Bartels-Stewart procedure via Schur decompositions to tensor Sylvester equations, followed by an explicit sequential sweep over entries. This sweep is defined directly from the factored form without re-using the target solution as an input or fitting any parameter. Pseudocode, MATLAB source, and cost formulas are supplied for independent reproduction. No equation reduces to a prior result by definition, no fitted input is relabeled as a prediction, and the central claim (N-independence via single sweep) is tested by direct implementation rather than by self-citation chains. The comparison to Chen-Kressner is external benchmarking, not load-bearing justification. Score remains 0 because the derivation chain contains no self-referential step.

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 0 invented entities

The algorithm rests on the existence of Schur decompositions for the coefficient matrices and on the correctness of the classical Bartels-Stewart reduction; no new free parameters, invented entities, or ad-hoc axioms are introduced.

axioms (1)
  • standard math Every square matrix over the complex numbers admits a Schur decomposition.
    The method relies only on Schur decompositions of the coefficient matrices as stated in the abstract.

pith-pipeline@v0.9.0 · 5845 in / 1292 out tokens · 44875 ms · 2026-05-23T07:17:20.922191+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Forward citations

Cited by 1 Pith paper

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score.

  1. A matrix-based spectral method for the numerical approximation of the fractional Laplacian and the fractional $p$-Laplacian of functions defined on $\mathbb R^n$

    math.NA 2026-05 unverdicted novelty 6.0

    A spectral method using diagonalizable differentiation matrices evaluates the Balakrishnan integral analytically on eigenvalues to approximate fractional Laplacians and p-Laplacians, demonstrated on 1D and 2D fraction...

Reference graph

Works this paper leans on

16 extracted references · 16 canonical work pages · cited by 1 Pith paper

  1. [1]

    R. H. Bartels and G. W. Stewart. Solution of the matrix equation AX +XB =C. Communications of the ACM , 15(9):820–826, 1972

    work page 1972

  2. [2]

    How and Why to Solve the O perator EquationAX−XB =Y

    Rajendra Bhatia and Peter Rosenthal. How and Why to Solve the O perator EquationAX−XB =Y . Bulletin of the London Mathematical Society , 29(1):1–21, 1997

    work page 1997

  3. [3]

    Recursive blocked algorithms for linear systems with Kronecker product structure

    Minhong Chen and Daniel Kressner. Recursive blocked algorithms for linear systems with Kronecker product structure. Numerical Algorithms, 84:1199–1216, 2020

    work page 2020

  4. [4]

    A Tensor Multigrid Method for Solving S ylvester Tensor Equations

    Yuhan Chen and Chenliang Li. A Tensor Multigrid Method for Solving S ylvester Tensor Equations. IEEE Transactions on Automation Science and Engineering , 21(3):4397–4405, 2024

    work page 2024

  5. [5]

    A Sylvester-Based IMEX Method via Differentiation Matrices for Solving Nonlinear Parabolic Equations

    Francisco de la Hoz and Fernando Vadillo. A Sylvester-Based IMEX Method via Differentiation Matrices for Solving Nonlinear Parabolic Equations. Communications in Computational Physics , 14(4):1001–1026, 2013

    work page 2013

  6. [6]

    The solution of two-dimen sional advection-diffusion equations via operational matrices

    Francisco de la Hoz and Fernando Vadillo. The solution of two-dimen sional advection-diffusion equations via operational matrices. Applied Numerical Mathematics , 72:172–187, 2013. 15

    work page 2013

  7. [7]

    Grasedyck

    L. Grasedyck. Existence and Computation of Low Kronecker-R ank Approximations for Large Linear Systems of Tensor Product Structure. Computing, 72:247–265, 2004

    work page 2004

  8. [8]

    Nicholas J. Higham. Functions of Matrices. Theory and Computation . Society for Industrial and Applied Mathematics, 2008

    work page 2008

  9. [9]

    Schur -decomposition for 3D matrix equations and its application in solving radiative discrete ordinates eq uations discretized by Cheby- shev collocation spectral method

    Ben-Wen Li, Shuai Tian, Ya-Song Sun, and Zhang-Mao Hu. Schur -decomposition for 3D matrix equations and its application in solving radiative discrete ordinates eq uations discretized by Cheby- shev collocation spectral method. Journal of Computational Physics , 229(4):1198–1212, 2010

    work page 2010

  10. [10]

    Massei and L

    S. Massei and L. Robol. A nested divide-and-conquer method f or tensor Sylvester equations with positive definite hierarchically semiseparable coefficients. IMA Journal of Numerical Analysis , 44(6):3482–3519, 2024

    work page 2024

  11. [11]

    B. N. Parlett. A recurrence among the elements of functions o f triangular matrices. Linear Algebra and its Applications , 14(2):117–121, 1976

    work page 1976

  12. [12]

    G. W. Stewart. Stochastic Automata, Tensors Operation, an d Matrix Equations, 1996. Technical Report, University of Maryland, UMIACS TR-96-11, CMSC TR-3598

    work page 1996

  13. [13]

    Sur l’´ equation en matrices px = xq

    James Joseph Sylvester. Sur l’´ equation en matrices px = xq. Comptes Rendus de l’Acad´ emie des Sciences, 99(2):67–71, 115–116, 1884. in French

    work page

  14. [14]

    Approximated Tensor Sum Preconditione r for Stochastic Automata Networks

    Abderezak Touzene. Approximated Tensor Sum Preconditione r for Stochastic Automata Networks. In IPDPS’06: Proceedings of the 20th international conferenc e on Parallel and distributed processing , 2006

    work page 2006

  15. [15]

    J. A. C. Weideman and S. C. Reddy. A MATLAB differentiation matr ix suite. ACM Transactions on Mathematical Software , 26(4):465–519, 2000. The codes are available at https://appliedmaths.sun.ac.za/∼ weideman/research/differ.html

    work page 2000

  16. [16]

    Developing iterative algorith ms to solve Sylvester tensor equations

    Xin-Fang Zhang and Qing-Wen Wang. Developing iterative algorith ms to solve Sylvester tensor equations. Applied Mathematics and Computation , 409:126403, 2021. 16

    work page 2021