An optimization approach for dynamical Tucker tensor approximation

An optimization-based approach for the Tucker tensor approximation of parameter-dependent data tensors and solutions of tensor differential equations with low Tucker rank is presented. The problem of updating the tensor decomposition is reformulated as fitting problem subject to the tangent space without relying on an orthogonality gauge condition. A discrete Euler scheme is established in an alternating least squares framework, where the quadratic subproblems reduce to trace optimization problems, that are shown to be explicitly solvable and accessible using SVD of small size. In the presence of small singular values, instability for larger ranks is reduced, since the method does not need the (pseudo) inverse of matricizations of the core tensor. Regularization of Tikhonov type can be used to compensate for the lack of uniqueness in the tangent space. The method is validated numerically and shown to be stable also for larger ranks in the case of small singular values of the core unfoldings. Higher order explicit integrators of Runge-Kutta type can be composed.


Introduction
Low-rank matrix and tensor approximation occurs in a variety of applications as model reduction approach [7].This paper concerns the problem of fitting parameter or time-dependent matrices and tensors with small (Tucker) rank within the Tucker tensor format [15], that is, updating the factors and core tensor of the model instead of working with the tensor itself.In a discretetime setting for tensor streams, incremental fitting procedures for dimensionality reduction were introduced as dynamic higher order generalizations of principal component analysis [25].While this approach seeks to minimize a reconstruction error by a sequence of Tucker tensor approximations sharing same factors and is discrete in time, a different time continuous setting was recently discussed by Koch and Lubich [12,13], where differential equations for the factor matrices and core tensor of the Tucker approximation were derived.For a time-varying family of tensors Aptq P R I 1 ˆI2 ˆ¨¨¨ˆI d and the approximation manifold M r (Tucker tensors with rank r) an alternative to the best-approximation problem min XptqPMr }Xptq ´Aptq} (1) is to considered a dynamical tensor approximation problem based on derivatives 9 Y which lie in the tangent space T Y M r , i.e., min 9 YptqPT Y Mr } 9 Yptq ´9 Aptq} (2) with an initial condition Yp0q, typically a best-approximation in M r of Ap0q.Note that for the choice 9 Aptq " F pt, Yptqq this becomes the defect approximation of a tensor differential equation 9 Yptq " F pt, Yptqq.A prominent example of a tensor differential equation on M r is the multiconfiguration time-dependent Hartree approach (MCTDH) of molecular quantum dynamics [3,17], where a Schrödinger equation subject to M r is solved.More precisely, the MCTDH ansatz takes the form i 9 ψ " H ψ subject to ψpx, tq " where the discretized constraint on the many-particle wave function ψ leads to a Tucker tensor manifold.A time-dependent variational principle is used to extract differential equations for the coefficient tensor `aj 1 ¨¨¨j d ˘and the factors ϕ pκq jκ .This procedure is equivalent to solving problem (2) for 9 Aptq " H Yptq, with Y denoting the discretized wave function in Tucker format, by applying the Galerkin condition xi 9 Y ´H Y, Vy " 0 for all V P Ă T Y M r , for a gauged tangent space Ă T Y M r .In this sense, problem (2) yields an initial value problem on M r in form of a system of nonlinear differential equations for the factors and core tensor in the Tucker decomposition.This system is explicit if an orthogonality gauge condition in T Y M r is imposed.However, the arising equations involve the (pseudo) inverse of the core unfoldings, which might be close to singular.Consider for instance the two dimensional (low-rank matrix) case, where M r consists of matrices Y " UCV T P R I 1 ˆI2 with U P R I 1 ˆr, V P R I 2 ˆr and C P R rˆr nonsingular.The mentioned equations for the factors of the model take the form (omitting the argument t) [12] 9 U " pI ´UU T q 9 AVC ´1,

9
V " pI ´VV T q 9 AUC ´T ,

9
C " U T 9 AV. ( Nonzero singular values of Y are those of the core C, which has to be inverted in the first two equations of (3).Hence, very small singular values lead to severe problems concerning stability of numerical explicit or implicit integrators [10].Unfortunately, this is a very realistic case, since in practice the prescribed initial ranks should be chosen larger than the unknown effective ranks (overestimation) because underestimating necessarily leads to large model errors.Therefore, a numerical scheme for (2) has to manage the rank deficient case.Recently, Lubich and Oseledets [18] introduced an operator splitting approach for the low-rank matrix approximation, which behaves robust in the mentioned problematic case.In the meanwhile the method was extended to the tensor train format [19,23] and applied to MCTDH [17] with related analysis presented in [21].Recently, also projected time-stepping methods for tensor differential equations on a low-rank tensor train manifold where established [5] and analyzed [11].Moreover, we stress that in applications the equations (3) (resp.higher dimensional versions) are often not well-defined due to initial states, which exhibit rank(s) strictly smaller than r, leading to singular values of core unfoldings that are exactly zero.As a consequence, the factor derivatives or core unfoldings have to be initially modified appropriately or, if feasible, the full tensor equations are evolved over a short time to generate a proper initial state [3,22].
In this paper we break out in a different direction.We will develop an integration scheme using the optimization framework (2) itself rather than a system of nonlinear differential equations.
As in previous approaches the approximation manifold is equipped with a full rank condition for the factors while the rank deficient case of core unfoldings is tackled by regularizing a fitting problem which exhibits the resulting (local) error.For that purpose, it turns out to be beneficial to give up the gauge condition in the tangent space.Remarkably, the proposed scheme does not need any inversion or (pseudo) inverse of the core unfoldings (density matrices in MCTDH).The approach is motivated from static fitting problems of data tensors, where quasi best-approximation in the Tucker format can be derived very efficiently from alternating least squares (ALS) or general nonlinear minimization of the norm of the residual [8,15].For instance, the key of the higher order orthogonal iteration (HOOI) [4] is the reduction of the quadratic subproblems in the ALS scheme to the task of computing an orthonormal basis of a dominant eigenspace [14].We stress that alternating linear schemes, where ALS is a special case, were successfully applied to optimization tasks subject to a low-rank tensor train manifold [9] and local convergence was proven in the case of convex functionals [24].Here, we establish a discrete Euler scheme for the dynamical problem by alternating least squares, where we impose orthonormal columns for the updated factor matrices as a constraint.As a consequence, the originally quadratic subproblems get linear with explicit solutions via SVD.The method does not require any pseudo inverse of unfoldings of the core.Hence, the approach is far less affected by small singular values, although the stability of the SVD solutions of the involved trace optimization problems can still be negatively affected in severe cases.It turns out that Tikhonov regularization reduces this problem while also being a compensation for the lack of uniqueness of the tangent space representation.The regularization only acts on the column space of the factor matrices and would be inactive if the conventional orthogonality gauge condition would have been imposed.
In the following section we introduce tensor notation and give some details to Tucker decompositions and the approximation manifold.In section 3 we review the approach of Koch and Lubich [13].Section 4 formulates the optimization framework and establishes the Euler scheme.Finally, we validate our approach by means of numerical examples.

Prerequisites
In the following a brief description of tensors, Tucker decomposition and the manifold of rankpr 1 , . . ., r d q tensors is given.The reader is referred to the review article [15] or the book [8] for further details and references therein.We treat real-valued matrices and tensors, but emphasize that the complex case (as in quantum mechanics applications) does not pose any problems.Specifically, the derivation of the Euler scheme of Sec. 4 with its Theorem 3 generalizes to this case.

Notation
We denote matrices by boldface capital letters, e.g., A and vectors by boldface lowercase letters, e.g., a.A tensor is denoted by boldface Euler script letters, e.g., X. Scalars are denoted by lowercase letters, e.g., α.
We use the Frobenius inner product and norm for matrices, denoted by x., .yand }.}, respectively.The norm of a tensor X P R I 1 ˆI2 ˆ¨¨¨ˆI d , d ě 2, is the square root of the sum of the squares of all elements, or equivalently, the Euclidean norm of the vectorized tensor The inner product of two tensors X, Y P R I 1 ˆI2 ˆ¨¨¨ˆI d is the sum of the product of their entries, i.e., xX, Yy " vecpXq T vecpYq. ( There holds }X} 2 " xX, Xy.Fibers of a tensor X P R I 1 ˆI2 ˆ¨¨¨ˆI d are the generalizations of matrix rows and columns.A mode-n fiber of X consists of the elements with all indices fixed except the nth index.The nth unfolding of a tensor X P R I 1 ˆI2 ˆ¨¨¨ˆI d is the matrix X pnq P R Inˆś j‰n I j that arranges the moden fibers to be the columns of the resulting matrix.We denote the corresponding matricization as rXs pnq " X pnq .
The n-rank of a tensor and the vector r " pr 1 , . . ., r d q T is called the Tucker rank of X.
The n-mode product of a tensor X P R I 1 ˆI2 ˆ¨¨¨ˆI d with a matrix U P R JnˆIn results in a tensor which is defined as There holds (with appropriate dimensions) and for n ‰ m We will further denote the range of a matrix U P R JnˆIn (column space) with RpUq and the orthogonal complement with RpUq K .
2.2 Tucker decomposition and manifold of rank-pr 1 , . . ., r d q tensors A Tucker decomposition [26] of a tensor X P R I 1 ˆI2 ˆ¨¨¨ˆI d has the form where the core tensor C P R r 1 ˆr2 ˆ¨¨¨ˆr d and the factor matrices In the following we further assume the factor matrices having orthonormal columns, that is Note that this does not lead to a unique Tucker decomposition, since r A Tucker decomposition can be written as a linear combination of rank-1 tensors that are formed as the outer products of the columns of the factor matrices, that is where C " pc j 1 ,...,j d q.The nth unfolding of a Tucker tensor is given by where b denotes the Kronecker product.
For given Tucker rank r " pr 1 , . . ., r d q T the set is the manifold of rank-pr 1 , . . ., r d q tensors.Every element X P M r has a representation as a Tucker decomposition, where w.l.o.g we assume the factor matrices having orthonormal columns, that is U T k U k " I.We will use M r as an approximation manifold with (typically) r n !I n .In the special case of d " 2 elements of M r are usually referred to as low-rank matrices, where r " pr, rq T .In this case, the conventional notation is X " UCV T P R I 1 ˆI2 with U P R I 1 ˆr, V P R I 2 ˆr and C P R rˆr nonsingular, but not necessarily diagonal as in the case of a (truncated) singular value decomposition (SVD).For Y P M r the corresponding tangent space is denoted with T Y M r , where 9 with 9 C P R r 1 ˆr2 ˆ¨¨¨ˆr d and 9 U n P T Un V In,rn , where is the tangent space of U P V I,r , the Stiefel manifold of real I ˆr matrices with orthonormal columns.If the gauge condition (orthogonality constraint) is imposed, we denote the tangent space with Ă T Y M r .We stress that the gauge leads to a unique parametrization of the tangent space, in fact, the RpUq-component of the variation can be absorbed by the first term in ( 16) [8].

Related approach of Koch and Lubich
Dynamical approximation of time-dependent tensors Aptq were previously considered in [13] with help of an orthogonal projection on the tangent space Ă T Y M r .This can be accomplished by a Galerkin condition, also known as Dirac-Frenkel variational principle [3], that is (omitting the argument t) The core and the factors of the decomposition of Y fulfilling (19) are explicitly given by a system of differential equations.We state here the associated theorem from [13].
Theorem 1 ( [13]).For a tensor Y P M r with n-mode factors having orthonormal columns the condition (19) is equivalent to with the core and the factors satisfying where P K n " I ´Un U T n is the orthogonal projection onto RpU n q K .The equations ( 21) can be integrated for a given initial tensor Yp0q P M r to obtain the Tucker decompositions Yptq P M r , t ą 0, which are approximations to Aptq P R I 1 ˆI2 ˆ¨¨¨ˆI d in the sense of the Galerkin condition (19), or equivalently, the optimization problem min The error of ( 19) can be bounded in terms of the best-approximation error (1), which can be extended to the case of tensor differential equations if (one-sided) Lipschitz conditions for F are assumed [13].
4 Optimization approach and Euler scheme via ALS

Dynamical approximation without orthogonality constraint in the tangent space
We can derive generalizations of (21) if we give up the orthogonality constraint U T n 9 U n " 0. We will then numerically tackle these ungauged equations by an optimization scheme for the defect approximation (23) circumventing the explicit form of the factor equations.This will allow us to overcome numerical problems coming from the orthogonal projections and the necessary inversion of core unfoldings.
Theorem 2. For a tensor Y P M r with n-mode factors having orthonormal columns the optimization problem with the core and the factors satisfying which reduces to (21 Proof.We look at the optimization problem min Since a tangent tensor in T Y M r has the form (24) we will minimize for given Yptq P M r the objective function with respect to 9 C and 9 U k .Rewriting yields This is a convex quadratic function in 9 C, hence, setting the gradient with respect to 9 C equal to zero gives the first equation in (25).We now define the part of the objective function involving only terms in 9 U n with fixed index n (utilizing nth unfolding) as the convex quadratic function The derivative of f n set to zero gives the second equation of (25).
Note that a discretization scheme using the derivatives (21) or (25) requires the solution of a least squares problem for the factor matrices by calculating the pseudo inverse C T pnq `Cpnq C T pnq ˘´1 of the nth-unfolding of the core tensor.Small singular values in C pnq will affect the solution and make the method unstable [10].However, the advantage of (25) lies in the fact that the derivative 9 U n is not orthogonal to the range of U n .This becomes apparent in a discretization scheme that we will establish in the forthcoming section, where we will establish a numerical time-stepping scheme that, in particular, circumvents the explicit formulation of the factor differential equations of (25).

Euler method under orthogonal columns constraint
We will consider Euler's method for solving (23) via an alternating least squares (ALS) scheme.For that purpose, all but one factor matrix derivative will be kept fixed consecutively, where the core derivative is assumed to be already computed according to the first equation in (25).Consider an Euler step at t for the second equation of (25) with length h ą 0, that is Note that the orthogonality gauge U T n 9 U n " 0 prevents the updated factor matrices from having orthonormal columns, that would be U n pt `hq T U n pt `hq " I. On the other hand, as will become apparent from Theorem 3 below, the function f n (26) with respect to the variable U h n :" U n pt `hq gets linear if U h n T U h n " I is imposed.We will therefore seek for an Euler step under the orthogonality constraint for the updated factor matrices.More precisely, for fixed t we define for V n P R Inˆrn (omitting the argument t) r f n pV n q :" f n `1 h pV n ´Un q ˘, ( and define U h n as the solution of the constraint optimization problem min where the involved derivative 9 Aptq is treated explicitly.The following theorem characterizes the updated factor matrices U h n from one Euler step according to (29) as the solution to a trace optimization problem, which is computable from an SVD of a matrix of size I n ˆrn .Theorem 3. The optimization problem min is equivalent to the trace optimization problem max with Moreover, given the (economy size) SVD B n " R n Σ n T T n with R n P R Inˆrn , T n P R rnˆrn and the diagonal matrix Σ n P R rnˆrn , the solution U h n to (31) is given by where we used V T n V n " U T n U n " I. Thus, the optimization problem (30) is equivalent to max which is the trace optimization problem (31).Consider now the singular value decomposition T n P R rnˆrn .Note that we assume r n ď I n .Cyclic permutation of the arguments in the trace yields where we define the matrix W n :" r R T n V n T n P R Inˆrn .Owing to W T n W n " I, we have for the diagonal elements |w jj | ď 1.Hence, the maximum of the objective in (31) is attained for w jj " 1, which implies I " R T n V n T n and consequently V n " R n T T n .
Note that the matrices B n do not require inversion of the potentially problematic matrices C pnq C T pnq .Theorem 3 readily generalizes to the case of complex valued dynamical tensor approximation.More precisely, in the subproblems of the ALS scheme the functions f n boil down to the real part of the inner product, that is Rerx 9 U n , B 1 n ys.As a consequence, the equivalent trace optimization problem takes the form max with the solution where B " RΣT H the economy sized SVD of B.

Algorithm for Euler step via ALS
The computation of an Euler step for (25) via alternating least squares for ( 23) is summarized in Algorithm 1.
Algorithm 1 Euler step from t Ñ t `h via ALS for ( 23) An Euler method for the dynamical approximation (23) using Algorithm 1 starts from a rank-r decomposition or (best-) approximation of Ap0q, e.g., via higher order orthogonal iteration (HOOI), see [15].The termination criterion uses the discrete version of the norm of the defect (resp.residual) (23), that is Algorithm 1 computes one time step of the underlying explicit Euler scheme for the initial value problem associated with the defect approximation (23).Higher order schemes of Runge-Kutta type can be composed based on the above Euler scheme, see section 5.2.The derivative at time t, i.e., 9 Aptq or F pt, Yptqq respectively, is treated explicitly and assumed to be efficiently computable.The overall computational effort depends heavily on this derivative evaluation, compare for instance with the numerical example in section 5.4.One inner loop evaluation in the Euler-ALS scheme described in Algorithm 1 involves the computation of the rectangular matrix B n of size I n ˆrn at the costs of O `pś d k"1 r k q ř d k"1 I k ˘operations plus the costs for the term r 9 A Ś k‰n U T k s pnq .If we assume 9 A P M r 1 (e.g. after projection/compression), the latter term takes O `In r 1 n ś k‰n r k `pś k‰n r 1 k q ř k‰n r k ˘operations.The economy sized SVD scales with OpI n r 2 n q.

Regularization
In the course of the derivation of Algorithm 1 the orthogonality gauge U T n 9 U n " 0 of [13] was given up.As a consequence, elements of the tangent space T Y M r have no unique representation any more.Although, we have not observed any related problems in the numerical tests, we found that including a regularization of Tikhonov type helps to stabilize the SVD in cases where B n is effectively rank deficient and where C pnq C T pnq is numerically close to singular, respectively.More precisely, we can treat the problem (α ą 0) compare with (26).Note that this positively affects definiteness of the quadratic term in (26), which becomes x 9 U n , 9 U n pC pnq C T pnq `αIqy.The discrete scheme changes as follows.We treat the modified optimization problem min which is equivalent to max according to Theorem 3. Observe that the regularization term is only effective because U T n 9 U n ‰ 0. In fact, if 9 U n " 1 h pV n ´Un q P RpU n q K we get where P n is the orthogonal projection onto RpU n q.Hence, the regularization term would degenerate to TrrU T n V n s " The choice α 9 h 2 makes the regularization term in (38) independent of the step size h.If regularization was considered in the numerical examples below we have always chosen α " h 2 .

Numerical validation
We use the Tensor Toolbox Version 2.6 [2], which implements efficient tensor arithmetics and fitting algorithms for different tensor classes in MATLAB [1].Results were achieved by using the Vienna Scientific Cluster 3 (VSC3).

Koch/Lubich Example
We consider the numerical example as in [13].A time-dependent data tensor A P R 15ˆ15ˆ15ˆ15 is constructed as where U k P R 15ˆ10 are random matrices with orthonormal columns, B P R 10ˆ10ˆ10ˆ10 a random core tensor, and C P R 15ˆ15ˆ15ˆ15 a random perturbation.The dynamical Tucker tensor approximations are computed with r " p10, 10, 10, 10q T , fixed h " 10 ´4 and different ε P t10 ´5, 10 ´4, 10 ´3u.We use Euler steps described in Algorithm 1 without regularization.
For initialization we use ∆U n " 0, n " 1 . . ., d.Moreover, a quasi best-approximation for Ap0q « A 0 is used, which is derived from HOOI.Note that this example solves the tensor initial value problem Results are shown in Figure 1, where relative errors are given, i.e., }Aptq ´Yptq}{}Aptq}.Alg. 1 was used with a tolerance of 10 ´5 for the change in the relative defect } 9 Aptq ´9 Yptq}{} 9 Aptq} (cf.(36)), which was basically reached after two loops.The norm of the relative defect itself was constantly in the order of ε (about 5ε).

Second Order Method
We repeat Example 1 with r " p8, 9, 10, 11q T and adjusted sizes in (40) for U k and B. The second order improved Euler method (extrapolated Euler) is now used, i.e., for an initial value problem y 1 " f pt, yq yp0q " η 0 the new iterate η ν at time t ν " t ν´1 `h is calculated from the previous iterate η ν´1 at time t ν´1 according to η ν " η ν´1 `hf pt ν´1 `h 2 , η ν´1 `h 2 f pt ν´1 , η ν´1 qq, which is a 2-stage explicit Runge-Kutta formula.The half and full Euler steps, respectively, are calculated using Alg. 1. Figure 2 shows relative errors for different step sizes h and ε " 10 ´10 on the interval r0, 1s (left) and compares these errors among each other (right).We can observe second order convergence in h.Again the norm of the relative defect for the time steps was in the order of about ε independent of the step size.

Small singular values in higher dimensions
We generalize the example from [10] to higher dimensions.A time-dependent tensor Aptq P R I 1 ˆI2 ˆ¨¨¨ˆI d with I k " I is constructed as with skew-symmetric random matrices W k P R IˆI and a super-diagonal tensor C P R I 1 ˆI2 ˆ¨¨¨ˆI d with diagonal entries c j " 1{2 pd´1qj , j " 1, . . ., I. We now use the Euler method of Alg. 1 and regularization with α " h 2 .Figure 3 shows relative errors at t " 0.3 against step length h for different ranks in the two dimensional case with I " 100.We repeat these computations for d " 3 and I " 50 and give associated results in Figure 4.
The tolerance for the change of the fit was 10 ´6 and basically reached after two loops in Algorithm 1.We observe in both cases the desired behavior of linearly decreasing error with respect to h towards the model accuracy, and stability also for higher ranks and larger step size.
Next we take a reaction-diffusion example in d " 2, 3 space dimensions similar to [22]: Bu Bt " 0.01 ∆u `0.1 u 3 , x P Ω " p0, 1q d , t ą 0, upx, tq " 0, x P BΩ, t ą 0, upx, 0q " 10d d ź n"1 e ´100 px pnq ´0.5q 2 , x P Ω. (44) We discretize the Laplace operator uniformly by second order centered finite differences leading to a d-term Kronecker sum with K 1 , I P R IˆI and K 1 " 1{k 2 tridiagp1, ´2, 1q, k " 1{I.Note that such a discretization of an operator as a sum of Kronecker product matrices applied to a grid function represented in Tucker format results in a sum of Tucker tensors of same rank, e.g., [6].This makes computations on large grids possible.The nonlinearity in (44) leads to elementwise (Hadamard ) products of tensors, which results in cubed Tucker ranks.The Hadamard products are recompressed via HOOI, however, more sophisticated efficient techniques for Hadamard products in the Tucker format were recently introduced [16] but not considered here.We solve (44) subject to the lowrank manifold M r for different ranks r " r P t3, 4, 5u, where we use the improved Euler method of example 5.2 with basic time steps realized by Alg. 1 using regularization.The computations of Figure 6 were carried out for the case d " 2 with I " 400 and step size h " 10 ´4, where the norm of the relative defect as a measure for local (time-stepping and projection) errors for different ranks is shown during time integration towards the blow-up.One observes decreasing 0 0.002 0.004 0.006 0.008 0.01 0.012 local defect in the initial stabilization phase and linearly increasing defect afterwards.The latter is due to effective rank growth caused by the nonlinear term in (44).Fig. 7 shows the initial data and the numerical solution at time t " 0.012 corresponding to problem (44).Fig. 8 shows the difference between a reference solution obtained from an adaptive fourth order Runge-Kutta scheme applied to the full 2d system and the rank 3 and 5 low-rank numerical solutions.We repeat the computations for the case of d " 3, where treatment of the full equations gets to expensive both in storage and computational effort.Fig. 9 shows local relative defects for different ranks, where we can observe a similar behavior as for the 2d case.Fig. 10 shows average timings for one iteration in Alg. 1 for varying mode size I and ranks r in the case d " 3. The function evaluation involves compression (via HOOI) of Hadamard products with squared ranks, which is the most expensive part.For larger ranks the part of B n involving the already computed right hand side scales with r 3 .The rest, especially scalings with respect to the mode size, scales constant or at most linearly in the range of consideration.

Conclusion
An optimization-based approach for dynamical low-rank matrix and Tucker tensor approximation of parameter-dependent data tensors and solutions of tensor differential equations such as MCTDH of molecular quantum dynamics was presented.The approach is motivated by the effectiveness of alternating schemes for static tensor fitting problems, which reduce quadratic optimization subproblems to eigenproblems accessible via SVD.The presented method uses alternating least squares to establish an Euler scheme for the dynamical fitting problem on the tangent space without a gauge constraint.The quadratic subproblems are found to be equivalent to trace optimization with solutions being explicitly computable via SVD of small size.Remarkably, the method needs no pseudo inverse of matrix unfoldings of the core tensor (or inverse of density matrices in MCTDH), which enhances stability in the presence of small singular values.This makes the method robust and eliminates numerical stability problems coming from "overestimated" initial values, which is confirmed by means of computational results.Regularization of Tikhonov type is suggested to compensate for the lack of uniqueness of the tangent space representations, having the positive side effect to further stabilize in rank deficient cases.Higher order methods can be composed, where a second order example is given.
In forthcoming work the approach can be extended to the hierarchical Tucker format and tensor trains (matrix product states) [20].

Figure 2 :
Figure2: Errors of the dynamical approximation of (40) with r " p8, 9, 10, 11q T and ε " 10 ´10 using the improved Euler method (left) and comparison of the errors among each other (right).

Figure 3 :
Figure 3: Errors at t " 0.3 of the dynamical approximation of (41) in the case d " 2 and I " 100 for different step sizes h and ranks r " pr, rq.

Figure 6 :
Figure 6: Local relative defect of the dynamical approximation of (44) in the case d " 2 and I " 400 for different ranks r " pr, rq.

Figure 7 :Figure 8 :
Figure 7: Initial data (left) and the numerical solution at time t " 0.012 (right) corresponding to (44) in the case d " 2.

Figure 9 :
Figure 9: Local relative defect of the dynamical approximation of (44) in the case d " 3 and I " 400 for different ranks r " pr, r, rq.

Figure 10 :
Figure 10: Timings for one iteration in Alg. 1 for the discretization of problem (44) in the case d " 3. Left: Varying rank r for fixed mode size I " 300.Right: Varying mode size I for fixed rank r " 5.In each case timings is given for the right hand side evaluation (tF), the part of the B n computation involving the already evaluated right hand side (tBF), as well as the part without right hand side and SVD (tB F).