A CUR Krylov Solver for Large-Scale Linear Matrix Equations

TL;DR

The paper addresses solving very large-scale, multi-term linear matrix equations (LMEs) that arise from implicit time integration and nonlinear PDE discretizations. It develops a CUR-based low-rank framework on the manifold $\mathcal{M}_r$ combined with a Krylov solver to solve two small, decoupled subproblems corresponding to selected columns and rows, using a fixed-point iteration to remove cross-dependencies. Key contributions include CUR-DEIM sampling for near-optimal low-rank approximations, a Krylov-based solver operating on compressed forms, rank adaptivity, and extensions to non-Sylvester LMEs and Newton iterations for nonlinear problems. Demonstrations cover implicit time integration for MDEs, very large Lyapunov equations with up to $\sim 10^{13}$ unknowns, and non-Sylvester LMEs, revealing substantial speedups and scalable accuracy for industrial-scale PDE discretizations.

Abstract

Developing efficient solvers for large-scale multi-term linear matrix equations remains a central challenge in numerical linear algebra and is still largely unresolved. This paper introduces a methodology leveraging CUR decomposition for solving large-scale generalized Sylvester as well as non-Sylvester multi-term equations on low-rank matrix manifolds. The approach decomposes the original equation into two smaller subproblems: one involving all columns with a small subset of rows, and the other involving all rows with a small subset of columns. The rows and columns are strategically selected using the discrete empirical interpolation method. We further utilize the CUR properties and propose a novel iterative scheme that removes the dependencies between selected and unselected rows (and likewise for columns), thereby enabling the subset problems to be solved independently. We present a Krylov-based scheme for solving the resulting subproblems, which scales effectively to large problems and does not rely on a Sylvester structure. The method incorporates rank adaptivity, dynamically adjusting computational rank to reach the desired accuracy. The methodology is demonstrated in three representative settings: (i) implicit time integration of matrix differential equations on low-rank manifolds, leading to multi-term linear matrix equations; (ii) large-scale steady-state generalized Lyapunov equations including cases of size up to $10^{13}$ unknown entries; and (iii) non-Sylvester linear matrix equations with Hadamard product terms, such as those arising in nonlinear partial differential equations.

Figures (9)

• Figure 1: Computational domains used in the demonstrations: (a) two-dimensional plate with a hole denoted with $\Omega^{(1)}_{2D}$; (b) a square denoted with $\Omega^{(2)}_{2D}$; and (c) a three-dimensional domain denoted with $\Omega_{2D}$ and obtained by extruding $\Omega^{(1)}_{2D}$ in the $z$ direction.
• Figure 2: 3D FEM heat conduction on $\Omega_{3D}$: All the errors are computed against explicit Runge Kutta fourth order with time step $\Delta t= 1e-5$. The number of grid points $n_{xy} = 462$ in the $x$ and $y$ directions combined and $n_{z} = 61$ in the $z$ direction. Final time is $T_f=0.6$. Panel (a) shows the relative error over time-step size $\Delta t$ for FOM and TDB–CUR for Euler, BDF2 and BDF3 integrators at ranks $r = 8, 15$, and panel (b) shows the relative error versus rank size $r$ for Euler, BDF2, and BDF3 integrators at time step $\Delta t = 0.001$.
• Figure 3: 3D FEM heat conduction on $\Omega_{3D}$: the figure shows elapsed time for FOM and TDB-CUR methods at rank $r=8, 15$. Elapsed time, measured using MATLAB’s tic/toc, reflects total wall-clock time. The horizontal axis $n$ is the total number of grid points, i.e., $n=n_{xy} n_z$.
• Figure 4: First three modes from low-rank decomposition of the solution. The top row shows the first three left singular vectors $\mathbf{u}_1, \mathbf{u}_2, \mathbf{u}_3$, representing spatial modes in the $x$-$y$ plane. The bottom row shows the corresponding right singular vectors $\mathbf{y}_1, \mathbf{y}_2, \mathbf{y}_3$, which capture the spatial variation along the $z$-direction. All the plots are for $t=0.6$s.
• Figure 5: Steady-state heat conduction-radiation on $\Omega^{(1)}_{2D}$ with $n = 18620$. Panel (a) shows the singular values of the solution over the number of Newton iterations for TDB-CUR with ($r=20$) and all the singular values of FOM. The right panel shows the Frobenius norm of the Newton discrepancy matrix of the TDB-CUR solver for multiple values of rank and FOM.

Theorems & Definitions (1)

• Definition 1