Projector-Splitting Integrator (PSI)

• PSI is a time-integration method for dynamical low-rank approximation that avoids inverting nearly singular core factors.
• It employs a Lie–Trotter splitting of the tangent-space projector into three subproblems, ensuring robust, exact updates on fixed-rank manifolds.
• PSI extends to Tucker tensors, kinetic equations, and neural network training, offering stability and efficiency in high-dimensional applications.

The projector-splitting integrator (PSI) is a time-integration method for dynamical low-rank approximation in which the evolution is constrained to a fixed-rank manifold and the tangent-space projector is split into simpler subprojections that are advanced sequentially. In its original matrix form, PSI was introduced for approximating a time-dependent matrix $A(t)$ by a rank-$r$ factorization $Y(t)=U(t)S(t)V(t)^T$, with $U$ and $V$ orthonormal and $S$ invertible, and it was designed to avoid the numerical pathologies caused by nearly singular core factors in direct factor-ODE formulations (Lubich et al., 2013). Subsequent work established PSI as a reference robust integrator for matrix and tensor dynamical low-rank approximation, extended it to Tucker tensors and multiconfiguration wavefunction methods, and used it as a building block in kinetic plasma simulation, random parabolic PDE solvers, and low-rank training algorithms (Lubich et al., 2017).

1. Geometric formulation on low-rank manifolds

In dynamical low-rank approximation, one replaces a full evolution equation by a projected evolution on a rank-constrained manifold. For matrices, the basic setting is

$\dot A(t)=F(t,A(t)), \qquad A(t_0)=A_0,$

or, when $A(t)$ is explicitly given,

$$\dot Y(t)=\arg\min_{Z\in T_{Y(t)}\mathcal M_r}\|\dot A(t)-Z\|_F, \qquad \dot Y(t)=P(Y(t))\,\dot A(t),$$

with $\mathcal M_r$ the manifold of rank-$r$0 matrices and $r$1 the orthogonal projector onto the tangent space $r$2. In the matrix factorization

$r$3

the tangent projector has the explicit form

$r$4

Equivalently, using range projectors,

$r$5

This projector identity is the algebraic foundation of PSI (Lubich et al., 2013).

If one imposes the gauge conditions

$r$6

then the factor equations become

$r$7

For $r$8, the same formulas hold with $r$9 replaced by $Y(t)=U(t)S(t)V(t)^T$0. The appearance of $Y(t)=U(t)S(t)V(t)^T$1 and $Y(t)=U(t)S(t)V(t)^T$2 is the central numerical obstruction: standard integrators such as Runge–Kutta or implicit midpoint can become unstable or inaccurate when $Y(t)=U(t)S(t)V(t)^T$3 is nearly singular, which occurs in over-approximation or when singular values approach a truncation threshold (Lubich et al., 2013).

This geometric formulation generalizes beyond matrices. For Tucker tensors, the projected evolution is

$Y(t)=U(t)S(t)V(t)^T$4

with fixed multilinear rank. In MCTDH and ML-MCTDH, the same tangent-space projection viewpoint is used to reinterpret variational wavefunction propagation and to expose the role of projector splitting in avoiding reduced-density-matrix inversions (Bonfanti et al., 2018).

2. Canonical $Y(t)=U(t)S(t)V(t)^T$5-$Y(t)=U(t)S(t)V(t)^T$6-$Y(t)=U(t)S(t)V(t)^T$7 construction

The defining idea of PSI is to apply a Lie–Trotter splitting to the projected low-rank evolution by decomposing the tangent projector into three terms,

$Y(t)=U(t)S(t)V(t)^T$8

and then advancing successively under the corresponding subproblems. For a step from $Y(t)=U(t)S(t)V(t)^T$9 to $U$0, starting from

$U$1

the matrix PSI can be written in factorized form as follows (Lubich et al., 2013).

In the prescribed-matrix case, with $U$2, the practical one-step algorithm is: $U$3 followed by a QR or SVD factorization

$U$4

Then one updates the core by

$U$5

Finally one updates the right factor via

$U$6

followed by

$U$7

and sets

$U$8

The same method may be expressed at the ODE level by the three substeps recalled in later work: $U$9 QR factorization $V$0, then

$V$1

and finally

$V$2

with $V$3 and $V$4 (Ceruti et al., 2020).

Two structural features are essential. First, each subproblem is exactly solvable in the abstract splitting description and remains on the rank-$V$5 manifold. Second, the method does not invert $V$6, so it avoids the direct singular-factor mechanism present in the gauge-fixed factor ODEs (Lubich et al., 2013). This is the basis for the method’s robustness under over-approximation.

The basic PSI is first order. A symmetric second-order scheme is obtained by composition. One formulation uses midpoint data $V$7 and the sequence

$V$8

then a half-step core correction, a full $V$9-update,

$S$0

followed by the reverse half-step sequence and final reconstruction $S$1. Further compositions yield higher order (Lubich et al., 2013).

3. Exactness, robustness, and rank-adaptive interpretation

The canonical theoretical results for PSI are exactness on rank-preserving evolutions and robustness with respect to small singular values. If the exact matrix solution remains of rank at most $S$2 and the initial value satisfies

$S$3

then the projector-splitting integrator is exact in one step: $S$4 The proof depends on the specific order $S$5; a different ordering such as $S$6 does not preserve this exactness property (Lubich et al., 2013).

The second hallmark is robustness under over-approximation. For near-rank-deficient problems, if

$S$7

where $S$8 has rank $S$9, then the rank-$\dot A(t)=F(t,A(t)), \qquad A(t_0)=A_0,$0 PSI solution stays close to the rank-$\dot A(t)=F(t,A(t)), \qquad A(t_0)=A_0,$1 PSI solution: $\dot A(t)=F(t,A(t)), \qquad A(t_0)=A_0,$2 In later formulations of the robustness theory, the error estimate is stated as

$\dot A(t)=F(t,A(t)), \qquad A(t_0)=A_0,$3

with constants depending on Lipschitz and boundedness data and the final time, but not on the singular values of the exact or numerical solution (Ceruti et al., 2020). This independence from small singular values is the central robustness property associated with PSI.

The practical consequences are direct. PSI is explicit and inexpensive in the original formulation; it requires matrix–matrix products with increments or right-hand sides and QR or SVD factorizations of skinny matrices rather than large full decompositions. Because it avoids $\dot A(t)=F(t,A(t)), \qquad A(t_0)=A_0,$4, it remains usable when the chosen rank exceeds the effective rank. This also makes adaptive rank changes natural: lowering the rank is trivial, and increasing the rank can be done without the singular-factor breakdown that obstructs direct factor-equation integrators (Lubich et al., 2013).

The same structure supports an optimization interpretation. In low-rank matrix optimization, when one has an updated matrix $\dot A(t)=F(t,A(t)), \qquad A(t_0)=A_0,$5 and needs an efficient rank-$\dot A(t)=F(t,A(t)), \qquad A(t_0)=A_0,$6 truncation or retraction, a single PSI step applied to $\dot A(t)=F(t,A(t)), \qquad A(t_0)=A_0,$7 at $\dot A(t)=F(t,A(t)), \qquad A(t_0)=A_0,$8 yields a low-rank update without a full SVD of the ambient matrix (Lubich et al., 2013). A Tucker-tensor analogue of this retraction-like role was later demonstrated by applying one nested Tucker PSI step to

$\dot A(t)=F(t,A(t)), \qquad A(t_0)=A_0,$9

for approximate addition on a low-rank tensor manifold (Lubich et al., 2017).

4. Generalizations to Tucker tensors and multiconfiguration wavefunctions

The tensor generalization of PSI in Tucker format is obtained by applying the matrix projector-splitting idea recursively to tensor unfoldings. For a tensor

$A(t)$0

one seeks a Tucker approximation

$A(t)$1

of fixed multilinear rank $A(t)$2, satisfying

$A(t)$3

The nested Tucker integrator proceeds mode by mode: it matricizes in mode $A(t)$4, applies the matrix PSI to the unfolding, and replaces the large $A(t)$5-substep by a recursively reduced tensor ODE. This yields a mode-recursive $A(t)$6-$A(t)$7-$A(t)$8 scheme whose final step is a core update $A(t)$9 and reconstruction

$$\dot Y(t)=\arg\min_{Z\in T_{Y(t)}\mathcal M_r}\|\dot A(t)-Z\|_F, \qquad \dot Y(t)=P(Y(t))\,\dot A(t),$$0

The method is exact on tensors of the prescribed multilinear rank, provided certain overlap matrices are invertible, and satisfies the error bound

$$\dot Y(t)=\arg\min_{Z\in T_{Y(t)}\mathcal M_r}\|\dot A(t)-Z\|_F, \qquad \dot Y(t)=P(Y(t))\,\dot A(t),$$1

with constants independent of small singular values of the unfoldings (Lubich et al., 2017).

In MCTDH, projector splitting was revisited from a tangent-space projection standpoint by decomposing the full variational projector as

$$\dot Y(t)=\arg\min_{Z\in T_{Y(t)}\mathcal M_r}\|\dot A(t)-Z\|_F, \qquad \dot Y(t)=P(Y(t))\,\dot A(t),$$2

and then each mode projector as

$$\dot Y(t)=\arg\min_{Z\in T_{Y(t)}\mathcal M_r}\|\dot A(t)-Z\|_F, \qquad \dot Y(t)=P(Y(t))\,\dot A(t),$$3

The central device is a QR decomposition of the mode-$$\dot Y(t)=\arg\min_{Z\in T_{Y(t)}\mathcal M_r}\|\dot A(t)-Z\|_F, \qquad \dot Y(t)=P(Y(t))\,\dot A(t),$$4 coefficient matricization,

$$\dot Y(t)=\arg\min_{Z\in T_{Y(t)}\mathcal M_r}\|\dot A(t)-Z\|_F, \qquad \dot Y(t)=P(Y(t))\,\dot A(t),$$5

which induces transformed single-particle functions

$$\dot Y(t)=\arg\min_{Z\in T_{Y(t)}\mathcal M_r}\|\dot A(t)-Z\|_F, \qquad \dot Y(t)=P(Y(t))\,\dot A(t),$$6

and orthonormal transformed single-hole functions. In this representation, the reduced density matrix factorizes as

$$\dot Y(t)=\arg\min_{Z\in T_{Y(t)}\mathcal M_r}\|\dot A(t)-Z\|_F, \qquad \dot Y(t)=P(Y(t))\,\dot A(t),$$7

so the problematic inverse $$\dot Y(t)=\arg\min_{Z\in T_{Y(t)}\mathcal M_r}\|\dot A(t)-Z\|_F, \qquad \dot Y(t)=P(Y(t))\,\dot A(t),$$8 disappears from the propagation of the transformed SPFs. The resulting projector-splitting algorithm is a Strang-splitting method that alternates QR-based orthogonalization, forward propagation of the transformed SPFs, and a backward gauge-correction step for $$\dot Y(t)=\arg\min_{Z\in T_{Y(t)}\mathcal M_r}\|\dot A(t)-Z\|_F, \qquad \dot Y(t)=P(Y(t))\,\dot A(t),$$9 (Bonfanti et al., 2018).

In ML-MCTDH, PSI was further developed as a singularity-free time integrator for tree tensor network wavefunctions. Standard ML-MCTDH equations contain inverses of mean-field density matrices,

$\mathcal M_r$0

which become numerically unstable when the wavefunction is weakly entangled. PSI avoids this by switching to a representation in which the inverse is the identity for the selected node, then performing an Euler-tour traversal of the tree. At each node, a transformed coefficient tensor $\mathcal M_r$1 is propagated forward, orthogonally decomposed,

$\mathcal M_r$2

and followed by backward propagation of an auxiliary $\mathcal M_r$3. This implementation was reported to be stable for large ML-MCTDH wavefunctions containing up to hundreds of thousands of nodes, with roughly $\mathcal M_r$4–$\mathcal M_r$5 orders of magnitude fewer Hamiltonian evaluations and $\mathcal M_r$6–$\mathcal M_r$7 orders of magnitude fewer Hamiltonian applications than standard ML-MCTDH on the spin-boson benchmarks considered there (Lindoy et al., 2021).

5. Kinetic and PDE realizations

PSI has been adapted to high-dimensional kinetic equations by representing the phase-space density in separable low-rank form and splitting the projected dynamics into lower-dimensional subproblems. For the Vlasov–Poisson equation, one uses

$\mathcal M_r$8

with orthonormal spatial and velocity bases. The tangent-space projector splits as

$\mathcal M_r$9

which yields a three-substep low-rank Lie splitting. With

$r$00

the $r$01- and $r$02-substeps become systems of $r$03-dimensional advection equations, while the $r$04-substep becomes a finite-dimensional matrix ODE. In the hierarchical extension, the $r$05 and $r$06 are themselves given low-rank tensor structure, reducing each time step further to sets of one-dimensional advection problems. The paper emphasizes that storage drops from $r$07 to $r$08 for the basic low-rank method, and that no CFL restriction is imposed when semi-Lagrangian or FFT-based subsolvers are used (Einkemmer et al., 2018).

For the Vlasov–Maxwell system, the low-rank ansatz

$r$09

again leads to projected evolution on a low-rank manifold. The projector is written as

$r$10

and the three substeps produce reduced PDEs in $r$11 or $r$12, plus an $r$13-dimensional ODE for $r$14. A complication absent from Vlasov–Poisson is that the low-rank and time-splitting approximations can break the consistency between the computed charge density and Gauss’ law for the electric field. To restore compatibility, a correction potential $r$15 is introduced through

$r$16

with a Poisson equation for $r$17 chosen so that the updated electric field satisfies the discrete Gauss constraint. The correction enforces Gauss’ law up to machine precision in the reported experiments, although the paper also notes that in strongly nonlinear regimes it can slightly worsen qualitative field behavior when the rank is too small (Einkemmer et al., 2019).

Random parabolic equations provide a different PDE realization. In the Dual DO formulation, a rank-$r$18 approximation is written as

$r$19

and the discrete low-rank scheme updates the mean, the deterministic modes, and the stochastic modes in a Gauss–Seidel-like staggered manner. The zero-mean part is projected with the tangent-space projector

$r$20

and the authors show that the fully discrete scheme is, in the full-rank case, equivalent in exact arithmetic to the first-order $r$21 projector-splitting method of Lubich-type. The discrete variational structure yields stability bounds independent of the smallest singular value: the implicit scheme is unconditionally stable, while explicit and semi-implicit variants satisfy parabolic CFL-type conditions (Kazashi et al., 2020).

6. Stability, structural limitations, and post-2020 variants

Although PSI is robust with respect to small singular values, its stability for PDE discretizations depends on how the low-rank projection is combined with spatial and temporal discretization. For linear hyperbolic systems discretized first and then projected, a detailed $r$22-stability analysis showed that the backward-in-time $r$23-step can reverse the damping introduced by numerical fluxes. In the Lax–Friedrichs setting, this leads to amplification of certain Fourier modes and $r$24-instability that is independent of the time-step size. The same work advocates applying DLRA to the continuous PDE first and only then discretizing the resulting low-rank subproblems. The resulting stabilized PSI recovers the classic CFL condition

$r$25

It also compares PSI with the unconventional integrator of Ceruti–Lubich, which orders the substeps $r$26-$r$27-$r$28 and avoids the destabilizing backward substep; that method is reported to be generally superior for first-order moments, while stabilized PSI performs better for higher moments in the reported hyperbolic tests (Kusch et al., 2021).

A later von Neumann-type analysis of linear hyperbolic and parabolic model equations studied PSI under both discretize-then-project and project-then-discretize formulations. For hyperbolic transport with Lie–Trotter splitting and forward Euler in the substeps, both DtP and PtD were shown to satisfy the same stability condition

$r$29

even though the amplification factors differ. For parabolic problems, the negative $r$30-step does not preclude unconditional stability if the substeps are discretized appropriately: Crank–Nicolson yields unconditional stability, and a hybrid scheme using backward Euler in the $r$31- and $r$32-steps and forward Euler in the $r$33-step has amplification factor

$r$34

which matches backward Euler. Strang splitting substantially enlarges the hyperbolic stability region and remains unconditionally stable with Crank–Nicolson in the parabolic case (Zhang et al., 21 Jul 2025).

The backward $r$35-step is also the point of departure for alternative robust integrators. The unconventional robust integrator of Ceruti and Lubich retains the exactness and singular-value-independent error bounds associated with PSI, but updates the left and right basis matrices in parallel and advances the core with a forward-in-time equation,

$r$36

thus eliminating the backward middle substep. The paper emphasizes potential advantages for strongly dissipative problems, more parallelism, and preservation of symmetry or skew-symmetry in settings where PSI does not preserve these structures (Ceruti et al., 2020).

Rank adaptation and parallelism were later pursued more aggressively. A parallel rank-adaptive DLRA integrator was constructed to inherit the robust error bound associated with PSI and BUG while solving the left basis update, right basis update, and reduced coefficient evolution in parallel. Starting from $r$37, it computes enlarged bases $r$38 and $r$39, forms an augmented coefficient matrix

$r$40

and truncates by SVD. Its global error estimate has the PSI-like form

$r$41

with constants independent of singular values, but the exactness and energy or dissipation preservation properties of BUG are not inherited (Ceruti et al., 2023).

In low-rank neural network training, PSI has been applied to the projected gradient flow

$r$42

There, the middle substep can increase the loss because of its backward-gradient character. This motivated backward-corrected and augmented variants. The augmented backward-corrected projector-splitting integrator forms an enlarged basis

$r$43

uses it in the right-factor update, and truncates by SVD. The reported theory includes the robust bound

$r$44

again independent of singular values, and a discrete descent estimate

$r$45

which is nonincreasing for $r$46. In the reported benchmarks, the method reduces the per-iteration decomposition cost to one QR decomposition and one SVD, compared with two QR decompositions and one SVD for standard BUG-style rank-adaptive schemes (Kusch et al., 5 Feb 2025).

Across these developments, PSI remains the canonical robust integrator for dynamical low-rank approximation: it established the small-singular-value-independent error mechanism, provided exactness on rank-preserving flows, and furnished a unifying template from which tensor, wavefunction, kinetic, and machine-learning variants have been derived. At the same time, the backward $r$47-step remains the focal point of both its mathematical distinctiveness and its principal limitations, especially for dissipative, stiff, or discretized hyperbolic problems (Lubich et al., 2013).