Randomized methods for dynamical low-rank approximation

TL;DR

The paper tackles solving large-scale matrix differential equations by marrying dynamical low-rank approximation with randomized linear algebra. It introduces a dynamical rangefinder to estimate the evolving solution range and two post-processing time-stepping schemes, DR SVD and DGN, plus rank-adaptive variants, to efficiently evolve low-rank factors with robustness to stiffness. The methods exhibit exactness under reasonable range-inclusion assumptions and achieve favorable accuracy-cost trade-offs across stiff Lyapunov, Allen–Cahn, stochastic Burgers, and Vlasov–Poisson applications, while enabling potential conservation of physical quantities through range augmentation. This approach provides scalable, parallelizable tools for high-dimensional dynamics with practical impact on physics-informed simulations and uncertainty quantification.

Abstract

We introduce novel dynamical low-rank methods for solving large-scale matrix differential equations, motivated by algorithms from randomized numerical linear algebra. In terms of performance (cost and accuracy), our methods overperform existing dynamical low-rank techniques. Several applications to stiff differential equations demonstrate the robustness, accuracy and low variance of the new methods, despite their inherent randomness. Allowing augmentation of the range and corange, the new methods have a good potential for preserving critical physical quantities such as the energy, mass and momentum. Numerical experiments on the Vlasov-Poisson equation are particularly encouraging. The new methods comprise two essential steps: a range estimation step followed by a post-processing step. The range estimation is achieved through a novel dynamical rangefinder method. Subsequently, we propose two methods for post-processing, leading to two time-stepping methods: dynamical randomized singular value decomposition (DRSVD) and dynamical generalized Nyström (DGN). The new methods naturally extend to the rank-adaptive framework by estimating the error via Gaussian sampling.

Figures (7)

• Figure 1: Comparison of the dynamical rangefinder algorithm \ref{['algo: dynamical rangefinder']} with the rangefinder halko2011finding applied to the toy problem \ref{['eq: toy problem']}. The target rank of approximation is $r=5$ and the figure shows the relative error (at time $h=0.1$) made by the methods as we increase the oversampling parameter. The boxes show the median, first and third quartiles and extreme values of the $100$ simulations.
• Figure 2: Accuracy and computation time (in seconds) of the methods applied to the Lyapunov equation \ref{['eq: lyapunov differential equation']}. The integers next to the markers correspond to the target rank. The reference solver is a second-order exponential Runge--Kutta method with stepsize $\delta t = 10^{-4}$, and the substeps of the low-rank methods were solved with this reference solver.
• Figure 3: Rank-adaptive methods applied to the Allen--Cahn equation \ref{['eq: Allen--Cahn discrete equation']} with stepsize $h = 0.25$ and $30$ simulations. The prescribed tolerance is $\tau = 10^{-8}$. The reference is computed with a second order exponential Runge--Kutta method with time step $\delta t = 0.005$, and the substeps in the low-rank methods are also solved with the reference solver.
• Figure 4: Relative error over time of several methods applied to the stochastic Burgers equation \ref{['eq: matrix burgers equation']} with stepsize $h=0.01$ and $30$ simulations. The reference solver is scipy RK45 method with absolute and relative tolerance set to $10^{-12}$, and the substeps of the low-rank methods are solved with this reference solver.
• Figure 5: Numerical simulations of the linear Landau damping \ref{['eq: linear landau damping']} with stepsize $h=0.04$. The projector-splitting is performed with a Strang splitting ($\mathrm{KSL2}$) and the substeps are solved with the reference solver. The DRSVD is performed with oversampling parameter $p=5$ and $q=1$ power iteration. The DGN is performed with oversampling parameters $p=5$, $\ell=0$ and $q=1$ power iteration.

Theorems & Definitions (4)

• proof
• proof
• proof
• proof