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

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: given 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.


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 timedependent 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 [19].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 [1,4].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 [17] and [16,18], respectively.The projectorsplitting integrators have been shown to possess remarkable robustness to the typical presence of small singular values [10,18], as opposed to applying standard integrators to the differential equations of dynamical low-rank approximation that are given in [12,13].However, the projector-splitting integrators do not preserve symmetry or antisymmetry.
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 [19]: 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 projectorsplitting 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.

General matrices: recap of the projector-splitting integrator for dynamical low-rank approximation
The objective is to approximate large time-dependent matrices A(t) ∈ R m×n for 0 ≤ t ≤ T by rank-r matrices Y(t) with comparatively low rank r ≪ m, n, which require much less storage than A(t) when they are available in a factorized, SVD-like form.
The large, or often too large matrices A(t) may be given explicitly or they are the unknown solution to a matrix differential equation (with right-hand side function F : Dynamical low-rank approximation as presented in [12] determines Y(t) as the solution of a projected matrix differential equation, with a projection P(Y) onto the tangent space T Y Mr of the manifold of rank-r matrices at Y ∈ Mr, .
where Y 0 is a rank-r approximation to A 0 , typically obtained by a truncated singular value decomposition.(Here, F(t, Y) = .
The solution Y(t) to this projected matrix differential equation then stays in the rank-r manifold Mr.
To make this abstract formulation practically useful, rank-r matrices Y(t) are written (non-uniquely) in factored form where the slim matrices U(t) ∈ R m×r and V(t) ∈ R n×r each have r orthonormal columns, and the small square matrix S(t) ∈ R r×r is invertible.We choose the tangent space projection P(Y) as the orthogonal projection onto T Y (Mr) with respect to the Euclidean or Frobenius inner product A, B = vec(A) ⊤ vec(B), where vec(A) ∈ R mn is a vectorization of A. Then, P(Y) is given as an alternating sum of three subprojections [12], The projector-splitting integrator of [17] 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, K := US is updated, in the second substep S is updated, and in the third substep L := VS ⊤ .The algorithm alternates between the numerical solution of matrix differential equations (of dimensions m × r, r × r, n × r) and orthogonal decompositions of slim matrices (of dimensions m × r and n × r).One time step of integration from time t 0 to t 1 = t 0 + h starting from a factored rank-r matrix Y 0 = U 0 S 0 V ⊤ 0 proceeds as follows: and set S0 = S(t 1 ).

L-step : Update
Then, the approximation after one time step is given by To proceed, we iterate the procedure taking 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 [17] for the detailed description.
The projector-splitting integrator has remarkable properties.First, it reproduces rank-r matrices without error.
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 U(t), S(t), V(t), which contain a factor S(t) −1 on the right-hand sides [12,Prop. 2.1].Moreover, the local Lipschitz constant of the tangent space projection P(•) is proportional to the inverse of the smallest nonzero singular value [12,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 (Robust error bound, [10, Theorem 2.1]) Let A(t) denote the solution of the matrix differential equation (1).Assume that the following conditions hold in the Frobenius norm • = • F : 1. F is Lipschitz-continuous and bounded: for all Y, Y ∈ R m×n and 0 ≤ t ≤ T , 2. The non-tangential part of F(t, Y) is ε-small: for all Y ∈ M in a neighbourhood of A(t) and 0 ≤ t ≤ T .3. The error in the initial value is δ-small: Let Yn denote the rank-r approximation to A(tn) at tn = nh obtained after n steps of the projector-splitting integrator with step-size h > 0.Then, the error satisfies for all n with tn = nh ≤ T 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 [10,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 [17,10].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 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 V = U in the decomposition (3) of the (skew-)symmetric rank-r matrix Y = USU ⊤ .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 A(t) to (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 K = US in the decomposition Y = USU ⊤ .The second substep is a substantially modified update of S, which can be viewed as a Galerkin approximation in the basis provided by the first substep. Given 1 with a (skew-)symmetric r × r-matrix S 1 at time t 1 = t 0 + h by the following algorithm: Algorithm 1: One time step of the (skew-)symmetry preserving integrator Integrate from t = t 0 to t 1 the r × r differential equation Set S 1 = S(t 1 ).
To continue in time, we take 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.

Exactness property of the (skew-)symmetric integrator
The exactness result Theorem 1 extends in the following way.
Theorem 3 (Exactness property) Let A(t) ∈ R n×n be (skew-)symmetric and of rank r for t 0 ≤ t ≤ t 1 , so that A(t) has a factorization (3) with equal left and right factors, A(t) = U(t)S(t)U(t) ⊤ .Moreover, assume that the r×r matrix U(t 1 ) ⊤ U(t 0 ) is invertible.With Y 0 = A(t 0 ), the (skew-)symmetric integrator for .
Proof We note that the projector-splitting integrator and the (skew-)symmetric integrator have the same first step.Let U 1 ∈ R n×r be the matrix with orthonormal columns computed in the first substep.Due to the exactness of the matrix projectorsplitting integrator as given by Theorem 1 we know that U 1 and A(t 1 ) have the same range and therefore Denoting ∆A = A(t 1 )−A(t 0 ), the (skew-)symmetric integrator provides for the second substep the solution since Y 0 = A(t 0 ).The result after a time step of the (skew-)symmetric integrator is where the last equality holds because of ( 6) and the (skew-)symmetry of A(t 1 ).⊓ ⊔

Robustness to small singular values
The error bound of Theorem 2 extends in the following way.
Let Yn denote the rank-r approximation to A(tn) at tn = 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 tn = nh ≤ T 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 [10, 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 [8,Section II.4].However, like for the projectorsplitting integrator of [17], it is not known if an error bound of higher order in the stepsize 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 [21].
We now prepare for the proof of Theorem 4, which views the (skew-)symmetric integrator as a perturbed variant of the projector-splitting integrator.
In the following, we denote by U 1 ∈ R n×r the matrix with orthonormal columns obtained in the first substep of the integrator.We recall that the matrix projectorsplitting and the (skew-)symmetric integrator have the first substep in common.We denote by 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 Y 0 .For the local error analysis, the following lemma is needed.
Lemma 1 Let U 1 , A 1 be defined as above.The following estimate holds: Proof The local error analysis in [10] shows that the r ×n matrix Z = S ps 1 V ps,⊤

1
, where S ps 1 and V ps 1 are the matrices computed in the third substep of the projector-splitting algorithm, satisfies The square of the left-hand side can be split into two terms: Hence, From the second term it follows that By the (skew-)symmetry of A 1 , this implies which yields the result.⊓ ⊔ In the following lemma, we show that the approximation given after one time step is O(h(h + ε)) close to the solution of system (1) when the starting values coincide.

Lemma 2 (Local Error)
The following local error bound holds: 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 Y 1 = U 1 S 1 U ⊤ 1 and Lemma 1 we have that The analysis of the local error thus reduces to estimating To this end, we introduce the following quantity: for We observe that where R(t) is defined as the term in big brackets.From Lemma 1 and from the bound B of F, which yields for t 0 ≤ t ≤ t 1 we conclude that the remainder term is bounded by This yields that F(t, A(t)) can be written as where the defect is bounded via the Lipschitz continuity of F as We now compare the two differential equations By construction, the solution of the first differential equation at time The solution of the second differential equation is 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

General tensors: recap of the projector-splitting integrator for the dynamical low-rank approximation by Tucker tensors
The objective is to approximate time-dependent tensors1 A(t) ∈ C n1×•••×n d for 0 ≤ t ≤ T by tensors Y (t) of multilinear rank r = (r 1 , . . ., r d ), with r i ≪ n i .(We recall that r i is the rank of the ith matricization Mat i (Y ) ∈ C ni×n ′ i with n ′ i = j =i n j , which aligns all entries of Y with ith index k in the kth row.The retensorization is denoted by Ten i (•), such that Ten i (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 . Dynamical low-rank approximation as presented in [13] determines Y (t) as the solution of the projected matrix differential equation, with a projection P(Y ) onto the tangent space T Y Mr of the manifold of tensors of multilinear rank r at Y ∈ Mr, .
where Y 0 is a rank-r approximation to A 0 .(Here, F (t, Y ) = .
A(t) if A(t) is given explicitly.)Tensors Y (t) of multilinear rank r are represented non-uniquely in the Tucker form [5] (using here the multilinear notation of [14]) where the core tensor C(t) ∈ C r1×•••×r d is of full multi-linear rank and the basis ma- trices U i ∈ C n×ri have orthonormal columns.We choose the tangent space projection P(Y ) as the orthogonal projection onto T Y (Mr) with respect to the Euclidean inner product A, B = vec(A) * vec(B), where vec(A) is a vectorization of A. Then, P(Y) is given as an alternating sum of 2d − 1 subprojections [16], 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 [16] and [18].It is shown in [18,Section 6] that the proposed Tucker integrators are mathematically equivalent.The algorithm runs through the modes i = 1, . . ., 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 [11,3] for the formulation and implementation of this algorithm in the context of the MCTDH method [20] of molecular quantum dynamics in the chemical physics literature.Moreover, the Tucker integrator has been proved in [18] 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 [18, Theorems 4.1 and 5.1] for the precise formulation, which is very similar to the matrix case.

Symmetric and anti-symmetric tensors: a structure-preserving integrator for dynamical low-rank approximation
and A is anti-symmetric if for every permutation σ ∈ S(d), It follows from [5] and [7] that a symmetric/anti-symmetric tensor Y ∈ C n×•••×n of multi-linear rank r = (r, . . ., r) admits a Tucker decomposition where the core tensor C ∈ C r×•••×r is symmetric/anti-symmetric of full rank r and the basis matrix U ∈ C n×r is the same for all indices.
We assume that the right-hand side function in (7) is such that 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.

(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 U, is identical to the first substep of the general Tucker integrator in [18,16].The second substep is a Galerkin method with the updated basis and determines the updated (anti-)symmetric core tensor.

Given the (anti-)symmetric tensor
for some (anti-)symmetric core tensor With Y 0 = A(t 0 ) we obtain from the second substep of the algorithm which proves the exactness.⊓ ⊔

Robustness to small singular values
The robust error bounds from Theorem 4 and [18, 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 6 (Robust error bound) Let A(t) denote the (anti-)symmetric solution of the tensor differential equation (7) with F satisfying (10).Assume the following: 1. F is Lipschitz-continuous and bounded.
2. The non-tangential part of F (t, Y ) is ε-small: for all Y of multilinear rank (r, . . ., r) in a neighbourhood of A(t) and 0 ≤ t ≤ T .3. The error in the initial value is δ-small: Let Yn denote the (anti-)symmetric approximation of multinear rank (r, . . ., r) to A(tn) at tn = nh obtained after n steps of the (anti-)symmetric Tucker integrator with stepsize h > 0.Then, the error satisfies for all n with tn = nh ≤ T 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 6 proceeds similar to the proof of Theorem 4 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 U 1 .

Lemma 3
The following estimate holds: where c only depends on d and a bound for hL.
Proof The error bound of [18,Theorem 5.1] shows that there exists We observe that .
As in the matrix case we obtain Thanks to (anti-)symmetry we have To conclude, we observe 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 where ĉ 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 The problem reduces to estimating C We introduce the tensor In the same way as in the proof of Lemma 2 (replacing 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 6 from the local to the global errors by the standard argument of Lady Windermere's fan [9,Section II.3].

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 [2] and TensorLab package v3.0 [23].

Addition of symmetric tensors: a computationally inexpensive retraction
Let A ∈ C n×•••×n be a symmetric tensor of multi-linear rank r = (r, . . ., r) and let B ∈ C n×•••×n .We consider the addition of two given tensors, where C ∈ C n×•••×n is not necessarily of low rank and we want to compute a symmetric rank-(r, . . ., 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: We will compare the solution obtained by computing the full addition and retracting to the manifold of symmetric tensors of multilinear rank (r, . . ., 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×100×100 and multi-linear rank 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 [6,22] 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 .

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 A(t) = e tW e t D e tW ⊤ , 0 ≤ t ≤ 1 .
The matrix D ∈ R N ×N is diagonal with elements d j = 2 −j and the matrix W ∈ R N ×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 [12].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 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 [1,4], 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 We consider the Hamiltonian given by where x l ∈ R represents the position of the l-th particle and we choose the torsion potential 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 ddimensional tensor Y (t) ∈ C K×•••×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 [15] while the two-particle interaction is treated with the anti-symmetric low-rank integrator of Section 5. We introduce the space discretization Let Y (t) be a time-dependent tensor defined element-wise by Y j1,j2,...,j d (t) = ψ(t, x j1 , . . ., x j d ) .
The fermionic property of the system implies that the tensor Y (t) is anti-symmetric.Denoting F K the Fourier matrix, we define The Fourier collocation space discretization of ( 11) is equivalent to the system 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 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, Afterwards, we start considering the equation We matricize in the first mode, The solution can now be computed with the 1-dimensional split-step Fourier method.Denoting and tensorizing back in the first mode we have that Taking this as initial condition and iterating the same process for all the successive modes we obtain the updated anti-symmetric tensor We now apply the anti-symmetric low-rank integrator to the multi-particle interaction, Ẋ = −W [X], X(t 0 ) = X 0 , where 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 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 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 roundoff 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.

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 [4] for the physical description of the problem and its MCTDHF formulation.We consider the time-dependent Schrödinger equation i∂ t ψ = H(t)ψ, ψ(t 0 ) = ψ 0 , with the Hamiltonian with the torsion potential V as before and with parameters ω(t) := A 0 e −t 2 /τ 2 sin(Ωt), A 0 = 100, Ω = 100, τ = 0.2π .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.
t1−s) D(s) ds ≤ e Lh 2L(Bh + ϑ)h.The result now follows using the definition of ϑ. ⊓ ⊔ Thanks to the Lipschitz continuity of the function F, we conclude the proof of Theorem 4 from the local to the global errors by the standard argument of Lady Windermere's fan [9, Section II.3].

Fig. 3
Fig.3Evolution to the fermionic ground state energy computed with rank 5.