---
title: Projector-Splitting Integrator (PSI)
url: https://www.emergentmind.com/topics/projector-splitting-integrator-psi
type: topic
---

# Projector-Splitting Integrator (PSI)

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 [1301.1058]. 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 [1709.02594].

## 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\) matrices and \(P(Y)\) the orthogonal projector onto the tangent space \(T_Y\mathcal M_r\). In the matrix factorization
\[
Y(t)=U(t)S(t)V(t)^T,
\]
the tangent projector has the explicit form
\[
P(Y)Z = ZVV^T - UU^T ZVV^T + UU^T Z.
\]
Equivalently, using range projectors,
\[
P(Y)Z = ZP_{\mathcal R(Y^T)} - P_{\mathcal R(Y)}ZP_{\mathcal R(Y^T)} + P_{\mathcal R(Y)}Z.
\]
This projector identity is the algebraic foundation of PSI [1301.1058].

If one imposes the gauge conditions
\[
U(t)^T\dot U(t)=0,\qquad V(t)^T\dot V(t)=0,
\]
then the factor equations become
\[
\dot U = (I-UU^T)\,\dot A\,V\,S^{-1}, \qquad
\dot V = (I-VV^T)\,\dot A^T\,U\,S^{-T}, \qquad
\dot S = U^T\dot A\,V.
\]
For \(\dot A=F(A)\), the same formulas hold with \(\dot A\) replaced by \(F(Y)\). The appearance of \(S^{-1}\) and \(S^{-T}\) is the central numerical obstruction: standard integrators such as Runge–Kutta or implicit midpoint can become unstable or inaccurate when \(S\) is nearly singular, which occurs in over-approximation or when singular values approach a truncation threshold [1301.1058].

This geometric formulation generalizes beyond matrices. For Tucker tensors, the projected evolution is
\[
\dot Y(t)=P(Y(t))F(t,Y(t)),
\qquad
Y(t)=C(t)\times_1 U_1(t)\times_2\cdots\times_d U_d(t),
\]
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 [1802.01058].

## 2. Canonical \(K\)-\(S\)-\(L\) 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,
\[
P(Y)Z = ZVV^T - UU^T ZVV^T + UU^T Z,
\]
and then advancing successively under the corresponding subproblems. For a step from \(t_0\) to \(t_1=t_0+h\), starting from
\[
Y_0=U_0S_0V_0^T,
\]
the matrix PSI can be written in factorized form as follows [1301.1058].

In the prescribed-matrix case, with \(\Delta A=A(t_1)-A(t_0)\), the practical one-step algorithm is:
\[
K_1 = U_0 S_0 + \Delta A\,V_0,
\]
followed by a QR or SVD factorization
\[
K_1 = U_1 S_1.
\]
Then one updates the core by
\[
\widetilde S_0 = S_1 - U_1^T \Delta A\,V_0.
\]
Finally one updates the right factor via
\[
L_1 = V_0 \widetilde S_0^T + \Delta A^T U_1,
\]
followed by
\[
L_1 = V_1 S_1^T,
\]
and sets
\[
Y_1 = U_1 S_1 V_1^T.
\]
The same method may be expressed at the ODE level by the three substeps recalled in later work:
\[
\dot{K}(t)=F(t,K(t)V_0^\top)V_0,\qquad K(t_0)=U_0S_0,
\]
QR factorization \(K(t_1)=U_1\widehat S_1\), then
\[
\dot{S}(t)=-U_1^\top F(t,U_1S(t)V_0^\top)V_0,\qquad S(t_0)=\widehat S_1,
\]
and finally
\[
\dot{L}(t)=F(t,U_1L(t)^\top)^\top U_1,\qquad L(t_0)=V_0\widetilde S_0^\top,
\]
with \(L(t_1)=V_1S_1^\top\) and \(Y_1=U_1S_1V_1^\top\) [2010.02022].

Two structural features are essential. First, each subproblem is exactly solvable in the abstract splitting description and remains on the rank-\(r\) manifold. Second, the method does not invert \(S\), so it avoids the direct singular-factor mechanism present in the gauge-fixed factor ODEs [1301.1058]. 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 \(A_{1/2}=A(t_0+h/2)\) and the sequence
\[
K_{1/2} = U_0S_0 + (A_{1/2}-A_0)V_0,\qquad (U_{1/2},S_{1/2})=\mathrm{QR}(K_{1/2}),
\]
then a half-step core correction, a full \(L\)-update,
\[
L_1 = V_0 \widetilde S_0^T + (A_1-A_0)^T U_{1/2},\qquad (V_1,S_1)=\mathrm{QR}(L_1),
\]
followed by the reverse half-step sequence and final reconstruction \(Y_1=U_1S_1V_1^T\). Further compositions yield higher order [1301.1058].

## 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 \(r\) and the initial value satisfies
\[
Y_0 = A(t_0),
\]
then the projector-splitting integrator is exact in one step:
\[
Y_1 = A(t_1).
\]
The proof depends on the specific order \(K\to S\to L\); a different ordering such as \(K\to L\to S\) does not preserve this exactness property [1301.1058].

The second hallmark is robustness under over-approximation. For near-rank-deficient problems, if
\[
A(t)=A'(t)+\varepsilon A''(t),
\]
where \(A'(t)\) has rank \(q<r\), then the rank-\(r\) PSI solution stays close to the rank-\(q\) PSI solution:
\[
\|Y_n-\widehat Y_n\| \le C(\varepsilon+h).
\]
In later formulations of the robustness theory, the error estimate is stated as
\[
\|Y_n-A(t_n)\|\le c_0\delta+c_1\varepsilon+c_2 h,
\]
with constants depending on Lipschitz and boundedness data and the final time, but not on the singular values of the exact or numerical solution [2010.02022]. 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 \(S^{-1}\), 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 [1301.1058].

The same structure supports an optimization interpretation. In low-rank matrix optimization, when one has an updated matrix \(A+\Delta A\) and needs an efficient rank-\(r\) truncation or retraction, a single PSI step applied to \(A+t\Delta A\) at \(t=1\) yields a low-rank update without a full SVD of the ambient matrix [1301.1058]. A Tucker-tensor analogue of this retraction-like role was later demonstrated by applying one nested Tucker PSI step to
\[
\dot Y = P(Y)B,\qquad Y(t_0)=A,
\]
for approximate addition on a low-rank tensor manifold [1709.02594].

## 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)\in \mathbb{R}^{n_1\times \cdots \times n_d}, \qquad \dot A(t)=F(t,A(t)),
\]
one seeks a Tucker approximation
\[
Y(t)=C(t)\times_1 U_1(t)\times_2 \cdots \times_d U_d(t),
\]
of fixed multilinear rank \((r_1,\dots,r_d)\), satisfying
\[
\dot Y(t)=P(Y(t))\,F(t,Y(t)).
\]
The nested Tucker integrator proceeds mode by mode: it matricizes in mode \(i\), applies the matrix PSI to the unfolding, and replaces the large \(L\)-substep by a recursively reduced tensor ODE. This yields a mode-recursive \(K\)-\(S\)-\(L\) scheme whose final step is a core update \(C^1=\ten_d(R_d^T(t_1))\) and reconstruction
\[
Y^1 = C^1\times_1 U_1^1 \times_2 \cdots \times_d U_d^1.
\]
The method is exact on tensors of the prescribed multilinear rank, provided certain overlap matrices are invertible, and satisfies the error bound
\[
\|Y_n-A(t_n)\|\le c_1 h + c_2 \varepsilon,
\]
with constants independent of small singular values of the unfoldings [1709.02594].

In MCTDH, projector splitting was revisited from a tangent-space projection standpoint by decomposing the full variational projector as
\[
\mathcal{P}(\Psi)=\mathcal{P}_0(\Psi)+\sum_\kappa \mathcal{P}_\kappa(\Psi),
\]
and then each mode projector as
\[
\mathcal{P}_\kappa(\Psi)=\mathcal{P}_\kappa^+(\Psi)-\mathcal{P}_\kappa^-(\Psi).
\]
The central device is a QR decomposition of the mode-\(\kappa\) coefficient matricization,
\[
A^{(\kappa)} = Q^{(\kappa)} S^{(\kappa)T},
\]
which induces transformed single-particle functions
\[
\tilde{\boldsymbol{\varphi}}^{(\kappa)} = S^{(\kappa)T}\boldsymbol{\varphi}^{(\kappa)}
\]
and orthonormal transformed single-hole functions. In this representation, the reduced density matrix factorizes as
\[
\rho^{(\kappa)} = S^{(\kappa)} S^{(\kappa)\dagger},
\]
so the problematic inverse \((\rho^{(\kappa)})^{-1}\) 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 \(S^{(\kappa)}\) [1802.01058].

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,
\[
\left[(\rho^{z_l})^{-1}\right]_{ki_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 \(\widetilde A^{z_l}\) is propagated forward, orthogonally decomposed,
\[
\widetilde{A}_{I^{z_l}i_{0}}^{z_l}=A_{I^{z_l}j_{0}}^{z_l}R_{j_{0}i_{0}}^{z_l},
\]
and followed by backward propagation of an auxiliary \(R^{z_l}\). This implementation was reported to be stable for large ML-MCTDH wavefunctions containing up to hundreds of thousands of nodes, with roughly \(3\)–\(4\) orders of magnitude fewer Hamiltonian evaluations and \(2\)–\(3\) orders of magnitude fewer Hamiltonian applications than standard ML-MCTDH on the spin-boson benchmarks considered there [2109.03134].

## 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
\[
f(t,x,v)\approx \sum_{i,j=1}^r X_i(t,x)\,S_{ij}(t)\,V_j(t,v),
\]
with orthonormal spatial and velocity bases. The tangent-space projector splits as
\[
P(f)g = P_{\overline X}g + P_{\overline V}g - P_{\overline X}P_{\overline V}g,
\]
which yields a three-substep low-rank Lie splitting. With
\[
K_j(t,x)=\sum_i X_i(t,x)S_{ij}(t), \qquad L_i(t,v)=\sum_j S_{ij}(t)V_j(t,v),
\]
the \(K\)- and \(L\)-substeps become systems of \(d\)-dimensional advection equations, while the \(S\)-substep becomes a finite-dimensional matrix ODE. In the hierarchical extension, the \(X_i\) and \(V_j\) 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 \(\mathcal O(n^{2d})\) to \(\mathcal O(n^d r)\) for the basic low-rank method, and that no CFL restriction is imposed when semi-Lagrangian or FFT-based subsolvers are used [1801.01103].

For the Vlasov–Maxwell system, the low-rank ansatz
\[
f(t,x,v)=\sum_{i,j=1}^r X_i(t,x)\,S_{ij}(t)\,V_j(t,v)
\]
again leads to projected evolution on a low-rank manifold. The projector is written as
\[
P(f)g = P_{\bar V}g - P_{\bar V}P_{\bar X}g + P_{\bar X}g,
\]
and the three substeps produce reduced PDEs in \(x\) or \(v\), plus an \(r^2\)-dimensional ODE for \(S\). 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 \(\phi\) is introduced through
\[
\partial_t E = \nabla\times B - j - \nabla \phi,
\]
with a Poisson equation for \(\phi\) 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 [1902.00424].

Random parabolic equations provide a different PDE realization. In the Dual DO formulation, a rank-\(R\) approximation is written as
\[
u(t)=\bar u(t)+U(t)Y(t)^\top,
\]
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
\[
\Pi_{UY^\top}[v] = P_{\mathcal Y}[v] + P_{\mathcal Y^\perp}\big[P_{\mathcal U}[v]\big],
\]
and the authors show that the fully discrete scheme is, in the full-rank case, equivalent in exact arithmetic to the first-order \(K\!S\!L\) 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 [2006.05211].

## 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 \(L^2\)-stability analysis showed that the backward-in-time \(S\)-step can reverse the damping introduced by numerical fluxes. In the Lax–Friedrichs setting, this leads to amplification of certain Fourier modes and \(L^2\)-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
\[
\lambda_{\max}(\tilde{\mathbf A})\frac{\Delta t}{\Delta x}\le 1
\quad\Longrightarrow\quad
\|\mathbf u^{n+1}\|_F \le \|\mathbf u^n\|_F.
\]
It also compares PSI with the unconventional integrator of Ceruti–Lubich, which orders the substeps \(K\)-\(L\)-\(S\) 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 [2107.07282].

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
\[
\nu=\lambda_{\max}\frac{\Delta t}{\Delta x}\le \frac13,
\]
even though the amplification factors differ. For parabolic problems, the negative \(S\)-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 \(K\)- and \(L\)-steps and forward Euler in the \(S\)-step has amplification factor
\[
G_{\text{hybrid}}(x)=\frac{1}{1+x},
\]
which matches backward Euler. Strang splitting substantially enlarges the hyperbolic stability region and remains unconditionally stable with Crank–Nicolson in the parabolic case [2507.15192].

The backward \(S\)-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,
\[
\dot{S}(t)=U_1^\top F(t,U_1S(t)V_1^\top)V_1,
\]
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 [2010.02022].

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 \(X_0=U_0S_0V_0^\top\), it computes enlarged bases \(\widehat U_1\) and \(\widehat V_1\), forms an augmented coefficient matrix
\[
\widehat S_1 =
\begin{pmatrix}
\bar S(t_1) & \widetilde S_1^L \\
\widetilde S_1^K & \mathbf 0
\end{pmatrix},
\]
and truncates by SVD. Its global error estimate has the PSI-like form
\[
\|X_n-X(t_n)\| \le C_0\delta + C_1 h + C_2 \varepsilon + C_3 \vartheta/h,
\]
with constants independent of singular values, but the exactness and energy or dissipation preservation properties of BUG are not inherited [2304.05660].

In low-rank neural network training, PSI has been applied to the projected gradient flow
\[
\dot Y(t) = -P(Y(t))\nabla \ell(Y(t)).
\]
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
\[
\widehat U_{k+1} = \operatorname{QR}\big([K_k \mid K_{k+1}]\big),
\]
uses it in the right-factor update, and truncates by SVD. The reported theory includes the robust bound
\[
\|Y(t_n)-W(t_n)\| \le \epsilon + c_1 h + c_2\delta + \frac{\vartheta}{h},
\]
again independent of singular values, and a discrete descent estimate
\[
\ell(\widehat Y_1)\le \ell(Y_0) -\left(1-\frac{hc_l}{2}\right)h\,\|P_{\widehat U_1}\nabla\ell(Y_0)\|^2,
\]
which is nonincreasing for \(h\le 2/c_l\). 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 [2502.03006].

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 \(S\)-step remains the focal point of both its mathematical distinctiveness and its principal limitations, especially for dissipative, stiff, or discretized hyperbolic problems [1301.1058].

Source: https://www.emergentmind.com/topics/projector-splitting-integrator-psi