# Iterative thresholding low-rank time integration

Matthieu Dolbeault, Markus Bachmayr, Polina Sachsenmaier

► To cite this version:

Matthieu Dolbeault, Markus Bachmayr, Polina Sachsenmaier. Iterative thresholding low-rank time integration. 2025. hal-05180198

HAL Id: hal-05180198

https://hal.science/hal-05180198v1

Preprint submitted on 23 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ITERATIVE THRESHOLDING LOW-RANK TIME INTEGRATION

MARKUS BACHMAYR, MATTHIEU DOLBEAULT, POLINA SACHSENMAIER

**ABSTRACT.** We develop time integration methods in low-rank representation that can adaptively adjust approximation ranks to achieve a prescribed accuracy, while ensuring that these ranks remain proportional to the corresponding best approximation ranks. Our approach relies on an iterative scheme combined with soft thresholding of the iterates. A model case of a time-dependent Schrödinger equation with low-rank matrix approximation is analyzed in detail, and the required modifications for second-order parabolic problems are described. Numerical tests illustrate the results for both cases.

**Keywords.** low-rank approximation, evolution equations, quasi-optimal ranks

**Mathematics Subject Classification.** Primary 65F55, 65M12; Secondary 65Y20, 65L70

## 1. INTRODUCTION

Standard numerical methods for solving partial differential equations (PDEs), when applied to problems on domains of large dimension, typically suffer from the *curse of dimensionality*: namely, their computational cost, as well as the space required to store the solution, scale exponentially in the dimension, making the methods impractical even at low resolution. In many cases of interest, however, such limitations can be overcome by appropriate compressed representations of approximate solutions, in particular by low-rank tensor representations.

Here we consider both dispersive and parabolic PDEs. The former arise as high-dimensional time-dependent problems in quantum mechanics, with the linear Schrödinger equation

$$ (1.1) \qquad i\partial_t u = -\Delta u + Vu $$

posed over $\mathbb{R}^{3N}$ when describing a system of $N$ particles in $\mathbb{R}^3$ subject to the potential function $V$. When working with low-rank approximations, such problems are commonly formulated in the framework of second quantization in terms of occupation numbers of orbitals. Another motivating class of examples, this time of parabolic type, are Fokker-Planck equations of the form

$$ (1.2) \qquad \partial_t u = \operatorname{div}(D\nabla u + uP), $$

which describe probability densities $u$ for trajectories of stochastic processes, as in the overdamped Langevin equation in molecular dynamics. Here, the dimension of the underlying space corresponds to the total number of degrees of freedom of the considered molecular system, the matrix $D$ is a diffusion coefficient, and the vector field $P$ corresponds to the acting forces.

Approximations based on low-rank tensor representations have a long history for both of these classes of problems. In many cases of interest, such representations can avoid the curse of dimensionality by a reduction to sums of tensor products of lower-dimensional factors. However, their efficiency depends crucially on the sizes of certain rank parameters in the

<footnote>INSTITUT FÜR GEOMETRIE UND PRAKTISCHE MATHEMATIK, RWTH AACHEN UNIVERSITY, TEMPLER-GRABEN 55, 52062 AACHEN, GERMANY</footnote>

<footnote>E-mail address: {bachmayr,dolbeault,sachsenmaier}@igpm.rwth-aachen.de.</footnote>

<footnote>M.B. acknowledges funding by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project number 442047500/SFB 1481 *Sparsity and Singular Structures* and by the European Union (ERC, COCOA, 101170147). M.D. was supported by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project number 442047500/SFB 1481 *Sparsity and Singular Structures*. P.S. was supported by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project number 320021702/GRK2326 *Energy, Entropy, and Dissipative Dynamics*.</footnote>
tensor formats, which reduce to matrix ranks in the case of tensor order two. In this regard, a central issue in working with such low-rank representations is that they do not form linear spaces: under vector space operations and other mappings, their rank parameters generally increase. For reasonable computational costs, methods operating on low-rank tensor representations thus need to ensure that rank parameters remain under control.

For time-dependent PDEs, several approaches for low-rank approximation exist that use different mechanisms for controlling ranks. These range from methods that keep the ranks fixed, such as dynamical low-rank approximation, which may lead to uncontrolled errors, to methods that approximate standard time stepping schemes up to any desired accuracy, but may instead lead to unnecessarily large ranks. While for methods based on space-time variational formulations, first results on how to achieve controlled ranks for given guaranteed errors exist, the availability of suitable such formulations can be a demanding condition.

Our main result in this work is a method for time integration that does not require a space-time variational formulation, but still yields provably convergent approximations while at the same time ensuring control of the ranks of the produced approximations and all intermediate quantities. More precisely, we show that these ranks remain comparable to their natural benchmark quantities, that is, to the best approximation ranks for the sought solutions at the achieved accuracy.

We focus here on representations in terms of matrices in low-rank representation, which involve fewer technicalities than higher-order tensors. However, for our purposes, they exhibit all main difficulties of the higher-order case, since the target tensor formats for higher-dimensional problems, such as tensor trains and hierarchical tensors, are themselves based on matrix ranks of certain matricizations (see, for example, the monograph [23] and the surveys [6,22]). In particular, the standard methods for error-controlled rank reduction that are available for these higher-order tensor formats reduce to truncation or thresholding of the singular value decomposition in the matrix case. The schemes proposed here can thus be applied to general hierarchical tensor representations without essential modifications.

1.1. **Rationale and main result.** We consider here time-dependent problems in two spatial variables on a product domain $\Omega = \Omega_1 \times \Omega_2$ with $\Omega_i \subseteq \mathbb{R}^{d_i}$ for $i = 1, 2$, using low-rank approximations to obtain a reduction to operations on lower-dimensional factors; this is also the main aim in more general higher-dimensional problems.

The low-rank approximations that we aim to construct for solutions $u \in C([0, T]; L^2(\Omega))$ are of the basic form

$$ (1.3) \qquad u(t, x_1, x_2) \approx \sum_{k=1}^{r(t)} u_k^{(1)}(t, x_1) u_k^{(2)}(t, x_2), \quad \text{for } t \in [0, T], \quad x_1 \in \Omega_1, \ x_2 \in \Omega_2, $$

where $u_k^{(i)} \in L^2(\Omega_i)$ for $i = 1, 2$ and $k = 1, \dots, r(t)$. Here a main difficulty is to adapt the rank parameter $r(t)$ to a prescribed error.

For each fixed $t$, low-rank best approximations of $(x_1, x_2) \mapsto u(t, x_1, x_2)$ as an element of $L^2(\Omega) = L^2(\Omega_1) \otimes L^2(\Omega_2)$ for each given rank can be obtained by singular value decomposition, where the error of low-rank approximation is given by the $\ell^2$ norm of the sequence of omitted singular values. For each given $\varepsilon > 0$, these tail norms determine a minimum rank $\hat{r}_\varepsilon(t)$ such that the best approximation of this rank has $L^2$-error at most $\varepsilon$.

Our main objective is to obtain, for any prescribed accuracy $\varepsilon$, low-rank approximations as in (1.3) with quasi-optimal ranks $r(t)$ – that is, we aim at a method that produces approximations with $r(t)$ comparable, up to constants, to the ranks $\hat{r}_\varepsilon(t)$ of best low-rank approximations of the same accuracy. At the same time, the ranks arising in the computational method should also remain controlled in terms of the ranks of the final result.

This requires a balance between error accumulation and applying sufficient rank truncations. Such a balance is difficult to achieve in step-truncation methods based on classical time stepping schemes: to ensure sufficiently small cumulative errors, also the rank truncation tolerances for each time step need to tighten accordingly, and a comparison to best
approximation ranks is difficult to achieve – indeed, no result of this type appears to be known for such methods.

In view of the extension of our approach to higher dimensions, what we mainly aim for is a reduction to operations on lower-dimensional factors in low-rank representations. At the same time, the methods that we consider here do not rely on direct manipulations of lower-dimensional components of tensor representations and are thus directly generalizable to various low-rank formats for higher-order tensors. Especially in the high-dimensional case, the overall computational costs are then driven mainly by the ranks of tensor representations.

To obtain ranks that are appropriate for the achieved approximation error, we use iterative refinement of approximate solutions in basis representations in time, applied on subintervals of length $h$ of the given total time interval $[0, T]$. Here, unlike standard time stepping methods, we do not necessarily let $h$ tend to zero to enforce convergence, but rather consider fixed (or slowly decreasing) $h$ in combination with a sufficiently refined temporal basis in each subinterval. Unlike space-time variational formulations, we rely only on standard integral formulations of evolution problems with uniform approximation in time.

To ensure that $h$ can be chosen to be relatively large and remains in particular independent of the spatial discretization, we use techniques similar to exponential integrators [28], but with particular adaptations to the low-rank setting. The basic idea was presented by Lawson in [34] for Runge-Kutta methods and later used in exponential integrators [30], see also [27].

For the Laplacian $\Delta$ on $\Omega$, we have $\Delta = \Delta_1 \otimes I + I \otimes \Delta_2$, where $\Delta_i$ is the Laplacian on $\Omega_i$ for $i = 1, 2$. For $u_i \in L_2(\Omega_i)$ and any $z \in \mathbb{C}$ with $\mathrm{Re} z \ge 0$,

$$ (1.4) \qquad e^{z\Delta}(u_1 \otimes u_2) = e^{z\Delta_1}u_1 \otimes e^{z\Delta_2}u_2. $$

As a consequence, the action of $e^{z\Delta}$ leaves ranks unchanged and thus can directly be reduced to lower-dimensional operations. Our starting point to make use of this fact is Duhamel's formula, which in the case of the time-dependent Schrödinger equation (1.1) with initial data $u_0$ reads

$$ u(t) = e^{it\Delta}u_0 - i \int_0^t e^{i(t-s)\Delta}Vu(s)ds. $$

For this class of problems, we additionally use the (rank-preserving) transition to the twisted variable $v(t) = e^{-it\Delta}u(t)$ to arrive at the formulation

$$ (1.5) \qquad v(t) = u_0 - i \int_0^t e^{-is\Delta}Ve^{is\Delta}v(s)ds. $$

We now combine this formulation with a choice of temporal basis functions, where we use polynomials of degree $J$ on each subinterval; although other choices could be of interest, we focus on this particular case in this work.

Under appropriate assumptions on $V$, (1.5) can serve as the basis for fixed-point iterations for refining approximations of $v$ on each subinterval. These converge for subinterval sizes $h$ that depend only on $V$, but not on the spatial discretizations. The iterations are performed by considering (1.5) in suitable sets of collocation points. As a consequence, the corresponding fixed points are steps of high-order implicit Runge-Kutta schemes on each subinterval. For the particular choice of Gauss-Legendre points, we obtain the Gauss-type methods, which for the unitary evolution of (1.1) have favorable norm and energy preservation properties.

We arrive at an iterative scheme that successively solves on each subinterval using Picard iteration or a variant known as spectral deferred correction, where the approximate solutions in each collocation point are stored in low-rank form. The corresponding ranks increase in each iteration by a factor depending on $J$ and representation rank of $V$. We subsequently reduce them in each step by a soft thresholding of singular values, which is a non-expansive operation in $L^2(\Omega)$, with a new technique for controlling the thresholds inspired by [9].

In our main result, Theorem 27, we show that with the strategy for adjusting error tolerances and thresholding parameters that we propose for the numerical scheme, we achieve a guaranteed error

$$ \varepsilon = C_{J,V}h^{J+1}\exp(c_{J,V}T) $$
for some positive constants $C_{J,V}$ and $c_{J,V}$, and obtain quasi-optimal bounds on the ranks $r(t)$ of all computed iterates. That is, if for an $s > 0$,

$$ \max_t \hat{r}_\eta(t) \lesssim \eta^{-\frac{1}{s}} $$

with the best approximation rank for $\eta > 0$ as defined above, then also

$$ \max_t r(t) \lesssim h^{-3}(1 + J \operatorname{rank}(V)) \eta^{-\frac{1}{s}}, $$

where $\operatorname{rank}(V)$ is the rank of the function $V$ (for simplicity, we here assume that the potential is sufficiently smooth and of bounded ranks). In other words, the largest ranks produced by the numerical scheme remain bounded up to a constant by the ranks of the result of applying the fixed point mapping to a low-rank best approximation of the solution. Since the constant depends also on $h^{-3}$, the subinterval size should not be chosen too small, and we rather aim to enforce convergence by refining the approximation in each subinterval, that is, by increasing $J$. In the more favorable case of exponential-type decrease of singular values,

$$ \max_t \hat{r}_\eta(t) \lesssim (1 + |\log \eta|)^\beta $$

with $\beta > 0$, then accordingly

$$ \max_t r(t) \lesssim h^{-3}(1 + J \operatorname{rank}(V)) (1 + |\log \eta|)^\beta. $$

We focus on the case of time-dependent Schrödinger equations, for which we give a detailed analysis, but also describe the modifications required for the case of second-order parabolic problems motivated by (1.2), where we use an adapted strategy based on Gauss-Radau points.

## 1.2. Novelty and relation to previous work.
While the basic building blocks of the scheme considered here are known, they are combined in a particular manner suitable for low-rank approximations, and the main novelty lies in a careful choice of error tolerances and rank truncation thresholds so that guaranteed convergence is achieved while maintaining estimates in terms of the ranks of best low-rank approximations of comparable accuracy.

Such bounds have been obtained previously for a specific type of method based on a space-time variational formulation. In such methods, instead of proceeding by time steps to update a spatial approximation, a joint approximation in spatial and temporal variables is gradually refined towards a sufficiently accurate approximation of the entire evolution. In the context of low-rank approximations, such methods were considered in [3, 10, 19], though without estimates for the computed ranks. However, such approaches are very well suited to ensure an appropriate balance between approximation errors and tensor ranks, and a method ensuring near-optimal ranks for high-dimensional parabolic problems has been proposed and analyzed in [8]. However, suitable well-posed space-time formulations do not exist for every problem of interest, and they often strongly constrain the norms in which approximation errors can be controlled.

Here we are thus interested in methods that do not require a well-posed space-time variational formulation as a starting point. With this restriction, to the best of our knowledge, all previously known methods either ensure convergence or control ranks, but not both.

Dynamical low-rank approximation [31] falls into the latter category of controlling ranks but not convergence. It is based on approximations by substitute problems on manifolds of matrices (or tensors) of fixed ranks. In such schemes, the right-hand side is projected onto the tangent space of the current base point on the manifold, and error bounds can only be obtained under the assumption that the error in this projection is small, which need not be the case (see, for example, [6, Rem. 7.2] and [33, Ex. 1]). The Basis Update & Galerkin (BUG) method [15] (see also [39] for high order extensions) is closely related to splitting integrators for dynamical low-rank approximation [37], but can be naturally modified to include the adaptation of ranks [14]. However, the issue that errors of tangent space projections need to be assumed to be small is still present in these methods.
There also exists a variety of methods that directly use low-rank representations of the right-hand side, rather than only its projections. Step-truncation methods are based on performing time steps of a standard scheme (such as a Runge-Kutta method) in tensor format and then performing rank truncation, as in [18, 33, 40]. A variation with conceptual similarities to BUG integrators are methods that use the right-hand side for basis enrichment and thus avoid the restrictions imposed by tangent space projections, as in [4, 44]. While in principle, such methods can be ensured to converge, the rank truncation tolerances need to be tightened with decreasing time step length, and the resulting ranks are unclear.

An important conceptual ingredient for the feasibility of our approach for PDE problems are concepts related to exponential integrators [30], in particular in a form originally suggested by [34]. Modern exponential integrators, however, generally rely on inverses of operators in the approximation of exponentials. Such methods have been used, for example, in combination with approximations on fixed-rank manifolds in [12, 42]. Inverses of operators such as the Laplacian, however, are generally not easy to realize in low-rank format – in fact, one of the most reliable techniques for obtaining low-rank approximations of such inverses consists in a reduction to sums of exponentials, which in view of (1.4) is easier to handle in low-rank format. For this reason, unlike standard exponential integrators, the methods considered here rely only on operator exponentials of Kronecker rank one.

While we prove our main results for a scheme based on Picard iteration, we also consider a variant based on spectral deferred correction (SDC), an approach that was originally proposed in [20], see also [11, 25, 35]. Connections to Runge-Kutta methods and Picard iteration are highlighted in particular in [13, 43]. SDC has also recently been used in the context of low-rank methods for the construction of higher-order step-truncation methods [36], considering also soft thresholding following [9], but without an analysis of approximation ranks.

The approach can be interpreted as a hybrid between time-stepping schemes and methods based on space-time variational formulations. It retains flexibility in its application to different classes of PDEs and achieves error control pointwise in time.

### 1.3. Outline.
In Section 2, we begin by devising the integral formulations and their discretizations on subintervals that form the basis of our schemes. In Section 3, we describe rank reduction techniques that we use and their basic properties. The proof of our main result is conducted in Section 4. In Section 5, we consider the adaptation of our approach to parabolic problems. Finally, we present results of numerical experiments in Section 6 and conclude with some open problems.

## 2. ITERATIVE FIXED-POINT METHODS

Our aim is to construct high-order methods for low-rank representations of evolution problems

$$ \partial_t u(t) = \mathcal{A}_t(u(t)) $$

on $L^2(\Omega)$, with a potentially unbounded operator $\mathcal{A}_t$. In particular, we focus on the case of affine linear differential equations of the form

$$ \mathcal{A}_t(v(t)) = (\mathcal{L} + \mathcal{G}_t)v(t) + f(t) $$

where $\mathcal{L}$ and $\mathcal{G}_t$ each denote linear operators on $L^2(\Omega)$ and $f(t)$ might encode some additional nonlinear information. We further assume that the operator $\mathcal{L}$ captures the majority of the stiffness of the resulting differential operator. In order to simplify the notation, we denote by $|| \cdot ||$ the $L^2(\Omega)$ norm and by $\langle \cdot, \cdot \rangle$ the corresponding inner product. Additionally, we define

$$ W(s, t) = L^\infty(s, t; L^2(\Omega)), \quad \text{for } 0 \le s \le t \le T. $$

While it is rather complicated to provide both rank bounds and a sufficient accuracy, for general time-stepping methods, works like [9] and [7] suggest that iterative solvers such as fixed-point iterations in low-rank format are able to retain quasi-optimal rank control on the iterates. The idea would thus be to subdivide the time interval $[0, T]$, on which we would like to evolve the above problem, into $n \in \mathbb{N}$ smaller intervals of size $h > 0$ and to
approximate solutions on those iteratively using a fixed-point iteration combined with low-rank techniques. Once a sufficiently accurate solution has been obtained on a subinterval, one can move to the next one and repeat the procedure.

A classical approach for bounded $A_t$ would consist in using the integral formulation

$$ (2.1) \qquad u(t) = u(t_0) + \int_{t_0}^{t} A_t(u(s)) ds $$

for $t_0 = 0$ and Picard iteration to successively obtain approximations at the discrete points $h, 2h, \dots, Nh \in [0, T]$. This basic scheme is then often combined with approximating the exact solution $u(t)$ on a resulting subinterval of the form $t \in [nh, (n+1)h]$ using Chebyshev polynomial series, see for example [16] and [17] for early contributions.

It is well-known that stability constraints for the convergence of fixed-point iterations depend strongly on the Lipschitz constant of the operator $A_t$. It is thus natural to try to reformulate the integral formulation, which is used as the base of iterative solvers in order to reduce the imposed step size restriction, since the ranks of computed low-rank approximations can potentially increase exponentially in the number of time intervals.

**2.1. Modified fixed-point formulation.** In physical applications, such as the time-dependent Schrödinger equation, it turns out that the stiffness of the time-independent linear operator $\mathcal{L}: L^2(\Omega) \to L^2(\Omega)$ is mainly responsible for the magnitude of the Lipschitz constant. Applying the discrete variation of constants formula to (2.1), one obtains

$$ u(t) = e^{(t-t_0)\mathcal{L}}u(t_0) + \int_{t_0}^{t} e^{(t-s)\mathcal{L}}(G_s u(s) + f(s)) ds. $$

The favorability of this formulation can be seen by comparing the resulting Lipschitz constants of our main application.

**Example 1** (Schrödinger equation with time-dependent potential). In the setting of the Schrödinger equation (1.1) with a potentially time-dependent potential $V_t$ we have

$$ \mathcal{L} = i\Delta, \quad G_t = -iV_t \quad \text{and} \quad f = 0. $$

The resulting Lipschitz constant

$$ \| -\Delta + V_t \|_{W(0,T) \to W(0,T)} $$

is unbounded due to the Laplacian and in computations, it will depend strongly on the spatial discretization: if we use for example an $L^2$-orthonormal basis expansion with sine product basis functions, the Laplacian in matrix form becomes a diagonal operator with entries being proportional to $\Delta_{k,k} \approx k^2$, which is mainly responsible for the resulting stiffness. Thus, we would need to progressively decrease the stepsize $h$ in order to obtain a better spatial resolution.

This issue can be circumvented by instead considering Duhamel's formula

$$ u(t) = e^{i(t-t_0)\Delta}u(t_0) - i \int_{t_0}^{t} e^{i(t-s)\Delta}V_s u(s) ds, $$

which is equivalent to the Picard integral formula

$$ (2.2) \qquad v(t) = v(t_0) + \int_{t_0}^{t} \mathcal{F}_s v(s) ds, \qquad \mathcal{F}_s = -ie^{-is\Delta}V_s e^{is\Delta}, $$

for the twisted variable $v(t) = e^{-it\Delta}u(t)$ which formally satisfies

$$ (2.3) \qquad \partial_t v(t) = \mathcal{F}_t v(t) = -ie^{-it\Delta}V_t e^{it\Delta}v(t). $$

This change of variable, also known as Lawson method for Runge-Kutta methods [34], reduces the stiffness of the PDE by explicitly integrating the term in $-i\Delta$. Here, $V_t$ should be understood as an operator and not as a function, hence $e^{-it\Delta}V_t g(t) = e^{-it\Delta}(V_t g)(t)$ for $g \in W(0, T)$. In this setting, we apply the resulting scheme on the twisted variable and after having obtained an approximation on the whole time interval $[0, T]$, we can revert the applied exponential transformations and 'untwist' the approximations.
Due to the norm preservation property

$$ \|e^{it\Delta}u(t)\| = \|u(t)\|, \quad \text{for all } t \in \mathbb{R}, $$

the Lipschitz constant reduces to $(t_0 - t)$ multiplied with

$$ (2.4) \qquad C_V = \|V_t\|_{W(0,T) \to W(0,T)}. $$

In particular, it is robust with respect to the spatial discretization for a bounded potential. Note that more realistic potentials with Coulomb singularities can be treated using Strichartz estimates if we consider $L^2(\Omega)$ norms in time, see [19]. In the following, we restrict our considerations to the case of a time-independent potential for notational convenience, even though the framework can also be extended to the case of a time-dependent potential.

**Example 2** (Second-order parabolic equation with anisotropic diffusion). We also consider a parabolic working example on $L^2(\Omega_1 \times \Omega_2)$ with anisotropic diffusion, given by

$$ (2.5) \qquad \partial_t u(t) = \operatorname{div} \begin{pmatrix} a & b \\ b & a \end{pmatrix} \nabla u(t) + f $$

for constants $a > 2b > 0$ and a source term $f \in L^2(\Omega_1 \times \Omega_2)$, which can be decomposed as

$$ \mathcal{L} = a\Delta, \quad \mathcal{G}_t = \mathcal{G} = 2b\partial_{x_1}\partial_{x_2} \quad \text{and} \quad f(t) = f. $$

We can again use Duhamel's formula, here in the form

$$ (2.6) \qquad u(t) = e^{a(t-t_0)\Delta}u(t_0) + \int_{t_0}^{t} e^{a(t-s)\Delta} (2b\partial_{x_1}\partial_{x_2}u(s) + f) ds, $$

as our starting point. However, an approach using twisted variables as in (2.2) is not applicable due to the unboundedness of $e^{-t\Delta}$ for $t > 0$, and for the same reason, additional care is needed to use $e^{a(t-s)\Delta}$ only when $t > s$.

Since parabolic problems thus require a separate treatment, we postpone the discussion of this case to Section 5. In the following, we first focus on the time-dependent Schrödinger equation.

**2.2. Fixed-point formulations.** Given a subinterval of the form $[t_0, t_0 + h]$, a (complex-valued) function $u \in W(t_0, t_0 + h)$, and a strictly increasing sequence of nodes

$$ t_0 \le t_1 < \cdots < t_J \le t_0 + h $$

for some $J \in \mathbb{N}$, we can define the Lagrange polynomial interpolation of $u$ as

$$ t \mapsto \sum_{j=1}^{J} u(t_j)\ell_j(t), \qquad \ell_j(t) = \prod_{m \neq j} \frac{t - t_m}{t_j - t_m}. $$

**2.2.1. Gauss-Legendre-Picard iteration.** A classical fixed-point iteration following directly from the Picard integral formulation in the twisted variables (2.2) is then given by

$$ (2.7) \qquad \Phi_{\text{Picard}}(\mathbf{v})(t) = v(t_0) + \int_{t_0}^{t} \sum_{j=1}^{J} \mathcal{F}_{t_j} v_j \ell_j(s) ds, $$

for $\mathbf{v} = (v_j)_{1 \le j \le J}$ an approximation of $(v(t_j))_{1 \le j \le J}$. This allows us to formulate a fixed-point application $\Phi_{\text{Picard}} : L^2(\Omega)^J \to L^2(\Omega)^J$ as

$$ (2.8) \qquad \Phi_{\text{Picard}}(\mathbf{v})_j = \Phi_{\text{Picard}}(\mathbf{v})(t_j) = v(t_0) + \sum_{m=1}^{J} \omega_{j,m} \mathcal{F}_{t_m} v_m, $$

where the quadrature weights are given by

$$ \omega_{j,m} = \int_{t_0}^{t_j} \ell_m(s) ds, \quad \text{for } 1 \le j, m \le J. $$
Naturally, one would like to maximize the degree of accuracy of the quadrature formula

$$ (2.9) \qquad \int_{t_0}^{t_0+h} g(s) \, ds \approx \sum_{m=1}^{J} \omega_m g(t_m), \quad \text{for } g \in W(t_0, t_0+h), $$

when choosing the quadrature nodes $t_1, t_2, \dots, t_J$ and corresponding weights

$$ \omega_m = \int_{t_0}^{t_0+h} \ell_m(s) \, ds, \quad \text{for } 1 \le m \le J. $$

It is well-known that the maximal degree of accuracy, which can be attained by such a quadrature formula is $2J - 1$, and that this bound can be reached only by the so-called *Gaussian quadrature formulas*, see for example [41] for scalar-valued functions.

**Example 3 (Gauss-Legendre quadrature).** For this family of quadrature formulas, the $J$ nodes are defined as the roots of the $J$th Legendre polynomial, which are contained in the open interval $(t_0, t_0+h)$.

The Gauss-Legendre quadrature nodes and weights can be computed efficiently using the Golub-Welsch algorithm [21]. Defining the tridiagonal matrix

$$ (2.10) \qquad T_J = \begin{pmatrix} t_0 + h/2 & \beta_1 & & & & 0 \\ \beta_1 & t_0 + h/2 & \beta_2 & & & \\ & \beta_2 & t_0 + h/2 & \ddots & & \\ & & \ddots & \ddots & \beta_{J-2} & \\ & & & \beta_{J-2} & t_0 + h/2 & \beta_{J-1} \\ & 0 & & \beta_{J-1} & t_0 + h/2 & \end{pmatrix}, $$

where $\beta_m = \frac{mh}{2\sqrt{4m^2-1}}$, it can be shown that the $J$ eigenvalues of $T_J$ correspond exactly to the Gauss-Legendre nodes $t_1, t_2, \dots, t_J$ and that the weights can be expressed as $\omega_m = 2v_{m,1}^2/\|v_m\|_2^2$, where $v_{1,1}, v_{2,1}, \dots, v_{J,1}$ denote the first entries of the eigenvectors $v_1, v_2, \dots, v_J$, respectively.

In the following, we thus restrict our considerations to the use of the Gauss-Legendre quadrature for the approximation of the time-dependent Schrödinger equation. The advantages of these particular nodes will be discussed in the subsequent sections.

**Remark 4 (Spectral integration matrix).** We refer to [20] and the references therein for a discussion on the numerical stability properties of the spectral integration matrix $(\omega_{j,m})_{1\le j,m\le J}$. By translation and scaling symmetries, its entries do not depend on $t_0$ and are proportional to $h$; they can be precomputed exactly using Gauss-Legendre quadrature in the smaller intervals $[t_0, t_j]$, $1 \le j \le J$.

**Remark 5 (Extended fixed-point iteration and transition to boundary value).** To transition to the next interval, following the approach described in [25], we approximate the boundary value by

$$ (2.11) \qquad \Phi_{\text{Picard}}(\mathbf{v})(t_0+h) = v(t_0) + \sum_{m=1}^{J} \omega_m \mathcal{F}_{t_m} v_m. $$

One could also define the fixed-point formulation on the $J+1$ nodes $t_1, \dots, t_J, t_0+h$, however the last value is not actively used inside the fixed-point iterations.

### 2.2.2. *Spectral Deferred Correction.*
Another possible fixed-point iteration, first introduced in [20], is spectral deferred correction. Observe that (2.8) can be rewritten inductively in $j$ as

$$ \mathbf{\Phi}_{\text{Picard}}(\mathbf{v})_j = \mathbf{\Phi}_{\text{Picard}}(\mathbf{v})_{j-1} + \sum_{m=1}^{J} \tilde{\omega}_{j,m} \mathcal{F}_{t_m} v_m, $$
where $\Phi_{\text{Picard}}(\mathbf{v})_0 = v_0$ and

$$ \tilde{\omega}_{j,m} = \omega_{j,m} - \omega_{j-1,m} = \int_{t_{j-1}}^{t_j} \ell_m(s) \, ds. $$

Given a low-order integrator $\Psi$ such that $\Psi_{s,t}(u(s)) \approx u(t)$, SDC is constructed by adding to the above formula a correction, obtained by integrating the residual error from $t_{j-1}$ to $t_j$, yielding an approximation of higher order [11]. In order to keep control on the ranks of the updates, we restrict our considerations to explicit low-order integrators $\Psi$. The SDC iterations can then be written as

$$ (2.12) \qquad \Phi_{\text{SDC}}(\mathbf{v})_j = \Phi_{\text{SDC}}(\mathbf{v})_{j-1} + \sum_{m=1}^{J} \tilde{\omega}_{j,m} \mathcal{F}_{t_m} v_m + \Psi_{t_{j-1},t_j} (\Phi_{\text{SDC}}(\mathbf{v})_{j-1} - v_{j-1}), $$

where $\Phi_{\text{Picard}}(\mathbf{v})_0 = v_0$. While the explicit Euler method

$$ \Psi_{t_{j-1},t_j} = (t_j - t_{j-1}) \mathcal{F}_{t_{j-1}} $$

seems quite standard [13], in the case of the Schrödinger equation, or more generally Hamiltonian systems, we prefer the symmetric integrator

$$ \Psi_{t_{j-1},t_j} = (t_j - t_{j-1}) \mathcal{F}_{(t_{j-1}+t_j)/2} $$

motivated by [32], since the midpoint approximation of $\mathcal{F}$ is expected to be more accurate. Note that we do not exactly preserve the energy, time-reversibility or symplecticity, since the integrator remains explicit. However, we only need to apply the potential operator $V$ once, which is advantageous for the rank control.

The transition to the subinterval boundary is then carried out analogously to (2.11).

### 2.3. Reformulation in matrix form.
In order to better compare the Picard and SDC iterations and understand how the corresponding fixed-point formulations are defined on the Gauss-Legendre nodes, following [29], we derive matrix versions of the corresponding fixed-point formulations.

For the Picard iteration, we simply have

$$ (2.13) \qquad \Phi_{\text{Picard}}\mathbf{v} = \mathbf{1} \otimes v_0 + \mathbf{\Sigma}\mathbf{v}, $$

where $\mathbf{\Sigma} = (\omega_{j,m}\mathcal{F}_{t_m})_{1\le j,m\le J}$. Regarding SDC, we start by iterating (2.12) to obtain

$$ \Phi_{\text{SDC}}(\mathbf{v})_j = \Phi_{\text{Picard}}(\mathbf{v})_j + \sum_{i=1}^{j} \Psi_{t_{i-1},t_i} (\Phi_{\text{SDC}}(\mathbf{v})_{i-1} - v_{i-1}). $$

One can then rearrange the above equation as

$$ (2.14) \qquad (\mathbf{I} - \mathbf{\Psi}) \Phi_{\text{SDC}}(\mathbf{v}) = \mathbf{1} \otimes v_0 + (\mathbf{\Sigma} - \mathbf{\Psi})\mathbf{v}, $$

where the matrix $\mathbf{\Psi} = (\mathbf{1}_{j>i} \Psi_{t_i, t_{i+1}})_{1\le j,i\le J}$ is lower triangular, allowing to explicitly compute $\Phi_{\text{SDC}}(\mathbf{v})$ one coordinate at a time.

These matrix equations will allow us to adapt the soft thresholded fixed point method from [9], originally designed for the Richardson iteration, to our setting of time integration, as will be seen in particular in Lemma 18.

### 2.4. Collocation method limit of the fixed-point iterations.
Considering the formulations (2.13) and (2.14), it becomes evident that the solutions to the fixed-point problems $\Phi_{\text{Picard}}(\mathbf{v}^*) = \mathbf{v}^*$ and $\Phi_{\text{SDC}}(\mathbf{v}^*) = \mathbf{v}^*$ coincide and satisfy

$$ (2.15) \qquad (\mathbf{I} - \mathbf{\Sigma})\mathbf{v}^* = \mathbf{1} \otimes v_0. $$

It turns out that this equation is equivalent to an implicit collocation formulation, see also [29]. Indeed, it holds that $\mathbf{v}^* = (\Phi_{\text{Picard}}(\mathbf{v}^*)(t_j))_{1\le j\le J}$, where $\Phi_{\text{Picard}}(\mathbf{v}^*)$ is a polynomial of degree $J$ in time satisfying

$$ (2.16) \qquad \partial_t \Phi_{\text{Picard}}(\mathbf{v}^*)(t_m) = \sum_{j=1}^{J} \mathcal{F}_{t_j} v_j^* \ell_j(t_m) = \mathcal{F}_{t_m} v_m^* = \mathcal{F}_{t_m} (\Phi_{\text{Picard}}(\mathbf{v}^*)(t_m)), $$
that is, $\Phi_{\text{Picard}}(\mathbf{v}^*)$ is a solution to (2.3) at times $t_m$, for all $1 \le m \le J$. As a consequence, it is also the solution to an implicit Runge-Kutta method with $J$ stages, whose Butcher tableau is given by

$$ (2.17) \qquad \frac{\mathbf{c}}{\mathbf{b}^T} \qquad \text{where} \quad \mathbf{A}_{m,j} = \frac{\omega_{j,m}}{h}, \quad \mathbf{b}_j = \frac{\omega_j}{h} \quad \text{and} \quad c_m = \frac{t_m - t_0}{h}, $$

see [26, Theorem 7.7]. Such methods are of particular interest, since they can reach, depending on the underlying quadrature, up to double the order of convergence of the number of interpolating conditions or stages, that is, $2J$. The upper bound can only be reached using Gaussian quadrature formulas, in which case the resulting Runge-Kutta methods are denoted as *Gauss methods*.

As for any Runge-Kutta scheme, we can analyse the linear stability by considering the scalar test problem

$$ (2.18) \qquad \dot{y}(t) = \lambda y(t), \quad \lambda \in \mathbb{C}, $$

for which the approximation generated by the scheme after one time step $h$ is given by $R(h\lambda)y(0)$ for some rational function $R$. It can be shown that the Gauss collocation method defined by (2.17) is *A-stable*, that is, the stability function $R$ satisfies

$$ |R(z)| \le 1 \quad \text{for all } z \in \mathbb{C} \text{ with } \text{Re}(z) \le 0. $$

It is in addition *isometry-preserving* in the sense that

$$ |R(z)| = 1 \quad \text{for } \text{Re}(z) = 0. $$

As the Schrödinger equation preserves the $L^2(\Omega)$ norm, the eigenvalues of its evolution operator are purely imaginary, and the above property translates into preservation of norm for the Gauss method. This will be addressed further in Lemma 15. Altogether, these properties make the Gauss-Legendre nodes an adequate choice in our setting.

**Remark 6** (Original formulation of SDC). In the original paper [20, Rem. 4.1] the transition to the subinterval boundary is carried out by evaluation of the interpolation polynomial, that is, $\sum_{m=1}^J v_m^k \ell_m(t_0 + h)$. However, as observed in several numerical experiments, the resulting method is not isometry-preserving, which is particularly important when computing approximations to the time-dependent Schrödinger equation. Moreover, the resulting scheme is then also no longer equivalent to a Gauss method and loses the superconvergence property at the subinterval boundaries.

**Remark 7** (Advantages of using SDC over Picard iteration). It was shown in [13] that both Picard iteration and SDC (using the explicit Euler as the low-order integrator) pick up at least one order of accuracy in each iteration. However, SDC is expected to converge even faster in the neighborhood of the fixed point, allowing one to reduce the number of iterations. Additionally, [13] suggest that using higher-order base solvers can indeed lead to the SDC iteration being able to pick up more than one order of accuracy with each iteration. Moreover, the used solver also has a direct influence on the resulting stability properties of the scheme. It was also proven in [25] that if the fixed-point method is iterated until convergence, the maximal order of convergence $2J$ for Gauss-Legendre nodes, can in fact be reached, see also Lemma 16. In that perspective, SDC can also be viewed as a preconditioner for the Picard iteration, as can be deduced from the matrix formulation (2.14).

**2.5. Convergence analysis of the fixed-point iterations.** We are now in the position to show the contractivity of the constructed fixed-point iterations, for small enough step sizes $h$. For $\mathbf{v} \in L^2(\Omega)^J$, define the norm

$$ \|\mathbf{v}\|_J = \max_{1 \le j \le J} \|v_j\|, $$

where we recall that $\| \cdot \| = \| \cdot \|_{L^2(\Omega)}$.
**Proposition 8** (Contraction constant for the Schrödinger equation). Let

$$ (2.19) \qquad \Lambda_J = \frac{1}{h} \int_0^h \sum_{m=1}^J |\ell_m(s)| ds \quad \text{and} \quad \rho = h \Lambda_J C_V. $$

Then the fixed-point iterations $\mathbf{\Phi}_{\text{Picard}}$ and $\mathbf{\Phi}_{\text{SDC}}$ defined by (2.8) and (2.12) satisfy

$$ \|\mathbf{\Phi}_{\text{Picard}}\mathbf{v} - \mathbf{\Phi}_{\text{Picard}}\mathbf{w}\|_{W(0,h)} \leq \rho \|\mathbf{v} - \mathbf{w}\|_J $$

and

$$ \|\mathbf{\Phi}_{\text{SDC}}\mathbf{v} - \mathbf{\Phi}_{\text{SDC}}\mathbf{w}\|_J \leq 2\rho e^{\rho} \|\mathbf{v} - \mathbf{w}\|_J $$

for any $\mathbf{v}, \mathbf{w} \in L^2(\Omega)^J$.

*Proof.* For the Picard iteration, we see that for any time $t \in [0, h]$, the expression

$$ \|\mathbf{\Phi}_{\text{Picard}}(\mathbf{v})(t) - \mathbf{\Phi}_{\text{Picard}}(\mathbf{w})(t)\| = \left\| \int_0^t \sum_{m=1}^J \ell_m(s) \mathcal{F}_{t_m}(v_m - w_m) ds \right\| $$

is bounded by

$$ \int_0^t \sum_{m=1}^J |\ell_m(s)| ds \ C_V \|v_m - w_m\| \leq \rho \ \|\mathbf{v} - \mathbf{w}\|. $$

For the SDC iteration, letting $\phi_j = (\mathbf{\Phi}_{\text{SDC}}\mathbf{v})_j - (\mathbf{\Phi}_{\text{SDC}}\mathbf{w})_j$, we have

$$ \begin{aligned} \|\phi_j\| &\leq \|\phi_{j-1}\| + \|\Psi_{t_{j-1},t_j}(\phi_{j-1} - v_{j-1} + w_{j-1})\| + \sum_{m=1}^J |\tilde{\omega}_{j,m}| \|\mathcal{F}_{t_m}(v_m - w_m)\| \\ &\leq (1 + \tau_j C_V) \|\phi_{j-1}\| + \tau_j C_V \|v_{j-1} - w_{j-1}\| + \sum_{m=1}^J |\tilde{\omega}_{j,m}| C_V \|v_m - w_m\|. \end{aligned} $$

By induction, using the upper bounds $1 + c \leq e^c$ for all $c \in \mathbb{R}$ and $\sum_{j=1}^J \tau_j \leq h$, we obtain

$$ \max_{1 \leq j \leq J} \|\phi_j\| \leq e^{hC_V} \left( \sum_{j,m=1}^{J} |\tilde{\omega}_{j,m}| + h \right) C_V \|\mathbf{v} - \mathbf{w}\|_J, $$

and we conclude by observing that

$$ \sum_{j,m=1}^{J} |\tilde{\omega}_{j,m}| = \sum_{j=1}^{J} \int_{t_{j-1}}^{t_j} \sum_{m=1}^{J} |\ell_m(s)| ds \leq h \Lambda_J \quad \text{and} \quad 1 = \frac{1}{h} \int_{0}^{h} \sum_{m=1}^{J} \ell_m(s) ds \leq \Lambda_J. \quad \square $$

*Remark 9* (Choice of quadrature nodes). The constant $\Lambda_J$ defined in (2.19) can be bounded from above by the Lebesgue constant. The asymptotic behavior of the latter is in the order $O(\sqrt{J})$ in the case of the Gauss-Legendre and Radau-Legendre quadrature points, see [24]. Since this constant is a limiting factor for the step size $h$, one could also think of using Chebyshev nodes, whose asymptotic behavior is in $O(\log J)$. Unfortunately, using Chebyshev nodes does not result in a scheme of order $2J$ at the subinterval boundaries. However, the averaged constant $\Lambda_J$ differs much less between the Gauss-Legendre and Chebyshev points, and can in both cases be bounded by $2$ for $J \leq 10$, which is sufficient for our purposes.

Moreover, the numerical experiments performed in [5] for the particular application of celestial mechanics suggest that the usage of an iteration based on Gauss-Legendre nodes leads to improved convergence properties over the Chebyshev nodes.

### 3. THRESHOLDING STRATEGY AND RESULTING ALGORITHMS

We now turn to our strategy for reducing the ranks of the intermediate approximations. Since we use a fixed-point iteration in order to successively approximate the true solution, it seems natural to combine the iterations with a rank control strategy.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Illustration of the soft thresholding operation applied to some element $v \in L^2(\Omega_1 \times \Omega_2)$.</figcaption></figure>

**3.1. Rank reduction.** For every $v \in L^2(\Omega_1 \times \Omega_2)$, there exists a singular value decomposition

$$ (3.1) \qquad v = \sum_{k=1}^{\infty} \sigma_k v_k^{(1)} \otimes v_k^{(2)}, $$

with the singular values $\sigma_1 \ge \sigma_2 \ge \dots \ge 0$ and orthonormal families of singular vectors $(v_k^{(1)})_{k \in \mathbb{N}}$ in $L^2(\Omega_1)$ and $(v_k^{(2)})_{k \in \mathbb{N}}$ in $L^2(\Omega_2)$. We have already noted the well-known fact that the truncation of the singular value decomposition to rank $r$ yields the best approximation of rank $r$ in $L^2(\Omega_1 \times \Omega_2)$, with

$$ \left\| v - \sum_{k=1}^{r} \sigma_k v_k^{(1)} \otimes v_k^{(2)} \right\|^2 = \sum_{k>r} \sigma_k^2. $$

We define the nonlinear *hard thresholding* operator on $L^2(\Omega_1 \times \Omega_2)$ by

$$ (3.2) \qquad \mathcal{H}_{\alpha}(v) = \sum_{\substack{k \in \mathbb{N} \\ \sigma_k > \alpha}} \sigma_k v_k^{(1)} \otimes v_k^{(2)}, $$

and the corresponding *recompression* operator by

$$ (3.3) \qquad \mathcal{R}_{\delta}(v) = \sum_{1 \le k \le r_{\delta}} \sigma_k v_k^{(1)} \otimes v_k^{(2)}, \quad \text{where} \quad r_{\delta} = \min\{r \in \mathbb{N}_0: \sum_{k>r} \sigma_k^2 \le \delta^2\}. $$

In general, the addition of two low-rank tensor representations leads to an increase in the rank parameter: the sum of two elements in $L^2(\Omega_1 \times \Omega_2)$ with respective ranks $r_1$ and $r_2$ has rank up to $r_1 + r_2$. Even though the used low-rank format is closed under addition, the resulting representation may be redundant. This redundancy can be removed in practice by applying a recompression operator (3.3) with a sufficiently small tolerance $\delta$ to the sum.

Another possibility to reduce the rank of an element $v \in L^2(\Omega_1 \times \Omega_2)$ would be to apply the *soft thresholding* operator defined by

$$ (3.4) \qquad \mathcal{S}_{\alpha}(v) = \sum_{k \in \mathbb{N}} \max\{\sigma_k - \alpha, 0\} v_k^{(1)} \otimes v_k^{(2)}. $$

A visualization of the application of the operator to a given sequence of singular values can be found in Figure 1. The errors resulting from the application of the truncation operators are related by

$$ \|v - \mathcal{S}_{\alpha}(v)\|^2 = \|v - \mathcal{H}_{\alpha}(v)\|^2 + \alpha^2 \operatorname{rank}(\mathcal{H}_{\alpha}(v)). $$

Note that $\mathcal{S}_{\alpha}(v)$ and $\mathcal{H}_{\alpha}(v)$ have the same rank $r = |\{k \in \mathbb{N}: \sigma_k > \alpha\}|$.

The following lemma allows to relate the errors caused by soft thresholding an element in $L^2(\Omega_1 \times \Omega_2)$ with different thresholds.

**Lemma 10.** For $0 < \alpha \le \beta$ it holds that

$$ \|v - \mathcal{S}_{\alpha}(v)\| \le \|v - \mathcal{S}_{\beta}(v)\| \le \frac{\beta}{\alpha} \|v - \mathcal{S}_{\alpha}(v)\|, $$
for all $v \in L^2(\Omega_1 \times \Omega_2)$.

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Visualization of Lemma 10: deviations between the singular values of some element $v \in L^2(\Omega_1 \times \Omega_2)$ and the output of a soft thresholding operator applied to $v$. The black dots represent coinciding values of the red and blue dots.</figcaption></figure>

*Proof.* For $v$ of the form (3.1), we have

$$ \|v - S_{\alpha}(v)\|^2 = \sum_{k \in \mathbb{N}} \min\{\alpha, \sigma_k\}^2, $$

hence the result follows from the inequalities $\min\{\alpha, \sigma_k\} \le \min\{\beta, \sigma_k\} \le \frac{\beta}{\alpha} \min\{\alpha, \sigma_k\}$. $\square$

3.2. **Iterative schemes.** An additional feature of the soft thresholding operator is its non-expansiveness property, that is, for all $\alpha \ge 0$ and $v, w \in L^2(\Omega)$, it holds that

$$ \|S_{\alpha}(v) - S_{\alpha}(w)\| \le \|v - w\| \quad \text{and} \quad \|(I - S_{\alpha})(v) - (I - S_{\alpha})(w)\| \le \|v - w\|, $$

where $\mathcal{I}$ denotes the identity operator. We refer to [9] and [6] for a proof. This feature can now be exploited for the rank control within our iterative methods: using the previously introduced notation, we set for the Gauss-Legendre-Picard iteration

$$ (3.5) \qquad \mathbf{v}^{k+1} = S_{\alpha_k}(\mathbf{\Phi}_{\text{Picard}}(\mathbf{v}^k)) = (S_{\alpha_k}(\mathbf{\Phi}\mathbf{v}^k(t_j)))_{1 \le j \le J}, $$

that is, we apply the soft thresholding operator with threshold $\alpha_k$ component-wise to each of the updated approximations at the Gauss-Legendre nodes. Soft thresholding the fixed-point iterations as in (3.5) does not interfere with the contractivity property, that is,

$$ \|\mathbf{v}^{k+1} - \mathbf{w}^{k+1}\| = \|S_{\alpha_k}(\mathbf{\Phi}_{\text{Picard}}(\mathbf{v}^k)) - S_{\alpha_k}(\mathbf{\Phi}_{\text{Picard}}(\mathbf{w}^k))\| \le \rho \|\mathbf{v}^k - \mathbf{w}^k\|, $$

for all $\mathbf{v}^k, \mathbf{w}^k \in L^2(\Omega)^J$, where $\rho$ is contraction factor from (2.19).

For the SDC iteration, we instead apply the soft thresholding operator after updating each of the Gauss-Legendre nodes. Then (2.12) modifies according to

$$ (3.6) \qquad v_j^{k+1} = S_{\alpha_k}(\phi_j^k), \quad \phi_j^k = v_{j-1}^{k+1} + \sum_{m=1}^{J} \tilde{\omega}_{j,m} \mathcal{F}_{t_m} v_m^k + \Psi_{t_{j-1}, t_j} (v_{j-1}^{k+1} - v_{j-1}^k). $$

We observe in the numerical experiments that applying $S_{\alpha_k}$ after each of the $J$ updates results in much smaller ranks of the intermediate approximations than thresholding only at the end of the iteration, as we would do with $S_{\alpha_k}(\mathbf{\Phi}_{\text{SDC}}(\mathbf{w}^k))$. However this makes the analysis more intricate, which is also the reason why we will focus on the Picard iteration in the following.

Using the soft thresholding operator as in (3.5) with fixed $\alpha_k = \alpha$ at all iterations, results in a modified fixed-point limit $\mathbf{v}_\alpha \in L^2(\Omega)$ defined by $\mathbf{v}_\alpha = S_\alpha(\mathbf{\Phi}_{\text{Picard}}(\mathbf{v}_\alpha))$, as illustrated on the left of Figure 3. In particular, it does not coincide with $\mathbf{v}^*$ from (2.15), which we are however interested in. This issue can be circumvented by choosing the soft thresholds such that $\lim_{k \to \infty} \alpha_k = 0$. An adaptive strategy presented and analyzed in [9] consists in setting

$$ \alpha_{k+1} = \begin{cases} \theta \alpha_k, & \text{if } \|S_{\alpha_k}(\mathbf{\Phi}_{\text{Picard}}(\mathbf{v}^k)) - \mathbf{v}^k\|_J \le c \|\mathbf{\Phi}_{\text{Picard}}(\mathbf{v}^k) - \mathbf{v}^k\|_J, \\ \alpha_k, & \text{otherwise,} \end{cases} $$
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Left: illustration of the fixed-point $v^*$ of $\Phi$, and of the modified version $v_\alpha$ due to the soft thresholding operator $S_\alpha$. Right: convergence of the fixed-point iteration coupled with soft thresholding, $v^k$, when the parameter $\alpha$ is decreased every second iteration.</figcaption></figure>

for some $c, \theta \in (0, 1)$. This amounts to decrease $\alpha$ whenever $v^k$ gets too close to $v_\alpha$, compared to its distance with $v^*$, see also Figure 3 on the right for an illustration. The initial threshold $\alpha_0$ is taken large enough to ensure that $S_{\alpha_0}(\Phi_{\text{Picard}}(v^0)) = S_{\alpha_0}(\Phi_{\text{Picard}}(0)) = 0$.

In this paper, however, we show that it is also possible to use a simpler strategy, consisting in multiplying $\alpha$ by a fixed factor $\theta \in (0, 1)$ at each iteration. Using the same initialization, the thresholds are then given by $\alpha_k = \theta^k \alpha_0$. As we will see in the following, this procedure yields both a convergence result and rank bounds on the iterates of the resulting fixed-point iteration. The pseudo-codes for the Picard iteration with constant decrease of $\alpha$ and SDC with adaptive decrease of $\alpha$ can be found in Algorithms 3.1 and 3.2. The other options, namely Picard iteration with adaptive decrease of $\alpha$ and SDC with constant decrease of $\alpha$, are obtained in a similar way.

<table><caption>**Algorithm 3.1** Gauss-Legendre-Picard iteration with constant decrease of $\alpha$</caption><tr><td>1:</td><td>for n = 1, ..., N do</td><td>▷ Time stepping</td></tr><tr><td>2:</td><td>Adjust the nodes t<sub>1</sub>, ..., t<sub>J</sub> to the current interval [(n − 1)h, nh]</td><td></td></tr><tr><td>3:</td><td>Initialize k ← −1, α<sub>0</sub> ← σ<sub>max</sub>(v<sub>0</sub>), v<sup>0</sup> ← 0, res ← 1</td><td></td></tr><tr><td>4:</td><td>while res &gt; ε<sub>n</sub> do</td><td>▷ Fixed-point iteration</td></tr><tr><td>5:</td><td>Increment k ← k + 1</td><td></td></tr><tr><td>6:</td><td>Compute φ<sup>k</sup> = Φ<sub>Picard</sub>(v<sup>k</sup>) from (2.8)</td><td></td></tr><tr><td>7:</td><td>Update res ← ||φ<sup>k</sup> − v<sup>k</sup>||<sub>J</sub></td><td></td></tr><tr><td>8:</td><td>Compute v<sup>k+1</sup> ← S<sub>α<sub>k</sub></sub>(φ<sup>k</sup>)</td><td>▷ Soft thresholding</td></tr><tr><td>9:</td><td>Update α<sub>k+1</sub> ← θα<sub>k</sub></td><td>▷ Fixed decrease of α</td></tr><tr><td>10:</td><td>end while</td><td></td></tr><tr><td>11:</td><td>v<sub>0</sub> ← v<sub>0</sub> + &Sigma;<sup>J</sup><sub>j=1</sub> ω<sub>j</sub> F<sub>t<sub>j</sub></sub> v<sub>j</sub><sup>k</sup></td><td>▷ Approximation at subinterval boundary</td></tr><tr><td>12:</td><td>v<sub>0</sub> ← R<sub>δ<sub>n</sub></sub>(v<sub>0</sub>)</td><td>▷ Recompression to a tolerance δ<sub>n</sub></td></tr><tr><td>13:</td><td>end for</td><td></td></tr></table>

In practice, we apply a recompression operator $\mathcal{R}_\delta$ of the form (3.3), with $\delta = \delta_{\text{rel}}$ res, after every addition of two elements in $L^2(\Omega_1 \times \Omega_2)$. Here res denotes the current residual error estimator, and $\delta_{\text{rel}}$ is a small parameter ensuring a bound on the relative error. In the case of SDC, we need to use an additional recompression operator with a smaller relative tolerance in order to approximate the residual occurring implicitly in line 7 of the algorithm accurately enough. Here, the recompression tolerance is chosen relative to the current residual norm. See also [20] for a derivation of SDC based on residual and error functions using the spectral integration matrix.
<caption>**Algorithm 3.2 SDC with adaptive decrease of $\alpha$**</caption>

1: **for** $n = 1, \dots, N$ **do**
2: &nbsp;&nbsp;&nbsp;&nbsp;Adjust the nodes $t_1, \dots, t_J$ to the current interval $[(n-1)h, nh]$
3: &nbsp;&nbsp;&nbsp;&nbsp;Initialize $k \leftarrow -1$, $\alpha_0 \leftarrow \sigma_{\max}(v_0)$, $\mathbf{v}^0 \leftarrow \mathbf{0}$, $\text{res} \leftarrow 1$
4: &nbsp;&nbsp;&nbsp;&nbsp;**while** $\text{res} > \varepsilon_n$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Increment $k \leftarrow k + 1$
6: &nbsp;&nbsp;&nbsp;&nbsp;**for** $j = 1, \dots, J$ **do**
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute $\phi_j^k$ from (3.6)
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Threshold $v_j^{k+1} \leftarrow S_{\alpha_k}(\phi_j^k)$
9: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
10: &nbsp;&nbsp;&nbsp;&nbsp;Update res $\leftarrow \|\phi^k - \mathbf{v}^k\|_J$
11: &nbsp;&nbsp;&nbsp;&nbsp;Compute err $\leftarrow \|\mathbf{v}^{k+1} - \mathbf{v}^k\|_J$
12: &nbsp;&nbsp;&nbsp;&nbsp;**if** err $\le c$ res **then**
13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\alpha_{k+1} \leftarrow \theta\alpha_k$
14: &nbsp;&nbsp;&nbsp;&nbsp;**else**
15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\alpha_{k+1} \leftarrow \alpha_k$
16: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
17: &nbsp;&nbsp;&nbsp;&nbsp;**end while**
18: &nbsp;&nbsp;&nbsp;&nbsp;$v_0 \leftarrow v_0 + \sum_{j=1}^J \omega_j \mathcal{F}_{t_j} v_j^k$
19: &nbsp;&nbsp;&nbsp;&nbsp;$v_0 \leftarrow \mathcal{R}_{\delta_n}(v_0)$
20: **end for**

$\triangleright$ Time stepping

$\triangleright$ Fixed point iteration

$\triangleright$ Soft thresholding

$\triangleright$ Adaptive decrease of $\alpha$

$\triangleright$ Approximation at subinterval boundary

$\triangleright$ Recompression to a tolerance $\delta_n$

**3.3. Quasi-optimal ranks.** Given $v \in L^2(\Omega)$ and $\varepsilon > 0$, the low-rank best approximation of $v$ satisfies

$$ (3.7) \qquad \mathcal{R}_{\varepsilon}(v) \in \underset{\|v-w\| \le \varepsilon}{\arg\min} \operatorname{rank}(w). $$

However, the functions we want to approximate are not exactly known, hence we cannot hope to compute this optimal low-rank representation. We thus aim for a relaxed notion of ranks, where we replace the constraint $\|v - w\| \le \varepsilon$ by a penalization of $\|v - w\|$. Soft thresholding is naturally suited for this purpose, as illustrated by the following result.

**Proposition 11.** Let $v \in L^2(\Omega)$ and $\alpha > 0$. It then holds that

$$ \|v - S_{\alpha}v\|^2 = \min_{w \in L^2(\Omega)} \{\|v - w\|^2 + \alpha^2 \operatorname{rank}(w)\}. $$

*Proof.* For any $w \in L^2(\Omega)$, using the Mirsky inequality from Lemma 34, one can estimate

$$ \begin{aligned} \|v - w\|^2 + \alpha^2 \operatorname{rank}(w) &\ge \sum_{k \in \mathbb{N}} |\sigma_k(v) - \sigma_k(w)|^2 + \alpha^2 \mathbb{1}_{\sigma_k(w) \ne 0} \\ &\ge \sum_{k \in \mathbb{N}} \min(\sigma_k(v)^2, \alpha^2) = \|v - S_{\alpha}v\|^2, \end{aligned} $$

where we distinguished the cases $\sigma_k(w) = 0$ and $\sigma_k(w) \ne 0$ to obtain the second line. Equality is attained for $w = \mathcal{H}_{\alpha}(u)$. $\square$

Given $v \in L^2(\Omega)$ and $\varepsilon \in (0, \|v\|)$, there exists $\alpha > 0$ such that $\varepsilon = \|v - S_{\alpha}(v)\|$. In view of Proposition 11 we will use the quantity

$$ (3.8) \qquad \frac{\|v - S_{\alpha}(v)\|^2}{\alpha^2} = \frac{\varepsilon^2}{\alpha^2} $$

as a benchmark for the ranks of an approximation of $v$.

**Remark 12.** When the singular values of $v$ decay regularly, these quantities correspond to quasi-optimal ranks, that is, they are bounded by the optimal ranks up to a constant. Indeed,
it is shown in [9, Prop. 3.6] that

$$ \mathrm{rank}(\mathcal{R}_{\varepsilon}(v)) \leq C\varepsilon^{-\frac{1}{s}} \quad \text{implies} \quad \frac{\|v - \mathcal{S}_{\alpha}(v)\|^2}{\alpha^2} \leq C_s\varepsilon^{-\frac{1}{s}}, $$

where $C_s > 0$ depends on $C > 0$ and $s > 0$. Moreover

$$ \mathrm{rank}(\mathcal{R}_{\varepsilon}(v)) \leq c(1 + |\log \varepsilon|)^{\beta} \quad \text{implies} \quad \frac{\|v - \mathcal{S}_{\alpha}(v)\|^2}{\alpha^2} \leq C_{c,\beta}(1 + |\log \varepsilon|)^{\beta}, $$

with $C_{c,\beta} > 0$ depending on $c$ and $\beta$. Thus in both cases, (3.8) yields quasi-optimal ranks.

### 4. ANALYSIS OF THE PICARD ITERATION

In this section, we analyze the accuracy of the proposed Gauss-Legendre-Picard iteration scheme and provide bounds on the ranks of the approximations at both intermediate time steps $t_1, t_2, \dots, t_J \in ((n-1)h, nh)$ and subintervals boundaries $nh$ for $1 \le n \le N$. The resulting estimates can then be combined together to a global quasi-optimality result.

For notational convenience, we introduce the following simplified notation: since we restrict our considerations to the case of the Gauss-Legendre-Picard iteration, we denote by $\Phi$ and $\mathbf{\Phi}$ the equivalent fixed-point iterations $\Phi_{\mathrm{Picard}}$ and $\mathbf{\Phi}_{\mathrm{Picard}}$ but formulated in untwisted variables, that is,

$$ \mathbf{\Phi}\mathbf{u}(t) = e^{it\Delta}u_0 - i \int_0^t \sum_{m=1}^J e^{i(t-t_m)\Delta}V u_m \ell_m(s) ds $$

and

$$ (4.1) \qquad (\mathbf{\Phi}\mathbf{u})_j = \mathbf{\Phi}\mathbf{u}(t_j) = e^{it_j\Delta}u_0 - i \sum_{m=1}^J \omega_{j,m} e^{i(t_j-t_m)\Delta}V u_m $$

for $\mathbf{u} \in L^2(\Omega)^J$. Note that the fixed-point iterations on the original and untwisted variable have both the same contraction constant $\rho$. We also omit the parentheses when applying the (nonlinear) operators $\Phi$, $S_\alpha$, $H_\alpha$ and $R_\delta$.

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Visualization of the decomposition of the local error in Proposition 13.</figcaption></figure>

Let $u^*(t)$ be the exact solution of the collocation formulation (2.16) on $[0, h]$ with initial data $u_0^*$, that is,

$$ u^*(t) = e^{it\Delta}u_0^* - i \int_0^t \sum_{m=1}^J e^{i(t-t_m)\Delta}V u_m^* \ell_m(s) ds, \quad \text{for } t \in [0, h], $$

with $u_m^* = u^*(t_m)$, and let $u^\diamond(t)$ be the solution of the collocation formulation on $[0, h]$ with initial data $u_0^\diamond$. We further write $u(t)$ for the exact solution of the Schrödinger equation with initial data $u_0^\diamond$ and $\tilde{u}(t)$ for the approximate solution produced by the scheme initialized at $u_0^*$, see Figure 4.
The following analysis requires some smoothness assumptions on the potential $V$. We therefore assume that the potential is uniformly bounded, that is, $C_V < \infty$. Moreover, we assume that the quantities

$$ \| (\partial_t - i\Delta + iV)^{2J+1} u^{\diamond} \|_{W(0,h)} \quad \text{and} \quad \| \partial_t^J (e^{-it\Delta}Vu) \|_{W(0,h)} $$

can also be bounded from above. In order to derive rank bounds for the iterates of the considered algorithm, we also assume that the potential function has finite rank

$$ V = \sum_{\ell=1}^{\text{rank}(V)} V_1^{(\ell)} \otimes V_2^{(\ell)}, \quad \text{for} \quad V_i^{(\ell)} \in L^2(\Omega_i), \quad i=1,2, $$

with $\text{rank}(V) \in \mathbb{N}$. We emphasize that this is the rank of the function $(x_1, x_2) \mapsto V(x_1, x_2)$, which differs from the (potentially infinite) rank of the multiplication operator $u \mapsto Vu$.

**4.1. Local accuracy bound.** Since we aim at approximating the solution of the Schrödinger equation in the untwisted variable, we formulate the accuracy results accordingly in terms of the original variable. However, we still make use of the iteration (2.2) in the proofs. Note that the accuracy results in this section hold independently of the particular strategy for decreasing the soft thresholding parameter $\alpha$.

The following theorem provides a local accuracy bound for the analyzed scheme. Note that here, if iterated over the subintervals, the bound will depend only linearly on the number of needed subintervals, which would correspond to a linear error accumulation and therefore achieve the optimal rate one can hope for in this setting.

**Proposition 13.** Let $\delta$ be the recompression parameter for the approximation at the subinterval boundary and $\varepsilon$ the tolerance of the Gauss-Legendre-Picard iteration. The scheme then achieves an accuracy of

$$ \|u(h) - \tilde{u}(h)\| \le \delta + \frac{\rho \varepsilon}{1 - \rho} + \|u_0^{\diamond} - u_0^*\| + \kappa_{2J} h^{2J+1}, $$

where we recall that $\rho = \Lambda_J C_V h$ denotes the contraction constant and

$$ (4.2) \qquad \kappa_{2J} = \frac{(J!)^4}{(2J+1)!(2J)!^2} \|(\partial_t - i\Delta + iV)^{2J+1} u^{\diamond}\|_{W(0,T)}. $$

*Remark 14.* By developing the term $(\partial_t - i\Delta + iV)^{2J+1}u^{\diamond}$, and using commutator identities, one can see that the constant $\kappa_{2J}$ involves $4J$ derivatives of $V$, but only $2J$ derivatives of $u^{\diamond}$, which translate into $2J$ derivatives on $u$. This gain in regularity requirements based on commutator estimates is already observed in [1], where it is also shown to hold when the equation contains the nonlinearity $|u|^2u$.

The proof of the above proposition will make use of two auxiliary lemmas. For simplicity of notation, we formulate them for the first interval $[0, h]$. However, they also analogously hold on each of the other subintervals $[nh, (n+1)h]$. First of all, we note that the proposed scheme is norm preserving.

**Lemma 15 (Isometry preservation).** It holds that $\|u^*(h) - u^{\diamond}(h)\| = \|u_0^* - u_0^{\diamond}\|$.

Note that this result in fact only holds for the fixed-point solutions, and not for the intermediate iterations.

*Proof.* According to (2.16), in the twisted variables, $p(t) = e^{-it\Delta}(u^*(t) - u^{\diamond}(t))$ is a polynomial of degree $J$, which satisfies the collocation identity

$$ \partial_t p(t_j) = \mathcal{F}_{t_j}(p(t_j)) = -ie^{-it_j\Delta}V e^{it_j\Delta}p(t_j), \quad 1 \le j \le J. $$
As $\langle \partial_t p(t), p(t) \rangle$ is a polynomial of degree $2J-1$, its Gauss quadrature is exact, hence

$$ 
\begin{aligned}
\|u^*(h) - u^\diamond(h)\|^2 - \|u_0^* - u_0^\diamond\|^2 &= \|p(h)\|^2 - \|p(0)\|^2 = \int_0^h \frac{d}{dt} \|p(s)\|^2 ds \\
&= 2 \operatorname{Re} \int_0^h \langle \partial_t p(s), p(s) \rangle ds = 2 \operatorname{Re} \sum_{j=1}^J \omega_j \langle \partial_t p(t_j), p(t_j) \rangle \\
&= 2 \operatorname{Re} \left( -i \sum_{j=1}^J \|e^{it_j \Delta} p(t_j)\|^2 \right) = 0. \quad \square
\end{aligned}
$$

Next, we prove that the fixed-point is able to reach a local accuracy of order $2J+1$ at the subinterval boundary, which can then be translated into a global error bound of order $2J$.

**Lemma 16 (Local error at subinterval boundary).** It holds that

$$ \|u^\diamond(h) - u(h)\| \le \kappa_{2J} h^{2J+1}, $$

where the constant $\kappa_{2J}$ is defined as in (4.2).

*Proof.* Let $e^{it(\Delta-V)}$ be the flow map of the Schrödinger equation, which is well defined on $L^2(\Omega)$. For $t \in [0, h]$, setting

$$ f(t) = e^{i(h-t)(\Delta-V)} u^\diamond(t) $$

yields

$$ (4.3) \qquad u^\diamond(h) - u(h) = f(h) - f(0) = \int_0^h \partial_t f(t) dt $$

where

$$ \partial_t f(t) = e^{i(h-t)(\Delta-V)} (\partial_t - i\Delta + iV) u^\diamond(t). $$

Observe in particular that $\partial_t f(t_j) = 0$ since $u^*(t)$ is the solution of a collocation method. We can thus bound (4.3) by

$$ \|u^\diamond(h) - u(h)\| = \left\| \int_0^h \partial_t f(t) dt - \sum_{j=1}^J \omega_j \partial_t f(t_j) \right\| \le \frac{(J!)^4}{(2J+1)!(2J)!^2} \|\partial_t^{2J+1} f\|_{W(0,h)} h^{2J+1}. $$

The last bound is a standard accuracy result for Gauss-Legendre quadrature, and can be found for example in [41]. To make sure that it also works in the Hilbert-space-valued case, we provide a proof with explicit constants in Lemma 32 of the appendix.

The statement of the lemma then follows from

$$ 
\begin{aligned}
\|\partial_t^{2J+1} f\|_{W(0,h)} &= \max_{t \in [0,h]} \|e^{i(h-t)(\Delta-V)} (\partial_t - i\Delta + iV)^{2J+1} u^\diamond(t)\| \\
&= \|(\partial_t - i\Delta + iV)^{2J+1} u^\diamond\|_{W(0,h)}.
\end{aligned}
\quad \square
$$

Inside the interval $[0, h]$, we obtain a local accuracy bound of order $J+1$.

**Lemma 17 (Local error bound inside subintervals).** It holds that

$$ \|u^*(t) - u(t)\| \le \frac{\kappa_J}{1-\rho} h^{J+1} + \frac{1}{1-\rho} \|u_0^* - u_0^\diamond\|, $$

for all $t \in [0, h]$, where

$$ \kappa_J = \frac{J!}{\sqrt{2J+1}(2J)!} \|\partial_t^J (e^{-it\Delta} V u)\|_{W(0,h)}. $$
Proof. For the twisted variables $v^*(t) = e^{-it\Delta}u^*(t)$ and $v(t) = e^{-it\Delta}u(t)$ it holds that

$$ \|u^*(t) - u(t)\| \le \|u_0^* - u_0^\diamond\| + \left\| \int_0^t \sum_{m=1}^J \mathcal{F}_{t_m}(v_m^* - v(t_m)) \ell_m(s) ds \right\| \\ + \left\| \int_0^t \left( \sum_{m=1}^J \mathcal{F}_{t_m} v(t_m) \ell_m(s) - \mathcal{F}_s v(s) \right) ds \right\|. $$

By Proposition 8, the second term is bounded by $\rho\|v^* - v\|_{W(0,h)} = \rho\|u^* - u\|_{W(0,h)}$.

Applying Lemma 33 to $g(s) = \mathcal{F}_s v(s) = -ie^{is\Delta}V u(s)$, we bound the third term by $\kappa_J h^{J+1}$.

We conclude by taking the supremum over all $t \in [0, h]$, and passing the second term to the left hand side. $\square$

We are finally in the position to prove our local accuracy bound.

Proof of Proposition 13. Let $K$ denote the value of $k$ at the end of Algorithm 3.1, for which we know that

$$ \|\mathbf{\Phi}\mathbf{u}^K - \mathbf{u}^K\|_J \le \varepsilon. $$

We can then decompose the error according to

$$ \|u(h) - \tilde{u}(h)\| \le \|u(h) - u^{\diamond}(h)\| + \|u^{\diamond}(h) - u^{*}(h)\| + \|u^{*}(h) - \mathbf{\Phi}\mathbf{u}^{K}(h)\| + \|\mathbf{\Phi}\mathbf{u}^{K}(h) - \tilde{u}(h)\|. $$

The first term can be bounded by $\kappa_2 J h^{2J+1}$ using Lemma 15, and the second one is equal to $\|u_0^\diamond - u_0^*\|$ according to Lemma 16. For the third term, Proposition 8 yields

$$ \|u^* - \mathbf{\Phi}\mathbf{u}^K\|_{W(0,h)} \le \rho\|\mathbf{u}^* - \mathbf{u}^K\|_J \le \rho(\|u^* - \mathbf{\Phi}\mathbf{u}^K\|_{W(0,h)} + \|\mathbf{u}^K - \mathbf{\Phi}\mathbf{u}^K\|_J), $$

and thus

$$ \|u^*(h) - \mathbf{\Phi}\mathbf{u}^K(h)\| \le \|u^* - \mathbf{\Phi}\mathbf{u}^K\|_{W(0,h)} \le \frac{\rho}{1-\rho}\|\mathbf{u}^K - \mathbf{\Phi}\mathbf{u}^K\|_J \le \frac{\rho\varepsilon}{1-\rho}. $$

Finally, by the characterisation (3.7) of $\mathcal{R}_\delta$, the last term is bounded by the recompression tolerance $\delta$. $\square$

4.2. Local rank bounds. Next, we proceed to the derivation of local rank bounds for the iterates of Algorithm 3.1 with constant decrease of $\alpha_k$, motivated by the notion of quasi-optimality discussed in Section 3.3. As introduced in the previous proof, let $K$ denote the value of $k$ at the end of Algorithm 3.1. Additionally, we set

$$ \mathbf{u}^* = (u^*(t_m))_{1 \le m \le J} \quad \text{and} \quad \mathbf{u} = (u(t_m))_{1 \le m \le J}. $$

The following lemma shows that the Gauss-Legendre-Picard iterations produce quasi-optimal low-rank approximations. It simplifies the analysis from [9] by allowing a constant decrease of $\alpha$.

**Lemma 18.** Starting from the initial data $u_0^*$, the approximations generated by Algorithm 3.1 in untwisted variables satisfy

$$ \|\mathbf{u}^{k+1} - \mathbf{u}^*\|_J \le \frac{2}{1-\rho}\|\mathbf{u}^* - \mathbf{S}_{\alpha_k}\mathbf{u}^*\|_J, \quad \text{for all } 0 \le k \le K. $$

Proof. Since $\mathbf{S}_{\alpha_0} \mathbf{\Phi}$ is contractive with fixed point at 0, it holds that

$$ \|\mathbf{u}^*\|_J - \|\mathbf{u}^* - \mathbf{S}_{\alpha_0} \mathbf{u}^*\|_J \le \|\mathbf{S}_{\alpha_0} \mathbf{u}^*\|_J = \|\mathbf{S}_{\alpha_0} \mathbf{\Phi} \mathbf{u}^*\|_J \le \rho \|\mathbf{u}^*\|_J $$

and thus

$$ \|\mathbf{u}^0 - \mathbf{u}^*\|_J = \|\mathbf{u}^*\|_J \le \frac{1}{1-\rho} \|\mathbf{u}^* - \mathbf{S}_{\alpha_0} \mathbf{u}^*\|_J. $$

We then proceed by induction on $k$. For $k \ge 0$, we have

$$ \|\mathbf{u}^{k+1} - \mathbf{u}^*\|_J \le \|\mathbf{S}_{\alpha_k} \mathbf{\Phi} \mathbf{u}^k - \mathbf{S}_{\alpha_k} \mathbf{\Phi} \mathbf{u}^*\|_J + \|\mathbf{u}^* - \mathbf{S}_{\alpha_k} \mathbf{\Phi} \mathbf{u}^*\|_J \le \rho \|\mathbf{u}^k - \mathbf{u}^*\|_J + \|\mathbf{u}^* - \mathbf{S}_{\alpha_k} \mathbf{u}^*\|_J. $$

By Lemma 10, for $0 \le l \le k$, it holds that

$$ \theta^l \|\mathbf{u}^* - \mathbf{S}_{\alpha_{k-l}} \mathbf{u}^*\|_J \le \|\mathbf{u}^* - \mathbf{S}_{\alpha_k} \mathbf{u}^*\|_J. $$
By induction, we then obtain

$$
\begin{aligned}
\|\mathbf{u}^{k+1} - \mathbf{u}^*\|_J &\le \sum_{\ell=0}^k \rho^\ell \|\mathbf{u}^* - \mathbf{S}_{\alpha_{k-\ell}} \mathbf{u}^*\|_J + \rho^{k+1} \|\mathbf{u}^0 - \mathbf{u}^*\|_J \\
&\le \|\mathbf{u}^* - \mathbf{S}_{\alpha_k} \mathbf{u}^*\|_J \left( \sum_{\ell=0}^k \frac{\rho^\ell}{\theta^\ell} + \frac{\rho^{k+1}}{\theta^k} \frac{1}{1-\rho} \right).
\end{aligned}
$$

Setting $\theta = \sqrt{\rho}$, the last constant can be bounded by

$$ \sum_{\ell=0}^{k-1} \frac{\rho^\ell}{\theta^\ell} + \frac{\rho^k}{\theta^k} \left(1 + \frac{\rho}{1-\rho}\right) = \frac{1-\sqrt{\rho^k}}{1-\sqrt{\rho}} + \frac{\sqrt{\rho^k}}{1-\rho} \le \frac{2}{1-\rho}. \quad \square $$

Remark 19. It is in fact possible to take an arbitrary $\theta \in [\rho, 1)$ in the proof. For $\theta = \rho$, we achieve the optimal rate of decay of $(\alpha_k)_{k \ge 0}$, and the constant in the error bound is at most $k + \frac{1}{1-\rho}$. For $\theta = \sqrt{\rho}$, the algorithm needs twice as many iterations in order to reach a desired value of the threshold, but one obtains a constant independent of the iteration number $k$.

We also need the following lemma, which was originally proposed in [9, Lemma 4.3] and allows to compare the ranks of the scheme to the ranks of the exact solution.

**Lemma 20.** For any $u, v \in L^2(\Omega)$ and $\alpha, \beta > 0$, it holds that

$$ \mathrm{rank}(\mathcal{S}_{\alpha+\beta}(u)) \le \mathrm{rank}(\mathcal{S}_{\alpha}v) + \frac{1}{\beta^2} \|u-v\|^2. $$

With this, we obtain the quasi-optimality of the ranks of the iterates $\mathbf{u}^k$.

**Lemma 21.** For $\varepsilon \ge 4 \frac{1+\rho}{1-\rho} \|\mathbf{u} - \mathbf{u}^*\|_J$, let $\bar{\alpha}$ be the value of the soft threshold such that

$$ \|u - S_{\bar{\alpha}}u\|_{W(0,h)} = \frac{1 - \rho}{1 + \rho} \frac{\varepsilon}{4}. $$

Then for $0 \le k \le K$,

$$ (4.4) \quad \mathrm{rank}(\mathbf{u}^k) := \max_{1 \le j \le J} \mathrm{rank}(u_j^k) \le \bar{r} := \frac{\varepsilon^2}{\rho \bar{\alpha}^2} = \frac{16 (1+\rho)^2}{\rho (1-\rho)^2} \frac{\|u - S_{\bar{\alpha}} u\|_{W(0,h)}^2}{\bar{\alpha}^2}. $$

*Proof.* The value of $\bar{\alpha}(t)$ is well-defined since the soft thresholding error is continuous and monotonic. The ranks of the first two iterates $\mathbf{u}^0$ and $\mathbf{u}^1$ are simply zero. For the remaining iterations, we let $0 \le k \le K-2$ and apply Lemma 20 to obtain

$$ \mathrm{rank}(\mathbf{u}^{k+2}) = \mathrm{rank}(\mathcal{S}_{\alpha_{k+1}} \mathbf{\Phi} \mathbf{u}^{k+1}) \le \mathrm{rank}(\mathcal{S}_{\alpha_{k+1}/2} \mathbf{u}) + \frac{4}{\alpha_{k+1}^2} \|\mathbf{\Phi} \mathbf{u}^{k+1} - \mathbf{u}\|_J^2. $$

The first term is controlled by

$$ (4.5) \quad \mathrm{rank}(\mathcal{S}_{\alpha_{k+1}/2} \mathbf{u}) \le \frac{4}{\alpha_{k+1}^2} \|\mathbf{u} - \mathcal{S}_{\alpha_{k+1}/2} \mathbf{u}\|_J^2 \le \frac{4}{\alpha_{k+1}^2} \|\mathbf{u} - \mathcal{S}_{\alpha_k} \mathbf{u}\|_J^2. $$

For the second term, notice that

$$ \|\mathbf{\Phi}\mathbf{u}^{k+1} - \mathbf{u}\|_J \le \|\mathbf{\Phi}\mathbf{u}^{k+1} - \mathbf{u}^*\|_J + \|\mathbf{u}^* - \mathbf{u}\|_J \le \rho\|\mathbf{u}^{k+1} - \mathbf{u}^*\|_J + \|\mathbf{u}^* - \mathbf{u}\|_J. $$

Using Lemma 18 and the non-expansivity of $\mathcal{I} - \mathcal{S}_{\alpha_k}$, we obtain that

$$ \|\mathbf{u}^{k+1} - \mathbf{u}^*\|_J \le \frac{2}{1-\rho} \|\mathbf{u}^* - \mathcal{S}_{\alpha_k} \mathbf{u}^*\|_J \le \frac{2}{1-\rho} (\|\mathbf{u} - \mathcal{S}_{\alpha_k} \mathbf{u}\|_J + \|\mathbf{u}^* - \mathbf{u}\|_J), $$

which yields

$$ \begin{aligned} \|\mathbf{\Phi}\mathbf{u}^{k+1} - \mathbf{u}\|_J^2 &\le \left( \frac{2\rho}{1-\rho} \|\mathbf{u} - \mathcal{S}_{\alpha_k} \mathbf{u}\|_J + \frac{1+\rho}{1-\rho} \|\mathbf{u}^* - \mathbf{u}\|_J \right)^2 \\ &\le \frac{8\rho^2}{(1-\rho)^2} \|\mathbf{u} - \mathcal{S}_{\alpha_k} \mathbf{u}\|_J^2 + 2 \frac{(1+\rho)^2}{(1-\rho)^2} \|\mathbf{u}^* - \mathbf{u}\|_J^2, \end{aligned} \quad (4.6) $$
where we used the inequality $(a+b)^2 \le 2a^2 + 2b^2$ to obtain the second line.

Combining equations (4.5) and (4.6), we get

$$ (4.7) \qquad \begin{aligned} \mathrm{rank}(\mathbf{u}^{k+2}) &\le \frac{4}{\alpha_{k+1}^2} \left( \frac{(1-\rho)^2 + 8\rho^2}{(1-\rho)^2} \|\mathbf{u} - \mathbf{S}_{\alpha_k} \mathbf{u}\|_J^2 + 2 \frac{(1+\rho)^2}{(1-\rho)^2} \|\mathbf{u}^* - \mathbf{u}\|_J^2 \right) \\ &\le \frac{8}{\rho} \frac{(1+\rho)^2}{(1-\rho)^2} \frac{\|\mathbf{u} - \mathbf{S}_{\alpha_k} \mathbf{u}\|_J^2 + \|\mathbf{u}^* - \mathbf{u}\|_J^2}{\alpha_k^2}, \end{aligned} $$

where the second line comes from the fact that $\alpha_{k+1}^2 = \rho\alpha_k^2$, together with the inequality $(1-\rho)^2 + 8\rho^2 \le 2(1+\rho)^2$, which is true since $\rho^2 \le \rho \le 1$.

It remains to replace $\alpha_k$ by $\bar{\alpha}$. As $k+1 < K$, the stopping criterion implies that

$$ \frac{\varepsilon}{1+\rho} \le \|\mathbf{u}^{k+1} - \mathbf{u}^*\|_J \le \frac{2}{1-\rho} \|\mathbf{u}^* - \mathbf{S}_{\alpha_k} \mathbf{u}^*\|_J, $$

where we again make use of Lemma 18. Using the condition on $\varepsilon$ and the non-expansivity of $\mathcal{I} - S_{\alpha_k}$, we obtain that

$$ \begin{aligned} \|u - S_{\alpha_k} u\|_{W(0,h)} &\ge \|\mathbf{u} - \mathbf{S}_{\alpha_k} \mathbf{u}\|_J \ge \|\mathbf{u}^* - \mathbf{S}_{\alpha_k} \mathbf{u}^*\|_J - \|\mathbf{u} - \mathbf{u}^*\|_J \\ &\ge \frac{1-\rho}{1+\rho} \frac{\varepsilon}{2} - \|\mathbf{u} - \mathbf{u}^*\|_J \ge \frac{1-\rho}{1+\rho} \frac{\varepsilon}{4} = \|u - S_{\bar{\alpha}} u\|_{W(0,h)}. \end{aligned} $$

This shows that $\alpha_k \ge \bar{\alpha}$, thus

$$ \frac{\|\mathbf{u} - \mathbf{S}_{\alpha_k} \mathbf{u}\|_J^2}{\alpha_k^2} \le \frac{\|\mathbf{u} - \mathbf{S}_{\bar{\alpha}} \mathbf{u}\|_J^2}{\bar{\alpha}^2} \quad \text{and} \quad \frac{\|\mathbf{u}^* - \mathbf{u}\|_J^2}{\alpha_k^2} \le \frac{\|\mathbf{u}^* - \mathbf{u}\|_J^2}{\bar{\alpha}^2}. $$

As a consequence, we can replace $\alpha_k$ by $\bar{\alpha}$ in equation (4.7), and conclude by observing that

$$ \|\mathbf{u}^* - \mathbf{u}\|_J \le \frac{1-\rho}{1+\rho} \frac{\varepsilon}{4} = \|u - S_{\bar{\alpha}} u\|_{W(0,h)} \quad \text{and} \quad \|\mathbf{u} - \mathbf{S}_{\bar{\alpha}} \mathbf{u}\|_J \le \|u - S_{\bar{\alpha}} u\|_{W(0,h)}. \quad \square $$

As discussed in Section 3.3, the quantity $\|u - S_{\bar{\alpha}} u\|_{W(0,h)}^2 / \bar{\alpha}^2$ can indeed be considered to describe a quasi-optimal rank for the corresponding error. Note that the factor $1/\rho$ in the bound (4.4) can be improved by choosing $\theta$ closer to 1.

We proceed in the derivation of a local rank bound by considering how the application of the fixed-point operator modifies the rank of the current iterate.

**Lemma 22.** For $0 \le k \le K$ the ranks of the fixed-point iteration $\mathbf{\Phi}$ applied to the current iterate $\mathbf{u}^k$ can be bounded in terms of

$$ \mathrm{rank}(\mathbf{\Phi}\mathbf{u}^k) \le \mathrm{rank}(u_0^*) + J \mathrm{rank}(V) \mathrm{rank}(\mathbf{u}^k). $$

*Proof.* It suffices to use the estimates $\mathrm{rank}(f+g) \le \mathrm{rank}(f) + \mathrm{rank}(g)$ and $\mathrm{rank}(fg) \le \mathrm{rank}(f)\mathrm{rank}(g)$ together with the definition of the fixed-point formulation. $\square$

Finally, we need an estimate for the ranks after the final recompression step carried out by applying $\mathcal{R}_\delta$ as defined in (3.3). This estimate is illustrated in Figure 5.

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Visualization of Lemma 23: the squares of the red singular values sum up to at least $\delta^2$, and the blue ones to at most $\|v\|^2 - \delta^2$, thus $r_\delta$ cannot be too close to $\mathrm{rank}(v)$.</figcaption></figure>
**Lemma 23.** For any $v \in L^2(\Omega)$ of finite rank, it holds that

$$r_{\delta} = \mathrm{rank}(\mathcal{R}_{\delta}v) \leq 1 + \frac{\|v\|^2 - \delta^2}{\|v\|^2} \mathrm{rank}(v).$$

*Proof.* Recall the definitions from (3.3):

$$\mathcal{R}_{\delta}(v) = \sum_{1 \leq k \leq r_{\delta}} \sigma_k v_k^{(1)} \otimes v_k^{(2)}, \quad \text{where} \quad r_{\delta} = \min\{r \in \mathbb{N}_0: \sum_{k>r} \sigma_k^2 \leq \delta^2\}.$$

By the optimality of $r_{\delta}$, we know that

$$\delta^2 < \sum_{k \geq r_{\delta}} \sigma_k^2, \quad \text{and} \quad \|v\|^2 = \sum_{k=1}^{\mathrm{rank}(v)} \sigma_k^2.$$

Therefore, as the singular values are non-increasing,

$$\frac{\|v\|^2}{\|v\|^2 - \delta^2} < \frac{\|v\|^2}{\|v\|^2 - \sum_{k \geq r_{\delta}} \sigma_k^2} = 1 + \frac{\sum_{k \geq r_{\delta}} \sigma_k^2}{\sum_{k < r_{\delta}} \sigma_k^2} \leq 1 + \frac{(\mathrm{rank}(v) - r_{\delta} + 1)\sigma_{r_{\delta}}^2}{(r_{\delta} - 1)\sigma_{r_{\delta}}^2} = \frac{\mathrm{rank}(v)}{r_{\delta} - 1}. \quad \square$$

In order to apply Lemma 23, we need $v$ to be of bounded rank and have a small $L^2(\Omega)$ norm. Instead of the natural choice $v = \Phi\mathbf{u}^K(t)$, we will take $v = \Phi\mathbf{u}^K(h) - S_{\tilde{\alpha}}\Phi\mathbf{u}^K(h)$, which has smaller norm, but the same singular values below the threshold $\tilde{\alpha}$.

**Lemma 24.** For $\varepsilon \geq 4\frac{1+\rho}{1-\rho}\|u - u^*\|_{W(0,h)}$, the rank of $\tilde{u}(h) = \mathcal{R}_{\delta}\Phi\mathbf{u}^K(h)$ can be bounded by

$$\mathrm{rank}(\tilde{u}(h)) \leq \max\left\{\frac{\bar{r}\rho}{2(1-\rho)^2}, 1 + \frac{\varepsilon^2 - \delta^2(1-\rho)^2}{\varepsilon^2} \mathrm{rank}(\Phi\mathbf{u}^K(h))\right\},$$

where $\bar{r}$ is defined as in (4.4).

*Proof.* Let $\tilde{\alpha}$ be the threshold value such that

$$\|\Phi\mathbf{u}^K(h) - S_{\tilde{\alpha}}\Phi\mathbf{u}^K(h)\| = \frac{\varepsilon}{1-\rho}.$$

As

$$\|\Phi\mathbf{u}^K(h) - u(h)\| \leq \frac{\rho\varepsilon}{1-\rho} + \|u(h) - u^*(h)\| \leq \varepsilon \left( \frac{\rho}{1-\rho} + \frac{1-\rho}{4(1+\rho)} \right),$$

the non-expansiveness of $\mathcal{I} - S_{\tilde{\alpha}}$ ensures that

$$\|u(h) - S_{\tilde{\alpha}}u(h)\| \geq \varepsilon \left( \frac{1}{1-\rho} - \frac{\rho}{1-\rho} - \frac{1-\rho}{4(1+\rho)} \right) \geq \frac{3}{4} \frac{1-\rho}{1+\rho} \varepsilon \geq 3 \|u(h) - S_{\tilde{\alpha}}u(h)\|.$$

This implies that $\tilde{\alpha} \geq 3\tilde{\alpha}$ and thus

$$\begin{aligned} \mathrm{rank}(S_{\tilde{\alpha}}\Phi\mathbf{u}^K(h)) &\leq \mathrm{rank}(S_{3\tilde{\alpha}}\Phi\mathbf{u}^K(h)) \\ &\leq \mathrm{rank}(S_{\tilde{\alpha}}u(h)) + \frac{\|\Phi\mathbf{u}^K(h) - u(h)\|^2}{4\tilde{\alpha}^2} \\ &\leq \frac{\|u - S_{\tilde{\alpha}}u\|_{W(0,h)}^2}{\tilde{\alpha}^2} + \frac{\varepsilon^2}{4\tilde{\alpha}^2(1-\rho)^2} \leq \frac{\bar{r}\rho}{2(1-\rho)^2}. \end{aligned}$$

Now, consider the recompressed version $\tilde{u}(h) = \mathcal{R}_{\delta}(\Phi\mathbf{u}^K(h))$. If $\sigma_{\min}(\tilde{u}(h)) > \tilde{\alpha}$ (which could happen even when $\delta < \varepsilon$, since hard thresholding makes a smaller error than soft thresholding), then

$$\mathrm{rank}(\tilde{u}(h)) \leq \mathrm{rank}(S_{\tilde{\alpha}}\Phi\mathbf{u}^K) \leq \frac{\bar{r}\rho}{2(1-\rho)^2}.$$

Otherwise, for $v = \Phi\mathbf{u}^K(h) - S_{\tilde{\alpha}}\Phi\mathbf{u}^K(h)$, it holds that

$$\mathrm{rank}(\tilde{u}(h)) = \mathrm{rank}(\mathcal{R}_{\delta}(\Phi\mathbf{u}^K(h))) = \mathrm{rank}(\mathcal{R}_{\delta}(v)),$$

and we conclude by applying Lemma 23 to $v$:

$$\mathrm{rank}(\mathcal{R}_{\delta}(v)) \leq 1 + \frac{\varepsilon^2/(1-\rho)^2 - \delta^2}{\varepsilon^2/(1-\rho)^2} \mathrm{rank}(v) = 1 + \frac{\varepsilon^2 - \delta^2(1-\rho)^2}{\varepsilon^2} \mathrm{rank}(\Phi\mathbf{u}^K(h)). \quad \square$$
We finish this section by combining the previous lemmas.

**Proposition 25.** Assume that $\rho \le 1/2$, and let

$$ \varepsilon \ge 4 \frac{1 + \rho}{1 - \rho} \|u - u^*\|_{W(0,h)} \quad \text{and} \quad \delta = \frac{\rho \varepsilon}{1 - \rho}. $$

If the initial rank can be bounded by

$$ \mathrm{rank}(u_0^*) \le \frac{1 + J \mathrm{rank}(V)\bar{r}}{\rho^2}, $$

then the rank of the approximation at the right subinterval boundary can be estimated by

$$ \mathrm{rank}(\tilde{u}(h)) \le \frac{1 + J \mathrm{rank}(V)\bar{r}}{\rho^2}, $$

and all other intermediate ranks, that is, the ranks of the iterates $\mathbf{u}^k$ and $\mathbf{\Phi}\mathbf{u}^k$ produced by the algorithm for $0 \le k \le K$, are bounded by $\frac{1+\rho^2}{\rho^2}(1 + J \mathrm{rank}(V)\bar{r})$.

*Proof.* Combining Lemmas 24, 22 and 21, we either have $\mathrm{rank}(\tilde{u}(h)) \le \frac{\bar{r}\rho}{2(1-\rho)^2} \le \bar{r}$, or

$$ \begin{aligned} \mathrm{rank}(\tilde{u}(h)) &\le 1 + \frac{\varepsilon^2 - (1-\rho)^2\delta^2}{\varepsilon^2} \left( \frac{1 + J \mathrm{rank}(V)\bar{r}}{\rho^2} + J \mathrm{rank}(V)\bar{r} \right) \\ &\le 1 + \frac{1-\rho^2}{\rho^2} + (1-\rho^2)J \mathrm{rank}(V)\bar{r} \left( \frac{1}{\rho^2} + 1 \right) \\ &= \frac{1}{\rho^2} + J \mathrm{rank}(V)\bar{r} \left( \frac{1}{\rho^2} - \rho^2 \right). \end{aligned} $$

The bound on intermediate ranks follows by another application of Lemmas 22 and 21. $\square$

*Remark 26.* The assumption $\rho \le 1/2$ is easily satisfied, by taking a time step $h$ twice smaller than what is already needed for the contractivity in Proposition 8. It is included only to avoid further technicalities, by allowing to bound the first option in the maximum of Lemma 24 by $\bar{r}$. If it was removed, the bounds on the ranks in Proposition 25 would have to be replaced by

$$ \max \left\{ \frac{\bar{r}\rho}{2(1-\rho)^2}, \frac{1 + J \mathrm{rank}(V)\bar{r}}{\rho^2} \right\}. $$

**4.3. Global bounds.** We have now all necessary tools to prove the main theorem. Assume that we know some upper bounds on the constants $\kappa_J$ and $\kappa_{2J}$, and denote by $\tilde{u}_n$ the approximation generated by the scheme at time $nh$, $1 \le n \le N = T/h$. Let $\eta$ denote the target accuracy, then our target for the soft threshold is

$$ \alpha_n = \inf_{(n-1)h \le t \le nh} \min_{\|u(t)-v\| \le \eta} \sigma_{\max}(v) = \max\{\alpha > 0: \forall t \in [(n-1)h, nh], \|u(t) - S_{\alpha}u(t)\| \le \eta\}, $$

since the minimizer is of the form $v = S_{\alpha(t)}u(t)$ for some $\alpha(t) > 0$. Our objective for the ranks is then given by

$$ r_n = \frac{\eta^2}{\alpha_n^2} = \frac{\|u - S_{\alpha_n}u\|_{W((n-1)h, nh)}^2}{\alpha_n^2}. $$

We further define

$$ \xi_n = (\eta + \kappa_J h^{J+1} + n\kappa_{2J} h^{2J+1}) \exp \left( 8\rho n \frac{1+\rho}{(1-\rho)^3} \right) $$

and make the choice

$$ \varepsilon_n = 4\xi_{n-1} \frac{1+\rho}{(1-\rho)^2} \quad \text{and} \quad \delta_n = \frac{\rho\varepsilon_n}{1-\rho}. $$
**Theorem 27.** The proposed scheme achieves a global error bounded by

$$ \sup_{1 \le n \le N} \|\tilde{u}_n - u(nh)\| \le (\eta + \kappa_J h^{J+1} + \kappa_{2J} h^{2J} T) \exp \left( 8\Lambda_J C_V \frac{1+\rho}{(1-\rho)^3} T \right). $$

The ranks of the scheme are bounded by

$$ \mathrm{rank}(\tilde{u}_n) \le \frac{16}{\rho^3} \left(\frac{1+\rho}{1-\rho}\right)^2 \left(1 + J \mathrm{rank}(V) \max\{r_1, r_2, \dots, r_n\}\right), $$

for all $1 \le n \le N$, and the intermediate ranks are at most twice as large.

*Proof.* For the global error bound, we prove by induction that

$$ (4.8) \qquad \|\tilde{u}_n - u(nh)\| \le \xi_n - \kappa_J h^{J+1}. $$

This is sufficient since $\rho = \Lambda_J C_V h$ and $nh \le T$. We initialise the scheme with $\tilde{u}_0 = \mathcal{H}_{\alpha_0} u(0)$, where $\alpha_0 = \min_{\|u(0)-v\| \le \eta} \sigma_{\max}(v) \ge \alpha_1$ and $r_0 = \eta^2 / \alpha_0^2 \le r_1$, which ensures that

$$ \|\tilde{u}_0 - u(0)\| \le \eta = \xi_0 - \kappa_J h^{J+1} \quad \text{and} \quad \mathrm{rank}(\tilde{u}_0) \le \frac{\|u(0) - S_{\alpha_0} u(0)\|^2}{\alpha_0^2} = r_0. $$

For $1 \le n \le N$, assume that (4.8) holds at index $n-1$, and denote $u_0^\dagger = u((n-1)h)$ and $u_0^* = \tilde{u}_{n-1}$ the values of the solution and the scheme at time $(n-1)h$. By Lemma 17,

$$ \|u^* - u\|_{W((n-1)h, nh)} \le \frac{\kappa_J h^{J+1} + \|\tilde{u}_{n-1} - u((n-1)h)\|}{1-\rho} \le \frac{\xi_{n-1}}{1-\rho} = \frac{1-\rho}{1+\rho} \frac{\varepsilon_n}{4}. $$

hence we satisfy the condition of Lemma 21. In addition, by Proposition 13,

$$ \begin{aligned} \|\tilde{u}_n - u(nh)\| & \le \delta_n + \frac{\rho \varepsilon_n}{1-\rho} + \kappa_{2J} h^{2J+1} + \xi_{n-1} - \kappa_J h^{J+1} \\ & \le \left(1 + 8\rho \frac{1+\rho}{(1-\rho)^3}\right) \xi_{n-1} + \kappa_{2J} h^{2J+1} - \kappa_J h^{J+1} \le \xi_n - \kappa_J h^{J+1}, \end{aligned} $$

where we used the fact that $1+x \le e^x$, and this completes the induction.

Next, we define $\bar{\alpha}_n$ for each interval $[(n-1)h, nh]$ such that

$$ \|u - S_{\bar{\alpha}_n} u\|_{W((n-1)h, nh)} = \frac{1-\rho}{1+\rho} \frac{\varepsilon_n}{4}, $$

and again set $\bar{r}_n = \varepsilon_n^2 / (\rho \bar{\alpha}_n^2)$. For the ranks, observe that for all $n \le T/h$, we have $\frac{1-\rho}{1+\rho} \frac{\varepsilon_n}{4} \ge \eta$, thus $\bar{\alpha}_n \ge \alpha_n$ and therefore

$$ \frac{\|u(t) - S_{\bar{\alpha}_n} u(t)\|^2}{\bar{\alpha}_n^2} \le \frac{\|u(t) - S_{\alpha_n} u(t)\|^2}{\alpha_n^2} \le r_n. $$

In particular, for any $1 \le n \le N$, the value of $\bar{r}_n$ satisfies

$$ \bar{r}_n \le \frac{16 (1+\rho)^2}{\rho (1-\rho)^2} r_n, $$

and we conclude by induction on Proposition 25. $\square$

**Remark 28.** We actually defined a scheme on the whole interval $[0, T]$, which satisfies

$$ \sup_{t \in [0, T]} \|\tilde{u}(t) - u(t)\| \le \frac{1}{1-\rho} (\eta + \kappa_J h^{J+1} + \kappa_{2J} h^{2J} T) \exp \left( 8\Lambda_J C_V \frac{1+\rho}{(1-\rho)^3} T \right). $$

The proof is essentially the same, with one additional use of Lemma 17 in the interval containing $t$, combined with (4.8) at $n = \lfloor t/h \rfloor$. The ranks of $\tilde{u}(t)$ are bounded in the same way as above, which justifies the quasi-optimality statements made in the introduction.
Remark 29. In Theorem 27, the bound on the ranks can only grow as n grows. If the ranks of the solution decrease over time, a closer inspection at the proof of Proposition 25 reveals that the bound on the ranks of the scheme can actually decrease, by a factor at best $1 - \rho^2$, from one interval to the next. The rate of decay $(1 - \rho^2)^n$ is however quite pessimistic, and the ranks computed numerically can track a faster decay, as shown in Section 6 (see Figure 8) for the parabolic case.

Remark 30 (Analysis of the SDC scheme). A central difference between Algorithms 3.1 and 3.2 is the fact that the soft thresholding operator is applied to the intermediate approximations at different points in the algorithm. While for the Gauss-Legendre-Picard iteration this is done after each application of the fixed-point iteration $\Phi_{\text{Picard}}$ in matrix form, the same can not be said about the SDC iterations. Here, the soft thresholding operator is applied after the component-wise application of $\Phi_{\text{SDC}}$, making it more difficult to use the non-expansiveness property. One possibility to fully analyze Algorithm 3.2 would thus be to decouple $\Phi_{\text{SDC}}$ into fixed-point iterations on each Gauss-Legendre node. Combined with the general contractivity of the SDC scheme, one could think about using different soft thresholding parameters (and strategies to decrease them) for the $J$ decoupled fixed-point iterations.

## 5. ADAPTATION TO PARABOLIC PROBLEMS

An essential ingredient of the described iterative algorithms is certainly the choice of quadrature points $t_1, t_2, \dots, t_J$. While for the time-dependent Schrödinger equation, we had numerous reasons to take Gauss-Legendre nodes, they turn out to be rather unfavorable in the parabolic setting: because of the isometry-preservation, the damping of high frequencies is correctly reproduced only for very small time steps $h$.

Recall the scalar test problem (2.18) used to analyze the stability properties of Runge-Kutta schemes. A somewhat stronger notion of stability can be defined by the so-called *L-stable* methods, for which the stability function $R$ needs to additionally satisfy

$$|R(z)| \to 0 \text{ as } \mathrm{Re}(z) \to -\infty,$$

that is, the numerical method reproduces the damping of modes with negative singular values. Note that this rules out preservation of isometries, since for rational functions $R$ we have $\lim_{x \to -\infty} |R(x)| = \lim_{y \to \infty} |R(iy)|$. Such a method can be obtained by instead using Radau-Legendre quadrature nodes on the respective subintervals.

**Example 31 (Radau-Legendre quadrature).** In contrast to the Gauss-Legendre quadrature, we force the last node to coincide with the right interval boundary $t_J = t_0 + h$. Optimizing over the remaining degrees of freedom, the quadrature formula (2.9) can at most attain a degree of accuracy $2J - 2$, that is, one less than the Gauss-Legendre quadrature formula. The optimal choice of points can again be computed by considering a tridiagonal matrix: it turns out it suffices to modify the bottom right entry of matrix (2.10), in such a way that the largest eigenvalue of the resulting matrix is exactly $t_0 + h$, see for example [2]. The Radau-Legendre quadrature weights and nodes can then be analogously derived from the eigenvalues and eigenvectors of the modified matrix. The resulting Gauss method (2.17) has then the order of convergence $2J - 1$.

The modified fixed-point formulation for the parabolic problem (2.5) was discussed in Example 2. Unfortunately, as already noted, we cannot directly translate the application of the Radau-Legendre-Picard iteration to this formulation, because the equivalent of formula (4.1) for the definition of $\Phi$ would include exponential operators of the form $e^{a(t_j - t_m)\Delta}$, which are unbounded when $t_j - t_m < 0$ and would thus lead to an extremely unfavorable coupling to the resolution of the spatial discretization. As before, we aim to avoid such a coupling between time steps size and spatial discretization in order to allow larger time steps.
In order to avoid taking exponentials of negative Laplacians, we instead only interpolate the second part of the integrand in (2.6), namely $G u(s) + f(s)$, by a Lagrange polynomial

$$p(s) = \sum_{m=1}^{J} (G u(t_m) + f(t_m)) \ell_m(s) = \sum_{m=1}^{J} 2b \partial_{x_1} \partial_{x_2} u(t_m) \ell_m(s) + f.$$

To compute the fixed-point application $\mathbf{\Phi}_{\text{Picard}}(\mathbf{u})_j$, for $1 \le j \le J$, we thus need to approximate integrals of the type

$$\int_{t_{j-1}}^{t_j} e^{a(t_j-s)\Delta} p(s) ds.$$

As the above integrand is still singular at $s = t_j$, we divide the integration domain by applying $N_j \in \mathbb{N}$ bisections that concentrate near $t_j$, that is,

$$[t_{j-1}, t_j] = [t_j - \tau_j, t_j - \frac{\tau_j}{2}] \cup \dots \cup [t_j - \frac{\tau_j}{2^{N_j-1}}, t_j - \frac{\tau_j}{2^{N_j}}] \cup [t_j - \frac{\tau_j}{2^{N_j}}, t_j],$$

where $\tau_j = t_j - t_{j-1}$, and use a quadrature formula $Q_a^b(g) \approx \int_a^b g(s) ds$ of moderate degree on each of these intervals. In the end, Duhamel's formula (2.6) is approximated by

$$\mathbf{\Phi}_{\text{Picard}}(\mathbf{u})_j = e^{a\tau_j\Delta} \mathbf{\Phi}_{\text{Picard}}(\mathbf{u})_{j-1} + \sum_{n=1}^{N_j} Q_{t_j - \frac{\tau_j}{2^{n-1}}}^{t_j - \frac{\tau_j}{2^n}} (e^{a(t_j-\cdot)\Delta} p(\cdot)) + Q_{t_j - \frac{\tau_j}{2^{N_j}}}^{t_j} (e^{a(t_j-\cdot)\Delta} p(\cdot)).$$

Note that we only need to compute the fixed-point application at the nodes $t_j$, since the interval boundary $t_0 + h$ is equal to $t_J$.

Concerning $\mathbf{\Phi}_{\text{SDC}}$, we use the same construction, with an additional explicit low-order integrator $\Psi_{t_{j-1},t_j}$ as in the dispersive setting (2.12). In the numerical tests of Section 6, we take the explicit Euler integrator

$$\Psi_{t_{j-1},t_j} = \tau_j e^{a\tau_j\Delta} G$$

applied to $\mathbf{\Phi}_{\text{SDC}}(\mathbf{u})_{j-1} - u_{j-1}$, and do not need to consider the source term $f$ since the correction is a difference of two solutions.

Due to the introduced secondary quadrature, the resulting (modified) Radau-Legendre-Picard iteration is no longer equivalent to the collocation formulation and corresponding Runge-Kutta method (2.17). However, as the numerical tests in the following section show, the resulting schemes still perform as expected, which makes them promising candidates for further analysis.

## 6. NUMERICAL EXPERIMENTS

In this section, we present numerical experiments for the time-dependent Schrödinger and parabolic model problems (1.1) and (2.5). All tests were performed using Julia version 1.10.0. Recall that we are considering evolution problems on product domains $\Omega_1 \times \Omega_2$, of the form

$$\partial_t u(t) = \mathcal{A}_t(u(t)), \quad t \in [0, T].$$

Here, we choose $\Omega_1 = \Omega_2 = [0, 1]$ and impose homogeneous Dirichlet boundary conditions. We thus have the $L^2(\Omega)$ orthonormal basis expansion of the solution

$$u(t)(x_1, x_2) = \sum_{k_1, k_2=1}^{\infty} 2u_{k_1, k_2}(t) \sin(\pi k_1 x_1) \sin(\pi k_2 x_2),$$

with $u_{k_1, k_2}(t) \in \mathbb{C}$ (or $\mathbb{R}$ in the parabolic setting) for $k_1, k_2 \in \mathbb{N}$. Galerkin projection onto the first $K$ frequencies in each coordinate leads to a system of ordinary differential equations

$$\partial_t \mathbf{u}(t) = \mathbf{A}_t(\mathbf{u}(t)), \quad t \in [0, T],$$

for some discrete operator $\mathbf{A}_t: \mathbb{C}^{K \times K} \to \mathbb{C}^{K \times K}$, and $\mathbf{u}(t)[k_1, k_2] \approx u_{k_1, k_2}(t)$. Note that using this Fourier basis expansion, we can compute the involved matrix exponentials explicitly: the Laplacian operator in matrix form is diagonal and thus its exponential can be computed by
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Left: decrease of the residual error estimator over the iterations for varying values of $\theta$ on a single interval $[t_0, t_0 + h]$. Right: evolution of the ranks using different strategies for the decrease of $\alpha$ and different iteration tolerances ($10^{-3}$ and $10^{-6}$) compared to the corresponding optimal ranks.</figcaption></figure>

applying the exponentials to its eigenvalues. Our aim is to obtain low-rank approximations of the discrete solution $\mathbf{u}(t)$.

6.1. **Schrödinger equation.** First, we consider the time-dependent Schrödinger equation with potential

$$V(x_1, x_2) = \cos(n\pi x_1) \cos(m\pi x_2)$$

with some $m, n \in \mathbb{N}$. The spatial semidiscretization yields

$$\partial_t \mathbf{u}(t) = -i (\mathbf{S} \otimes \mathbf{I} + \mathbf{I} \otimes \mathbf{S} + \mathbf{V}_n \otimes \mathbf{V}_m) \mathbf{u}(t),$$

with the matrices

$$\mathbf{S}[k_1, k_2] = \int_0^1 2 (\partial_{x_1} \sin(k_1 \pi x_1)) (\partial_{x_1} \sin(k_2 \pi x_1)) dx_1,$$

$$\mathbf{V}_n[k_1, k_2] = \int_0^1 2 \cos(n\pi x_1) \sin(k_1\pi x_1) \sin(k_2\pi x_1) dx_1.$$

As initial data, we choose

$$\mathbf{u}_0 = \mathbf{v}_K \otimes \mathbf{v}_K, \text{ where } \mathbf{v}_K \in \mathbb{C}^K \text{ with } v_K[k] = \frac{1}{k}$$

and $K = 300$ frequencies. In the following, we use $n = m = 1$ for the cosine potential and a time step $h = 0.1$ with $N = 5$ intervals and $J = 11$ Gauss-Legendre quadrature nodes. Note that here it is important that with the fixed-point formulation and iterative schemes that we use, $h$ can be chosen independently of the spatial discretization determined by $K$.

For the implementations of the two iterative solvers, we follow Algorithms 3.1 and 3.2. For target accuracies of $\eta = 10^{-3}$ and $\eta = 10^{-6}$, the recompression parameter at subintervals boundaries is taken as $\delta_n = 10^{-4}$ and $\delta_n = 10^{-7}$, respectively. In practice, when computing the updated approximations $\mathbf{v}^{k+1}$, we use additional recompressions up to a relative error of $10^{-3}$ after every operation on low-rank functions. Note that in case of SDC, we follow the approach described in [20] for the computation of the $\phi_j^k$ in line 7 of the algorithm from the error and residual functions and hence need here an additional smaller recompression tolerance of $10^{-6}$ (relative to the current approximation of the residual function), as already mentioned in Section 3.2, due to an observed accumulation of recompression errors.

Another important parameter is the rate by which $\alpha$ is decreased. As can be deduced from the left part of Figure 6, smaller values of $\theta \in (0, 1)$ allow for a faster convergence of the fixed-point iteration, even though we only have rank bounds for $\theta \ge \rho$, and they are tighter for larger values of $\theta$. One can also observe that the adaptive decrease of $\alpha$ requires
<figure><img src="image_8.png" /><figcaption>FIGURE 7. Left: evolution of the absolute error with respect to a reference solution. Right: deviation of the norm of the computed approximations from $||u_0|| = 1$ over time.</figcaption></figure>

about twice as many iterations as the constant decrease of $\alpha$. This can be explained by the following consideration: in the adaptive setting, $\alpha_{k+1}$ is smaller than $\alpha_k$ whenever

$$ \|\mathbf{u}^{k+1} - \mathbf{u}^k\|_J \le \frac{1}{6} \|\mathbf{u}^k - \mathbf{\Phi}\mathbf{u}^k\|_J, $$

is satisfied, where $\mathbf{\Phi}$ now denotes the discretized fixed-point application on $(\mathbb{C}^{K \times K})^J$. This condition only seems to occur when $\mathbf{u}^{k+1}$ and $\mathbf{u}^k$ have been thresholded with the same parameter $\alpha_k = \alpha_{k-1}$, because in that case they are both close to the modified fixed point $\mathbf{u}_{\alpha_k} = \mathbf{u}_{\alpha_{k-1}}$. Hence, $\alpha$ can only be decreased every second iteration. Moreover, the distance to the fixed point is dominated by the soft thresholding error, which is almost proportional to $\alpha$, that is,

$$ (6.1) \qquad \|\mathbf{u}^{k+1} - \mathbf{u}^*\|_J \approx \|(\mathcal{I} - S_{\alpha_k})(\mathbf{\Phi}\mathbf{u}^k)\|_J \approx \alpha_k \sqrt{\mathrm{rank}(\mathbf{u}^{k+1})}, $$

justifying the dependence of the accuracy in the number of executed iterations. In the following, we thus use $\theta = \frac{1}{2}$ and $\theta = \frac{1}{5}$ for the constant and adaptive decrease of $\alpha$, respectively, in order to maintain comparable numbers of iterations. Note that the illustration on the left of Figure 6 does not change visibly if we replace the SDC scheme by the simpler Picard iteration, which is another indicator of the fact that the iteration error is dominated by the soft thresholding error.

We continue our numerical investigations by considering the performance of the Picard iterative solver combined with different iteration tolerances. The obtained evolution of ranks is visualized on the right of Figure 6: for both iteration tolerances the Picard method combined with the two possible reduction strategies of $\alpha$ follows, up to a constant, the trajectories of the optimal ranks. The fixed-point reference solution was computed in full matrix format using the Picard iteration with iteration tolerance $10^{-12}$ and without any recompression or truncation.

As can be deduced from the left part of Figure 7, the absolute error accumulates only slowly over time and remains below the tolerances used for the iterative solvers over the entire interval $[0, 0.5]$. In particular, we have a significantly smaller error at the subinterval boundaries: while on the inner nodes, the error is in general expected to be of order $O(h^{J+1})$, we expect an accuracy of $O(h^{2J})$ at the subinterval boundaries.

In addition, the recompression at the boundary nodes is expected to preserve the $L^2(\Omega)$ norm of the scheme much better than soft thresholding. Indeed, for a function $u$ with singular values ($\sigma_k$) and a threshold $\alpha \ll 1$,

$$ \|u\|^2 - \|S_{\alpha}(u)\|^2 = \sum_{k \in \mathbb{N}} \sigma_k^2 - \max(\sigma_k - \alpha, 0)^2 = \sum_{\sigma_k \ge \alpha} (2\alpha\sigma_k - \alpha^2) + \sum_{\sigma_k < \alpha} \sigma_k^2 = O(\alpha). $$
On the other hand, the recompression behaves similarly to hard thresholding, for which the norm deviation is much smaller

$$ \|u\|^2 - \|\mathcal{H}_\alpha(u)\|^2 = \sum_{\sigma_k > \alpha} \sigma_k^2 = \mathcal{O}(\alpha^2). $$

This can be seen on the right of Figure 7. In fact, we observe a preservation of the norm at the subinterval boundaries up to iteration and recompression tolerances.

6.2. **Parabolic problem.** The second running example that we consider is the heat equation with anisotropic diffusion (2.5). The problem is approximated, analogously to the first example, by an ODE of the form

$$ \partial_t \mathbf{u} = (-a(\mathbf{S} \otimes \mathbf{I} + \mathbf{I} \otimes \mathbf{S}) + 2b\mathbf{B} \otimes \mathbf{B}) \mathbf{u} + \mathbf{f} \otimes \mathbf{f}, \quad \mathbf{u} \in \mathbb{R}^{K \times K}, $$

with parameters $K = 500$, $a = 1$ and $b = -\frac{1}{2}$, and where the matrix $\mathbf{B}$ is given by

$$ \mathbf{B}[k_1, k_2] = \int_0^1 2(\partial_{x_1} \sin(k_1 \pi x_1)) \sin(k_2 \pi x_1) dx_1. $$

The initial data $\mathbf{u}_0$ and source term $\mathbf{f}$ are chosen in order to first observe an increase in ranks, followed by a quick decrease due to the smoothing properties of the equation, and finally a slight long-term rise of the ranks caused by the source term. More precisely, we set

$$ \mathbf{u}_0[k_1, k_2] = \begin{cases} \delta_{k_1, k_2}, & \text{for } 21 \le k_1, k_2 \le 30, \\ 0, & \text{otherwise,} \end{cases} \quad \text{and} \quad \mathbf{f}[k] = \begin{cases} k, & \text{for } 1 \le k \le 10, \\ 0, & \text{otherwise.} \end{cases} $$

In order to best capture the initial rank evolution, we use a stepsize of $h = 0.001$ with $N = 10$ intervals and $J = 11$ Radau-Legendre quadrature nodes. Note that with this choice of points, since $t_J = t_0+h$, we also apply the soft thresholding operator to the approximations computed at the right subinterval boundaries instead of the recompression operator. Larger step sizes would be possible for this problem without any stability issues. However, resolving the initial rank fluctuations then requires a correspondingly higher polynomial degree.

For the secondary quadrature, we take $N_j = 5$ refinements by bisection of the intervals $[t_0, t_j]$ for $1 \le j \le J$ and use Radau-Legendre quadrature $Q_a^b$ with 5 nodes on each of the resulting subintervals. Without this special form of integral approximation, a much larger number of interpolation points (at least $7J$) is needed to reach tolerances below an accuracy of $10^{-9}$. Otherwise, the fixed-point iterations only converge up to a certain error, and then rapidly diverge.

We again use Algorithms 3.1 and 3.2, respectively, and make use of additional low-rank recompressions in intermediate steps. For recompressing intermediate solutions after the summation of two low-rank matrices, we use a relative tolerance of $10^{-4}$, see also the previous example, except when we approximate the residual function within the SDC iteration, where a lower relative recompression tolerance of $10^{-6}$ is required.

The fixed-point reference solution was computed in full matrix format using the Picard iteration with an iteration tolerance of $10^{-9}$, and 10 dyadic refinements with 5 quadrature nodes on each subinterval for the secondary quadrature.

The quasi-optimality of the proposed schemes is highlighted in Figure 8: here, we compare the resulting ranks and accuracies using different iterative solvers, which are iterated up to an accuracy of $10^{-4}$. The parameter $\theta$ is chosen again such that the numbers of needed iterations per interval are comparable, that is, $\theta = 0.5$ in the case of a constant decrease of the soft threshold and $\theta = 0.2$ in the adaptive setting, in which $\alpha$ is decreased whenever

$$ \|\mathbf{u}^{k+1} - \mathbf{u}^k\|_J \le c \|\mathbf{u}^k - \mathbf{\Phi}(\mathbf{u}^k)\|_J $$

is satisfied with $c = \frac{2}{5}$. Note that $c$ was chosen such that the iterations on an interval were carried out at least twice (in our example: exactly twice), before being reduced. It becomes evident that for all four possible combinations of iterative thresholding schemes the obtained approximation ranks exhibit a quasi-optimal behavior by following, up to a constant only slightly larger than one, the trajectory of the optimal ranks, which are obtained by hard
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Left: evolution of the ranks using different fixed-point iterations and strategies for decreasing of $\alpha$, with an iteration tolerance of $10^{-4}$, compared to the optimal ranks of a truncated reference solution. Right: corresponding evolution of the absolute errors with respect to the reference solution.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9. Comparison of the ranks of the computed approximation using Picard iterations and a constant decrease of $\alpha$ with $\theta = 0.5$, its truncation to a tolerance of $10^{-3}$, the arising intermediate ranks in the computation of the approximation and the truncated reference solution.</figcaption></figure>

thresholding the full-rank reference solution to a tolerance of $10^{-4}$. We also observe a linear increase in the absolute error over time, as one would expect from the analysis. Note that the errors are larger when $\alpha$ is decreased adaptively, or when using the Picard iteration, which is compensated by ranks slightly closer to optimal for these methods.

Regarding these numerical experiments, it is not completely clear how much of an advantage the SDC scheme yields over the Picard iteration, since the resulting error seems dominated by the value of the soft thresholding parameter. This observation also explains the much smaller absolute error in the first interval: as the ranks are higher, the fixed-point iteration reaches a smaller value of $\alpha$, as suggested by (6.1), leading to a better accuracy.

In Figure 9, we plot the minimal and maximal values of the ranks of all intermediate low-rank matrices appearing in the computation of the fixed-point iterations, at any of the times $t_j$. These ranks remain comparable to the ranks of the scheme, which only involves the last iteration in each interval, up to a factor about 2, corroborating the last statement of Theorem 27. In particular, this shows that the secondary quadrature used in the approximation of the integrals over the exponentials does not cause an unbounded growth of the intermediate ranks. The figure is obtained for the Picard iterative solver with truncation tolerance $10^{-4}$ and a constant decrease of $\alpha$; using instead an adaptive decrease of $\alpha$ produces almost identical results. Moreover, the scheme is able to recover the optimal ranks almost exactly, when truncating to a slightly lower tolerance of $10^{-3}$ after the whole evolution has
been computed. At two points in time, we even observe smaller ranks for the truncated approximation, pointing out to the fact that the error of numerical integration does not necessarily result in an increase of the ranks.

## 7. CONCLUSION

In this work, we have proposed and analyzed a high-order time integration scheme that can be ensured to achieve the desired accuracy, while maintaining ranks that are controlled in terms of best approximation ranks. In particular, we avoid systematic errors that can occur in schemes based on tangent space projections. To the best of our knowledge, this is the first instance of a time integration scheme that achieves such a balance between accuracy and approximation ranks, aside from the method based on a space-time variational formulation in [8].

While the present method does not require such a variational formulation, in aiming for the iterative refinement of solutions on longer time intervals it can be regarded as a compromise between space-time methods and truncated time integration schemes. As such, a question for future work is how to choose different – and potentially, adaptive – basis functions for the time variable. Note, however, that the present choice of polynomials in time has certain advantages concerning, for example, the well-understood properties of mass and energy preservation. There also remain some open questions regarding the analysis of the present high-order polynomial approximations, in particular the asymptotics of averaged Lebesgue constants noted in Remark 9.

We will address two generalizations as subjects for further work: in the context of Schrödinger-type problems, applications to more general classes of potentials and the application to nonlinear problems, such as nonlinear Schrödinger equations; and more generally, to further types of different PDEs. In particular, a more detailed analysis of the more involved scheme for the parabolic case remains open at this point. A second major direction of generalization concerns the extension of the method to high-dimensional problems via higher-order low-rank tensor representations. Although the basic scheme can be extended to this case directly based on results for hierarchical tensors in [7,9], the analysis becomes more subtle due to the need for tracking dependencies on the dimensionality.

## ACKNOWLEDGEMENTS

Co-funded by the European Union (ERC, COCOA, 101170147). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

## REFERENCES

[1] Yvonne Alama Bronsard. Error analysis of a class of semi-discrete schemes for solving the Gross-Pitaevskii equation at low regularity. *Journal of Computational and Applied Mathematics*, 418:114632, 2023.

[2] Hessah Alqahtani, Carlos F. Borges, Dušan Lj. Djukić, Rada M. Mutavdžić Djukić, Lothar Reichel, and Miodrag M. Spalević. Computation of pairs of related Gauss-type quadrature rules. *Applied Numerical Mathematics*, 208:32–42, 2025.

[3] Roman Andreev and Christine Tobler. Multilevel preconditioning and low-rank tensor iteration for space-time simultaneous discretizations of parabolic PDEs. *Numerical Linear Algebra with Applications*, 22(2):317–337, 2015.

[4] Daniel Appelö and Yingda Cheng. Robust implicit adaptive low rank time-stepping methods for matrix differential equations. *Journal of Scientific Computing*, 102(3):81, 2025.

[5] Jeffrey M. Aristoff, Joshua T. Horwood, and Aubrey B. Poore. Orbit and uncertainty propagation: A comparison of Gauss-Legendre-, Dormand-Prince-, and Chebyshev-Picard-based approaches. *Celestial Mechanics and Dynamical Astronomy*, 118(1):13–28, 2014.

[6] Markus Bachmayr. Low-rank tensor methods for partial differential equations. *Acta Numerica*, 32:1–121, 2023.

[7] Markus Bachmayr and Wolfgang Dahmen. Adaptive near-optimal rank tensor approximation for high-dimensional operator equations. *Foundations of Computational Mathematics*, 15(4):839–898, 2015.
[8] Markus Bachmayr and Manfred Faldum. A space-time adaptive low-rank method for high-dimensional parabolic partial differential equations. *Journal of Complexity*, 82:Paper No. 101839, 70, 2024.

[9] Markus Bachmayr and Reinhold Schneider. Iterative methods based on soft thresholding of hierarchical tensors. *Foundations of Computational Mathematics*, 17:1037–1083, 2017.

[10] Thomas Boiveau, Virginie Ehrlacher, Alexandre Ern, and Anthony Nouy. Low-rank approximation of linear parabolic equations by space-time tensor Galerkin methods. *ESAIM Math. Model. Numer. Anal.*, 53(2):635–658, 2019.

[11] Tommaso Buvoli. A class of exponential integrators based on spectral deferred correction. *SIAM Journal on Scientific Computing*, 42(1):A1–A27, 2020.

[12] Benjamin Carrel and Bart Vandereycken. Projected exponential methods for stiff dynamical low-rank approximation problems. *arXiv preprint arXiv:2312.00172*, 2023.

[13] Mathew Causley and David Seal. On the convergence of spectral deferred correction methods. *Communications in Applied Mathematics and Computational Science*, 14(1):33–64, 2019.

[14] Gianluca Ceruti, Jonas Kusch, and Christian Lubich. A rank-adaptive robust integrator for dynamical low-rank approximation. *BIT Numerical Mathematics*, 62(4):1149–1174, 2022.

[15] Gianluca Ceruti and Christian Lubich. An unconventional robust integrator for dynamical low-rank approximation. *BIT Numerical Mathematics*, 62(1):23–44, 2022.

[16] Charles W. Clenshaw. The numerical solution of linear differential equations in Chebyshev series. *Mathematical Proceedings of the Cambridge Philosophical Society*, 53(1):134–149, 1957.

[17] Charles W. Clenshaw and Harry J. Norton. The solution of nonlinear ordinary differential equations in Chebyshev series. *The Computer Journal*, 6(1):88–92, 1963.

[18] Alec Dektor, Abram Rodgers, and Daniele Venturi. Rank-adaptive tensor methods for high-dimensional nonlinear PDEs. *Journal of Scientific Computing*, 88(2):36, 2021.

[19] Mi-Song Dupuy, Virginie Ehrlacher, and Clément Guillot. A space-time variational formulation for the many-body electronic Schrödinger evolution equation. *arXiv preprint arXiv:2405.18094*, 2024.

[20] Alok Dutt, Leslie Greengard, and Vladimir Rokhlin. Spectral deferred correction methods for ordinary differential equations. *BIT Numerical Mathematics*, 40:241–266, 2000.

[21] Gene H. Golub and John H. Welsch. Calculation of Gauss quadrature rules. *Math. Comp.*, 23:221–230; addendum, ibid. 23 (1969), no. 106, loose microfiche suppl. A1–A10, 1969.

[22] Lars Grasedyck, Daniel Kressner, and Christine Tobler. A literature survey of low-rank tensor approximation techniques. *GAMM-Mitteilungen*, 36(1):53–78, 2013.

[23] Wolfgang Hackbusch. *Tensor spaces and numerical tensor calculus*, volume 56 of Springer Series in Computational Mathematics. Springer, Cham, 2019. Second edition.

[24] William W. Hager, Hongyan Hou, and Anil V. Rao. Lebesgue constants arising in a class of collocation methods. *IMA Journal of Numerical Analysis*, 2016.

[25] Thomas Hagstrom and Ruhai Zhou. On the spectral deferred correction of splitting methods for initial value problems. *Communications in Applied Mathematics and Computational Science*, 1(1):169–205, 2007.

[26] Ernst Hairer, Gerhard Wanner, and Syvert P. Nørsett. *Solving Ordinary Differential Equations I*, volume 8 of Springer Series in Computational Mathematics. Springer Berlin Heidelberg, Berlin, Heidelberg, 1993.

[27] Marlis Hochbruck, Jan Leibold, and Alexander Ostermann. On the convergence of Lawson methods for semilinear stiff problems. *Numerische Mathematik*, 145:553–580, 2020.

[28] Marlis Hochbruck and Alexander Ostermann. Exponential integrators. *Acta Numerica*, 19:209–286, 2010.

[29] Jingfang Huang, Jun Jia, and Michael Minion. Accelerating the convergence of spectral deferred correction methods. *Journal of Computational Physics*, 214(2):633–656, 2006.

[30] Aly-Khan Kassam and Lloyd N. Trefethen. Fourth-order time-stepping for stiff PDEs. *SIAM Journal on Scientific Computing*, 26(4):1214–1233, 2005.

[31] Othmar Koch and Christian Lubich. Dynamical low-rank approximation. *SIAM Journal on Matrix Analysis and Applications*, 29(2):434–454, 2007.

[32] Karolina Kropielnicka and Juan Carlos Del Valle. Second order exponential splittings in the presence of unbounded and time-dependent operators: Construction and convergence. *SIAM Journal on Numerical Analysis*, 63(1):288–305, 2025.

[33] Hei Yin Lam, Gianluca Ceruti, and Daniel Kressner. Randomized low-rank Runge-Kutta methods. *SIAM Journal on Matrix Analysis and Applications*, 46(2):1587–1615, 2025.

[34] John Douglas Lawson. Generalized Runge-Kutta processes for stable systems with large Lipschitz constants. *SIAM Journal on Numerical Analysis*, 4(3):372–380, 1967.

[35] Anita T. Layton and Michael L. Minion. Implications of the choice of quadrature nodes for Picard integral deferred corrections methods for ordinary differential equations. *BIT*, 45(2):341–373, 2005.

[36] Shun Li, Yan Jiang, and Yingda Cheng. High-order implicit low-rank method with spectral deferred correction for matrix differential equations. *arXiv preprint arXiv:2412.09400*, 2024.

[37] Christian Lubich and Ivan V. Oseledets. A projector-splitting integrator for dynamical low-rank approximation. *BIT*, 54(1):171–188, 2014.
[38] Aleksandr Semenovich Markus. The eigen- and singular values of the sum and product of linear operators. *Russian Mathematical Surveys*, 19(4):91, 1964.

[39] Fabio Nobile and Sébastien Riffaud. Robust high-order low-rank BUG integrators based on explicit Runge-Kutta methods. *arXiv preprint arXiv:2502.07040*, 2025.

[40] Abram Rodgers and Daniele Venturi. Implicit integration of nonlinear evolution equations on tensor manifolds. *Journal of Scientific Computing*, 97(2):33, 2023.

[41] Josef Stoer and Roland Bulirsch. *Introduction to Numerical Analysis*. Springer New York, 1980.

[42] Marco Sutti and Bart Vandereycken. Implicit low-rank Riemannian schemes for the time integration of stiff partial differential equations. *Journal of Scientific Computing*, 101(1):3, 2024.

[43] Martin Weiser. Faster SDC convergence on non-equidistant grids by DIRK sweeps. *BIT Numerical Mathematics*, 55(4):1219-1241, 2015.

[44] Mingru Yang and Steven R. White. Time-dependent variational principle with ancillary Krylov subspace. *Physical Review B*, 102(9):094315, 2020.

## 8. APPENDIX

In the following, we include some auxiliary lemmas needed in the proof of Theorem 27, and recall their proofs for the sake of completeness.

**Lemma 32** (Gauss-Legendre quadrature error). Let $\omega_j$ and $t_j$ for $1 \le j \le J$ be the Gauss-Legendre quadrature weights and nodes on the interval $[0, h]$. For any function $g$ in $C^{2J}(0, h; L^2(\Omega))$, it holds that

$$ \left\| \sum_{j=1}^{J} \omega_j g(t_j) - \int_{0}^{h} g(s) ds \right\| \le \frac{(J!)^4}{(2J+1)!(2J)!^2} \max_{s \in [0,h]} \| \partial_t^{2J} g(s) \| h^{2J+1}. $$

*Proof.* Denote by $p \in \mathbb{P}_{2J-1}$ the polynomial that interpolates the function values $g(t_j)$ and the corresponding derivatives $g'(t_j)$ at all $J$ quadrature nodes $t_j$. As the Gauss-Legendre quadrature is exact for elements of $\mathbb{P}_{2J-1}$, we have

$$ \sum_{j=1}^{J} \omega_j g(t_j) = \sum_{j=1}^{J} \omega_j p(t_j) = \int_{0}^{h} p(s) ds, $$

hence

$$ \left\| \sum_{j=1}^{J} \omega_j g(t_j) - \int_{0}^{h} g(s) ds \right\| \le \int_{0}^{h} \|g(s) - p(s)\| ds. $$

Next, define the function $z(t) = \prod_{j=1}^{J}(t-t_j)$ and set, for some $s \in [0, h] \setminus \{t_1, \dots, t_J\}$,

$$ h(t) = g(t) - p(t) - \frac{z(t)^2}{z(s)^2} (g(s) - p(s)). $$

As $h(x) = h(t_j) = \partial_t h(t_j) = 0$ for all $1 \le j \le J$, the function $h$ has at least $2J+1$ zeros, counted with multiplicity, in the interval $[0, h]$. By a repeated application of Rolle's theorem, its derivative of order $2J$, given by

$$ \partial_t^{2J} h(t) = \partial_t^{2J} g(t) - (2J)! \frac{g(s) - p(s)}{z(s)^2}, $$

cancels at some point $\xi(s) \in [0, h]$. This implies that

$$ \|g(s) - p(s)\| = \left\| \frac{z(s)^2}{(2J)!} \partial_t^{2J} g(\xi(s)) \right\| \le \frac{z(s)^2}{(2J)!} \max_{\xi \in [0,h]} \|\partial_t^{2J} g(\xi)\|. $$

Recalling that the $L^2(\Omega)$ normalized Legendre polynomial $L_J$ of degree $J$ on the reference interval $[-1, 1]$ can be written as

$$ L_J(x) = \gamma_J \prod_{j=1}^{J} (x - x_j) = \gamma_J z\left(\frac{x+1}{2}h\right), \quad \text{where} \quad \gamma_J = \sqrt{\frac{2J+1}{2} \frac{(2J)!}{2^J J!^2}}, $$
since the Gauss-Legendre nodes $t_j = \frac{x_j+1}{2}h$ are rescaled versions of the roots $x_j$ of $L_J$. We conclude with

$$ \int_0^h \frac{z(s)^2}{(2J)!} ds = \frac{h/2}{(2J)!} \int_{-1}^1 \left( \prod_{j=1}^J \frac{x-x_j}{2} h \right)^2 dx = \frac{(h/2)^{2J+1}}{(2J)!\gamma_J^2} = \frac{(J!)^4 h^{2J+1}}{(2J+1)!(2J)!^2}. \quad \square $$

A similar, but weaker, result also holds within the interval $[0, h]$.

**Lemma 33 (Quadrature error inside the interval).** For $t \in [0, h]$ and $g \in C^J(0, h; L^2(\Omega))$ it holds that

$$ \left\| \int_0^t \left( g(s) - \sum_{m=1}^J \ell_m(s) g(t_m) \right) ds \right\| \le \frac{J!}{\sqrt{2J+1}(2J)!} \max_{s \in [0,h]} \| \partial_t^J g(s) \| h^{J+1}. $$

*Proof.* Define again $z(t) = \prod_{j=1}^J (t - t_j)$, and let

$$ p = \sum_{m=1}^{J} g(t_m) \ell_m \in \mathbb{P}_{J-1} $$

be the Lagrange interpolant of $g$. For $s \ne t_1, \dots, t_J$, the function

$$ h : t \mapsto g(t) - p(t) - \frac{z(t)}{z(s)}(g(s) - p(s)) $$

vanishes at $s$ and all nodes $t_1, \dots, t_J$, hence there exists some $\xi(s) \in [0, h]$ such that

$$ 0 = \partial_t^J h(\xi(s)) = \partial_t^J g(\xi(s)) - \frac{J!}{z(s)} (g(s) - p(s)), $$

from which we can conclude that

$$ \|g(s) - p(s)\| \le \frac{|z(s)|}{J!} \max_{\xi \in [0,h]} \|\partial_t^J g(\xi)\|. $$

We can therefore bound the quadrature error by

$$ \left\| \int_0^t (g(s) - p(s)) ds \right\| \le \int_0^t \|g(s) - p(s)\| ds \le \int_0^h \frac{|z(s)|}{J!} ds \max_{\xi \in [0,h]} \|\partial_t^J g(\xi)\|, $$

and conclude the proof by a Cauchy-Schwarz inequality

$$ \left( \int_0^h \frac{|z(s)|}{J!^2} ds \right)^2 \le \frac{h}{J!^2} \int_0^h z(s)^2 ds = \frac{h^{2J+2}}{2J+1} \frac{J!^2}{(2J)!^2}. \quad \square $$

Finally, we need the following comparison between Frobenius distance and $\ell^2$ distance of singular values.

**Lemma 34 (Mirsky inequality).** Let $u, v \in L^2(\Omega_1 \times \Omega_2)$, and denote their singular value decomposition as

$$ u = \sum_{k=1}^{\infty} \sigma_k(u) u_k^{(1)} \otimes u_k^{(2)} \quad \text{and} \quad v = \sum_{k=1}^{\infty} \sigma_k(v) v_k^{(1)} \otimes v_k^{(2)}, $$

with the singular values $\sigma_k$ in non-increasing order. It then holds that

$$ \|u - v\|_{L^2(\Omega_1 \times \Omega_2)}^2 \ge \sum_{k=1}^{\infty} |\sigma_k(u) - \sigma_k(v)|^2. $$

This result was originally shown to hold in a finite-dimensional setting. Its extension to Hilbert-Schmidt operators, or equivalently to $L^2$ functions on a separable tensor product domain, can be found in [38, Corollary 5.3]. We include a self-contained proof below.
*Proof.* First, we restrict ourselves to hermitian operators, by considering eigendecompositions

$$U := \begin{pmatrix} 0 & u \\ u^* & 0 \end{pmatrix} = \sum_{k \in \mathbb{Z} \setminus \{0\}} \lambda_k U_k U_k^* \quad \text{and} \quad V := \begin{pmatrix} 0 & v \\ v^* & 0 \end{pmatrix} = \sum_{k \in \mathbb{Z} \setminus \{0\}} \mu_k V_k V_k^*$$

with

$$\lambda_{\pm k} = \pm \sigma_k(u), \quad U_{\pm k} = \frac{1}{\sqrt{2}} \begin{pmatrix} u_k^{(1)} \\ \pm u_k^{(2)} \end{pmatrix}, \quad \mu_{\pm k} = \pm \sigma_k(v) \quad \text{and} \quad V_{\pm k} = \frac{1}{\sqrt{2}} \begin{pmatrix} v_k^{(1)} \\ \pm v_k^{(2)} \end{pmatrix}.$$

Using $(U_k)_{k \in \mathbb{Z} \setminus \{0\}}$ and $(V_\ell^*)_{\ell \in \mathbb{Z} \setminus \{0\}}$ as orthonormal bases for the columns and rows, we obtain in the Frobenius norm

$$\|U - V\|_F^2 = \left\| \sum_{k, \ell \in \mathbb{Z} \setminus \{0\}} (\lambda_k - \mu_\ell) U_k (U_k^* V_\ell) V_\ell^* \right\|_F^2 = \sum_{k, \ell \in \mathbb{Z} \setminus \{0\}} p_{k, \ell} (\lambda_k - \mu_\ell)^2,$$

where $p_{k,\ell} = |U_k^* V_\ell|^2$. The matrix $(p_{k,\ell})_{k,\ell \in \mathbb{Z} \setminus \{0\}}$ is bistochastic, that is, it has non-negative entries satisfying

$$\sum_{\ell \in \mathbb{Z} \setminus \{0\}} p_{k,\ell} = U_k^* \sum_{\ell \in \mathbb{Z} \setminus \{0\}} V_\ell V_\ell^* U_k = \|U_k\|^2 = 1, \quad \text{for all } k \in \mathbb{Z} \setminus \{0\},$$

and

$$\sum_{k \in \mathbb{Z} \setminus \{0\}} p_{k, \ell} = V_\ell^* \sum_{k \in \mathbb{Z} \setminus \{0\}} U_k U_k^* V_\ell = \|V_\ell\|^2 = 1, \quad \text{for all } \ell \in \mathbb{Z} \setminus \{0\}.$$

Therefore, it is in the closed convex hull of permutation matrices $p^\sigma = (1_{\ell=\sigma(k)})_{k,\ell}$, with $\sigma$ a permutation of $\mathbb{Z} \setminus \{0\}$. For such matrices, a simple reordering inequality shows that

$$\sum_{k, \ell \in \mathbb{Z} \setminus \{0\}} p_{k, \ell}^\sigma (\lambda_k - \mu_\ell)^2 = \sum_{k \in \mathbb{Z} \setminus \{0\}} (\lambda_k - \mu_{\sigma(k)})^2 = \sum_{k \in \mathbb{Z} \setminus \{0\}} \lambda_k^2 - 2\lambda_k \mu_{\sigma(k)} + \mu_{\sigma(k)}^2 \geq \sum_{k \in \mathbb{Z} \setminus \{0\}} (\lambda_k - \mu_k)^2.$$

Taking convex combinations, this remains true with $p$ instead of $p^\sigma$, and we conclude the inequality with

$$2\|u-v\|_{L^2(\Omega_1 \times \Omega_2)}^2 = \|U-V\|_F^2 \geq \sum_{k \in \mathbb{Z} \setminus \{0\}} (\lambda_k - \mu_k)^2 = 2 \sum_{k=1}^{\infty} |\sigma_k(u) - \sigma_k(v)|^2. \quad \square$$