Time integration of symmetric and anti-symmetric low-rank matrices and Tucker tensors

TL;DR

This paper introduces a novel numerical integrator for symmetric and anti-symmetric low-rank matrices and Tucker tensors that preserves symmetry properties and exhibits robust error behavior, improving upon existing methods.

Contribution

It presents a new integrator that maintains symmetry or anti-symmetry in low-rank matrix and tensor approximations, unlike previous projector-splitting methods.

Findings

Exact reproduction of low-rank matrices and tensors.

Robust error behavior in the presence of small singular values.

Numerical experiments demonstrating the integrator's effectiveness.

Abstract

A numerical integrator is presented that computes a symmetric or skew-symmetric low-rank approximation to large symmetric or skew-symmetric time-dependent matrices that are either given explicitly or are the unknown solution to a matrix differential equation. A related algorithm is given for the approximation of symmetric or anti-symmetric time-dependent tensors by symmetric or anti-symmetric Tucker tensors of low multilinear rank. The proposed symmetric or anti-symmetric low-rank integrator is different from recently proposed projector-splitting integrators for dynamical low-rank approximation, which do not preserve symmetry or anti-symmetry. However, it is shown that the (anti-)symmetric low-rank integrators retain favourable properties of the projector-splitting integrators: low-rank time-dependent matrices and tensors are reproduced exactly, and the error behaviour is robust to the presence of small singular values, in contrast to standard integration methods applied to the differential equations of dynamical low-rank approximation. Numerical experiments illustrate the behaviour of the proposed integrators.

Equations

$$\text{accentset}\text{mbox}.\mathbf{A}(t)=\mathbf{F}(t,\mathbf{A}(t)),\mathbf{A}(t_0)={\mathbf{A}}_0.$$

$$\text{accentset}\text{mbox}.\mathbf{Y}(t)=\mathrm{P}(\mathbf{Y}(t))\mathbf{F}(t,\mathbf{Y}(t)),\mathbf{Y}(t_0)={\mathbf{Y}}_0,$$

$$\mathbf{Y}(t)=\mathbf{U}(t)\mathbf{S}(t)\mathbf{V}(t{)}^{\top },$$

$$\mathrm{P}(\mathbf{Y})\text{Z}={\text{ZVV}}^{\top }-{\text{UU}}^{\top }{\text{Z VV}}^{\top }+{\text{UU}}^{\top }\text{Z}.$$

$$\dot{\text{K}}(t)=\mathbf{F}(t,\text{K}(t){\mathbf{V}}_0^{\top }){\mathbf{V}}_0,\text{K}(t_0)={\mathbf{U}}_0{\mathbf{S}}_0.$$

$$\dot{\mathbf{S}}(t)=-{\mathbf{U}}_1^{\top }\mathbf{F}(t,{\mathbf{U}}_1\mathbf{S}(t){\mathbf{V}}_0^{\top }){\mathbf{V}}_0,\mathbf{S}(t_0)={\hat{\mathbf{S}}}_1,$$

$$\dot{\text{L}}(t)=\mathbf{F}(t,{\mathbf{U}}_1\text{L}(t{)}^{\top }{)}^{\top }{\mathbf{U}}_1,\text{L}(t_0)={\mathbf{V}}_0{\tilde{\mathbf{S}}}_0^{\top }.$$

$${\mathbf{Y}}_1={\mathbf{U}}_1{\mathbf{S}}_1{\mathbf{V}}_1^{\top }.$$

$$\Vert \mathbf{F}(t,\mathbf{Y})-\mathbf{F}(t,\tilde{\mathbf{Y}})\Vert \le L\Vert \mathbf{Y}-\tilde{\mathbf{Y}}\Vert ,\Vert \mathbf{F}(t,\mathbf{Y})\Vert \le B.$$

$$\Vert (\mathbf{I}-\mathrm{P}(\mathbf{Y}))\mathbf{F}(t,\mathbf{Y})\Vert \le \epsilon$$

$$\Vert {\mathbf{Y}}_0-{\mathbf{A}}_0\Vert \le \delta .$$

$$\Vert {\mathbf{Y}}_n-\mathbf{A}(t_n)\Vert \le c_0\delta +c_1\epsilon +c_{2}h,$$

$$\mathbf{F}(t,\mathbf{Y})$$

$$\dot{\text{K}}(t)=\mathbf{F}(t,\text{K}(t){\mathbf{U}}_0^{\top }){\mathbf{U}}_0,\text{K}(t_0)={\mathbf{U}}_0{\mathbf{S}}_0.$$

$$\dot{\mathbf{S}}(t)={\mathbf{U}}_1^{\top }\mathbf{F}(t,{\mathbf{U}}_1\mathbf{S}(t){\mathbf{U}}_1^{\top }){\mathbf{U}}_1,$$

$$\mathbf{S}(t_0)={\mathbf{U}}_1^{\top }{\mathbf{Y}}_0{\mathbf{U}}_1=({\mathbf{U}}_1^{\top }{\mathbf{U}}_0){\mathbf{S}}_0({\mathbf{U}}_1^{\top }{\mathbf{U}}_0{)}^{\top }.$$

$${\mathbf{U}}_1{\mathbf{U}}_1^{\top }\mathbf{A}(t_1)=\mathbf{A}(t_1).$$

$${\mathbf{S}}_1={\mathbf{U}}_1^{\top }{\mathbf{Y}}_0{\mathbf{U}}_1+{\mathbf{U}}_1^{\top }(\mathbf{A}(t_1)-\mathbf{A}(t_0)){\mathbf{U}}_1={\mathbf{U}}_1^{\top }\mathbf{A}(t_1){\mathbf{U}}_1,$$

$${\mathbf{Y}}_1={\mathbf{U}}_1{\mathbf{S}}_1{\mathbf{U}}_1^{\top }={\mathbf{U}}_1{\mathbf{U}}_1^{\top }\mathbf{A}(t_1)({\mathbf{U}}_1{\mathbf{U}}_1^{\top })=\mathbf{A}(t_1),$$

$$\Vert {\mathbf{Y}}_n-\mathbf{A}(t_n)\Vert \le c_0\delta +c_1\epsilon +c_{2}h,$$

$$\vartheta (h,\epsilon ):=(4e^{Lh}BL+9BL)h^2+(3e^{Lh}+4)\epsilon h,$$

$$\Vert {\mathbf{U}}_1{\mathbf{U}}_1^{\top }{\mathbf{A}}_1{\mathbf{U}}_1{\mathbf{U}}_1^{\top }-{\mathbf{A}}_1\Vert \le 2\vartheta (h,\epsilon ).$$

$$\Vert {\mathbf{U}}_1\text{Z}-{\mathbf{A}}_1\Vert \le \vartheta :=\vartheta (h,\epsilon ).$$

$$\Vert {\mathbf{U}}_1\text{Z}-{\mathbf{A}}_1{\Vert }^2$$

$$=\Vert {\mathbf{U}}_1{\mathbf{U}}_1^{\top }({\mathbf{U}}_1\text{Z}-{\mathbf{A}}_1)+(\text{I}-{\mathbf{U}}_1{\mathbf{U}}_1^{\top }){\mathbf{A}}_1{\Vert }^2$$

$$=\Vert {\mathbf{U}}_1{\mathbf{U}}_1^{\top }({\mathbf{U}}_1\text{Z}-{\mathbf{A}}_1){\Vert }^2+\Vert (\text{I}-{\mathbf{U}}_1{\mathbf{U}}_1^{\top }){\mathbf{A}}_1{\Vert }^2.$$

$$\Vert {\mathbf{U}}_1{\mathbf{U}}_1^{\top }({\mathbf{U}}_1\text{Z}-{\mathbf{A}}_1){\Vert }^2+\Vert (\text{I}-{\mathbf{U}}_1{\mathbf{U}}_1^{\top }){\mathbf{A}}_1{\Vert }^2\le {\vartheta }^2.$$

$$\Vert {\mathbf{U}}_1{\mathbf{U}}_1^{\top }{\mathbf{A}}_1-{\mathbf{A}}_1\Vert \le \vartheta .$$

$$\Vert {\mathbf{U}}_1{\mathbf{U}}_1^{\top }{\mathbf{A}}_1{\mathbf{U}}_1{\mathbf{U}}_1^{\top }-{\mathbf{A}}_1\Vert$$

$$\le \Vert {\mathbf{U}}_1{\mathbf{U}}_1^{\top }({\mathbf{A}}_1{\mathbf{U}}_1{\mathbf{U}}_1^{\top }-{\mathbf{A}}_1)\Vert +\Vert {\mathbf{U}}_1{\mathbf{U}}_1^{\top }{\mathbf{A}}_1-{\mathbf{A}}_1\Vert$$

$$=\Vert {\mathbf{U}}_1{\mathbf{U}}_1^{\top }({\mathbf{U}}_1{\mathbf{U}}_1^{\top }{\mathbf{A}}_1-{\mathbf{A}}_1{)}^{\top }\Vert +\Vert {\mathbf{U}}_1{\mathbf{U}}_1^{\top }{\mathbf{A}}_1-{\mathbf{A}}_1\Vert$$

$$\le 2\Vert {\mathbf{U}}_1{\mathbf{U}}_1^{\top }{\mathbf{A}}_1-{\mathbf{A}}_1\Vert$$

$$\le 2\vartheta ,$$

$$\Vert {\mathbf{Y}}_1-{\mathbf{A}}_1\Vert \le h({\hat{c}}_1\epsilon +{\hat{c}}_{2}h),$$

Full text

1 Introduction

In this paper we propose and analyse an algorithm that computes a symmetric or skew-symmetric low-rank approximation to large symmetric or skew-symmetric time-dependent matrices that are either given explicitly or are the unknown solution to a matrix differential equation. A related algorithm is given for the approximation of symmetric or anti-symmetric time-dependent tensors by symmetric or anti-symmetric Tucker tensors of low multilinear rank.

In the matrix case, motivation for this work comes from Lyapunov and Riccati differential equations, which have large symmetric matrices as solutions, which can often be well approximated by low-rank matrices mena2018numerical . For tensors, our main motivation comes from the quantum dynamics of bosonic or fermionic systems, where the symmetric or anti-symmetric wave function is approximated by low-rank symmetric or anti-symmetric Tucker tensors in the MCTDHB and MCTDHF methods for bosons and fermions, respectively AlonSC08 ; Caillat:MCTDHF . An efficient integrator that preserves symmetry and anti-symmetry and uses them to reduce the computational complexity, is needed in these and other applications, such as using a step of the integrator as a computationally efficient retraction in optimization algorithms for (anti-)symmetric low-rank matrices and tensors.

The algorithms proposed in this paper are non-trivial modifications of the projector-splitting integrators for the dynamical low-rank approximation of matrices and Tucker tensors that were proposed in LubichOseledets and Lubich:MCTDH ; LubichVandWalach , respectively. The projector-splitting integrators have been shown to possess remarkable robustness to the typical presence of small singular values KieriLubichWalach ; LubichVandWalach , as opposed to applying standard integrators to the differential equations of dynamical low-rank approximation that are given in KochLubich07 ; KochLubich10 . However, the projector-splitting integrators do not preserve symmetry or anti-symmetry.

We will show that the (anti-)symmetry-preserving integrators proposed here retain the robustness with respect to small singular values of the projector-splitting algorithms. This relies on an exactness property, namely that time-dependent matrices and tensors of the approximation rank are reproduced exactly by the integrator. This exactness property will also be shown to be retained from the projector-splitting integrators. We note, however, that the integrators proposed here can no longer be interpreted as splitting integrators.

The new (anti-)symmetry-preserving integrators are favourable also from the computational viewpoint: compared with the projector-splitting integrator, the computational cost is halved in the (skew-)symmetric matrix case; in the case of $d$-dimensional (anti-)symmetric tensors, the computational cost for the core tensor is reduced by the factor $1/d!$, and that for the basis matrices by $1/d$.

A first attempt to modify the projector-splitting integrator and preserve the symmetry in the matrix setting, can be found in mena2018numerical : numerical examples show the correct behaviour of the approximate solution but no convergence analysis or extension to multi-dimensional arrays is provided, and no use of the symmetry is made to reduce the computational effort.

The outline of the paper is the following: in Section 2, we briefly restate the idea of dynamical low-rank approximation for matrices and we present the matrix projector-splitting integrator with some of its properties. In Section 3, we consider the case of (skew-)symmetric matrices; we present the (skew-)symmetry-preserving low-rank integrator and study its properties. In Section 4, we recapitulate the projector-splitting integrator for low-rank Tucker tensors. In Section 5, we present the integrator for (anti)-symmetric tensors of low multilinear rank and study its properties. In the final section, we present numerical experiments that illustrate the approximation properties and the robustness to small singular values.

Throughout the paper, we use the convention to denote matrices by boldface capital letters and tensors by italic capital letters.

2 General matrices: recap of the projector-splitting integrator for dynamical low-rank approximation

The objective is to approximate large time-dependent matrices ${\mathbf{A}}(t)\in{\mathbb{R}}^{m\times n}$ for $0\leq t\leq T$ by rank-$r$ matrices ${\mathbf{Y}}(t)$ with comparatively low rank $r\ll m,n$, which require much less storage than ${\mathbf{A}}(t)$ when they are available in a factorized, SVD-like form. The large, or often too large matrices ${\mathbf{A}}(t)$ may be given explicitly or they are the unknown solution to a matrix differential equation (with right-hand side function ${\mathbf{F}}:{\mathbb{R}}\times{\mathbb{R}}^{m\times n}\to{\mathbb{R}}^{m\times n}$)

[TABLE]

Dynamical low-rank approximation as presented in KochLubich07 determines ${\mathbf{Y}}(t)$ as the solution of a projected matrix differential equation, with a projection $\mathrm{P}({\mathbf{Y}})$ onto the tangent space $T_{\mathbf{Y}}\mathcal{M}_{r}$ of the manifold of rank-$r$ matrices at ${\mathbf{Y}}\in\mathcal{M}_{r}$,

[TABLE]

where ${\mathbf{Y}}_{0}$ is a rank-$r$ approximation to ${\mathbf{A}}_{0}$, typically obtained by a truncated singular value decomposition. (Here, ${\mathbf{F}}(t,{\mathbf{Y}})=\accentset{\mbox{\large.}}{{\mathbf{A}}}(t)$ if ${\mathbf{A}}(t)$ is given explicitly.) The solution ${\mathbf{Y}}(t)$ to this projected matrix differential equation then stays in the rank-$r$ manifold $\mathcal{M}_{r}$.

To make this abstract formulation practically useful, rank-$r$ matrices ${\mathbf{Y}}(t)$ are written (non-uniquely) in factored form

[TABLE]

where the slim matrices ${\mathbf{U}}(t)\in{\mathbb{R}}^{m\times r}$ and ${\mathbf{V}}(t)\in{\mathbb{R}}^{n\times r}$ each have $r$ orthonormal columns, and the small square matrix ${\mathbf{S}}(t)\in{\mathbb{R}}^{r\times r}$ is invertible. We choose the tangent space projection $\mathrm{P}({\mathbf{Y}})$ as the orthogonal projection onto $T_{\mathbf{Y}}(\mathcal{M}_{r})$ with respect to the Euclidean or Frobenius inner product $\langle{\mathbf{A}},{\mathbf{B}}\rangle=\textbf{vec}({\mathbf{A}})^{\top}\textbf{vec}({\mathbf{B}})$, where $\textbf{vec}({\mathbf{A}})\in{\mathbb{R}}^{mn}$ is a vectorization of ${\mathbf{A}}$. Then, $\mathrm{P}({\mathbf{Y}})$ is given as an alternating sum of three subprojections KochLubich07 ,

[TABLE]

The projector-splitting integrator of LubichOseledets is a Lie–Trotter or Strang splitting method that splits the right-hand side of (2) according to the three terms in (4). It turned out that such a splitting combines very well with the factorization (3). In the first substep of a Lie–Trotter splitting, $\textbf{K}:={\mathbf{U}}{\mathbf{S}}$ is updated, in the second substep ${\mathbf{S}}$ is updated, and in the third substep $\textbf{L}:={\mathbf{V}}{\mathbf{S}}^{\top}$. The algorithm alternates between the numerical solution of matrix differential equations (of dimensions $m\times r$, $r\times r$, $n\times r$) and orthogonal decompositions of slim matrices (of dimensions $m\times r$ and $n\times r$). One time step of integration from time $t_{0}$ to $t_{1}=t_{0}+h$ starting from a factored rank-$r$ matrix ${\mathbf{Y}}_{0}={\mathbf{U}}_{0}{\mathbf{S}}_{0}{\mathbf{V}}_{0}^{\top}$ proceeds as follows:

K-step : Update ${\mathbf{U}}_{0}\rightarrow{\mathbf{U}}_{1},{\mathbf{S}}_{0}\rightarrow\hat{{\mathbf{S}}}_{1}$

Integrate from $t=t_{0}$ to $t_{1}$ the $m\times r$ matrix differential equation

[TABLE]

Perform a QR factorization $\textbf{K}(t_{1})={\mathbf{U}}_{1}\hat{{\mathbf{S}}}_{1}$. 2. 2.

S-step : Update $\hat{{\mathbf{S}}}_{1}\rightarrow\tilde{{\mathbf{S}}}_{0}$

Integrate from $t=t_{0}$ to $t_{1}$ the $r\times r$ matrix differential equation

[TABLE]

and set $\tilde{{\mathbf{S}}}_{0}={\mathbf{S}}(t_{1})$. 3. 3.

L-step : Update ${\mathbf{V}}_{0}\rightarrow{\mathbf{V}}_{1},\tilde{{\mathbf{S}}}_{0}\rightarrow{\mathbf{S}}_{1}$

Integrate from $t=t_{0}$ to $t_{1}$ the $n\times r$ matrix differential equation

[TABLE]

Perform a QR factorization $\textbf{L}(t_{1})={\mathbf{V}}_{1}{\mathbf{S}}_{1}^{\top}$.

Then, the approximation after one time step is given by

[TABLE]

To proceed, we iterate the procedure taking ${\mathbf{Y}}_{1}$ as starting point for the next step.

The above algorithm describes the first-order Lie–Trotter splitting. The algorithm for the second-order Strang splitting is obtained by concatenating the above algorithm with the same algorithm in reversed order, each for half the step-size; see LubichOseledets for the detailed description.

The projector-splitting integrator has remarkable properties. First, it reproduces rank-$r$ matrices without error.

Theorem 2.1 (Exactness property, (LubichOseledets, , Theorem 4.1))

Let ${\mathbf{A}}(t)\in\mathbb{R}^{m\times n}$ be of rank $r$ for $t_{0}\leq t\leq t_{1}$, so that ${\mathbf{A}}(t)$ has a factorization (3), ${\mathbf{A}}(t)={\mathbf{U}}(t){\mathbf{S}}(t){\mathbf{V}}(t)^{\top}$. Moreover, assume that the $r\times r$ matrix ${\mathbf{V}}(t_{1})^{\top}{\mathbf{V}}(t_{0})$ is invertible. With ${\mathbf{Y}}_{0}={\mathbf{A}}(t_{0})$, the projector-splitting integrator for $\accentset{\mbox{\large.}}{{\mathbf{Y}}}(t)=\mathrm{P}({\mathbf{Y}}(t))\accentset{\mbox{\large.}}{{\mathbf{A}}}(t)$ is then exact: ${\mathbf{Y}}_{1}={\mathbf{A}}(t_{1})$.

The second remarkable property is the robustness of the algorithm to the presence of small singular values of the solution or its approximation. This is in contrast to standard integrators applied to (2) or the equivalent differential equations for the factors ${\mathbf{U}}(t)$, ${\mathbf{S}}(t)$, ${\mathbf{V}}(t)$, which contain a factor ${\mathbf{S}}(t)^{-1}$ on the right-hand sides (KochLubich07, , Prop. 2.1). Moreover, the local Lipschitz constant of the tangent space projection $\mathrm{P}(\cdot)$ is proportional to the inverse of the smallest nonzero singular value (KochLubich07, , Lemma 4.2). The appearance of small singular values is typical in applications, because the smallest singular value retained in the approximation cannot be expected to be much larger than the largest discarded singular value of the solution, which needs to be small to obtain good accuracy of the low-rank approximation.

Theorem 2.2 (Robust error bound, (KieriLubichWalach, , Theorem 2.1))

Let ${\mathbf{A}}(t)$ denote the solution of the matrix differential equation (1). Assume that the following conditions hold in the Frobenius norm $\|\cdot\|=\|\cdot\|_{F}$:

${\mathbf{F}}$* is Lipschitz-continuous and bounded: for all ${\mathbf{Y}},\widetilde{{\mathbf{Y}}}\in\mathbb{R}^{m\times n}$ and $0\leq t\leq T$,*

[TABLE] 2. 2.

The non-tangential part of ${\mathbf{F}}(t,{\mathbf{Y}})$ is $\varepsilon$-small:

[TABLE]

for all ${\mathbf{Y}}\in\mathcal{M}$ in a neighbourhood of ${\mathbf{A}}(t)$ and $0\leq t\leq T$. 3. 3.

The error in the initial value is $\delta$-small:

[TABLE]

Let ${\mathbf{Y}}_{n}$ denote the rank-$r$ approximation to ${\mathbf{A}}(t_{n})$ at $t_{n}=nh$ obtained after n steps of the projector-splitting integrator with step-size $h>0$. Then, the error satisfies for all $n$ with $t_{n}=nh\leq T$

[TABLE]

where the constants $c_{i}$ only depend on $L,B,$ and $T$. In particular, the constants are independent of singular values of the exact or approximate solution.

It is further shown in (KieriLubichWalach, , Section 2.6.3) that an inexact solution of the matrix differential equations in the projector-splitting integrator leads to an additional error that is bounded in terms of the local errors in the inexact substeps, again with constants that do not depend on small singular values.

Numerical experiments with the matrix projector-splitting integrator and comparisons with standard numerical integrators are reported in LubichOseledets ; KieriLubichWalach . These experiments show good behaviour also for spatially discretized partial differential equations where the Lipschitz constant becomes large, a case that as of now is not covered by the theory.

3 Symmetric and skew-symmetric matrices: a structure-preserving integrator for dynamical low-rank approximation

We now assume that the right-hand side function in (1) is such that

[TABLE]

This condition ensures that the solutions to the matrix differential equation (1) and the projected differential equation (2) are (skew-)symmetric provided the initial values are (skew-)symmetric. For (2), this is seen from formula (4) for the tangent space projection with equal left and right factors ${\mathbf{V}}={\mathbf{U}}$ in the decomposition (3) of the (skew-)symmetric rank-$r$ matrix ${\mathbf{Y}}={\mathbf{U}}{\mathbf{S}}{\mathbf{U}}^{\top}$.

While the projector-splitting integrator for dynamical low-rank approximation described in the previous section has favourable properties, it does not preserve symmetry or skew-symmetry of the solution ${\mathbf{A}}(t)$ to (1).

3.1 (Skew-)symmetry preserving integrator

We now propose a modified integrator that preserves symmetry and skew-symmetry and still retains the exactness and robustness properties of the projector-splitting integrator. A step with this integrator consists of two substeps. The first substep is identical to the first substep (K-step) of the projector-splitting integrator: it updates $\textbf{K}={\mathbf{U}}{\mathbf{S}}$ in the decomposition ${\mathbf{Y}}={\mathbf{U}}{\mathbf{S}}{\mathbf{U}}^{\top}$. The second substep is a substantially modified update of ${\mathbf{S}}$, which can be viewed as a Galerkin approximation in the basis provided by the first substep.

Given ${\mathbf{Y}}_{0}={\mathbf{U}}_{0}{\mathbf{S}}_{0}{\mathbf{U}}_{0}^{\top}$ with a (skew-)symmetric $r\times r$-matrix ${\mathbf{S}}_{0}$ at time $t_{0}$, we compute the factorization ${\mathbf{Y}}_{1}={\mathbf{U}}_{1}{\mathbf{S}}_{1}{\mathbf{U}}_{1}^{\top}$ with a (skew-)symmetric $r\times r$-matrix ${\mathbf{S}}_{1}$ at time $t_{1}=t_{0}+h$ by the following algorithm:

To continue in time, we take ${\mathbf{Y}}_{1}$ as starting value for the next step and perform another step of the integrator.

Note that in this integrator the factor R in the $QR$-decomposition of the first substep is not reused in the second substep, in contrast to the projector-splitting integrator. The computational cost is approximately halved, since the L-step is not needed here.

We will now show that the (skew-)symmetric integrator retains the exactness and robustness properties of the projector-splitting integrator, using these known results in the proof.

3.2 Exactness property of the (skew-)symmetric integrator

The exactness result Theorem 2.1 extends in the following way.

Theorem 3.1 (Exactness property)

Let ${\mathbf{A}}(t)\in\mathbb{R}^{n\times n}$ be (skew-)symmetric and of rank $r$ for $t_{0}\leq t\leq t_{1}$, so that ${\mathbf{A}}(t)$ has a factorization (3) with equal left and right factors, ${\mathbf{A}}(t)={\mathbf{U}}(t){\mathbf{S}}(t){\mathbf{U}}(t)^{\top}$. Moreover, assume that the $r\times r$ matrix ${\mathbf{U}}(t_{1})^{\top}{\mathbf{U}}(t_{0})$ is invertible. With ${\mathbf{Y}}_{0}={\mathbf{A}}(t_{0})$, the (skew-)symmetric integrator for $\accentset{\mbox{\large.}}{{\mathbf{Y}}}(t)=\mathrm{P}({\mathbf{Y}}(t))\accentset{\mbox{\large.}}{{\mathbf{A}}}(t)$ is then exact: ${\mathbf{Y}}_{1}={\mathbf{A}}(t_{1})$.

Proof

We note that the projector-splitting integrator and the (skew-)symmetric integrator have the same first step. Let ${\mathbf{U}}_{1}\in\mathbb{R}^{n\times r}$ be the matrix with orthonormal columns computed in the first substep. Due to the exactness of the matrix projector-splitting integrator as given by Theorem 2.1 we know that ${\mathbf{U}}_{1}$ and ${\mathbf{A}}(t_{1})$ have the same range and therefore

[TABLE]

Denoting $\Delta{\mathbf{A}}={\mathbf{A}}(t_{1})-{\mathbf{A}}(t_{0})$, the (skew-)symmetric integrator provides for the second substep the solution

[TABLE]

since ${\mathbf{Y}}_{0}={\mathbf{A}}(t_{0})$. The result after a time step of the (skew-)symmetric integrator is

[TABLE]

where the last equality holds because of (6) and the (skew-)symmetry of ${\mathbf{A}}(t_{1})$. ∎

3.3 Robustness to small singular values

The error bound of Theorem 2.2 extends in the following way.

Theorem 3.2 (Robust error bound)

Let ${\mathbf{A}}(t)$ denote the (skew-)symmetric solution of the matrix differential equation (1) with ${\mathbf{F}}$ satisfying (5). Assume that conditions $1$.-$\,3$. of Theorem 2.2 are fulfilled.

Let ${\mathbf{Y}}_{n}$ denote the rank-$r$ approximation to ${\mathbf{A}}(t_{n})$ at $t_{n}=nh$ obtained after n steps of the (skew-)symmetric integrator of Algorithm 1 with step-size $h>0$. Then, the error satisfies for all $n$ with $t_{n}=nh\leq T$

[TABLE]

where the constants $c_{i}$ only depend on $L,B,$ and $T$. In particular, the constants are independent of singular values of the exact or approximate solution.

As in (KieriLubichWalach, , Section 2.6.3), it can be further shown that an inexact solution of the matrix differential equations in the projector-splitting integrator leads to an additional error that is bounded in terms of the local errors in the inexact substeps, again with constants that do not depend on small singular values.

Remark 1

The method of Algorithm 1 is of order 1, and higher order can be obtained simply by composition as, e.g., in (HLW, , Section II.4). However, like for the projector-splitting integrator of LubichOseledets , it is not known if an error bound of higher order in the step-size $h$ can be obtained with constants that are independent of small singular values. Numerical experiments with the Strang version of the projector-splitting integrator, which is of order 2, indicate an order reduction in some examples with very small singular values ostermann2018convergence .

We now prepare for the proof of Theorem 3.2, which views the (skew-)symmetric integrator as a perturbed variant of the projector-splitting integrator.

Let us introduce the quantity

[TABLE]

which represents the local error bound after one time step of the projector-splitting integrator, as proved in (KieriLubichWalach, , Theorem 2.1).

In the following, we denote by ${\mathbf{U}}_{1}\in{\mathbb{R}}^{n\times r}$ the matrix with orthonormal columns obtained in the first substep of the integrator. We recall that the matrix projector-splitting and the (skew-)symmetric integrator have the first substep in common.

We denote by ${\mathbf{A}}_{1}$ the (skew-)symmetric solution at time $t_{1}$ of the full problem (1), where we consider the initial data to coincide with the (skew-)symmetric rank-$r$ matrix ${\mathbf{Y}}_{0}$. For the local error analysis, the following lemma is needed.

Lemma 1

Let ${\mathbf{U}}_{1},{\mathbf{A}}_{1}$ be defined as above. The following estimate holds:

[TABLE]

Proof

The local error analysis in KieriLubichWalach shows that the $r\times n$ matrix ${\mathbf{Z}}={\mathbf{S}}_{1}^{\mathrm{ps}}{\mathbf{V}}_{1}^{\mathrm{ps},\top}$, where ${\mathbf{S}}_{1}^{\mathrm{ps}}$ and ${\mathbf{V}}_{1}^{\mathrm{ps}}$ are the matrices computed in the third substep of the projector-splitting algorithm, satisfies

[TABLE]

The square of the left-hand side can be split into two terms:

[TABLE]

Hence,

[TABLE]

From the second term it follows that

[TABLE]

By the (skew-)symmetry of ${\mathbf{A}}_{1}$, this implies

[TABLE]

which yields the result. ∎

In the following lemma, we show that the approximation given after one time step is $O(h(h+\varepsilon))$ close to the solution of system (1) when the starting values coincide.

Lemma 2 (Local Error)

The following local error bound holds:

[TABLE]

where the constants only depend on $L$ and $B$ and a bound of the step size. In particular, the constants are independent of singular values of the exact or approximate solution.

Proof

By the identity ${\mathbf{Y}}_{1}={\mathbf{U}}_{1}{\mathbf{S}}_{1}{\mathbf{U}}_{1}^{\top}$ and Lemma 1 we have that

[TABLE]

The analysis of the local error thus reduces to estimating $\|{\mathbf{S}}_{1}-{\mathbf{U}}_{1}^{\top}{\mathbf{A}}_{1}{{\mathbf{U}}}_{1}\|$. To this end, we introduce the following quantity: for $t_{0}\leq t\leq t_{1}$,

[TABLE]

We observe that

[TABLE]

where $\textbf{R}(t)$ is defined as the term in big brackets. From Lemma 1 and from the bound $B$ of ${\mathbf{F}}$, which yields for $t_{0}\leq t\leq t_{1}$

[TABLE]

we conclude that the remainder term is bounded by

[TABLE]

This yields that ${\mathbf{F}}(t,{\mathbf{A}}(t))$ can be written as

[TABLE]

where the defect

[TABLE]

is bounded via the Lipschitz continuity of ${\mathbf{F}}$ as

[TABLE]

We now compare the two differential equations

[TABLE]

By construction, the solution of the first differential equation at time $t_{1}$ is $\widetilde{\mathbf{S}}(t_{1})={\mathbf{U}}_{1}^{\top}{\mathbf{A}}_{1}{{\mathbf{U}}}_{1}$. The solution of the second differential equation is ${\mathbf{S}}_{1}$ as given by the second substep of the (skew-)symmetric integrator. We now apply the Gronwall inequality to the previous system and obtain

[TABLE]

The result now follows using the definition of $\vartheta$. ∎

Thanks to the Lipschitz continuity of the function ${\mathbf{F}}$, we conclude the proof of Theorem 3.2 from the local to the global errors by the standard argument of Lady Windermere’s fan (HairerNorsettWanner:ODE_BOOK1, , Section II.3).

4 General tensors: recap of the projector-splitting integrator for the dynamical low-rank approximation by Tucker tensors

The objective is to approximate time-dependent tensors111In view of the applications in quantum dynamics, we here consider tensors with complex entries. $A(t)\in{\mathbb{C}}^{n_{1}\times\dots\times n_{d}}$ for $0\leq t\leq T$ by tensors $Y(t)$ of multilinear rank ${\mathbf{r}}=(r_{1},\dots,r_{d})$, with $r_{i}\ll n_{i}$. (We recall that $r_{i}$ is the rank of the $i$th matricization $\mathbf{Mat}_{i}(Y)\in{\mathbb{C}}^{n_{i}\times n_{i}^{\prime}}$ with $n_{i}^{\prime}=\prod_{j\neq i}n_{j}$, which aligns all entries of $Y$ with $i$th index $k$ in the $k$th row. The retensorization is denoted by Ten${}_{i}(\cdot)$, such that Ten${}_{i}(\mathbf{Mat}_{i}(Y))=Y$.)

The tensors $A(t)$ may be given explicitly or they are the unknown solution to a tensor differential equation (with right-hand side function $F:{\mathbb{R}}\times{\mathbb{C}}^{n_{1}\times\dots\times n_{d}}\to{\mathbb{C}}^{n_{1}\times\dots\times n_{d}}$)

[TABLE]

Dynamical low-rank approximation as presented in KochLubich10 determines $Y(t)$ as the solution of the projected matrix differential equation, with a projection $\mathrm{P}(Y)$ onto the tangent space $T_{Y}\mathcal{M}_{\mathbf{r}}$ of the manifold of tensors of multilinear rank ${\mathbf{r}}$ at $Y\in\mathcal{M}_{\mathbf{r}}$,

[TABLE]

where $Y_{0}$ is a rank-${\mathbf{r}}$ approximation to $A_{0}$. (Here, $F(t,Y)=\accentset{\mbox{\large.}}{A}(t)$ if $A(t)$ is given explicitly.) Tensors $Y(t)$ of multilinear rank ${\mathbf{r}}$ are represented non-uniquely in the Tucker form DeLauthawer:HOSVD (using here the multilinear notation of KoldaBader:TensorDec )

[TABLE]

where the core tensor $C(t)\in\mathbb{C}^{r_{1}\times\dots\times r_{d}}$ is of full multi-linear rank and the basis matrices ${\mathbf{U}}_{i}\in\mathbb{C}^{n\times r_{i}}$ have orthonormal columns. We choose the tangent space projection $\mathrm{P}(Y)$ as the orthogonal projection onto $T_{Y}(\mathcal{M}_{r})$ with respect to the Euclidean inner product $\langle A,B\rangle=\textbf{vec}(A)^{*}\textbf{vec}(B)$, where $\textbf{vec}(A)$ is a vectorization of $A$. Then, $\mathrm{P}({\mathbf{Y}})$ is given as an alternating sum of $2d-1$ subprojections Lubich:MCTDH , and like in the matrix case, a projector-splitting integrator with favourable properties can be formulated and efficiently implemented. The matrix projector-splitting integrator proposed in Section 2.1 has been successfully extended to the Tucker tensor format in different algorithmic versions in Lubich:MCTDH and LubichVandWalach . It is shown in (LubichVandWalach, , Section 6) that the proposed Tucker integrators are mathematically equivalent. The algorithm runs through the modes $i=1,\dots,d$ and solves differential equations for matrices of the dimension of the slim basis matrices and for the core tensor in alternation with orthogonal decompositions of slim matrices. We refer also to KlossBL17 ; BonfantiB18 for the formulation and implementation of this algorithm in the context of the MCTDH method MeyerGW09 of molecular quantum dynamics in the chemical physics literature.

Moreover, the Tucker integrator has been proved in LubichVandWalach to satisfy analogous properties to the matrix projector-splitting integrator: the exactness property and the robust convergence in the presence of small singular values of matricizations of the core tensor. We refer to (LubichVandWalach, , Theorems 4.1 and 5.1) for the precise formulation, which is very similar to the matrix case.

5 Symmetric and anti-symmetric tensors: a structure-preserving integrator for dynamical low-rank approximation

A tensor $A=(a_{i_{1},\dots,i_{d}})\in{\mathbb{C}}^{n\times\dots\times n}$ is symmetric if for every permutation $\sigma\in S(d)$,

[TABLE]

and $A$ is anti-symmetric if for every permutation $\sigma\in S(d)$,

[TABLE]

It follows from DeLauthawer:HOSVD and Hackbusch:SymTenRepr that a symmetric/anti-symmetric tensor $Y\in{\mathbb{C}}^{n\times\dots\times n}$ of multi-linear rank $\textbf{r}=(r,\dots,r)$ admits a Tucker decomposition

[TABLE]

where the core tensor $C\in{\mathbb{C}}^{r\times\dots\times r}$ is symmetric/anti-symmetric of full rank ${\mathbf{r}}$ and the basis matrix ${\mathbf{U}}\in{\mathbb{C}}^{n\times r}$ is the same for all indices.

We assume that the right-hand side function in (7) is such that

[TABLE]

Like (5) in the matrix case, this ensures that the solutions to the tensor differential equation (7) and the projected differential equation (8) are (anti-)symmetric provided the initial tensors are (anti-)symmetric. As we noted already in the matrix case, the projector-splitting Tucker integrator does not preserve (anti-)symmetry.

5.1 (Anti-)symmetry preserving Tucker integrator

The numerical integrator defined in Section 3 for the matrix case extends in a natural way to the Tucker tensor format. The first substep, which updates the basis matrix ${\mathbf{U}}$, is identical to the first substep of the general Tucker integrator in LubichVandWalach ; Lubich:MCTDH . The second substep is a Galerkin method with the updated basis and determines the updated (anti-)symmetric core tensor.

Given the (anti-)symmetric tensor $Y_{0}=C_{0}\operatorname{{\hbox{\large\sf X}}}_{i=1}^{d}{\mathbf{U}}_{0}$, we compute the (anti-)symmetric approximation $Y_{1}=C_{1}\operatorname{{\hbox{\large\sf X}}}_{i=1}^{d}{\mathbf{U}}_{1}$ at time $t_{1}=t_{0}+h$ as follows:

To continue, we take $Y_{1}$ as the starting value for the next step.

5.2 Exactness property of the (anti-)symmetric Tucker integrator

The following result extends the exactness results of Theorem 3.1 and (LubichVandWalach, , Theorem 4.1) to (anti-)symmetric tensors.

Theorem 5.1 (Exactness property)

Let $A(t)\in{\mathbb{C}}^{n\times\dots\times n}$ be (anti-)symmetric and of multilinear rank $(r,\dots,r)$ for $t_{0}\leq t\leq t_{1}$, so that $A(t)=C(t)\operatorname{{\hbox{\large\sf X}}}_{i=1}^{d}{\mathbf{U}}(t)$, where the $n\times r$ basis matrix ${\mathbf{U}}$ has orthonormal columns. Moreover, assume that the $r\times r$ matrix ${\mathbf{U}}(t_{1})^{*}{\mathbf{U}}(t_{0})$ is invertible. With $Y_{0}=A(t_{0})$, the (anti-)symmetric Tucker integrator for $\accentset{\mbox{\large.}}{Y}(t)=\mathrm{P}(Y(t))\accentset{\mbox{\large.}}{A}(t)$ is then exact: $Y_{1}=A(t_{1})$.

Proof

The projector-splitting Tucker integrator and the (anti-)symmetric integrator have the same first substep. Let ${\mathbf{U}}_{1}\in\mathbb{R}^{n\times r}$ be the basis matrix with orthonormal columns computed in the first substep. Due to the exactness of the projector-splitting Tucker integrator as shown by (LubichVandWalach, , Theorem 4.1) we have that $A(t_{1})$ has the (anti)-symmetric Tucker representation

[TABLE]

for some (anti-)symmetric core tensor $\widehat{C}_{1}\in{\mathbb{C}}^{r\times\dots\times r}$. Using the rule $A\times_{i}{\mathbf{V}}\times_{i}\textbf{W}=A\times_{i}(\textbf{WV})$, this implies that

[TABLE]

With $Y_{0}=A(t_{0})$ we obtain from the second substep of the algorithm

[TABLE]

which proves the exactness. ∎

5.3 Robustness to small singular values

The robust error bounds from Theorem 3.2 and (LubichVandWalach, , Theorem 5.1) extend to the (anti)symmetric Tucker integrator as follows. The norm $\|B\|$ of a tensor $B$ used here is the Euclidean norm of the entries of $B$.

Theorem 5.2 (Robust error bound)

Let $A(t)$ denote the (anti-)symmetric solution of the tensor differential equation (7) with $F$ satisfying (10). Assume the following:

$F$* is Lipschitz-continuous and bounded.* 2. 2.

The non-tangential part of $F(t,Y)$ is $\varepsilon$-small:

[TABLE]

for all $Y$ of multilinear rank $(r,\dots,r)$ in a neighbourhood of $A(t)$ and $0\leq t\leq T$. 3. 3.

The error in the initial value is $\delta$-small:

[TABLE]

Let $Y_{n}$ denote the (anti-)symmetric approximation of multinear rank $(r,\dots,r)$ to $A(t_{n})$ at $t_{n}=nh$ obtained after n steps of the (anti-)symmetric Tucker integrator with step-size $h>0$. Then, the error satisfies for all $n$ with $t_{n}=nh\leq T$

[TABLE]

where the constants $c_{i}$ only depend on the Lipschitz constant $L$ and bound $B$ of $F$, on $T$, and on the dimension $d$. In particular, the constants are independent of singular values of matricizations of the exact or approximate solution.

It can be further shown that an inexact solution of the matrix differential equations in the integrator leads to an additional error that is bounded in terms of the local errors in the inexact substeps, again with constants that do not depend on small singular values.

The proof of Theorem 5.2 proceeds similar to the proof of Theorem 3.2 for the (skew)-symmetric matrix case. We begin with a key lemma and then analyze the local error produced after one time step, comparing the numerical solution with the exact solution that starts from the same initial value $A_{0}=Y_{0}$. We denote the value of this solution at $t_{1}$ by $A_{1}$. The basis matrix computed in the first substep of the integrator is denoted by ${\mathbf{U}}_{1}$.

Lemma 3

The following estimate holds:

[TABLE]

where $c$ only depends on $d$ and a bound for $hL$.

Proof

The error bound of (LubichVandWalach, , Theorem 5.1) shows that there exists $Z\in\mathbb{C}^{r\times n\times\dots\times n}$ such that

[TABLE]

We observe that

[TABLE]

As in the matrix case we obtain

[TABLE]

Thanks to (anti-)symmetry we have

[TABLE]

which yields

[TABLE]

To conclude, we observe

[TABLE]

and the result follows by an iteration of this argument. ∎

We are now in the position to analyse the local error produced after one time step of the integrator.

Lemma 4 (Local error)

The error of the (anti-)symmetric Tucker integrator after one time step satisfies

[TABLE]

where $\hat{c}$ only depends on $d$ and a bound of $hL$. In particular, the constant is independent of singular values of the exact or approximate solution.

Proof

By construction of the algorithm and by Lemma 3 we have

[TABLE]

The problem reduces to estimating $\|C_{1}-A_{1}\operatorname{{\hbox{\large\sf X}}}_{i=1}^{d}{\mathbf{U}}_{1}^{*}\|$. We introduce the tensor

[TABLE]

which satisfies

[TABLE]

In the same way as in the proof of Lemma 2 (replacing ${\mathbf{S}}$ by $C$) this is compared with the differential equation for $C(t)$ in the second substep of the integrator. This yields the stated result. ∎

Using the Lipschitz continuity of the function $F$, we conclude the proof of Theorem 5.2 from the local to the global errors by the standard argument of Lady Windermere’s fan (HairerNorsettWanner:ODE_BOOK1, , Section II.3).

6 Numerical Experiments

In this section we show results of various numerical experiments. The computations were done using Matlab R2017a software with Tensor Toolbox package v2.6 TTB_Software and TensorLab package v3.0 vervliet2016tensorlab .

6.1 Addition of symmetric tensors: a computationally inexpensive retraction

Let $A\in\mathbb{C}^{n\times\dots\times n}$ be a symmetric tensor of multi-linear rank $\textbf{r}=(r,\dots,r)$ and let $B\in\mathbb{C}^{n\times\dots\times n}$. We consider the addition of two given tensors,

[TABLE]

where $C\in\mathbb{C}^{n\times\dots\times n}$ is not necessarily of low rank and we want to compute a symmetric rank-$(r,\dots,r)$ approximation. Such a retraction is typically required in optimization problems on low-rank manifolds and needs to be computed in each iterative step of a descent algorithm. The approach considered here consists of reformulating the addition problem as the solution of the following differential equation at time $t=1$:

[TABLE]

We will compare the solution obtained by computing the full addition and retracting to the manifold of symmetric tensors of multilinear rank $(r,\dots,r)$ with the one obtained from the application of the symmetric low-rank tensor integrator. The advantage of the latter method is that the approximation is built inside the manifold, so that no truncation to rank $r$ is needed.

For our numerical example, we initialize A as a symmetric random Tucker tensor of size $100\times 100\times 100$ and multi-linear rank $\textbf{r}=(10,10,10)$; we take B as an element in the tangent space of the symmetric rank-r tensor manifold at $A$. We compare the dynamical low rank approximation $Y_{1}$ generated by the algorithm introduced in Section 5 with a low rank symmetric retraction de2000best ; regalia2013monotonically of the full solution denoted by $X$. For the last part, we use the built-in tucker_als and tucker_sym functions of the Tensor Toolbox Package.

We observe that the approximation $Y_{1}$ shows the correct behavior, at reduced computational cost. Decreasing the norm of the tensor $B$ decreases the approximation error as expected, proportional to $\|B\|^{2}$.

6.2 Robustness with respect to small singular values

We present two numerical examples and show robustness of the proposed symmetric integrator in the presence of small singular values. For the sake of presentation we consider the matrix case. Analogous examples can be implemented for Tucker tensors, and similar results are obtained.

In the first example, the time-dependent matrix is given explicitly as

[TABLE]

The matrix $\textbf{D}\in{\mathbb{R}}^{N\times N}$ is diagonal with elements $d_{j}=2^{-j}$ and the matrix $\textbf{W}\in{\mathbb{R}}^{N\times N}$ is skew-symmetric and randomly generated. We choose $N=100$ and final time $T=1$. We compare the symmetric low-rank integrator presented in Section 3 with a numerical solution obtained with a 4-th order explicit Runge-Kutta method applied to the system of differential equations for dynamical low-rank approximation as derived in KochLubich07 .

The numerical results for different ranks are shown in Figure 1. In contrast to the Runge–Kutta method, the proposed symmetric low-rank integrator does not suffer of a step-size restriction in the presence of small singular values.

In the second example, we integrate the Lyapunov matrix differential equation (cf. mena2018numerical )

[TABLE]

Here, we choose $\textbf{A}=\texttt{tridiag}(-1,2,-1)\otimes\textbf{I}+\textbf{I}\otimes\texttt{tridiag}(-1,2,-1)\in{\mathbb{R}}^{N\times N}$ as a discrete Laplacian. The positive definite matrix $\textbf{Q}\in{\mathbb{R}}^{N\times N}$ has rank 5 and is randomly generated. The orthonormal matrix $\textbf{U}_{0}\in{\mathbb{R}}^{N\times N}$ is randomly generated and $\textbf{S}_{0}\in{\mathbb{R}}^{N\times N}$ is of rank 1 with only one non-zero element, $s_{11}=1$.

The reference solution and the linear subproblems appearing in the definition of the symmetric low-rank integrator have been solved with the Matlab solver ode45 and stringent tolerance parameters {’RelTol’, 1e-10, ’AbsTol’, 1e-14} . We choose $N=100$ and final time $T=0.1$. The singular values of the reference solution and the absolute errors $\|Y_{n}-A(t_{n})\|_{F}$ at time $t_{n}=T$ of the approximate solutions for different ranks calculated with different step-sizes are shown in Figure 2.

6.3 Ground state of a fermionic multi-particle system

A natural field of application of the (anti-)symmetric low-rank algorithms is in quantum dynamics of systems of fermions or bosons, which is described by anti-symmetric or symmetric multivariate wave functions, respectively. In the MCTDHF and MCTDHB methods AlonSC08 ; Caillat:MCTDHF , the approximate wave function is sought for in the form of a low-rank Tucker tensor that is anti-symmetric and symmetric, respectively. It approximates the huge tensor of coefficients with respect to some fixed spatial basis, such as a Fourier basis. The (anti-)symmetric low-rank time integrator proposed in this paper appears ideally suited as a numerical integrator for such problems.

As a first illustration of the approach, we apply the method for the calculation of the ground state of a system of $d$ fermions in 1 space dimension. We calculate the ground state as the solution for large time of the imaginary-time Schrödinger equation

[TABLE]

We consider the Hamiltonian given by

[TABLE]

where $x_{l}\in{\mathbb{R}}$ represents the position of the $l$-th particle and we choose the torsion potential

[TABLE]

Choosing a collocation method with a tensor Fourier basis set (with $K$ basis functions per particle) for approximating the anti-symmetric wave function leads to a huge tensor differential equation (7), where a low-rank approximation to the anti-symmetric $d$-dimensional tensor $Y(t)\in{\mathbb{C}}^{K\times\dots\times K}$ is to be computed. This is done with a variational splitting method. The stiffness introduced by the Laplacian will be handled with a split-step Fourier method lubich2008quantum while the two-particle interaction is treated with the anti-symmetric low-rank integrator of Section 5.

We introduce the space discretization

[TABLE]

Let $Y(t)$ be a time-dependent tensor defined element-wise by

[TABLE]

The fermionic property of the system implies that the tensor $Y(t)$ is anti-symmetric. Denoting $\mathcal{F}_{K}$ the Fourier matrix, we define

[TABLE]

The Fourier collocation space discretization of (11) is equivalent to the system

[TABLE]

Using the trigonometric equality $\cos(x-y)=\cos(x)\cos(y)+\sin(x)\sin(y)$, the linear operator $H$ can be written in a multi-linear product form as

[TABLE]

In order to remove the stiffness introduced by the Laplacian, we split (12) in $(d+2)$ sub-problems. The solution of the first sub-problem at time $t_{1}=t_{0}+h$ is obtained updating the core tensor $C_{0}$ in the initial data,

[TABLE]

Afterwards, we start considering the equation

[TABLE]

We matricize in the first mode,

[TABLE]

with initial data,

[TABLE]

The solution can now be computed with the 1-dimensional split-step Fourier method. Denoting

[TABLE]

and tensorizing back in the first mode we have that

[TABLE]

Taking this as initial condition and iterating the same process for all the successive modes we obtain the updated anti-symmetric tensor

[TABLE]

We now apply the anti-symmetric low-rank integrator to the multi-particle interaction,

[TABLE]

where

[TABLE]

The core $C_{0}$ is renormalized after each step, since the absolute size of the tensor is irrelevant. We emphasize the fact that all along the implementation, it is crucial to use the structure of the Tucker tensor and avoid to build the huge matrix $\textbf{V}_{0}$ appearing in the definition of the integrator. The K and $C$ problems are linear and can be solved with few iterations of the Arnoldi process.

In our numerical experiment we choose $d=3$ particles in 1 space dimension and fix the number of Fourier basis functions per particle at $K=128$, the step-size at $h=0.01$ and we propagate the system until $T=40$.

We introduce the discrete energy

[TABLE]

Although the integrator preserves the anti-symmetry in theory, in a straightforward implementation round-off errors will destroy the anti-symmetry and take the system to the lowest state of energy: the bosonic ground state - the one achieved starting from a symmetric initial value. This behavior can be corrected in the integrator by enforcing the anti-symmetry of the small core tensor (which is violated only by round-off errors) at each step or ever after a few steps. In this way the computation tends to the fermionic ground state. The energy levels generated by the approximation of rank 5 are shown in Figure 3 for the bosonic system and the fermionic system with and without enforced anti-symmetrization.

6.4 MCTDHF - Ultra fast laser dynamics

In the second example we consider the situation where an external pulsing laser field is introduced in the system. We refer to Caillat:MCTDHF for the physical description of the problem and its MCTDHF formulation. We consider the time-dependent Schrödinger equation

[TABLE]

with the Hamiltonian

[TABLE]

with the torsion potential $V$ as before and with parameters

[TABLE]

As initial value, we choose the ground-state calculated at the previous step. We fix the number of Fourier basis functions per particle at $K=128$, the step-size at $h=0.005$ and we propagate the system until $T=1$ by the same algorithm as in the previous subsection, but this time for the real-time evolution instead of the imaginary-time evolution.

The time evolution of the energy obtained by the approximation of rank 5 is shown in Figure 4.

Acknowledgements.

We thank Balázs Kovács and Hanna Walach for their constructive comments and suggestions. This work was supported by Deutsche Forschungsgemeinschaft, Graduiertenkolleg 1838 “Spectral Theory and Dynamics of Quantum Systems”.