In this work, we aim at approximating the solution \(A(t)\) to large-scale matrix differential equations of the form

\begin{equation} \label {eq: equation} \dot {A}(t)=F(A(t)),\quad A(0)=A_0\in \mathbb {R}^{m\times n}. \end{equation}

For large \(m\) and \(n\), the solution of (1) becomes expensive; in fact, it may not even be possible to store the entire matrix \(A(t)\) explicitly. To circumvent this limitation, model order reduction techniques based on exploiting (approximate) low-rank structure of \(A(t)\) can be employed. In particular, dynamical low-rank approximation [3] approximates \(A(t)\) by evolving matrices \(Y(t)\) on the manifold \(\mathcal {M}_r\) of rank-\(r\) matrices, reducing memory usage when \(r \ll m, n\). By the Dirac-Frenkel variational principle, the matrix \(Y(t)\) is obtained by solving the differential equation

\begin{equation} \label {dlra-ode} \dot {Y}(t)=P_r(Y(t))F(Y(t)),\quad Y(0)=Y_0\in \mathcal {M}_r, \end{equation}

where \(P_r(Y(t))\) denotes the orthogonal projection onto \(T_{Y(t)} {\mathcal M}_r\), the tangent space of \({\mathcal M}_r\) at \(Y(t)\). However, the stiffness of this equation leads to a severe step size restriction for standard explicit time integration methods. To address this issue, special integrators for this equation have been proposed [4, 2, 5]. Under the assumption

\begin{equation} \label {eq: assumption} \|F(Y)- P_r(Y)F(Y)\|_F\leq \tilde {\epsilon },\quad \text {for all }Y\in \mathcal {M}_r\cap \{\text {suitable neighbourhood of $A(t)$}\} \end{equation}

all these methods exhibit at least first-order convergence up to \(\mathcal {O}(\tilde {\epsilon })\).

Assumption (3), which states that \(F(Y)\) is nearly contained in the tangent space, is arguably a strong assumption. According to [2] and the examples shown, it is possible that \(A(t)\) can be well approximated by a rank-\(r\) matrix even if (3) is not satisified with small \(\tilde {\epsilon }\). When this assumption fails for small \(\tilde {\epsilon }\), using tangent space projections in numerical methods risks introducing significant errors.

In this work, we develop low-rank time integration methods for (1) that do not rely on (3) but only require \(A(t)\) to admit accurate low-rank approximations. Our approach is based on the notion of projected integrators, which first perform a standard time integration step and then project back to the manifold. For the manifold \(\mathcal M_r\), the efficiency of projected integrators is impaired by the occurrence of high-rank matrices, e.g., during the intermediate stages of a Runge-Kutta method. Previous work [2] mitigated this with repeated tangent space projections. Here, we propose a novel alternative using randomized low-rank approximation, employing random sketches instead of tangent projections to control rank growth efficiently.

To the best of our knowledge, this is the first work to propose and analyze randomized low-rank approximation methods for time integration. The randomized low-rank Runge-Kutta (RK) methods proposed in this work combine explicit RK methods with randomized low-rank approximation. Assuming that the dynamics generated by \(F\) preserve rank-\(r\) matrices approximately, we derive a probabilistic result that establishes a convergence order (up to the level of rank-\(r\) approximation error) based on the so-called stage order of the underlying RK method, which matches the order established in [2] for projected RK methods. However, unlike the results in [2], our numerical experiments indicate that randomized low-rank RK methods actually achieve the usual convergence order of the RK method, which can be significantly higher. For the randomized low-rank RK 4, we also establish order \(4\) theoretically when allowing for modest intermediate rank increases in the stages. This compares favorably to order \(2\) implied by the techniques in [2].

References

• [1]  Hei Yin Lam, Gianluca Ceruti, and Daniel Kressner. Randomized low-rank Runge-Kutta methods. arXiv preprint arXiv:2409.06384, 2024.

• [2]  Emil Kieri and Bart Vandereycken. Projection methods for dynamical low-rank approximation of high-dimensional problems. Comput. Methods Appl. Math., 19(1):73–92, 2019.

• [3]  Othmar Koch and Christian Lubich. Dynamical low-rank approximation. SIAM J. Matrix Anal. Appl., 29(2):434–454, 2007.

• [4]  Christian Lubich and Ivan V. Oseledets. A projector-splitting integrator for dynamical low-rank approximation. BIT, 54(1):171–188, 2014.

• [5]  Gianluca Ceruti, Lukas Einkemmer, Jonas Kusch, and Christian Lubich. A robust second-order low-rank BUG integrator based on the midpoint rule. BIT, 64(3):Paper No. 30, 2024.