1 Introduction↩︎

Large-scale matrix differential equations arise in a variety of scientific computing applications. In kinetic theory, the Vlasov–Poisson equation describes the evolution of particle distributions in collisionless plasmas. In materials science, the Allen–Cahn equation models phase separation processes. In stochastic analysis, the Fokker–Planck equation governs the evolution of probability density functions. When discretized on tensor-product grids, these equations lead to matrix differential equations of the form \[\label{eq:fom95intro} \dot{A}(t) = F(t, A(t)), \quad A(0) = A_0 \in \mathbb{R}^{m \times n},\tag{1}\] where \(m\) and \(n\) can be very large. We refer to 1 as the full order model (FOM). In the rest of the paper, we will consider autonomous problems, without loss of generality, in order to keep the notation light.

In the last decade, various techniques for compressing the dynamics have been proposed. One of the most popular and well-studied techniques is the dynamical low-rank approximation (DLRA) [1], which has since been applied to kinetic equations [2], [3], weakly compressible flows [4], and high-dimensional nonlinear PDEs [5]. The key idea is to evolve the lower-dimensional dynamics by projecting the vector field onto the tangent space of the manifold \(\mathcal{M}_r\) of fixed rank-\(r\) matrices. This yields a system of coupled differential equations for the low-rank factors, reducing the complexity from \(O(mn)\) to \(O((m+n)r)\). The DLRA is defined by the Dirac–Frenkel variational principle: \[\dot{Y}(t) = \mathcal{P}_{Y(t)}F(Y(t)), \quad Y(0) = Y_0 \in \mathcal{M}_r,\] where \(\mathcal{P}_Y\) denotes the \(\ell_2\)-orthogonal projection onto the tangent space \(\mathcal{T}_Y \mathcal{M}_r\).

Practical numerical integrators for the DLRA include the projector-splitting integrator [6] and the BUG integrator [7], together with its higher-order and rank-adaptive variants [8], [9] and related projection methods [10], [11], as well as interpolatory low-rank integrators [12]. These methods always require QR decompositions at each time step in order to orthogonalize the updated bases. In the large-scale setting on distributed hardware, the Householder QR of \(K \in \mathbb{R}^{m \times r}\) suffers from a synchronization bottleneck: it performs \(r\) sequential panel factorizations, each requiring a global reduction on the row-distributed factor [13]; the panel factorizations themselves are BLAS-2 operations. Randomized numerical linear algebra (RNLA) [14]–[17] provides a natural alternative. These sketch-based orthogonalizations replace the Householder QR by working on a small sketch \(\Theta K \in \mathbb{R}^{\ell \times r}\) of the tall factor \(K \in \mathbb{R}^{m \times r}\) (with \(\Theta \in \mathbb{R}^{\ell \times m}\)). In its simplest form, randomized Cholesky QR computes the QR decomposition of \(\Theta K\) and uses the resulting triangular factor to orthogonalize \(K\), collapsing the whole orthogonalization to a single BLAS-3 matrix multiply and a single global reduction. The more robust randomized Gram–Schmidt (RGS) builds the basis column by column instead; at the typical sketch size \(\ell = O(r)\) its leading-order flop count \(O(m \ell r)\) matches Householder QR, while it requires far fewer global synchronizations — the source of the speedup on modern parallel hardware (Section 5.2).

Beyond the orthogonalization step, randomization has also been used directly inside the time integration of low-rank matrix ODEs, most notably through the randomized Runge–Kutta methods of [18] — a direction complementary to ours, which keeps the time integrator classical and randomizes only the basis update. Closely related are the dynamical randomized methods of [19], which build the low-rank integrator itself from randomized linear algebra — randomizing the range-finding and projection steps rather than only the orthogonalization — and report accurate, low-variance, structure-preserving behavior on stiff problems, including the conservation of physical invariants on Vlasov–Poisson. The sketch orthogonalization developed here can serve as the basis update within those methods as well.

There are two natural ways to incorporate sketching into DLRA:

• Oblique approach: sketch the Galerkin condition that defines the DLRA. This leads to an oblique tangent space projection involving the sketch matrices \(\Theta\) and \(\Omega\). The resulting dynamics are different from the standard DLRA.

• Orthogonal approach: evolve sketch-orthogonal bases (e.g.via RGS) but use standard orthogonal projections for the dynamics. This preserves the geometric structure of the classical DLRA. The sketch only affects the runtime of the orthogonalization process, not the dynamics.

The key mathematical insight of this paper is that the oblique approach changes the dynamics and leads to errors on problems with a large perpendicular residual (see Definition 2), while the orthogonal approach preserves the dynamics exactly. Table 1 summarizes the qualitative differences; a detailed comparison is given in Table 2. This is illustrated in Figure 1, which shows the energy over time of a two-stream instability simulation for both approaches. The oblique sketch DLRA exhibits energy drift, while the orthogonal sketch DLRA matches the standard DLRA with a near-constant energy in the nonlinear regime.

Contributions and outline↩︎

The main contributions of this paper are:

1. We propose an orthogonal sketch DLRA framework that evolves sketch-orthogonal bases \((P, W)\) while using standard orthogonal projections. This is connected to the classical DLRA via a bijective map \(\varphi\) that relates sketch-orthogonal and orthogonal factorizations.

2. We show that sketching the Galerkin condition leads to an oblique projection that is equivalent to the classical DLRA if and only if the sketches make the perpendicular component \(\mathcal{P}_Y^\perp F(Y)\) orthogonal to the tangent space (Proposition 5). This condition holds up to the sketching error on low-rank compatible problems (Section 4.2), while on problems with a large perpendicular residual (e.g.Vlasov–Poisson) the oblique approach fails in general.

3. We derive sketch versions of the projector-splitting and BUG integrators, replacing Householder QR with RGS.

4. We demonstrate the approach on the Allen–Cahn, Fokker–Planck, and Vlasov–Poisson equations. The orthogonal sketch DLRA matches the accuracy of the standard DLRA, while the oblique approach loses physical invariants (e.g.the electric energy) as soon as the perpendicular residual becomes large.

Preliminaries are collected in Section 2; the orthogonal and oblique constructions are developed in Sections 3 and 4; sketch integrators and their cost are derived in Section 5; numerical experiments are reported in Section 6.

2 Preliminaries↩︎

2.1 Oblivious subspace embeddings↩︎

We recall the definition of an oblivious subspace embedding (OSE), introduced by [20] in a sketching context and subsequently surveyed in [14]–[16], [21].

Typical choices for oblivious subspace embeddings include Gaussian random matrices, the subsampled randomized Hadamard transform (SRHT) [22], and sparse sign embeddings [23]. The sampling dimension \(\ell\) depends on the target subspace dimension: for the rank-\(r\) subspaces that arise in this paper, a practical choice for large-scale problems is \(\ell \in [2r, 4r]\) [24], while for the moderate sizes of Section 6 we find that \(\ell = 2r\) is already sufficient. In both regimes \(\ell \ll m, n\). The OSE property extends from individual vectors to tall matrices, as recorded in the next proposition.

The proposition follows by applying the OSE property in its equivalent inner-product form (Definition 1) to the fixed subspace \(\mathrm{Range}([U \;\, V]) \subset \mathbb{R}^m\), whose dimension is at most \(2r\), and summing over the \(r\) column pairs.

The two-sided pairings that appear in Sections 4 and 4.2 are controlled by the following corollary, obtained by applying Proposition 1 on each side and collecting the cross terms.

2.2 Dynamical low-rank approximation↩︎

Recall the full order model (FOM) equation \(\dot{A}(t) = F(A(t))\) with initial condition \(A(0) = A_0 \in \mathbb{R}^{m \times n}\), introduced in 1 . The dynamical low-rank approximation (DLRA) approximates the FOM trajectory by a curve on the manifold of fixed-rank-\(r\) matrices \(\mathcal{M}_r= \{Y \in \mathbb{R}^{m \times n}: \mathop{\mathrm{rank}}(Y) = r\}\), parameterized as \(Y = U \Sigma V^T\) with \(U \in \mathrm{St}(m,r)\), \(V \in \mathrm{St}(n,r)\), and \(\Sigma \in \mathrm{GL}(r)\). Its tangent space at \(Y\) is the set \[\mathcal{T}_Y \mathcal{M}_r= \{\delta U \Sigma V^T + U \delta \Sigma V^T + U \Sigma \delta V^T \mid \delta U \in \mathbb{R}^{m \times r},\;\delta \Sigma \in \mathbb{R}^{r \times r},\;\delta V \in \mathbb{R}^{n \times r},\;\delta U^T U = 0,\;\delta V^T V = 0\},\] where \(\delta U^T U = 0\) and \(\delta V^T V = 0\) are gauge conditions ensuring a unique representation of tangent vectors (see [25], [26] for the geometric background of Stiefel and low-rank manifolds). The orthogonal projection onto \(\mathcal{T}_Y \mathcal{M}_r\) is given by \[\label{eq:tangent95proj} \mathcal{P}_Y Z = UU^T Z + Z VV^T - UU^T Z VV^T.\tag{2}\]

The DLRA trajectory \(Y(t)\) is defined by the Galerkin condition \[\label{eq:galerkin} \left\langle \dot{Y}(t) - F(Y(t)), \, \delta Y \right\rangle = 0 \qquad \forall \delta Y \in \mathcal{T}_{Y(t)} \mathcal{M}_r.\tag{3}\] Decomposing the vector field as \(F(Y) = \mathcal{P}_Y F(Y) + \mathcal{P}_Y^\perp F(Y)\), the Galerkin condition automatically annihilates the perpendicular component against any tangent direction, since \(\left\langle \mathcal{P}_Y^\perp F(Y),\, \delta Y \right\rangle = 0\) for every \(\delta Y \in \mathcal{T}_Y \mathcal{M}_r\). The size of this perpendicular component is what governs the modeling error of the DLRA. Using the projector \(\mathcal{P}_Y\) above, the Galerkin condition is equivalent to the projected matrix differential equation \[\label{eq:dlra} \dot{Y}(t) = \mathcal{P}_{Y(t)} F(Y(t)), \quad Y(0) = Y_0 \in \mathcal{M}_r,\tag{4}\] which is the standard DLRA equation found in the literature [1]. Throughout, the vector field in the low-rank dynamics is evaluated at the current approximation \(Y\) — we write \(F(Y)\), or \(F(U \Sigma V^T)\) (resp.\(F(P S W^T)\)) in factored form — and we reserve \(F(A)\) for the full-order model 1 ; when \(Y \approx A\) the two fields nearly coincide.

From that description of the tangent space and with the Galerkin condition, it is straightforward to derive the three coupled differential equations \[\label{eq:coupled95equations} \left\{ \begin{align} &\dot{U} = (I - UU^T) F(U \Sigma V^T) V \Sigma^{-1}, \\ &\dot{\Sigma} = U^T F(U \Sigma V^T) V, \\ &\dot{V} = (I - VV^T) F(U \Sigma V^T)^T U \Sigma^{-T}. \end{align} \right.\tag{5}\]

As shown in [1], solving these three coupled equations is ill-conditioned when \(\Sigma\) has small singular values, because of the \(\Sigma^{-1}\) factor. This was a major numerical issue for DLRA and motivated the development of the projector-splitting and BUG integrators [6], [7], which avoid inverting \(\Sigma\) by splitting the dynamics into well-conditioned substeps.

A central notion for the upcoming analysis of the oblique sketch DLRA is low-rank compatibility, which classifies how much of the vector field \(F\) escapes the tangent space of \(\mathcal{M}_r\). We use \(\eta\) for the low-rank compatibility tolerance throughout, to avoid collision with the OSE parameter \(\epsilon\) of Section 2.1.

The special case \(\eta = 0\) yields exact \((r,k)\)-low-rank compatibility (the perpendicular component has rank at most \(k\)), while the case \(k = 0\) recovers the norm-based condition of [27].

2.3 Sketch-orthogonality and randomized Gram–Schmidt↩︎

We now introduce the notion of sketch-orthogonality, which will be central to our approach. Consider two oblivious subspace embeddings \(\Theta \in \mathbb{R}^{\ell \times m}\) and \(\Omega \in \mathbb{R}^{\ell \times n}\) with sampling parameter \(\ell > 0\).

A matrix \(P \in \mathrm{St}_\Theta(m, r)\) is said to have sketch-orthogonal columns with respect to \(\Theta\). Note that \(P\) itself does not have orthogonal columns in general; only its sketch \(\Theta P \in \mathbb{R}^{\ell \times r}\) does.

Given \(K \in \mathbb{R}^{m \times r}\) with columns \(k_1, \dots, k_r\), the simplest approach to produce a sketch-orthogonal basis \(P\) is to sketch \(\Theta K \in \mathbb{R}^{\ell \times r}\), compute its thin QR factorization \(\Theta K = Q R\), and set \(P = K R^{-1}\). Then \(\Theta P = Q\) has orthonormal columns, so \(P \in \mathrm{St}_\Theta(m, r)\). This process is called randomized Cholesky QR (rCholQR) [28], [29]. However, it can be unstable when \(K\) is ill conditioned or breaks down if \(\Theta K\) is rank-deficient, since \(R\) would then be singular. The randomized Gram–Schmidt (RGS) process described in [24] avoids this issue by progressively building the basis \(P\), as in the Gram–Schmidt process, while replacing orthogonal projections with oblique projections to produce a \(P\) whose columns are orthonormal with respect to the sketched inner product \(\left\langle \Theta \cdot, \Theta \cdot \right\rangle\). In more detail, RGS builds \(P = [p_1, \dots, p_r]\) and its sketch \(\Theta P\) column by column as follows: at step \(j\), given the columns built previously \(P_{<j} = [p_1, \dots, p_{j-1}]\) and their sketch \(\Theta P_{<j}\),

• sketch the new column, \(s_j = \Theta k_j\);

• compute the projection coefficients in the sketch space, \(r_{<j,j} = (\Theta P_{<j})^T s_j\);

• form the residuals, \(\tilde{s}_j = s_j - (\Theta P_{<j})\, r_{<j,j}\) and \(\tilde{p}_j = k_j - P_{<j}\, r_{<j,j}\);

• normalize using the sketched norm, \(r_{jj} = \left\lVert\tilde{s}_j\right\rVert_2\), \(\Theta p_j = \tilde{s}_j / r_{jj}\), \(p_j = \tilde{p}_j / r_{jj}\).

By construction \(\Theta P\) has orthonormal columns. In our experiments, RGS keeps the process numerically stable even when \(\Theta K\) is close to rank-deficient, and we apply a single rCholQR refinement to the resulting basis to enforce sketch-orthogonality up to machine precision [24], [30].

For a structured random sketch such as the subsampled randomized Hadamard transform, computing all \(r\) sketches \(s_j = \Theta k_j\) costs \(O(m r \log \ell)\) in total (or \(O(m \ell r)\) for a dense Gaussian sketch). The remaining per-step work is dominated by the long-vector update \(\tilde{p}_j = k_j - P_{<j}\, r_{<j,j}\), costing \(2m(j-1)\) flops at step \(j\); the sketch-side operations are \(O(\ell r)\) per step. With the typical sketch size \(\ell = O(r)\), summing over \(j = 1, \dots, r\) and ignoring lower-order terms yields a total cost of \(mr^2\) flops for RGS. In addition, RGS can exploit mixed precision. It is shown in [24] that the operation \(P_{<j}\, r_{<j,j}\) can be performed with lower precision while still obtaining a well-conditioned basis \(P\). We analyze the cost in more detail in Section 5.2, and refer to [24] for the cost and stability analysis of RGS.

3 Orthogonal sketch DLRA↩︎

In this section, we propose the orthogonal sketch DLRA: we evolve sketch-orthogonal bases but keep standard orthogonal projections in the dynamics. The sketch is used only for the orthogonalization step (replacing Householder QR with RGS), not for the projection itself. As we shall see in Section 4, this is in contrast with the more natural approach of sketching the Galerkin condition, which leads to an oblique projection that can fail with a large perpendicular residual.

3.1 Sketch-orthogonal manifold and bijective map↩︎

We define the set of fixed rank \(r\) matrices with sketch-orthogonality conditions (w.r.t.\(\Theta\) and \(\Omega\)) as \[\mathcal{M}_r^{\Theta, \Omega} = \left\{ P S W^T \mid S \in \mathrm{GL}(r), P \in \mathrm{St}_{\Theta}(m, r), W \in \mathrm{St}_{\Omega}(n, r) \right\},\] and its tangent space at \(Y = PSW^T \in \mathcal{M}_r^{\Theta, \Omega}\) as \[\mathcal{T}_Y \mathcal{M}_r^{\Theta, \Omega} = \left\{ \delta P S W^T + P \delta S W^T + P S \delta W^T \;\middle|\; \begin{align} & \delta S \in \mathbb{R}^{r \times r},\;\delta P \in \mathbb{R}^{m \times r},\;\delta W \in \mathbb{R}^{n \times r} \\ & \delta P^T \Theta^T \Theta P =0, \;\delta W^T \Omega^T \Omega W =0 \end{align} \right\}.\] The special case \(\Theta = I_m\) and \(\Omega = I_n\) recover the standard manifold \(\mathcal{M}_r\) and its tangent space. For a rectangular sketch (\(\ell < m, n\)), the identification with \(\mathcal{M}_r\) is pointwise rather than global: the inclusion \(\mathcal{M}_r^{\Theta, \Omega} \subseteq \mathcal{M}_r\) always holds, and a matrix \(Y \in \mathcal{M}_r\) admits a sketch-orthogonal factorization exactly when the sketches are injective on its column and row spaces, that is, \[\mathcal{M}_r^{\Theta, \Omega} \;=\; \left\{ Y \in \mathcal{M}_r\;\middle|\; \Theta \text{ injective on } \mathrm{Range}(Y) \text{ and } \Omega \text{ injective on } \mathrm{Range}(Y^T) \right\}.\] By Assumption 1, every iterate arising in this paper belongs to \(\mathcal{M}_r^{\Theta, \Omega}\), and at every such point, the tangent spaces coincide as subsets of \(\mathbb{R}^{m \times n}\): \[\mathcal{T}_Y \mathcal{M}_r^{\Theta, \Omega} = \mathcal{T}_Y \mathcal{M}_r\qquad \text{for all } Y \in \mathcal{M}_r^{\Theta, \Omega}.\] The sketch-orthogonality conditions only affect the parameterization, not the underlying geometric objects. This equivalence is made precise by the following map \(\varphi\): \[\begin{align} \varphi : \quad &\mathrm{St}_{\Theta}(m, r) \times \mathrm{GL}(r) \times \mathrm{St}_{\Omega}(n, r) \rightarrow \mathrm{St}(m, r) \times \mathrm{GL}(r) \times \mathrm{St}(n, r) \\ &\left( P,\;S,\;W \right) \mapsto \left( P (P^T P)^{-1/2},\;(P^T P)^{1/2} S (W^T W)^{1/2},\;W (W^T W)^{-1/2} \right), \end{align}\] which is a bijection onto its image — the orthogonal triples whose column and row spaces are embedded by \(\Theta\) and \(\Omega\), with an inverse map \[\begin{align} \varphi^{-1} : \quad &\operatorname{Im}\varphi \rightarrow \mathrm{St}_{\Theta}(m, r) \times \mathrm{GL}(r) \times \mathrm{St}_{\Omega}(n, r) \\ &\left( U,\;\Sigma,\;V \right) \mapsto \left( U (U_{\Theta}^T U_{\Theta})^{-1/2},\;(U_{\Theta}^T U_{\Theta})^{1/2} \Sigma (V_{\Omega}^T V_{\Omega})^{1/2},\;V (V_{\Omega}^T V_{\Omega})^{-1/2} \right), \end{align}\] where \(U_{\Theta} = \Theta U\), \(V_{\Omega} = \Omega V\), and \((\cdot)^{1/2}\) denotes the symmetric positive square root of a symmetric positive definite matrix. The domain \(\operatorname{Im}\varphi = \{(U, \Sigma, V) \in \mathrm{St}(m, r) \times \mathrm{GL}(r) \times \mathrm{St}(n, r) : \Theta U \text{ and } \Omega V \text{ have full column rank}\}\) is exactly the set of orthogonal triples whose column and row spaces are embedded by \(\Theta\) and \(\Omega\); this full-column-rank condition is what makes \(U_{\Theta}^T U_{\Theta}\) and \(V_{\Omega}^T V_{\Omega}\) invertible, so that \(\varphi^{-1}\) is well defined. In particular, given a sketch-orthogonal factorization \(Y = PSW^T\), the corresponding orthogonal factorization is \(\varphi(Y) = U \Sigma V^T\) with \[U = P (P^T P)^{-1/2}, \quad \Sigma = (P^T P)^{1/2} S (W^T W)^{1/2}, \quad V = W (W^T W)^{-1/2}.\]

3.2 Orthogonal sketch DLRA↩︎

With the bijective map in hand, we can express the standard orthogonal projection 2 in terms of the sketch-orthogonal bases \(P\) and \(W\). Since \(\mathrm{Range}(Y) = \mathrm{Range}(P) = \mathrm{Range}(U)\) and \(\mathrm{Range}(Y^T) = \mathrm{Range}(W) = \mathrm{Range}(V)\), we have \(P(P^TP)^{-1}P^T = UU^T\) and \(W(W^TW)^{-1}W^T = VV^T\), so that \[\label{eq:orthogonal95proj95sketch} \mathcal{P}_{Y}Z = P (P^TP)^{-1} P^T Z - P (P^TP)^{-1} P^T Z W (W^T W)^{-1} W^T + Z W (W^T W)^{-1} W^T.\tag{6}\] This is exactly the same projection as 2 , only rewritten in terms of the sketch-orthogonal factors: the sketch matrices \(\Theta\) and \(\Omega\) enter only implicitly through the parameterization.

We then define the orthogonal sketch DLRA as the differential equation \[\label{eq:ortho95sketch95dlra} \dot{Y}(t) = \mathcal{P}_{Y(t)} F(Y(t)), \quad Y(0) = P_0 S_0 W_0^T \in \mathcal{M}_r^{\Theta, \Omega},\tag{7}\] where \(\mathcal{P}_Y\) is the standard \(\ell_2\)-orthogonal projection 6 . The sketch-orthogonality of \(P\) and \(W\) is maintained via RGS during time stepping, but the dynamics are governed by orthogonal projections. Equation 7 therefore coincides with the classical DLRA equation 4 up to a reparameterization via \(\varphi\), so the classical DLRA error bound of [1] applies verbatim, without any sketch-induced constant.

A remark on the gauge is in order. Since \(\mathcal{T}_Y\mathcal{M}_r^{\Theta,\Omega}=\mathcal{T}_Y\mathcal{M}_r\) at each \(Y\in\mathcal{M}_r^{\Theta,\Omega}\) (Section 3.1), the trajectory \(Y(t)\) of 7 is determined only by the orthogonal projection \(\mathcal{P}_Y\); a gauge choice only fixes how a tangent vector \(\dot{Y}\) is split into the factor velocities \((\dot{P},\dot{S},\dot{W})\) and leaves \(Y(t)\) unchanged. We therefore adopt in Proposition 3 the Euclidean gauge \(P^T\dot{P}=0\), \(W^T\dot{W}=0\) — rather than the sketch gauge \(\dot{P}^T\Theta^T\Theta P=0\), \(\dot{W}^T\Omega^T\Omega W=0\) built into \(\mathcal{T}_Y\mathcal{M}_r^{\Theta,\Omega}\) — precisely because it reduces the coupled equations ?? to the classical form 5 . This gauge does not keep \(P,W\) sketch-orthogonal along the continuous flow, but that is immaterial: the integrators of Section 5 re-orthogonalize the bases by RGS at every macro step, restoring the sketch-orthogonal parameterization at the discrete points where it is actually used.

4 Oblique tangent space projection↩︎

A natural approach to incorporate sketching into DLRA is to sketch the Galerkin condition 3 that defines the approximation, in the spirit of randomized Galerkin methods for model reduction [31], [32] or the Petrov–Galerkin condition in Krylov subspace methods [24], [33]. For \(Y = PSW^T \in \mathcal{M}_r^{\Theta, \Omega}\) (see Section 3.1), it is natural to impose the following sketch Galerkin condition: \[\begin{align} \label{eq:sketch95galerkin95oblique} &\left\langle \Theta (\dot{Y} - F(Y)) \Omega^T, \, \Theta \delta Y \Omega^T \right\rangle = 0, \qquad \delta Y \in \mathcal{T}_{Y(t)} \mathcal{M}_r^{\Theta, \Omega}. \end{align}\tag{8}\] In this section, we show that this sketch Galerkin condition is the classical DLRA transported to the sketch space (Lemma 1); from it we derive the resulting oblique tangent space projection as the full-space lift, show that it is equivalent to the classical DLRA exactly when a sketched orthogonality condition holds (Proposition 5), met up to the sketching error when the problem is low-rank compatible in the sense of Definition 2, and analyze how the dynamics deviate otherwise.

4.0.0.1 Standing conventions for this section.

Throughout Section 4 we keep \(P \in \mathrm{St}_\Theta(m, r) \subset \mathbb{R}^{m \times r}\) and \(W \in \mathrm{St}_\Omega(n, r) \subset \mathbb{R}^{n \times r}\), so that \(\Theta P \in \mathbb{R}^{\ell \times r}\) and \(\Omega W \in \mathbb{R}^{\ell \times r}\) are tall (\(\ell \geq r\)) with orthonormal columns. For these to be well defined, \(\Theta P\) and \(\Omega W\) must have full column rank \(r\), which holds by Assumption 1. Under sketch-orthogonality, Remark 2 gives \[\label{eq:pinv95convention} (\Theta P)^+ \;=\; (\Theta P)^T \;=\; P^T \Theta^T, \qquad (\Omega W)^+ \;=\; (\Omega W)^T \;=\; W^T \Omega^T,\tag{9}\] where \(^+\) denotes the Moore–Penrose pseudo-inverse (cf. Remark 2). Since \(\Theta P\) and \(\Omega W\) are tall (\(\ell \geq r\)) with full column rank, \((\Theta P)^+\) and \((\Omega W)^+\) are genuine left pseudo-inverses; under sketch-orthogonality, they coincide with \((\Theta P)^T\) and \((\Omega W)^T\) respectively.

The analysis of Section 4 is pointwise in time: each statement fixes a state \(Y \in \mathcal{M}_r^{\Theta, \Omega}\) together with a single pair \((\Theta, \Omega)\), and concerns the oblique sketch DLRA equation at that instant. In the discrete integrator of Section 5, a fresh pair \((\Theta, \Omega)\) is drawn at every macro step, and the results of this section apply independently at each step.

4.1 Oblique coupled equations↩︎

We now show that the sketch Galerkin condition 8 is the classical DLRA condition transported to the sketch space: it drives the sketched state \(\Theta Y \Omega^T\) by the sketched vector field \(\Theta F(Y) \Omega^T\), as the next lemma makes precise.

In sketch coordinates, this oblique projection is the ordinary orthogonal one: for every \(Z \in \mathbb{R}^{m \times n}\), \[\Theta\,(\mathcal{P}_Y^{\Theta, \Omega} Z)\,\Omega^T = \mathcal{P}_{\Theta Y \Omega^T}\!\left(\Theta Z \Omega^T\right),\] where \(\mathcal{P}_{\Theta Y \Omega^T}\) denotes the standard orthogonal tangent-space projection at the rank-\(r\) matrix \(\Theta Y \Omega^T \in \mathbb{R}^{\ell \times \ell}\). This is the projection-level form of Lemma 1.

In the rest of the paper, the oblique sketch DLRA refers to the projected differential equation \(\dot{Y}(t) = \mathcal{P}_{Y(t)}^{\Theta, \Omega} F(Y(t))\) with \(Y(0) = Y_0 \in \mathcal{M}_r^{\Theta, \Omega}\) — that is, the sketched Galerkin condition 8 closed by the canonical lift of Lemma 1.

4.2 Equivalence with the classical DLRA in the low-rank compatible regime↩︎

We now discuss the conditions under which the oblique sketch DLRA is (approximately) equivalent to the classical DLRA, and what happens when these conditions fail.

We now turn to the large-residual regime: the oblique sketch DLRA then solves a perturbed equation whose residual does not vanish, and is therefore no longer a faithful approximation of the classical DLRA.

Following the decomposition used in the proof of Proposition 5, we write \(F(Y(t)) = \mathcal{P}_{Y(t)} \left[ F(Y(t)) \right] + \mathcal{P}_{Y(t)}^{\perp} \left[ F(Y(t)) \right]\) and recast the sketch Galerkin condition 8 as \[\label{eq:oblique95residual} \left\langle \Theta (\dot{Y}(t) - \mathcal{P}_{Y(t)} \left[ F(Y(t)) \right]) \Omega^T, \, \Theta \delta Y \Omega^T \right\rangle = \left\langle \Theta \mathcal{P}_{Y(t)}^{\perp} \left[ F(Y(t)) \right] \Omega^T, \, \Theta \delta Y \Omega^T \right\rangle, \quad \forall \delta Y \in \mathcal{T}_{Y(t)} \mathcal{M}_r.\tag{10}\] For each fixed tangent direction \(\delta Y\), the left-hand side of 10 pairs two tangent vectors of rank at most \(2r\); whenever \(\Theta\) and \(\Omega\) embed their combined column and row spaces (of dimension at most \(4r\)), Corollary 1 (with \(d = 4r\)) shows that it reproduces the classical Galerkin pairing \(\left\langle \dot{Y}(t) - \mathcal{P}_{Y(t)} \left[ F(Y(t)) \right], \, \delta Y \right\rangle\) that drives the DLRA, up to an additive error \((2\epsilon + \epsilon^2) \left\lVert\dot{Y}(t) - \mathcal{P}_{Y(t)} \left[ F(Y(t)) \right]\right\rVert_F \left\lVert\delta Y\right\rVert_F\). The right-hand side is what we call the oblique residual: when the problem is \((r,k)\)-low-rank compatible, \(\mathcal{P}_{Y(t)}^{\perp} \left[ F(Y(t)) \right]\) has rank at most \(k\), so for each fixed tangent direction \(\delta Y\) the relevant column and row spaces have dimension at most \(2r+k\) and an \(\mathrm{OSE}(2r+k)\) controls this residual through the two-sided embedding (Corollary 1, \(d = 2r+k\)), with probability \(1 - 2\delta\): explicitly, since \(\left\langle \mathcal{P}_{Y(t)}^{\perp} \left[ F(Y(t)) \right], \, \delta Y \right\rangle = 0\), \[\bigl|\left\langle \Theta \mathcal{P}_{Y(t)}^{\perp} \left[ F(Y(t)) \right] \Omega^T, \, \Theta \delta Y \Omega^T \right\rangle\bigr| \;\leq\; (2\epsilon + \epsilon^2)\, \left\lVert\mathcal{P}_{Y(t)}^{\perp} \left[ F(Y(t)) \right]\right\rVert_F\, \left\lVert\delta Y\right\rVert_F.\] Otherwise, \(\mathcal{P}_{Y(t)}^{\perp} \left[ F(Y(t)) \right]\) is large relative to what a sketch of size \(\ell \ll m, n\) can control, so the oblique residual on the right-hand side of 10 is in general a non-vanishing perturbation of the classical Galerkin condition. This perturbation drives the dynamics away from the classical DLRA solution and is confirmed numerically in Section 6.3, where the oblique approach leads to a visible energy drift on the Vlasov–Poisson equation.

4.3 Recap: oblique vs.orthogonal↩︎

The key differences between the two approaches are summarized in Table 2.

In summary, the oblique approach is cheaper but only conditionally and approximately correct, while the orthogonal approach is always correct at the price of an extra randomized Gram–Schmidt per factor step. For problems known to be approximately low-rank compatible, the oblique approach is a reasonable choice; otherwise, the orthogonal approach should be preferred.

5 Sketch DLRA integrators↩︎

Classical DLRA integrators — projector-splitting (KSL) [6], BUG [7], and its augmented and rank-adaptive variants [8] — are built from K/L/S substeps composed with Householder QR extractions. Each lifts to a sketch integrator by (i) replacing the substep projectors as in Table 3 and (ii) replacing the Householder QR by the randomized Gram–Schmidt procedure of Section 2.3. The splitting order, Galerkin transfer, rank truncation, and adaptivity are inherited unchanged. Throughout, \(h > 0\) is the time step, the current iterate \(Y_0 = P_0 S_0 W_0^T\) is sketch-orthogonal, and we use the dot notation \(\dot{x}\) to denote the time derivative of \(x\). In general, a fresh independent pair \((\Theta, \Omega)\) may be drawn at each macro step, the state being re-expressed in the new sketch-orthogonal bases by RGS before the substeps are taken; in the experiments of Section 6 we instead hold a single pair \((\Theta, \Omega)\) fixed over the whole trajectory.

5.1 Sketch substeps and integrators↩︎

In each substep, one factor is evolved while the others are frozen; the projector acting on the right-hand side selects the orthogonal or oblique variant. The three substep ODEs and the RGS basis update are collected in Table 3; both columns collapse to the classical DLRA substeps when \(\Theta = I_m\) and \(\Omega = I_n\). The sketch-orthogonality of \(P\) and \(W\) makes \((\Theta P) \in \mathbb{R}^{\ell \times r}\) and \((\Omega W) \in \mathbb{R}^{\ell \times r}\) have orthonormal columns, so their left inverses \((\Theta P)^T = P^T \Theta^T\) and \((\Omega W)^T = W^T \Omega^T\) are well-defined (Remark 2) and the oblique formulas are unambiguous.

6pt

Every classical DLRA integrator lifts to the sketch setting by picking one column of Table 3 and using RGS for the basis extraction. Below, we report the two variants used in the experiments. The same recipe extends straightforwardly to the augmented BUG and rank-adaptive BUG integrators of [8], which we do not describe here.

5.1.0.1 Sketch projector-splitting integrator (sKSL).

A Lie–Trotter step performs the sequence K, \(-\)S, L on \([0, h]\): the K-substep starts from \(K(0) = P_0 S_0\) and produces \(P_1 = \mathop{\mathrm{RGS}}(K(h), \Theta)\); the S-substep uses the negated right-hand side; the L-substep ends with \(W_1 = \mathop{\mathrm{RGS}}(L(h), \Omega)\) [6]. A symmetric Strang splitting K/2, \(-\)S/2, L, \(-\)S/2, K/2 gives second order at the cost of one extra S-substep.

5.1.0.2 Sketch BUG integrator (sBUG).

The K- and L-substeps are solved in parallel from \((P_0, S_0, W_0)\), yielding \(P_1 = \mathop{\mathrm{RGS}}(K(h), \Theta)\) and \(W_1 = \mathop{\mathrm{RGS}}(L(h), \Omega)\) [7]. The S-substep is then run from the Galerkin transfer \(S_{\mathrm{init}} = M S_0 N^T\) with \[M = (P_1^T P_1)^{-1} P_1^T P_0, \quad N = (W_1^T W_1)^{-1} W_1^T W_0 \qquad \text{(orthogonal)},\] \[M = (\Theta P_1)^+ (\Theta P_0), \quad N = (\Omega W_1)^+ (\Omega W_0) \qquad \text{(oblique)}.\]

5.2 Advantages of sketch DLRA on modern hardware↩︎

At leading order, the two integrators have the same basis-update flop count, \(O(m r^2)\) vs.\(O(m \ell r)\) with \(\ell = O(r)\), so the practical speedup of the sketch approach is not explained by flop savings alone. It comes from three additional mechanisms that flop counts miss [13], [24], [34]: (i) fewer global synchronizations per basis update, (ii) suitability for exploiting mixed precision while maintaining numerical stability, and (iii) a shift from bandwidth-bound BLAS-2 panel factorizations to compute-bound BLAS-3 matrix multiplications for well-conditioned bases. We therefore report per-step cost along three axes: flops, global synchronizations among processors, and BLAS level of the dominant kernel. Throughout, the tall factors \(K \in \mathbb{R}^{m \times r}\) and \(L \in \mathbb{R}^{n \times r}\) are assumed row-distributed across processors, and \(F\) is assumed to factorize through the low-rank structure so that the cost of computing \(F(P S W^T) W\) and \(F(P S W^T)^T P\) stays contained, as for the test problems of Section 6.

5.2.0.1 Standard basis update.

Different algorithms can be used to compute an orthogonal basis, such as Householder QR, classical Gram–Schmidt, or modified Gram–Schmidt. For a row-distributed \(K\), the column-by-column Householder QR needs two global synchronizations per column to compute the factorization [13] and \(r\) synchronizations to explicitly compute the \(Q\) factor from its compact representation, for a total of \(3r\) synchronizations. The dominant kernel is a BLAS-2 operation on the full row-dimension \(m\). The flop count is \(4 m r^2\), higher than classical and modified Gram–Schmidt, as the basis needs to be explicitly constructed from its Householder-based representation. Communication-avoiding variants such as TSQR and CAQR collapse this to a single AllReduce, but require a bespoke reduction tree on Householder factors that most high-level scientific-computing stacks do not expose. Modified Gram–Schmidt requires \(j-1\) synchronizations to compute column \(j\) for a total of \(O(r^2)\) synchronizations: column \(j\) cannot start until the inner products and norm computed from column \(j-1\) have been combined between processors. Although classical Gram–Schmidt requires \(r\) synchronizations, it is known to be unstable in practice. Both classical and modified Gram–Schmidt require \(2m r^2\) flops.

5.2.0.2 Sketch basis update.

RGS has the same communication cost as classical Gram–Schmidt and, on the long-vector updates, half the flops: \(r\) synchronizations and \(mr^2\) flops, to which the sketch step adds \(O(mr\log\ell)\) for an SRHT sketch (sub-leading) or \(O(m\ell r) \approx 2mr^2\) for a dense Gaussian sketch (same order as Householder QR). However, it is shown to be numerically as stable as modified Gram–Schmidt and is also suitable for using mixed precision [24]. When \(K\) is well conditioned and randomized Cholesky QR can be used, the column-by-column process is replaced by a single sketch of the whole tall factor followed by a small QR:

1. \(\Theta K \in \mathbb{R}^{\ell \times r}\) — one BLAS-3 tall-and-skinny matmul; one global reduction of an \(\ell \times r\) block.

2. Thin QR of \(\Theta K\) — \(O(\ell r^2)\) flops, entirely local.

3. \(P = K R^{-1}\) — BLAS-3 matmul, embarrassingly parallel on the row-distributed \(K\).

The total per-basis synchronization count is \(1\), independent of \(r\), and every kernel is BLAS-3. Step 1 costs \(O(m \ell r)\) flops for a dense Gaussian sketch; subsampled randomized Hadamard transforms reduce this to \(O(m r \log \ell)\) [22], [24]. On GPUs, where BLAS-2 panel kernels are memory-bound, and on TPUs and tensor cores, which expose only batched matrix products, this BLAS-3 rCholQR path is the precondition for the basis update to approach peak hardware throughput.

5.2.0.3 Per-step totals.

The sKSL integrator runs the K- and L-basis updates sequentially (the S-substep needs \(P_1\) before L starts), so its per-step synchronization count is \(O(r^2)\) for modified Gram–Schmidt, \(O(r)\) for Householder QR, classical Gram–Schmidt, and RGS; and \(O(1)\) when randomized Cholesky QR is used in place of RGS. The sBUG variant dispatches the K- and L-updates to separate devices (Remark 8), which halves both counts on two devices; the rCholQR total stays \(O(1)\) independently of \(r\). On distributed or multi-GPU hardware where reduction latency is a material fraction of the step cost, this is a real advantage at the ranks \(r \in [10, 100]\) typical of DLRA applications [24], [34].

As a disclaimer, our experiments (Section 6) run on a single CPU on numpy/scipy at small problem sizes where reduction latency is negligible; they validate the theory of Sections 3 and 4 rather than measure wall-clock speedup. A speedup is expected only on large problems with \(m, n\) in the range \(10^5\)–\(10^7\), especially for the parallel BUG and augmented-BUG variants, or in a higher-dimensional tensor setting (see future directions in the conclusion section).

The bases \(P\) and \(W\) produced by RGS with \((r, \epsilon, \delta)\)-OSE sketches are also nearly Euclidean-orthonormal — \(\tfrac{1}{1+\epsilon} I \preceq P^T P,\, W^T W \preceq \tfrac{1}{1-\epsilon} I\) with probability at least \(1-\delta\) per basis (\(1-2\delta\) jointly, \(\Theta\) and \(\Omega\) being independent) — so the inverses \((P^T P)^{-1}\) and \((W^T W)^{-1}\) that appear in the orthogonal substeps of Table 3 are well-conditioned and can be computed stably via Cholesky; see [24] for the detailed RGS stability analysis.

6 Numerical experiments↩︎

We compare three families of integrators on three test problems that span the spectrum of low-rank compatibility. The integrators are: (i) the classical DLRA BUG method of [7] as an unsketched baseline; (ii) the oblique sketch DLRA variants of KSL and BUG analyzed in Section 4; and (iii) the orthogonal sketch DLRA variants of KSL and BUG of Section 3. All experiments are run on a MacBook Pro with an Apple M1 chip and \(16\,\)GB of RAM. The implementation is in Python and is publicly available on GitHub³. The sketch matrices are Gaussian \(\Theta \in \mathbb{R}^{\ell \times m}\) and \(\Omega \in \mathbb{R}^{\ell \times n}\), with sketch size \(\ell\) reported for each experiment. All sketch-based methods produce their sketch-orthogonal bases via the randomized Gram–Schmidt (RGS) process of [24], followed by a single randomized Cholesky QR refinement to enforce sketch-orthogonality up to machine precision [30] (Section 2.3). A single pair \((\Theta, \Omega)\) is drawn and held fixed over the whole trajectory (rather than resampled at each macro step as in Section 5), and the statistics below are taken over independent such draws. For each sketch method, the reported curve is the median relative Frobenius error \[\mathrm{err}(t) \;:=\; \frac{\left\lVert Y(t) - A_{\mathrm{ref}}(t)\right\rVert_F}{\left\lVert A_{\mathrm{ref}}(t)\right\rVert_F}\] over \(20\) independent draws of \((\Theta, \Omega)\) (\(5\) for Vlasov–Poisson, see Section 6.3); the shaded band is the interquartile range. The reference trajectory \(A_{\mathrm{ref}}(t)\) is obtained by integrating the full-order problem 1 with scipy’s adaptive RK45 at tight tolerance. As a benchmark for the best achievable rank-\(r\) approximation we also plot \(\left\lVert\sigma_r(A_{\mathrm{ref}}(t))\right\rVert_F / \left\lVert A_{\mathrm{ref}}(t)\right\rVert_F\), the relative error of the best rank-\(r\) SVD truncation. The experiments are organized from the most favorable case (Allen–Cahn, approximately LRC) to the most adversarial case (Vlasov–Poisson, large perpendicular residual).

6.1 Allen–Cahn equation↩︎

The Allen–Cahn equation [35] describes the process of phase separation in multi-component alloy systems. In its simplest form, it is expressed as the partial differential equation \[\label{eq:allen95cahn} \begin{align} &\dot{u}(t, x, y) = \varepsilon \Delta u(t, x, y) + u(t, x, y) - u(t, x, y)^3, &&(x,y) \in \mathcal{D} \subset \mathbb{R}^2, \quad t \in [0, T], \\ &u \text{ periodic on } \partial \mathcal{D}, &&t \in [0, T], \\ &u(0, x, y) = u_0(x, y) &&(x, y) \in \mathcal{D}, \end{align}\tag{11}\] where \(\Delta u\) is the diffusion term, \(u - u^3\) is the reaction term, and \(\varepsilon\) is a small parameter that controls the stiffness of the equation. By discretizing 11 using second-order centered finite differences on a uniform tensor-product grid with \(n\) points per direction, we obtain the Sylvester-like matrix differential equation \[\label{eq:allen95cahn95discrete} \dot{X}(t) \;=\; A\, X(t) + X(t)\, A^T + X(t) - X(t)^{*3},\tag{12}\] where \(A \in \mathbb{R}^{n \times n}\) is \(\varepsilon\) times the symmetric one-dimensional finite-difference Laplacian and the power is taken element-wise (\(X^{*3} = X * X * X\), with \(*\) the Hadamard product). The diffusion \(A X + X A^T\) preserves rank-\(r\) structure exactly, and for smooth initial data the reaction \(X - X^{*3}\) acts as a low-rank perturbation along the trajectory. This places Allen–Cahn squarely in the approximately \((r, k)\)-low-rank compatible regime for a small \(k\), and the low-rank-compatibility analysis of Section 4.2 predicts that the oblique sketch DLRA should perform comparably to its orthogonal counterpart on this approximately-LRC problem.

We run 12 on \(\mathcal{D} = [0, 2\pi]^2\) with \(n = 64\), \(\varepsilon = 0.01\), rank \(r = 8\), time horizon \(t \in [0, 10]\), a DLRA macro step \(\Delta t = 1/60\), and a solution snapshot saved every \(10\) DLRA steps (i.e., every \(1/6\) time units, for \(61\) snapshots in total). Within each DLRA step, the K, L, and S sub-problems are solved with scipy’s adaptive RK45 at relative and absolute tolerance \(10^{-12}\). Figure 2 shows the relative error of the oblique and orthogonal sketch variants at sketch size \(\ell = 2r = 16\), compared against the unsketched BUG baseline and the best rank-\(r\) SVD truncation of the reference trajectory. Both sketch families produce sensible low-rank approximations of the reference solution. The oblique variant tracks the orthogonal one in the median, confirming the LRC prediction — on a problem where the perpendicular residual is near low-rank, the equivalence regime of Section 4.2 applies — but with a visibly wider interquartile band, the first sign of the run-to-run sensitivity that becomes severe on the Vlasov–Poisson problem, where the perpendicular residual is large.

6.2 Fokker–Planck equation↩︎

The Fokker–Planck equation describes the evolution of the probability density function (PDF) of the \(d\)-dimensional state \(X_t\) of the Itô stochastic differential equation \[d X_t \;=\; \mu(X_t)\, dt + \sigma\, dW_t,\] where \(\mu : \mathbb{R}^d \to \mathbb{R}^d\) is the drift, \(\sigma > 0\) is a constant diffusion coefficient, and \(W_t\) is the \(d\)-dimensional standard Wiener process. The associated PDF \(f(t, x)\) satisfies \[\label{eq:fokker95planck} \partial_t f(t, x) \;=\; - \sum_{i=1}^d \partial_{x_i}\!\bigl(\mu_i(x)\, f(t, x)\bigr) + \frac{\sigma^2}{2} \sum_{i=1}^d \partial^2_{x_i} f(t, x),\tag{13}\] with initial condition \(f(0, x) = f_0(x)\).

Low-rank methods for the Fokker–Planck PDF are a natural setting: closely related are the dynamically orthogonal (DO) field equations [36], [37] and the rank-adaptive tensor approach of [5], whose drift we adopt below. Following [5], we take \(\sigma = 2\) and the drift \[\mu_i(x) \;=\; \bigl(\gamma(x_{i+1}) - \gamma(x_{i-2})\bigr)\, \xi(x_{i-1}) - \phi(x_i), \qquad i = 1, \ldots, d,\] with indices interpreted modulo \(d\) and \(\gamma, \xi, \phi\) the \(2\pi\)-periodic functions \(\gamma(x) = \sin(x)\), \(\xi(x) = \cos(x)\), \(\phi(x) = \exp(\sin(x)) + 1\). We specialize to \(d = 2\) on the flat torus \(\mathcal{D} = [0, 2\pi]^2\) with the initial condition \[f_0(x_1, x_2) \;=\; \frac{1}{m_0}\Bigl(e^{\sin(x_1 - x_2)^2} + \sin(x_1 + x_2)^2\Bigr),\] where \(m_0\) normalizes \(f_0\) to unit Frobenius norm on the discretization grid.

Discretizing 13 with \(n\) equispaced points per direction and centered finite differences, we obtain the Sylvester-like matrix differential equation \[\label{eq:fokker95planck95discrete} \dot{X}(t) \;=\; D_{xx}\, X(t) + X(t)\, D_{yy}^T - D_x\bigl(M_1 * X(t)\bigr) - \bigl(M_2 * X(t)\bigr)\, D_y^T, \qquad X(0) = X_0,\tag{14}\] where \(D_{xx}, D_{yy} \in \mathbb{R}^{n \times n}\) are the scaled periodic 1D Laplacians, \(D_x, D_y \in \mathbb{R}^{n \times n}\) are periodic centered first-difference operators, \(*\) is the Hadamard product, and \(M_1, M_2 \in \mathbb{R}^{n \times n}\) are the drift coefficients \(\mu_1, \mu_2\) evaluated on the grid.

Low-rank compatibility is less clear-cut here than for Allen–Cahn: the diffusion \(D_{xx} X + X D_{yy}^T\) is rank-preserving, but the drift matrices \(M_1, M_2\) are not separable in \((x, y)\), so the Hadamard products \(M_1 * X\) and \(M_2 * X\) generically raise the rank of \(X\). In practice, however, \(M_1\) and \(M_2\) built from the smooth periodic functions \(\gamma, \xi, \phi\) are exactly rank-\(3\) matrices⁴ This places Fokker–Planck in the approximate-LRC regime, and the analysis of Section 4.2 predicts that the oblique variant should produce sensible results, as for Allen–Cahn.

We run 14 with \(n = 64\), rank \(r = 8\), time horizon \(t \in [0, 1]\), a DLRA macro step \(\Delta t = 1/600\), and a solution snapshot saved every \(10\) DLRA steps (i.e., every \(1/60\) time units, for \(61\) snapshots in total). Within each DLRA step, the K, L, and S sub-problems are solved with scipy’s adaptive RK45 at relative and absolute tolerance \(10^{-12}\). Figure 3 shows the relative error of the oblique and orthogonal sketch variants at sketch size \(\ell = 2r = 16\), compared against the unsketched BUG baseline and the best rank-\(r\) SVD truncation of the reference trajectory. The orthogonal sketch variants (left panel) concentrate tightly around the unsketched BUG error — the interquartile band is essentially invisible over \(20\) Gaussian sketches — and the orthogonal KSL variant is actually slightly below BUG, reflecting the standard ordering between projector-splitting and BUG integrators rather than a sketching effect. The oblique sketch variants (right panel) also compute a good low-rank approximation, but exhibit a clearly visible interquartile band throughout the trajectory, including a wider spread on the final error than their orthogonal counterparts.

6.3 Vlasov–Poisson equation↩︎

The Vlasov–Poisson equation models the evolution of an electron distribution in a collisionless plasma, in the electrostatic regime with a uniform neutralizing ion background: \[\label{eq:vlasov95poisson} \begin{align} &\partial_t f + v \cdot \nabla_x f - E(f) \cdot \nabla_v f = 0, \\ &\nabla_x \cdot E(f) = 1 - \int f\, dv, \qquad \nabla_x \times E(f) = 0, \end{align}\tag{15}\] where \(f(t, x, v) \ge 0\) is the particle density and \(E(f)\) is the self-consistent electric field, with periodic boundary conditions in \(x\). DLRA methods for this equation have been developed in [2], [3], [38]. We specialize to \(d = 1\) and discretize on an \(n_x \times n_v\) grid using second-order centered finite differences; the unknown matrix \(X(t) \in \mathbb{R}^{n_x \times n_v}\) satisfies \[\label{eq:vlasov95poisson95discrete} \dot{X}(t) \;=\; -\, D_x\, X(t)\, \Lambda_v \;+\; E\bigl(t, X(t)\bigr)\, X(t)\, D_v^T, \qquad X(0) = X_0,\tag{16}\] with \(D_x, D_v\) the periodic centered first-difference operators, \(\Lambda_v = \operatorname{diag}(v_j)\) the velocity grid, and \(E(t, X) = \operatorname{diag}(E_1, \ldots, E_{n_x})\) the diagonal matrix of discrete electric-field values. In the nonlinear regime the solution is no longer well approximated on the rank-\(r\) manifold, so the perpendicular component \(\mathcal{P}_{Y}^{\perp} \left[ F(Y) \right]\) becomes large and the oblique residual of 10 amplifies quickly instead of (approximately) vanishing; the equivalence of Proposition 5 is then lost. This is the adversarial test case that motivates orthogonal sketch DLRA over oblique sketch DLRA.

We run 16 on the standard two-stream instability [39], [40] with \((\mathcal{D}_x, \mathcal{D}_v) = ([0, 10\pi], [-6, 6])\), \(n_x = n_v = 64\), rank \(r = 8\), time horizon \(t \in [0, 60]\), and initial condition \[f_0(x, v) \;=\; \frac{1}{2\sqrt{2\pi}}\Bigl(e^{-(v - v_0)^2 / 2} + e^{-(v + v_0)^2 / 2}\Bigr)\bigl(1 + \alpha \cos(k x)\bigr), \qquad v_0 = 2.4,\;k = 0.2,\;\alpha = 10^{-3},\] producing two counter-streaming beams perturbed by a small linear mode. We use a DLRA macro step \(\Delta t = 1/20\) with scipy adaptive RK45 substeps at tolerance \(10^{-12}\). Figure 4 reports the relative error and the electric energy \(\tfrac{1}{2}\int |E(t, x)|^2 dx\) at sketch size \(\ell = 2r = 16\), each curve a median over \(5\) Gaussian draws (fewer than for Allen–Cahn and Fokker–Planck because the failure is already unambiguous; individual oblique runs trigger inner-integrator failures and yield truncated curves). The orthogonal variants (blue) track the unsketched BUG baseline within a constant factor across the full window, including through the nonlinear saturation; the oblique variants (red) diverge exponentially, with the oblique KSL variant eventually causing the inner integrator to fail.

The failure is structural. Vlasov–Poisson is Hamiltonian and conserves total energy \(E_{\mathrm{kin}} + E_{\mathrm{el}}\) and total mass \(\iint f\,dx\,dv\), and the classical DLRA inherits both up to the rank-\(r\) truncation error; at our parameters the orthogonal sketch BUG preserves them to order \(10^{-2}\) over \([0, 60]\) (Figure 5). The oblique sketch DLRA, by contrast, solves a perturbed ODE \(\dot{Y} = \mathcal{P}^{\Theta, \Omega}_{Y} F(Y)\) whose invariants depend on the random sketches rather than on the physics: total energy drifts by more than \(100\%\) before \(t = 40\), and total mass goes negative at \(t \approx 38\). The mechanism is the one diagnosed in Section 4.2: in the nonlinear regime the perpendicular residual \(\mathcal{P}_{Y}^{\perp} \left[ F(Y) \right]\) is large and far from low-rank, so no sketch of size \(\ell \ll m, n\) can make the oblique residual of 10 vanish, and the equivalence of Proposition 5 is lost. Increasing the sketch size within the low-rank budget does not rescue the oblique scheme: the failure is intrinsic to the oblique projection on this problem, not an artifact of the sketch size \(\ell = 2r\). The orthogonal sketch DLRA sidesteps both issues. Its tangent projection is the standard orthogonal projector 6 regardless of \((\Theta, \Omega)\) (Section 3), so the trajectory coincides with the classical DLRA solution and inherits its conservation behavior. The sketches still enter through the RGS basis extraction, but there an OSE plays only a stability role — keeping \(P^TP\) and \(W^TW\) well-conditioned in the orthogonal substeps of Table 3 — rather than the correctness role it had to play in the oblique scheme through Proposition 5.

7 Conclusion↩︎

The orthogonal sketch DLRA we propose applies to any problem (LRC or not) and inherits the classical DLRA error bound, while reducing synchronizations — and, for structured sketches, flops — on the basis-update step. Sketching the Galerkin condition that defines the DLRA leads to an oblique tangent projection whose dynamics deviate from the classical DLRA and fail on problems with a large perpendicular residual. The orthogonal sketch DLRA instead evolves sketch-orthogonal bases under standard orthogonal projections, preserving the geometric structure of the classical DLRA and replacing the only BLAS-2 kernel of the pipeline (Householder QR) by randomized Gram–Schmidt and, when applicable, randomized Cholesky QR, which together let the basis update fully exploit GPU hardware (Section 5.2). We derived sketch versions of the projector-splitting and BUG integrators and validated them on the Allen–Cahn, Fokker–Planck, and Vlasov–Poisson equations. The contribution is exploratory: our small-CPU experiments validate the theory rather than measure wall-clock speedup, which is expected only at large \(m, n\).

Three directions follow. First, subspace iteration could be considered to address cases with slowly decaying singular values, where randomized methods are typically less accurate. Second, the parallelism of the sketch BUG integrators could be exploited: the K- and L-substeps are independent, each reducing to an ODE solve and a Level-3 BLAS sketch multiplication, well-suited to a two-GPU pipeline. Third, the framework could be extended to Tucker and tensor-train formats, where the QR bottleneck is more pronounced and RGS extends naturally to each mode. The sketch integrators developed here also combine directly with the low-rank Parareal parallel-in-time integrator of [41], adding parallelism across the time direction.

References↩︎