A nested divide-and-conquer method for tensor Sylvester equations with positive definite hierarchically semiseparable coefficients

TL;DR

This paper introduces a nested divide-and-conquer method for efficiently solving tensor Sylvester equations with positive definite hierarchically semiseparable coefficients, achieving near-optimal computational complexity.

Contribution

The paper develops a novel nested divide-and-conquer algorithm leveraging low-rank updates for tensor Sylvester equations with hierarchically semiseparable matrices, providing theoretical analysis and validation.

Findings

Achieves quasi-optimal computational cost of O(n^d (log(n) + log(κ)^2 + log(κ) log(ε^{-1})))

Provides worst-case residual norm amplification estimates

Validated on 2D and 3D numerical case studies

Abstract

Linear systems with a tensor product structure arise naturally when considering the discretization of Laplace type differential equations or, more generally, multidimensional operators with separable coefficients. In this work, we focus on the numerical solution of linear systems of the form $$ \left(I\otimes \dots\otimes I \otimes A_1+\dots + A_d\otimes I \otimes\dots \otimes I\right)x=b,$$ where the matrices $A_t\in\mathbb R^{n\times n}$ are symmetric positive definite and belong to the class of hierarchically semiseparable matrices. We propose and analyze a nested divide-and-conquer scheme, based on the technology of low-rank updates, that attains the quasi-optimal computational cost $\mathcal O(n^d (\log(n) + \log(\kappa)^2 + \log(\kappa) \log(\epsilon^{-1})))$ where $\kappa$ is the condition number of the linear system, and $\epsilon$ the target accuracy. Our theoretical analysis highlights the role of inexactness in the nested calls of our algorithm and provides worst case estimates for the amplification of the residual norm. The performances are validated on 2D and 3D case studies.