1 Introduction↩︎

Low-rank approximation of large-scale matrices has demonstrated significant value in various fields, including bioinformatics [1], machine learning [2], [3], mathematical finance [4], control theory [5], [6], image compression [7], and natural language processing [8]. In this work, we seek low-rank approximate solutions to large-scale semilinear stiff matrix differential equations of the following form: \[\label{eq1461} \dot{X}(t)=AX(t)+X(t)A^{\top}+F(t,X(t)), \quad X(t_{0})=X_{0},\tag{1}\] where \(X(t)\in\mathbb{R}^{m\times n},F:[t_{0},\infty)\times\mathbb{R}^{m\times n}\to\mathbb{R}^{m\times n}\) is nonlinear, and \(X_0\), \(A\in\mathbb{R}^{m\times n}\) are time-independent. Such equations arise in a variety of large-scale applications. For instance, spatial semi-discretizations of nonlinear PDEs, such as the Allen–Cahn equation, usually lead to stiff linear terms generated by discrete differential operators. Differential Riccati equations from control theory provide another important example, where large-scale linear dynamics are coupled with nonlinear quadratic matrix terms. These problems motivate the development of efficient low-rank methods that can handle stiffness, nonlinearity, and high dimensionality simultaneously.

In recent years, efficient numerical solutions for large-scale stiff matrix differential equations have garnered significant attention in the fields of physics, chemistry, and engineering [9], [10]. The primary challenge lies in balancing computational efficiency with robustness against stiffness and high-dimensional nonlinearity. The traditional Lie-Trotter splitting method addresses the linear part analytically by separating the stiff and nonstiff terms, with recent extensions demonstrating efficacy for nonlinear fractional PDEs [11]. Then an exponential integrator is employed [12], [13]. However, a dynamic low-rank approximation (DLRA) of the nonlinear terms still depends on the projective-splitting integrator [14]. It faces challenges related to sensitivity to small singular values and bottlenecks in computational efficiency [15]. Although the DLRA significantly reduces dimensionality through manifold constraints [16], its classical implementation is constrained by explicit step size limitations and fixed rank assumptions [17]. Subsequent improvements, such as the adaptive BUG method [18], [19], enhance the stability but have limited scalability for high-dimensional nonlinear problems. The emergence of randomized numerical linear algebra introduces new ideas for low-rank decomposition, including the randomized singular value decomposition (RSVD) proposed by Halko et al. [20] and the generalized Nyström method developed by Tropp et al. [21]. These methods accelerate matrix approximation through sketching techniques. Additionally, Carrel [18] extends these approaches to dynamic scenarios by developing the dynamic RSVD (DRSVD) and dynamic generalized Nyström (DGN) methods. These methods demonstrate low variance and high parallelism in time-dependent systems [22]. However, most existing studies concentrate on static or weakly stiff scenarios, and the deep integration of the splitting framework with randomized dynamic low-rank (RDLR) methods has not been thoroughly investigated.

In this work, inspired by [18], to address the stiffness, nonlinearity, and unbearably large size of Eq. 1 , splitting-based RDLR approximations are designed. Our strategy not only inherits the robustness of the splitting method in handling stiffness [23] but also enhances the computational efficiency of the nonlinear terms through randomization techniques. To balance the accuracy and computational efficiency of the proposed methods, especially when the optimal rank is unknown beforehand and varies greatly over time, we extend them to rank-adaptive settings. It also provides a novel solving paradigm for higher-dimensional stiff systems, which has both high efficiency and stability. The main contributions can be summarized as follows:

• A novel RDLR method based on an operator splitting strategy has been developed, specifically tailored to the structural characteristics of Eq. 1 . This effectively decouples and handles the linear stiff term and the nonlinear term, enabling efficient computation with longer time and steps simulations than [18]. For instance, the Allen-Cahn equation inherently favors these features.

• The linear stiff term is integrated exactly using matrix exponentials [24]. This overcomes explicit step size restrictions inherent in standard methods sizes [25], [26], ensuring unconditional stability.

• The nonlinear term is resolved efficiently via the RDLR method [18]. This approach bypasses the computational bottlenecks associated with high-dimensional matrix operations.

• The proposed framework unifies computational efficiency, numerical accuracy, and stiffness stability. It is particularly well-suited for large-scale scientific computations, including high-dimensional PDEs, control theory [5], [6] (e.g., Riccati equations), and dynamical systems [16].

The rest of this paper is organized as follows. Section 2 provides details of our proposed methods. We first establish the Lie-Trotter and Strang splitting formats for Eq. 1 , respectively. Then, based on these splittings, splitting-based RDLR approximations are designed. Section 3 extends these approaches to the rank-adaptive scenario. Numerical validation in Section 4 includes convergence tests and large-scale simulations, demonstrating the excellent numerical performance. Section 5 presents the main conclusions and outlines directions for future research. The pseudocode of some key algorithms discussed in this paper is provided in Appendix.

2 Splitting-based randomized dynamical low-rank approximations↩︎

By utilizing the variation-of-constants formula [27], the exact solution of Eq. 1 can be expressed as follows: \[X(t)=\mathrm{e}^{(t-t_0)A}X(t_0)\mathrm{e}^{(t-t_0)A^{\top}}+\int_{t_0}^t\mathrm{e}^{(t-s)A}F(s,X(s))\mathrm{e}^{(t-s)A^{\top}}\mathrm{d}s,\] for \(t_{0}\leq t\leq T\).

The stiffness of Eq. 1 limits the applicability of most standard numerical integrators [24]. Explicit methods must satisfy a CFL-like condition and require the use of tiny time step sizes [28], while implicit methods result in significant computational costs. Consequently, we propose RDLR approximations based on the splitting method.

2.1 Splitting into two subproblems↩︎

Taking the structure of 1 into account, the key idea is the use of splitting method. For foundational theory and modern developments in operator splitting, see [24], [29].

Now, we split 1 into two terms: the stiff (linear) part and the nonstiff (nonlinear) part. This yields the following two subproblems:

\[\label{eq2461} \dot{M}(t)=AM(t)+M(t)A^{\top},\quad\quad M(t_{0})=M_{0},\tag{2}\] and \[\label{eq2462} \dot{N}(t)=F(t,N(t)),\quad\quad N(t_0)=N_0,\tag{3}\] where \(t\in[t_{0},T]\). The principal advantage of this splitting strategy stems from its ability to precisely decouple stiffness dynamics: the stiff component is entirely encapsulated within the first subproblem 2 , whereas the second subproblem 3 is a nonstiff ODE. This separation achieves dual computational benefits: (i) it expands the admissible choices of numerical integrators for the nonstiff subsystem, and (ii) enables systematic derivation of splitting schemes with arbitrary convergence orders by reconfiguring the operator composition in 1 .

For a given positive integer \(\hat{M}\), the time step size \(\tau\) is defined as \(\tau=\frac{T-t_0}{\hat{M}}\). Then, the time interval \([t_{0},T]\) is decomposed as \([t_0,T]=\bigcup\limits_{k=1}^{\hat{M}} [t_0, t_0 + k\tau]\). We denote the solutions of the subproblems (2) and (3) with initial values \(M_0\) and \(V_0\) by \(\Phi_\tau^A(M_0)\) and \(\Phi_\tau^F(N_0)\), respectively. The simplest splitting method for solving Eq. 1 is the Lie-Trotter splitting scheme with step size \(\tau\): \[\mathcal{L}_\tau:=\Phi_\tau^A\circ\Phi_\tau^F.\] With this notation, we can obtain an approximation \(X_{\tau}^{L}\) of the solution to \(X(t)\) 1 at \(t = t_0+\tau\), i.e., \[X_{\tau}^{L}=\mathcal{L}_\tau(X_{0})=\Phi_\tau^A\circ\Phi_\tau^F(X_{0}),\] where \(N_{0} = X_{0}\). This is a first-order method, and by the symmetric formula, we can have a second-order method called Strang splitting: \[S_{\tau}:=\Phi_{\tau/2}^{A}\circ\Phi_{\tau}^{F}\circ\Phi_{\tau/2}^{A}.\] Similarly, the solution of 1 at \(t_0 + \tau\) is approached as: \[X_{\tau}^{S}=S_\tau(X_{0})=\Phi_{\tau/2}^A\circ\Phi_\tau^F\circ\Phi_{\tau/2}^A(X_{0}),\] where \(M_{0} = X_{0}\).

Based on these two splitting methods, a low-rank solution of 1 is investigated in the next.

2.2 The dynamical low-rank integration↩︎

Denoting the manifold of rank-\(r\) matrices by \[\mathcal{M}_r=\{M\in\mathbb{R}^{m\times n}\mid\text{rank}(M)=r\}.\] Our aim is to find a low-rank approximation \(Y \in \mathcal{M}_r\) \((r\ll\min\{m,n\})\) to the solution of 1 . Under the splitting strategy introduced in Section 2.1, our primary objective turns to seeking low-rank approximations \(W \in \mathcal{M}_r\) and \(Z \in \mathcal{M}_r\) for subproblems 2 and 3 , respectively.

We first concentrate on the stiff subproblem 2 . For any \(Y \in \mathcal{M}_r\), we observe that \(AY + YA^{\top} \in \mathcal{T}_Y\mathcal{M}_r\), where \(\mathcal{T}_Y\mathcal{M}_r\) is the tangent space of \(\mathcal{M}_r\) at the current state \(Y\) [30]. This implies that if the initial condition lies in \(\mathcal{M}_r\), the solution of 2 remains within \(\mathcal{M}_r\) for all times. In other words, subproblem 2 preserves the rank of its solution.

The exact solution of this homogeneous problem at \(t_0+\tau\) within initial value \(M(t_{0})=M_{0}\) is given by \[M(t_0+\tau)=\mathrm{e}^{\tau A}M_{0}\mathrm{e}^{\tau A^{\top}}.\] If \(M_{0}\) is given in rank-\(r\) matrix factorization \(USV^{\top}\), where \(U\in\mathbb{R}^{m\times r}\), \(V\in\mathbb{R}^{n\times r}\) have orthonormal columns, \(S\in\mathbb{R}^{r\times r}(r\ll\min\{m,n\})\) is a nonsingular matrix. Subsequently, we obtain a rank-\(r\) approximation \(W_{\tau}\) to the exact solution \(M(t_0+\tau)=\Phi_\tau^A(M_{0})\) by \[W_{\tau}=\mathrm{e}^{\tau A}USV(\mathrm{e}^{\tau A}V)^{\top}.\] For the matrix exponential-matrix multiplications \(\mathrm{e}^{\tau A}U\) and \(\mathrm{e}^{\tau A}V\), several high performance algorithms exist, such as Taylor series expansion [31], Leja point-based interpolation [32], and Krylov subspace approximations [33], [34]. In order to maintain orthogonality, the matrices \(\mathrm{e}^{\tau A}U\) and \(\mathrm{e}^{\tau A}V\) can undergo the QR decomposition to enforce orthogonality. Specifically, we factorize \(\mathrm{e}^{\tau A}U=\widetilde{U}R_1\) where \(\widetilde{U}\in\mathbb{R}^{m\times r}\) has orthonormal columns and \(R_1\in\mathbb{R}^{r\times r}\) is an upper triangular matrix. Analogously, we decompose \(\mathrm{e}^{\tau A}V=\widetilde{V}R_2\) yields another orthonormal basis \(\widetilde{V}\in\mathbb{R}^{n\times r}\) with upper triangular matrix \(R_2\in\mathbb{R}^{r\times r}\). To preserve the rank-\(r\) matrix structure, the core matrix becomes \(\widetilde{S}=R_1SR_2^\top\), resulting in the low-rank approximation \(W_\tau=\widetilde{U}\widetilde{S}\widetilde{V}^\top\) of the solution to 2 .

For the nonlinear subproblem 3 , we employ the RDLR methods proposed in [18]. These new methods primarily involve two steps: a postprocessing step and a range estimation step (see details in Section 2.3). Firstly, using a dynamical rangefinder method to get an orthonormal matrix \(Q\in\mathbb{R}^{m\times r}\) (\(r\leq\min(m,n)\)) such that \[\left\| N(t) - QQ^{\top} N(t) \right\|_F \stackrel{!}{=} \min.\] Then, we utilize the DRSVD (or DGN) method to obtain a low-rank approximation \(Z \in \mathcal{M}_r\) to 3 . The author [18] has demonstrated that the new method achieves high accuracy, low variance, and significant robustness. After applying the RDLR to 3 , a low-rank approximation of \(N(t)\) at \(t = t_0+\tau\) with a low-rank initial value \(N_0\) is \[Z_\tau = \tilde{\Phi}_\tau^F(N_0).\] Here, \(\tilde{\Phi}_\tau^F\) represents the low-rank solution for 3 .

By combining the low-rank solution flow \(\tilde{\Phi}_\tau^F\) of subproblem 3 and the exact flow \(\Phi_\tau^A\) of subproblem 2 , we can get a rank-\(r\) approximate solution \(X^L_1\) of Eq. 1 at \(t = t_0+\tau\). This is denoted as \[\label{eq:star1} \tilde{\mathcal{L}_\tau}:=\Phi_\tau^A\circ\tilde{\Phi}_\tau^F,\tag{4}\] for low-rank Lie-Trotter splitting. Starting with \(Y_0\), a rank-\(r\) approximation of the initial value \(X_0\). The rank-\(r\) approximation of 1 at time \(t_0+\tau\) is written as \[X^\mathcal{L}_1=\tilde{\mathcal{L}_\tau}(Y_0)=(\Phi_\tau^A\circ\tilde{\Phi}_\tau^F)(Y_0).\]

Similarly, we denote as \[\tilde{\mathcal{S}_\tau}:=\Phi_{\tau/2}^A\circ\tilde{\Phi}_\tau^F\circ\Phi_{\tau/2}^A,\] for low-rank Strang splitting. Given a rank-\(r\) approximation \(Y_0\) of the initial value \(X_0\), the low-rank approximation of 1 at \(t_{0}+\tau\) is: \[X^\mathcal{S}_1=\tilde{\mathcal{S}_\tau}(Y_0)=\left(\Phi_{\tau/2}^A\circ\tilde{\Phi}_\tau^F\circ\Phi_{\tau/2}^A\right)(Y_0).\] To ensure the comprehensiveness and readability of our work, in the next section, we will elaborate how to find a low-rank solution of the nonlinear subproblem 3 .

2.3 Two randomized dynamical low-rank methods for Eq. 3↩︎

The RDLR framework proposed in [18] addresses large-scale matrix differential equations through novel adaptations of randomized numerical linear algebra techniques. Its framework is: firstly, a dynamical randomized rangefinder estimates the evolving range subspace; subsequently, two distinct post-processing strategies yield specialized time integrators, i.e., DRSVD and DGN. The RDLR method exhibits distinctive advantages over conventional dynamical low-rank approximation approaches [35], [36]: they behave without step size restrictions, and demonstrate inherent robustness. Before propose our two RDLR methods for the nonlinear subproblem 3 , we firstly establish a dynamical randomized rangefinder.

To investigate the initial response characteristics of 3 , we focuses on analyzing the physical processes within the first time step \([t_0,t_0+\tau]\). The aim of the dynamical rangefinder is to compute an orthonormal basis \(Q_\tau\in \mathbb{R}^{m\times(r+p)}\), where \(p\geq0\) is the oversampling parameter, such that the column space of \(Q_\tau\) approximates the range of \(N(t_0+\tau)\). Equivalently, the projection error \[\left\lVert N(t_0+\tau) - Q_\tau Q_\tau^{\top} N(t_0+\tau) \right\rVert_F\] is expected to be small. Here \(||\cdot||_F\) is the Frobenius norm.

A solution to the above optimization problem can be obtained via the following two steps:

1. Let \(\Omega \in \mathbb{R}^{n\times(r+p)}\) be a Gaussian random matrix. Instead of explicitly forming \(N(t_0+\tau)\Omega\), we evolve the sketched variable \(B(t)\in\mathbb{R}^{m\times(r+p)}\) by solving the differential equation: \[\dot{B}(t) = F\left(t, B(t)(\Omega^\top\Omega)^{-1}\Omega^\top\right)\Omega, \qquad B(t_0)=N_0\Omega .\]

2. After obtaining \(B(t_0+\tau)\), extract the orthonormal basis by \[Q_\tau = \operatorname{orth}\bigl(B(t_0+\tau)\bigr).\] This can be done, for example, by the QR decomposition.

When higher accuracy is required, the estimated subspace can be further improved by performing \(q\geq0\) power iterations through alternating range and co-range updates. The detailed derivation and analysis can be found in [18], and the complete procedure is summarized in Algorithm 11. Now, we are in position to propose two RDLR methods for solving 3 .

2.3.1 The dynamical randomized SVD↩︎

A rank-\(r\) approximation of the initial data \(N_0\) is represented as \(N_0 \approx U_0S_0V_0^\top\), where \(U_0\) and \(V_0\) are two column orthogonal matrices, and \(S_0\in\mathbb{R}^{r\times r}\) is nonsingular. When this factorisation is computed by a truncated SVD, \(S_0\) is diagonal with positive entries. The dynamic rangefinder identifies an orthogonal matrix \(Q_\tau\) such that \(N(t_0 + \tau) \approx Q_\tau Q_\tau^{\top} N(t_0 + \tau)\). This step effectively captures the structure of the low-rank subspaces of \(N(t)\) within the first time interval \([t_0,t_0+\tau]\). In order to ensure that the matrix \(Q\) includes foundational information along with dynamically evolving orthogonal bases, we combine the initial basis \(U_0\) with the dynamic basis \(Q_\tau\) to construct the extended orthogonal basis \(Q\), defined as \[Q = \text{orth}([U_0, Q_\tau]).\] This leads to a novel approximation \(N(t) \approx Q Q^{\top} N(t)\) \((t\in[t_0,t_0+\tau])\). It reduces the computational complexity by projecting the original high-dimensional nonlinear equations into the lower-dimensional subspace defined by the matrix \(Q\). The projection coefficient matrix \(C(t)\) is defined below to derive \(n\times(2r+p)\) the differential equation \[\label{eq2467} \dot{C}(t) = F(t,QC(t)^{\top})^{\top}Q, \quad C(t_0) = N_0^{\top}Q.\tag{5}\]

If we apply a SVD to \(C(t_0 + \tau)^{\top}\): \[\tilde{U}_\tau\Sigma_\tau V_\tau^{\top}=\text{svd}(C(t_0+\tau)^{\top}).\] Then, we have \(N(t_0 + \tau) \approx N_1 = Q \tilde{U}_\tau \Sigma_\tau V_\tau^{\top}\). Finally, the resulting low-rank decomposition is truncated to the target rank \(N_1=\mathcal{T}_r(Q\widetilde{U}_\tau\Sigma_\tau V_\tau^{\top})\).

The computational complexity can be significantly reduced through the use of low-rank decomposition techniques. In particular, as established in [18], when the system satisfies the regularity condition \(\text{Range}(N(t))\subset\text{Range}([U_0,Q_\tau])=\text{Range}(Q)\), for all \(t\in[t_0,t_0+\tau]\), Eq. 5 becomes exact reconstruction of the solution.

The DRSVD combines stochastic projection with dynamic evolution by estimating the column space through a dynamic rangefinder, enables efficient low-rank approximation of time-dependent matrix differential equations. The detailed procedure is summarized in Algorithm 12.

2.3.2 The dynamical generalised Nyström method↩︎

The DGN method provides an alternative approach for solving the nonlinear subproblem 3 . In contrast to DRSVD, which mainly relies on an estimated column space, DGN constructs a low-rank approximation by dynamically estimating both the column space and the row space of \(N(t)\) through a dynamical rangefinder and a dynamical co-rangefinder.

Similar to Section 2.3.1, for the rank-r approximation of the initial value \(N_0 \approx U_0S_0V_0^\top\), the DGN method proceeds as follows with oversampling parameters \(p,\ell \geq 0\), and power iterations \(q \geq 0\). It worth mentioning that, in here, \(S_0\) is not required to be diagonal; if the initial factorisation is obtained from a truncated SVD, then \(S_0\) is diagonal with positive singular values.

1. Compute the dynamical rangefinder basis \(\tilde{Q}_1\) for the column space and the dynamical co-rangefinder basis \(\tilde{Q}_2\) for the row space: \[\begin{align} \tilde{Q}_1 &= \operatorname{DynamicalRangefinder}(N_0, \tau, F, r, p, q), \\ \tilde{Q}_2 &= \operatorname{DynamicalCoRangefinder}(N_0, \tau, F, r, p+\ell, q). \end{align}\]

2. Augment the initial bases and form the orthonormal matrices: \[Q_1 = \operatorname{orth}([U_0, \tilde{Q}_1]), \qquad Q_2 = \operatorname{orth}([V_0, \tilde{Q}_2]).\]

3. Solve the following three reduced matrix differential equations on the low-dimensional subspaces: \[\begin{align} \dot{B}(t) &= F\bigl(t,B(t)Q_2^\top\bigr)Q_2, & B(t_0)&=N_0Q_2, \\ \dot{C}(t) &= F\bigl(t,Q_1C(t)^\top\bigr)^\top Q_1, & C(t_0)&=N_0^\top Q_1, \\ \dot{D}(t) &= Q_1^\top F\bigl(t,Q_1D(t)Q_2^\top\bigr)Q_2, & D(t_0)&=Q_1^\top N_0Q_2. \end{align}\] These three reduced equations can be solved simultaneously.

4. At \(t=t_0+\tau\), assemble the rank-\(r\) approximation as \[N_1 = B(t_0+\tau) \bigl(D(t_0+\tau)\bigr)_r^\dagger C(t_0+\tau),\] where \((\cdot)_r^\dagger\) denotes the truncated pseudo-inverse to rank \(r\).

The detailed derivation and analysis of the DGN method can be found in [18], and the complete procedure is summarized in Algorithm 13.

3 Implementation and extension↩︎

In this section, we present the core implementation of our splitting-based RDLR methods, and extend them to the rank-adaptive scenario.

3.1 Implementation↩︎

For Eq. 1 , we focus on how to compute the low-rank solution at \(t_1 = t_0 + \tau\). The low-rank solutions at other time steps (\(t_n = t_0 + n\tau\) for \(n \geq 2\)) can be obtained analogously.

For the linear subproblem 2 , the core operation involves computing the action of the matrix exponential \(\mathrm{e}^{\tau A}\) (Lie-Trotter splitting) or \(\mathrm{e}^{\frac{\tau}{2} A}\) (Strang splitting) on the current low-rank solution factor. Although efficient algorithms exist for matrix exponential operations (e.g., Taylor series expansion [31], Leja point interpolation [32], and Krylov subspace approximation [33], [34]), we utilize MATLAB’s built-in expm function due to its robustness and numerical stability. Crucially, the matrix exponential \(\mathrm{e}^{\tau A}\) or \(\mathrm{e}^{\frac{\tau}{2} A}\) is computed only once, prior to the main time integration loop, and reused throughout the time integration, providing significant computational savings.

The nonlinear subproblem 3 is solved via using the DRSVD method (i.e., Algorithm 12) or the DGN method (i.e., Algorithm 13). Within this framework, the underlying matrix differential equation is numerically integrated using a fourth-order Runge-Kutta (RK4). To ensure sufficient temporal resolution and accuracy while maintaining computational efficiency, we set \(10\) uniform substeps within each main time step of size \(\tau\). Moreover, our proposed splitting-based RDLR approximations are denoted as DRSVD-LT, DRSVD-ST, DGN-LT, DGN-ST. For example, DRSVD-LT means the DRSVD method is used in the low-rank Lie-Trotter splitting 4 . Other notations are defined in a similar manner. In these notations, the distinguishing labels "LT" and "ST" are added to distinguish which splitting method is choose.

For clarification, the DRSVD-LT and DRSVD-ST methods are summarized in Algorithms 1 and 2, respectively. Algorithms for the DGN-LT and the DGN-ST methods can be obtained by replacing the DRSVD method with the DGN method in Algorithms 1 and 2, respectively. Next, we will extend these approximations to the rank-adaptive scenario.

3.2 Rank-adaptive extension of the proposed methods↩︎

In most applications, the target rank is typically unknown in advance [37], or varies greatly over time. Consequently, the adaptive randomized rangefinder based on tolerance to dynamically adjust the target rank-\(r\) is proposed in [20]. This method relies on the following lemma.

Let \(A=N(t_0+\tau)-QQ^{\top}N(t_0+\tau)\) and \(\alpha=10\) the following conclusion is obtained: \[\begin{align} &\|N(t_0+\tau)-QQ^{\top}N(t_0+\tau)\|_2 \\ &\quad \leq 10\sqrt{\frac{2}{\pi}} \cdot \max_{i=1,\ldots,\kappa} \biggl\| \begin{aligned} &N(t_0+\tau)\omega^{(i)} - QQ^{\top}N(t_0+\tau)\omega^{(i)}\biggr\|_2 \end{aligned} \end{align}\] with probability \(1-10^{-\kappa}\).

Recall that the dynamical rangefinder facilitates the approximation \(B(t_0+\tau)\approx N(t_0+\tau)\Omega\). Gaussian matrix \(\Omega\), which leads to \(N(t_0+\tau)\omega^{(i)}=N(t_0+\tau)\Omega_{[:,i]}\approx B(t_0+\tau)_{[:,i]}\). In other words, \(B(t_0+\tau)\) can be used as an accurate error estimate for the range at (almost) no extra cost. In addition, [18] opts for block computation (adding \(\kappa\) columns each time) to cut iterations, balancing accuracy and speed. This procedure can be implemented by Algorithm 14. This approach offers a flexible and efficient adaptive framework, making it highly suitable for large-scale problems with unknown or time-dependent target ranks.

Similar to the duality transformation from the dynamical rangefinder to the dynamical corangefinder, the adaptive dynamical corangefinder is straightforward. Specifically, by replacing the B-step in the adaptive dynamical rangefinder (Algorithm 14) with its dual C-step formulation 5 .

Within the rank-adaptive framework, we systematically extend DRSVD and DGN to the rank-adaptive scenario. For the adaptive DRSVD, the primary modification involves replacing fixed-rank truncation with tolerance-based truncation. This eliminates the need for prespecified oversampling parameters while preserving the rank-preserving properties for the stiff term integration. Analogously, the adaptive DGN implements rank optimization by replacing both rangefinder and corangefinder steps with adaptive processes. Moreover, our methods under the rank-adaptive framework are denoted as ADRSVD-LT, ADRSVD-ST, ADGN-LT and ADGN-ST.

We prioritize this strategy as it not only removes subjectivity in oversampling parameter selection but also permits greater rank growth flexibility to capture transient solution features, enhancing adaptability for high-dimensional dynamical systems. The numerical results of the proposed methods are presented below.

4 Numerical results↩︎

This section validates the performance of the proposed algorithm through numerical experiments. In Section 4.1, we quantitatively evaluate the convergence order and computational accuracy using benchmark cases: the Allen-Cahn equation and the differential Riccati equation (DRE), with comparative analysis against existing splitting low-rank techniques. Section 4.2 examines computational efficiency and stability in complex scenarios via long-time simulations, providing theoretical foundations for practical implementations.

4.1 Convergence tests↩︎

This section aims to validate the convergence behavior of the proposed DRLR methods on two typical stiff problems. In the following, we will conduct a comparative analysis against the low-rank techniques in [13], [14], [17], [19], [38]–[40]. By varying time steps and spatial grid points, we evaluate convergence orders and computational accuracy of our methods, while presenting essential quantitative metrics for each case.

In these examples, for convenience, we set the number of grid points in both the \(x\)- and \(y\)-directions to \(\hat{N}(\hat{N}\in \mathbb{Z}^+)\) and the grid spacing to \(h\). We define the relative error at time \(T = t_0 + {\hat{M}}\tau\) as: \[\text{relerr}(\tau, h) = \frac{\left\| X_{\hat{M}} - X(T) \right\|_{F}}{\left\| X(T) \right\|_{F}},\] where \(X_{\hat{M}}\) is the numerical approximation of \(X(T)\). Then, we denote \[\text{rate}_{\tau} = \log_{\tau_1/\tau_2} \frac{\text{relerr}(\tau_1,h)}{\text{relerr}(\tau_2,h)} \quad \text{and} \quad \text{rate}_{h} = \log_{h_1/h_2} \frac{\text{relerr}(\tau,h_1)}{\text{relerr}(\tau,h_2)}.\]

4.1.1 Allen-Cahn equation↩︎

The Allen-Cahn equation, a standard category of stiff matrix differential equations, serves as a classical phase-field modeling framework predominantly employed to elucidate phase separation, interface dynamics, and the evolution of multiphase systems [41], [42]. The most basic representation of this equation includes a nonlinear term in the form of a cubic polynomial (reaction term) alongside a diffusion term [25], specifically: \[\label{eq4461} \frac{\partial f}{\partial t}=\varepsilon\Delta f+f-f^3,\tag{6}\] in this context, we set \(\varepsilon=0.01\). We designate the initial value as follows: \[f_0(x,y)=\frac{\left[\mathrm{e}^{-\tan^2(x)}+\mathrm{e}^{-\tan^2(y)}\right]\sin(x)\sin(y)}{1+\mathrm{e}^{|\csc(-x/2)|}+\mathrm{e}^{|\csc(-y/2)|}}.\]

The governing equations are numerically solved on the computational domain \(\Omega=[0,2\pi]^{2}\), which is discretized using a uniformly spaced grid with \(1024\times1024\) nodal points. Periodic boundary conditions are imposed along all domain boundaries. Through implementation of a second-order centered finite difference discretization, Eq. 6 is transformed into a system of ODEs. Compared to Eq. 1 , the matrix \(A\) becomes the corresponding discrete Laplace matrix, and \(A = A^\top\). The nonlinear term \(F(t, X(t))\) is replaced by \(X(t) - X(t)^{*3}\), where \(^{*3}\) is the cubic operation on elementwise. More precisely, \[\dot{X}(t)=AX(t)+X(t)A+X(t)-X(t)^{*3}.\]

Before systematically investigating convergence orders, we first evaluate the robustness of various methods under a large step size. In Table ¿tbl:tab95new?, we compare the performance of our methods (more precisely, DRSVD-LT, DRSVD-ST, DGN-LT and DGN-ST) with some existing DLRA methods applied to 6 , where \(\dagger\) represents that the method did not converge. The spatial grid size is set to \(N = 1024\), the final time to \(T = 0.1\), and only two time steps are taken, corresponding to a step size of \(h = 0.05\). The target rank is set to \(r = 12\). The results show that existing DLRA methods perform poorly due to the stiffness of the problem and the large time step. In contrast, the proposed four methods (i.e., DRSVD-LT, DRSVD-ST, DGN-LT and DGN-ST) methods remain computationally stable and achieve higher accuracy.

Having confirmed the robustness of the proposed methods under a large step size, we now systematically investigate their convergence behavior. The time interval is set to \([0, T]\) with \(T=1\). Fig. 3 demonstrates that Eq. 6 solved at low rank \(r = 16\), DRSVD-LT, DRSVD-ST, DGN-LT and DGN-ST consistently achieve lower relative errors across various time step sizes \(\tau\) compared to the conventional projector-splitting method in [13]. This precision advantage holds for both Lie-Trotter and Strang splitting schemes. Fig. 4 further reveals that the four proposed methods exhibit superior stability, characterized by weaker error accumulation over time. Throughout the simulation, the error curves for both proposed methods remain below that of projector-splitting. This stability advantage is especially evident when using the Strang splitting scheme.

Tables ¿tbl:tab1? and ¿tbl:tab2? present the relative errors in Frobenius norm and temporal convergence orders for Eq. 6 \((\hat{N} = 1024)\) using DRSVD-LT, DGN-LT, DRSVD-ST and DGN-ST methods. For Lie-Trotter splitting (Table ¿tbl:tab1?), both DRSVD-LT and DGN-LT exhibit the first-order convergence rate with negligible differences in errors. On the other hand, relative errors decrease consistently with refining \(\hat{M}\). DRSVD-ST and DGN-ST (Table ¿tbl:tab2?) achieves near second-order convergence for \(r\geq16\). While DRSVD-type splitting or DGN-type splitting yield nearly identical results, minor discrepancies in convergence rates suggest slight method-dependent variations. Overall, Strang splitting (DRSVD-ST and DGN-ST) demonstrates superior accuracy and theoretical second-order behavior, whereas Lie-Trotter splitting (DRSVD-LT and DGN-LT) maintains robust first-order convergence.

Fig. 5 illustrates the temporal evolution of errors for both DRSVD-type splitting or DGN-type splitting methods, compared to the optimal rank approximation error, under the conditions of a fixed rank \(r = 14\) and the number of time steps \(\hat{M} = 128\). The left and right subplots represent the results obtained using the DRSVD and DGN techniques for addressing the nonlinear term, respectively. It is evident from the figure that, regardless of the nonlinear treatment employed (DRSVD or DGN), the error curve of the Strang splitting method consistently remains below that of Lie-Trotter splitting method, thereby confirming its inherent advantage in achieving the second-order accuracy.

Table ¿tbl:tab3? presents temporal convergence behaviors in ADRSVD-LT, ADGN-LT, ADRSVD-ST and ADGN-ST for Eq. 6 with \(\hat{N}= 1024\). For the truncation step, the relative tolerance is set to \(\varphi = 10^{-8}\) and the absolute tolerance to \(10^{-12}\). The rangefinder step employs a tolerance of \(10^{-8}\). For the Lie-Trotter splitting method, both ADRSVD and ADGN yield identical relative errors in the Frobenius norm, achieving the first-order convergence. In contrast, the Strang splitting method exhibits significantly higher accuracy and approaches the second-order convergence. These results highlight the superior performance of the Strang splitting method and the minimal impact of the nonlinear term approximation method (ADRSVD vs ADGN) on convergence rates. Moreover, Fig. 6 illustrates the evolution of the numerical rank over time when the number of time steps \(\hat{M}\) is 256.

Tables ¿tbl:tab4? and ¿tbl:tab5? present the relative errors in the Frobenius norm and the spatial convergence orders for Eq. 6 (\(\hat{M} = 1024\)) using DRSVD-LT, DGN-LT, DRSVD-ST and DGN-ST methods. For the Lie-Trotter splitting method (Table ¿tbl:tab4?), both DRSVD and DGN demonstrate a spatial convergence behavior where the relative errors decrease as \(\hat{N}\) increases. The observed spatial convergence order is approximately 2 for larger number of \(\hat{N}\), indicating a consistent and stable convergence pattern. The errors for DRSVD and DGN are quite similar, suggesting comparable performance in capturing the spatial dynamics, with convergence primarily influenced by the spatial discretization process, as expected. For the Strang splitting method (Table ¿tbl:tab5?), a spatial convergence order of around 2 is also achieved, with errors decreasing steadily. The convergence is robust, and the results for DRSVD and DGN are nearly identical, indicating that both methods effectively manage the spatial discretization for Strang splitting.

Table ¿tbl:tab6? illustrates that for Eq. 6 , both the Lie-Trotter splitting and Strang splitting methods, when coupled with ADRSVD and ADGN discretizations, exhibit the second-order spatial convergence. As the spatial grid numbers \(\hat{N}\) increases from \(16\) to \(256\), the relative errors systematically decrease with increasing \(\hat{N}\), and the convergence rates stabilize near the theoretical value of 2. The errors and convergence rates of ADRSVD and ADGN are indistinguishable within numerical precision, demonstrating equivalent performance in spatial discretization. Preliminary analysis suggests that minor fluctuations in convergence rates at coarse grids (\(\hat{N} \leq 64\)) are primarily due to the dominance of truncation errors, while the predominance of spatial discretization errors at finer grids ensures the recovery of theoretical convergence rates. These results validate the robustness of the methods in addressing nonlinear terms and confirm the theoretical accuracy of the spatial discretization schemes.

4.1.2 Differential Riccati equation↩︎

Another example of a stiff PDE is called DRE, which plays a crucial role in control theory, filtering and estimation, and financial mathematics by solving problems such as optimal control, state estimation, and portfolio optimization [43]–[45].

In this part, we focus on DRE in the optimal control problem of linear quadratic regulators with parabolic PDEs over a finite time interval \(T\). Thus, we construct the following linear control system: \[\label{eq4462} \dot{x}=Ax+Bu,\quad x(t_0)=x_0,\tag{7}\] where \(A\in\mathbb{R}^{d\times d}\) and \(B\in\mathbb{R}^{d\times m}\) are the system matrices, \(x\in\mathbb{R}^m\) is the state variable, \(u\in\mathbb{R}^m\) is the control. The function \(\mathcal{J}\) that need to be minimized is as follows: \[\mathcal{J}(u,x)=\frac{1}{2}\int_{t_0}^{T}\left(x(t)^\top C^\top QCx(t)+u(t)^\top Ru(t)\right)\mathrm{d}t,\] where \(C\in\mathbb{R}^{q\times d}\), \({Q}\in\mathbb{R}^{q\times q}\) is symmetric and positive semidefinite and \({R}\in\mathbb{R}^{m\times m}\) is symmetric and positive definite. Moreover, the optimal control is given in feedback form by \(u_{\mathrm{opt}}(t)=-R^{-1}B^{\top}X(t)x(t)\), where \(X(t)\) is the solution of the following DRE: \[\dot{X}(t)=A^\top X(t)+X(t)A+C^\top QC-X(t)BR^{-1}B^\top X(t).\]

To further demonstrate the applicability of our approaches, we utilize a test case based on the diffusion-convection equation, which was developed from [46]. Consider the following equation: \[\label{eq4463} \partial_tw=\Delta w-10x\partial_xw-100y\partial_yw,\quad w|_{\partial\Omega}=0,\tag{8}\] defined on the domain \(\Omega = (0, 1)^2\), subject to homogeneous Dirichlet boundary conditions. The matrix \(A\) arises from standard centered finite difference spatial discretization of (\(\ref{eq4463}\)) with \(\widetilde{d}\) uniformly spaced stiff points per dimension. Let \(x_{i}=i\delta\) (\(\delta = (\widetilde{d} + 1)^{-1}\), \(i = 1,\dots,\widetilde{d}\)) denote interior discretization points in the \(x\)-direction. For matrix \(B \in \mathbb{R}^{\widetilde{d} \times 1}\), its entries follow: \[B_i= \begin{cases} 1 & \mathrm{if}\quad0.1<x_i\leq0.3, \\ 0 & \mathrm{otherwise}. \end{cases}\] For  \(C \in \mathbb{R}^{1 \times \widetilde{d}}\), the entries are: \[C_i= \begin{cases} 1 & \mathrm{if}\quad0.7<x_i\leq0.9, \\ 0 & \mathrm{otherwise}. \end{cases}\] Additionally, \(R = I\) and \(Q = 100I\) are selected. Similarly for the \(y\)-direction.

In the following numerical experiments, we take the initial value \(X_0 = I\), the final time \(T = 0.1\), and the reference solution calculated with the Dormand-Prince 5(4) method (denoted as dop54) [47]. We show the results for \(\hat{N} = 400\) grid points of the spatial discretization.

Tables ¿tbl:tab7? and ¿tbl:tab8? illustrate the numerical performance of solving Eq. 7 using DRSVD-LT, DGN-LT, DRSVD-ST and DGN-ST methods. The experimental results indicate that for the Lie-Trotter splitting method (as shown in Table ¿tbl:tab7?), the error performance of both DRSVD and DGN is comparable. As the rank-\(r\) and the number of time steps \(\hat{M}\) increase, the relative error gradually decreases, with the Lie-Trotter splitting method typically demonstrating the first-order convergence. In Table ¿tbl:tab8?, we present the corresponding results for the Strang splitting method, which exhibits the anticipated second-order convergence for sufficiently high ranks.

Fig. 7 illustrates the temporal evolution of errors for Lie-Trotter splitting and Strang splitting methods under the conditions of a fixed rank \(r = 8\) and the number of time steps \(\hat{M} = 512\). Specifically, the left subplot displays the results obtained using DRSVD to handle the nonlinear term, while the right subplot presents the results using DGN to address the nonlinear term. As anticipated, regardless of which nonlinear treatment (DRSVD or DGN) is employed, the error curve of the Strang splitting method consistently remains below that of the Lie-Trotter splitting method, confirming its inherent advantage in achieving the second-order accuracy.

Table ¿tbl:tab9? presents the numerical results obtained using ADRSVD-LT, ADGN-LT, ADRSVD-ST and ADGN-ST methods. For the truncation step, the relative tolerance is set to \(\varphi = 10^{-8}\) and the absolute tolerance to \(10^{-12}\). The rangefinder step employs a tolerance of \(10^{-8}\). The datas in tables clearly indicate that, under the same splitting method, the relative errors (relerr) and convergence orders (rate) of ADRSVD and ADGN are nearly identical. This suggests that the differences in these low-rank methods do not significantly affect the results. On the other hand, the relative errors decrease steadily with increasing the number of time steps \(\hat{M}\) and strang splitting performs exceptionally well with fine time steps; its convergence order approaches the theoretical second-order, and the relative error is low, demonstrating a significant advantage over the Lie-Trotter splitting method. Furthermore, in Fig. 8, we show the change in rank over time with a fixed number of time steps \(\hat{M} = 512\).

4.2 Simulation↩︎

In this section, to validate the numerical stability of the ARSVD-LT, ADGN-LT, ADRSVD-ST and ADGN-ST methods in long-term simulations, we focus on the following Allen-Cahn equation with logarithmic Flory–Huggins potential [48]: \[\label{eq4464} \begin{cases} \partial_tu=\varepsilon^2\Delta u-f(u),\quad(t,\mathbf{x})\in(0,\infty)\times\Omega \\ F(u)=\theta\left[u\mathrm{ln}u+(1-u)\mathrm{ln}(1-u)\right]+2\theta_cu(1-u),\mathrm{~for}\quad0<u<1 \\ \left.u\right|_{t=0}=u^0, & \end{cases}\tag{9}\] where \(u\) represents a real-valued function, \({\mathbf{\varepsilon}}\) is defined as the mobility coefficient, \(F(u)\) signifies the logarithmic Flory-Huggins potential, with \(f(u) = F^\prime(u)\). \(\theta\) and \(\theta_{c}\) are the absolute temperature and critical temperature respectively. Here we set \({\mathbf{\varepsilon}} = 0.1\), \(\theta = 0.8\) and \(\theta_c = 1\).

Besides, we select two different initial values for featuring high-curvature interfaces:

(1) Star-shaped: This is borrowed from [49], where

\[u(x,y,0) = \tanh \left( \frac{0.25 + 0.1 \cos(6 \tilde{\theta}) - \sqrt{(x-.5)^2 + (y - 0.5^2)}}{\varepsilon \sqrt{2}} \right),\] with \((x,y) \in [0,1]^2\) and \[\tilde{\theta} = \begin{cases} \tan^{-1} \left( \frac{y - 0.5}{x - 0.5} \right), & x > 0.5, \\ \pi + \tan^{-1} \left( \frac{y - 0.5}{x - 0.5} \right), & \mathrm{others}. \end{cases}\]

(2) Butterfly-shaped: This employs the high-curvature, butterfly-shaped interface construction initial value proposed in [49], which is parameterized by the following equation: \[\begin{cases} x(\hat{\theta})=k(a+b\cos(c\hat{\theta})\sin(d\hat{\theta}))\cos(\hat{\theta}), \\ y(\hat{\theta})=k(a+b\cos(c\hat{\theta})\sin(d\hat{\theta}))\sin(\hat{\theta}), & \end{cases}\] where \(\hat{\theta}\in[0,2\pi]\). In this case, we set \(k = 1.8\), \(a = 0.40178\), \(b = 0.30178\), \(c = 6.0\), and \(d = 3.0\).

The temporal evolution of these two initial shapes are shown in Figs. 9 and 10, respectively. Fig. 9 illustrates the dynamic evolution of numerical solutions for Eq. 9 with the star-shaped initial condition under spatial grid numbers \(\hat{N} = 128\) and the time step size \(\tau\) is \(0.01\), where the first two rows present the results of ADRSVD-LT and ADGN-LT, and the last two rows show those of ADRSVD-ST and ADGN-ST. Each column corresponds to the phase-field distribution at \(T = 0,10,50,100\), respectively. The initially star-shaped interface exhibits a gradual morphological evolution: the tip regions (high-curvature areas) contract inward, while the gap regions (low-curvature areas) expand outward, ultimately approaching a circular morphology-a behavior consistent with curvature-driven flow dynamics. Numerical results demonstrate that both splitting schemes yield solutions with high consistency in interface positioning, curvature distribution, and energy dissipation characteristics.

Fig. 10 presents numerical simulations of the dynamic evolution for Eq. 9 with butterfly-shaped initial condition under spatial grid numbers \(\hat{N} = 128\), where the first two rows present the results of ADRSVD-LT and ADGN-LT, and the last two rows show those of ADRSVD-ST and ADGN-ST. Each column correspondsto phase-field distributions at \(T = 0, 1, 5, 10\). The initial butterfly-shaped interface evolves with high-curvature regions contracting inward and low-curvature regions expanding outward, leading to a smoothed morphology over time. Results demonstrate consistent accuracy and stability in capturing interface evolution for both splitting schemes, while the adaptive low-rank methods effectively preserve phase-field quality without significant numerical dissipation or interface distortion, validating their robustness in handling complex geometric initial-value problems.

5 Concluding remarks↩︎

This work proposes a splitting-based RDLR approximation that unifies computational efficiency and numerical stability. We thoroughly consider the characteristics of 1 , dividing the problem into a stiff linear component and a nonstiff nonlinear component. For the stiff part, exponential integrators overcome explicit step size constraints; for the nonlinear part, RDLR approximation enables efficient resolution while bypassing high-dimensional computational bottlenecks. We establish complete integration frameworks for Lie-Trotter and Strang splittings, respectively, with extending to the rank-adaptive schemes. Through numerical experiments on the Allen-Cahn equation and the DRE, we verify the convergence orders and computational advantages of the proposed method. Furthermore, we demonstrate its robustness through extensive simulations of high-curvature initial interfaces over extended time periods. Based on this work, three future research directions are suggested:

• Extension to tensor manifolds to tackle higher-dimensional stiff PDEs.

• Establishment of a systematic error analysis framework for the proposed low-rank numerical scheme.

• The current method, utilizing truncated SVD for low-rank approximation, is effective when initial value perturbations are sufficiently small. However, truncated SVD may become unstable or inaccurate when initial perturbations are large. Future work needs to investigate and develop more robust low-rank approximation techniques to effectively handle such more challenging scenarios.

Acknowledgments↩︎

This work is supported by the National Natural Science Foundation of China (12401536), Sichuan Science and Technology Program (2024NSFSC0441) and Key Laboratory of Large-Scale Electromagnetic Industrial Software, Ministry of Education (EMCAE202502).

Data availability↩︎

The datasets generated and analyzed during the current study are not publicly available but are available from the authors on reasonable request.

Conflict of interests↩︎

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix↩︎

For the clarity, this appendix collects the pseudocode of some algorithms in Sections 2 and 3.

References↩︎