PaperlayerView exact revision on arXiv

Source-linked AI summary

High-order robust basis-update & Galerkin integrators for dynamical low-rank approximation

Cory D. Hauck, Jonas Kusch, Steffen Schotthöfer

arXiv:2608.27749v1math.NA

TL;DR

Dynamical low-rank approximation requires integrators that remain accurate and stable for stiff dynamics, but existing high-order BUG methods rely on explicit schemes. This paper develops general-order augmented and parallel BUG integrators with flexible reduced-equation integration, proving order-p error bounds and demonstrating convergence and stability on stiff and non-stiff benchmarks.

  • Problem

    Existing high-order BUG integrators rely on explicit projected schemes, while stiffness can require specialized time integration and prohibitively small steps for stability.

  • Method

    The paper introduces augmented and parallel general-order BUG integrators that generate high-order bases independently of the reduced equation’s implicit or exponential integrator.

  • Results

    The integrators confirmed predicted convergence on non-stiff Schrödinger problems, while order-six augmented BUG remained stable and convergent at large steps for a stiff forced heat equation.

  • Takeaways & Limitations

    The construction supports high-order BUG integration in stiff settings and uses fewer basis vectors than the compared explicit RK-BUG integrator.

  • Takeaways & Limitations

    The framework leaves further BUG-p variants, including multistep constructions, for future development.

Abstract

from arXiv · show

Dynamical low-rank approximation has become a widely used numerical method in diverse disciplines. Its main idea is to represent the matrix or tensor-valued solution to a time-dependent differential equation as a low-rank factorization. The evolution of the factorization leads to highly stiff dynamics and requires the derivation of novel time integration methods that are not prone to this stiffness. A promising family of integrators are basis-update \& Galerkin (BUG) integrators as they enable implicit time integration and structure--preservation. However, current BUG integrators are limited to second--order accuracy while general-order BUG integrators are designed as projections of explicit time integration methods, thus severely limiting their use. In this work, we propose general-order BUG integrators that do not rely on an explicit time integration scheme while requiring a smaller number of basis functions to achieve high-order accuracy. We prove a general order error bound for the proposed augmented and parallel BUG integrators and demonstrate their behaviour for a series of stiff and non-stiff numerical benchmarks in which they significantly outperform previous BUG versions.

1 Introduction · 2 Low-rank Runge-Kutta integrators · 3 High-order BUG integrators

The paper develops arbitrary-order augmented and parallel BUG integrators that avoid explicit-scheme restrictions, use smaller bases than prior general-order approaches, and admit curvature-independent robust error bounds. These methods target stiffness in dynamical low-rank approximation while preserving the structural advantages of BUG integrators.

  • 1 Introduction: DLRA represents matrix or tensor solutions in low-rank factorizations, reducing computational complexity but introducing stiffness through the curvature of the low-rank manifold.The stiffness is particularly severe when the smallest singular value approaches zero, forcing prohibitively small time steps for standard schemes.
  • 2 Low-rank Runge-Kutta integrators: Existing general-order projected and Runge–Kutta BUG methods are limited by explicit time-step restrictions, large basis matrices, and inapplicability to stiff problems requiring implicit or exponential updates.The Runge–Kutta BUG basis may require 2rs* basis vectors, while its coefficient update is tied to an explicit projected Runge–Kutta method.
  • 3 High-order BUG integrators: The proposed high-order BUG integrators extend augmented and parallel BUG schemes to arbitrary order without relying on an underlying Runge–Kutta method, enabling implicit or exponential coefficient updates.They are designed for stiff problems and require smaller basis matrices than the approach in.
  • 3.1 Augmented BUG integrators: The augmented BUG construction builds bases first and then performs a sequential Galerkin coefficient update, with basis properties ensuring solution inclusion and consistency.The construction can use explicit low-storage schemes for augmentation, reducing the number of bases needed for order p.
  • 3.2 Parallel BUG integrators: The parallel BUG integrator pre-augments bases using an s-stage explicit Runge–Kutta method of order p and requires sr basis vectors before parallel coefficient solves.For p = s = 2, this choice resembles the second-order parallel integrator.
  • 3.3 Robust error bound: The integrators can also be made locally conservative for kinetic equations through conservative bases and truncation choices, including versions compatible with general time integration methods.For the parallel integrator, removing the coefficient update can directly yield a locally conservative scheme under the stated basis conditions.

4 Robust error bound for the high order augmented BUG integrator

Theorem 1 establishes a robust error bound for the order-p augmented BUG integrator by combining local-error estimates with stability of the exact problem. The analysis also discusses a refined bound when the solution is truncated to the original or higher rank.

  • Robust error bound: Theorem 1’s proof bounds the local error in three parts and derives the global error via Lady Windermere’s fan and exact-problem stability.The basis is assumed to satisfy properties P1 and P2, including a consistency relation with μ = h(h^p + hεr).
  • Robust error bound: R(t0) = 0 and R(t1) = O(μ) provide the residual estimates used in the subsequent error bound.The residual R is analyzed through its bounded p-th derivatives and the resulting defect estimate.
  • Robust error bound: The argument concludes that Y1 − A(t1) = O(μ), and then establishes the stated bound using Lipschitz continuity and Gronwall’s inequality.The defect is controlled using the Lipschitz constant L of F and the bound for R(t).
  • Rank truncation: When truncating to the original rank r or higher, the O(ϑ) term can be replaced by a term depending on the right-hand side’s normal components.The remark notes that adaptive rank is often required in practical applications.
  • Rank truncation: For the projected solution X(t1), the analysis uses ∥X(t1) − A(t1)∥ = O(hεr) to obtain the truncation-related bound.X(t) solves the projected differential equation with initial condition X(t0) = Y0.

5 Robust error bound for the high order parallel BUG integrator

The analysis establishes that the constructed bases approximate the full-rank and projected Runge–Kutta substeps accurately enough to yield a robust global error bound for the order-p parallel BUG integrator. It also shows that the S-step can be replaced by a cheaper locally conservative computation requiring only two parallel rank-(s−1)r solves.

  • Runge–Kutta substep approximation: Lemmas 3 and 4 show that the constructed bases accurately approximate the full-rank and projected Runge–Kutta substeps for j ≤ s − 1.The associated Runge–Kutta methods are introduced solely for numerical analysis and are not required in the actual K- and S-step computations.
  • Basis approximation: The constructed bases span the full-rank solution at t1 with sufficient accuracy under Assumptions A1–A3.This result is established in Lemma 5 and is a key ingredient in the robust error analysis.
  • Global error bound: The order-p parallel BUG integrator satisfies the one-step bound ∥bU bU^⊤ eA(t1)bV bV^⊤ − bY1∥ = O(h(h^p + ε^r)), with ∥bY1 − Y1∥ ≤ ϑ.The proof combines bounds for the S-step, K- and L-step contributions, and the remaining basis-projection term.
  • Reduced S-step: The S-step may be replaced by bS1 = bU0^⊤K(t1), requiring only two parallel rank-(s−1)r solves.Besides improving efficiency, this replacement yields a locally conservative integrator.

6 Implementation details: BUG-p integrator

The order-p augmented BUG integrator advances one step by assembling an augmented basis with an explicit s-stage Runge–Kutta constructor. The Galerkin S-step may use a different method, followed by rank truncations including a final tolerance-based truncation.

  • Time step: BUG-p advances the factorization from t_n to t_n+1 = t_n + h in one order-p augmented integrator step.The input is represented as Y_n = U_n S_n V_n^⊤ with U_n^⊤U_n = V_n^⊤V_n = I.
  • Basis construction: The augmented basis is assembled using an explicit s-stage Runge–Kutta method of order p solely to generate basis information.The Galerkin S-step is not required to use this Runge–Kutta method.
  • Truncation: The BUG-p step includes truncation to rank r_n and a final truncation with tolerance ϑ.The notation [[·]]_{r_n} denotes rank truncation, while [[·]]_ϑ denotes the final tolerance-based truncation.

1. RK basis construction. Compute augmented orthonormal bases

The RK basis construction uses coefficients from an explicit s-stage Runge–Kutta method of order p to form augmented bases. The method then projects the current approximation into the augmented space, whose tensor-product interpretation is an orthogonal projection.

  • RK basis construction: The basis construction uses coefficients from an explicit s-stage Runge–Kutta method of order p.
  • RK basis construction: The augmented bases have at most (s + 1)r_n columns without removing linear dependencies.The bound is written as r̄ ≤ (s + 1)r_n.
  • Projection step: The projection step projects the current approximation into the augmented space using r̄ × r_n matrices M = U.
  • Projection step: The resulting operation is the orthogonal projection of Y_n onto the tensor-product space spanned by U and V.

3. Galerkin S-step. Integrate the reduced matrix differential equation

The Galerkin S-step integrates the reduced matrix differential equation with a method chosen independently of the explicit Runge–Kutta basis constructor. To preserve order p, the reduced problem requires local error O(h^(p+1)), allowing explicit, implicit, exponential, or otherwise stiffly stable methods according to stiffness.

  • Galerkin S-step: The reduced-equation integrator is independent of the explicit Runge–Kutta method used in the basis constructor.This separates basis construction from reduced-problem time integration.
  • Galerkin S-step: Local error O(h^(p+1)) in the reduced equation is sufficient to retain order p.A method of order at least p for the reduced problem provides this requirement.
  • Galerkin S-step: Explicit methods suit non-stiff reduced problems, while implicit, exponential, or otherwise stiffly stable methods can address stiff reduced problems.The method choice is therefore adapted to the stiffness of the reduced equation.
  • RK basis constructor: Each Runge–Kutta basis-constructor stage augments the basis, applies a Galerkin projection to the stage combination, truncates to the current rank, and generates the next basis direction.The stage value produced by these operations is used to generate the next basis direction.
  • RK basis constructor: For ℓ < s, the next Runge–Kutta stage is computed by a BUG-type stage step.The constructor initializes the current factors and evaluates stage times and right-hand sides before computing stage directions.

1. Stage basis augmentation. Construct

The stage basis augmentation constructs orthonormal left and right spaces by combining the initial bases with Runge–Kutta stage directions, while previously computed stage bases are already contained in these spans. The construction supports low-rank implementation and can omit some stage directions only when the required stage-span property is preserved.

  • Stage basis augmentation: Previously computed stage bases U^(j) and V^(j), j ≤ ℓ, are already contained in the augmented spans and need not be appended explicitly.This containment is established by Lemma 1.
  • Stage basis augmentation: The RK basis constructor initializes the factors, evaluates s Runge–Kutta stages, updates temporary bases from stage directions, and returns the augmented bases U and V.It uses the inputs t_n, h, U_n, S_n, V_n and explicit RK coefficients a_ℓj and c_ℓ.
  • Low-rank implementation: All constructor operations can be performed in low-rank format without forming the full matrices Y^(ℓ) and Z̄_ℓ, provided the right-hand side admits the required projected products.The implementation uses products such as F^(ℓ)V^(ℓ) and F^(ℓ)ᵀU^(ℓ), together with small projected matrices.
  • Memory optimization: Zero Runge–Kutta coefficients can reduce the constructor, but omitted stages or directions must still preserve the stage-span property for every retained stage.A stage with b_ℓ = 0 may be omitted from the final basis only when it is not needed directly or indirectly for later retained stages.

7 Implementation details: parallel BUG-p integrator

The order-p parallel BUG integrator advances one step from t_n to t_{n+1}=t_n+h by constructing pre-augmented bases before performing its prescribed operations. An explicit s-stage Runge–Kutta method of order p supplies basis information and stops after stage s−1.

  • Time-step setup: Parallel BUG-p advances the solution from t_n to t_{n+1}=t_n+h in one time step.The step is described for Y_n=U_nS_nV_n^⊤ with V_n^⊤V_n=I.
  • Basis construction: The method first constructs pre-augmented bases bU_0 and bV_0 using an RK pre-basis constructor.The constructor is analogous to the one in Section 6.1.
  • Basis construction: The RK pre-basis constructor uses an explicit s-stage Runge–Kutta method of order p only for basis information and stops after stage s−1.The parallel BUG-p step then proceeds through its specified operations.

1. RK pre-basis construction. Compute pre-augmented orthonormal bases

The parallel BUG integrator constructs pre-augmented orthonormal bases, solves independent K-, L-, and S-step equations in parallel, then augments and truncates the resulting representation. Its projected evolutions use pre-augmented ranks, while repeated constructor truncations can bottleneck GPU computations.

  • RK pre-basis construction: Pre-augmented orthonormal bases are constructed from the current factors and Runge–Kutta-generated basis information.The basis constructor uses coefficients from an explicit s-stage Runge–Kutta method.
  • Parallel K-, L-, and S-steps: After basis construction, the K-, L-, and S-step matrix differential equations are independent and can be solved in parallel.This parallel structure is especially relevant for tensor-valued differential equations requiring multiple factors.
  • Parallel K-, L-, and S-steps: To retain order p, the subproblems must be integrated with a method whose local errors are O(hp+1).The computed endpoint values approximate K(tn+1), L(tn+1), and the Galerkin coefficient endpoint.
  • Final basis augmentation: The final bases append new information from the K- and L-steps while discarding numerically dependent columns.The expensive projected time evolutions are performed only at the pre-augmented ranks ˆr ≤ srn, even though the untruncated final representation may be larger.
  • Assembly and truncation: The augmented solution is assembled in the final bases and then compressed by truncation.For GPU computations, repeated truncations in the basis constructor can become infeasible and cause rank doubling at each step.

8 Numerical results

Numerical experiments evaluate the proposed BUG-p and parallel BUG-p integrators on non-stiff and stiff problems. BUG-p generally outperforms competing methods in the non-stiff setting and achieves order 6 with a rank 5r basis in the stiff test.

  • Non-stiff Schrödinger equation: For order p = 4, BUG-p has the best error-to-run-time inflection point among BUG-p, RK-BUG, and parallel BUG-p in both solver-order scenarios.The scenarios use order 8 Runge–Kutta evolution equations or Runge–Kutta schemes matching the BUG integrators’ order.
  • Non-stiff Schrödinger equation: On the non-stiff Schrödinger equation, BUG-p tends to outperform both RK-BUG and parallel BUG-p, although runtime comparisons are implementation- and hardware-dependent.The experiments compare fixed-order RK-8 and matching-order reduced solvers.
  • Stiff forced heat equation: The stiff benchmark uses a forced heat equation on an N × N grid with N = 128 and compares BUG-p against no explicit RK-BUG baseline because of parabolic stability restrictions.The reduced evolution equations use stiffly stable exponential time integrators.
  • Stiff forced heat equation: The stiff experiment constructs an order 6 BUG-p method with only a rank 5r basis and evolves the S matrix using a sixth-order exponential method.The reduced basis construction involves several substeps and singular value decompositions, minimizing the costs of the S-step.

9 Summary and conclusion

The paper introduces arbitrary-order BUG integrators whose coefficient updates are decoupled from the explicit RK method used to construct bases, enabling stiffly stable time integration. It derives augmented and parallel variants with robust order-p error bounds and validates their convergence, efficiency, and potential extensions numerically and theoretically.

  • Framework and novelty: The framework constructs BUG integrators of arbitrary order p without tying coefficient updates to the explicit Runge–Kutta scheme used for basis construction.The auxiliary RK method generates high-order basis directions, while reduced S-, K-, and L-step equations can use matching- or higher-order schemes.
  • Numerical validation: BUG-p=4 has the best error-to-runtime inflection point among BUG-p integrators of different order.Each marker in the efficiency experiment represents one time-step choice.
  • Theory: The augmented BUG-p and parallel BUG-p variants satisfy a curvature-independent error bound of order p under standard Lipschitz, boundedness, and small-curvature assumptions.The result is established in Theorem 1.
  • Theory: The required p-th-order basis information is contained in rank (s + 1)r for augmented BUG-p and sr for parallel BUG-p.Lemma 1 provides the rank characterization underlying the reduced basis construction.
  • Numerical validation: Numerical benchmarks confirmed predicted convergence rates and showed that coefficient-integrator decoupling can reduce error constants or increase observed convergence beyond nominal order.These effects were demonstrated on a non-stiff Schrödinger benchmark, alongside tests on a stiff forced heat equation.
  • Outlook: The general basis-construction conditions enable further BUG-p integrators, including multistep methods that may be especially relevant for parallel BUG-p.This identifies a direction for extending the proposed framework.

A Notation · B Recap: High-order BUG integrators

The paper establishes notation for low-rank solutions, BUG bases, projected Runge–Kutta stages, and truncation, then recaps midpoint and second-order parallel BUG integrators. Both recap methods augment bases with higher-order information before Galerkin or parallel updates and final truncation.

  • A Notation: The notation defines low-rank solutions Yn = UnSnVn^⊤, current rank r, truncation operators, orthonormal bases, and truncation tolerance ϑ.It also distinguishes augmented and parallel BUG bases, reduced coefficients, projected stage values, and full-rank quantities used in error analysis.
  • B Recap: High-order BUG integrators: Second-order BUG integrators extend augmented and parallel first-order methods by using augmentation steps to incorporate higher-order information.The recap presents these methods as predecessors to the proposed general-order constructions.
  • B.1 Midpoint BUG integrator: The midpoint BUG integrator first computes a half-step approximation with the first-order augmented BUG method without truncation, then performs a Galerkin update.The midpoint approximation has rank r ≤ br ≤ 2r before the subsequent augmented basis construction.
  • B.1 Midpoint BUG integrator: Truncating the midpoint approximation to rank r yields the rank 3r midpoint variant, whereas omitting truncation produces rank 4r basis matrices.The rank 4r formulation proceeds through basis augmentation, Galerkin integration of the reduced coefficient equation, and truncation.
  • B.2 Second–order parallel integrator: The second-order parallel integrator extends the first-order method by constructing pre-augmented bases from the old bases and right-hand-side products.The pre-augmented matrices have dimensions m × 2r and n × 2r.
  • B.2 Second–order parallel integrator: The parallel method independently advances K and L equations, appends their endpoint information to the pre-augmented bases, and orthonormalizes them into bases of rank up to 4r.It then forms the coefficient representation and truncates the factorized solution to a new rank r1.
Loading 2608.27749v1…