SIAM J. MATRIX ANAL. APPL.

c 2010 Society for Industrial and Applied Mathematics


Vol. 31, No. 5, pp. 2360–2375

DYNAMICAL TENSOR APPROXIMATION∗

Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

OTHMAR KOCH† AND CHRISTIAN LUBICH‡

Abstract. For the approximation of time-dependent data tensors and of solutions to tensor
differential equations by tensors of low Tucker rank, we study a computational approach that can
be viewed as a continuous-time updating procedure. This approach works with the increments
rather than the full tensor and avoids the computation of decompositions of large matrices. In this
method, the derivative is projected onto the tangent space of the manifold of tensors of Tucker
rank (r1 , . . . , rN ) at the current approximation. This yields nonlinear differential equations for the
factors in a Tucker decomposition, suitable for numerical integration. Approximation properties of
this approach are analyzed.

Key words. low-rank approximation, time-varying tensors, continuous updating, Tucker de-
composition, tensor differential equations

AMS subject classifications. 15A18, 15A69, 65F99, 65L05

DOI. 10.1137/09076578X

1. Introduction. Tensor approximation is an active research area with both

interesting theoretical questions and numerous applications in the compression and
retrieval of large structured data; see, e.g., [3, 5, 10, 15, 19] and the references therein.
While these papers deal with the approximation of a fixed tensor, we consider here
the approximation of time-varying (or parameter-dependent) tensors by tensors of low
Tucker rank. There arises the question of how to update a given tensor approximation,
a problem that has been addressed in a different, discrete-time setting in [16, 17]. The
approach considered in the present paper can be viewed as a continuous-time updating
technique, which works only with the increments in the tensors rather than the tensors
themselves and does not require computing any decompositions of large matrices. It
extends the dynamical low-rank approximation of matrices [7, 14] to the higher-order
tensor case.
Consider a time-varying family of tensors A(t) ∈ RI1 ×···×IN for 0 ≤ t ≤ t. Let
Mr denote the manifold of all order-N tensors of Tucker rank r = (r1 , . . . , rN ) (see
section 2 for the notion of Tucker rank or mode rank), where rn ≤ In (and typically
rn  In ) for n = 1, . . . , N . The best approximation to A(t) in Mr (with respect to
the Frobenius norm  · ) is

(1.1) X(t) ∈ Mr such that X(t) − A(t) = min .

Here, we consider instead the dynamical tensor approximation Y(t) ∈ Mr determined

from the condition that for every t the derivative Ẏ(t), which is in the tangent space
TY(t) Mr , be chosen as

(1.2) Ẏ(t) ∈ TY(t) Mr such that Ẏ(t) − Ȧ(t) = min .

∗ Received by the editors July 21, 2009; accepted for publication (in revised form) by L. De

Lathauwer April 23, 2010; published electronically July 15, 2010. This work was supported by the
Austrian Academy of Sciences’ APART program and by DFG, SPP 1324.
[Link]
† Institute for Analysis and Scientific Computing, Vienna University of Technology, Wiedner

Hauptstrasse 8-10, A-1040 Wien, Austria (othmar@[Link]).

‡ Mathematisches Institut, Universität Tübingen, Auf der Morgenstelle 10, D-72076 Tübingen,

Germany (lubich@[Link]).
2360

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 DYNAMICAL TENSOR APPROXIMATION 2361

This is complemented with an initial condition, ideally Y(0) = X(0). Note that for
given Y(t), the derivative Ẏ(t) is obtained in (1.2) by a linear projection, though
Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

onto a state-dependent vector space. Problem (1.2) yields an initial value problem of
nonlinear ordinary differential equations on Mr , which becomes numerically efficiently
accessible as differential equations for the factors in the Tucker decomposition of
tensors of Tucker rank (r1 , . . . , rN ).
In a different context, a closely related approach was developed in the multicon-
figuration time-dependent Hartree method of multiparticle quantum dynamics [2, 13],
where the multidimensional time-dependent wave function is approximated by a linear
combination of products of functions depending only on a single spatial variable.
We will study theoretical properties of the dynamical tensor approximation (1.2)
in this paper: in section 2 we derive the differential equations that are to be solved
numerically, and after an auxiliary section on the tangent space projection (section
3) we study approximation properties in section 4. Some numerical experiments are
given in section 5 (see also [14] for further numerical experiments with time-dependent
3-tensors that arise from a discretized reaction-diffusion PDE in 3 space dimensions).
The present paper extends the dynamical low-rank approximation of matrices,
studied in our paper [7], to tensors. Though there are many conceptual similarities
with the matrix case, the analysis of the tensor case is not simply a straightforward
extension and therefore requires a careful discussion. While we have kept the general
organization of the paper largely parallel to [7] to make similarities and differences
easily visible, we note that the material in section 2 is quite apart from the matrix
case in both notation and techniques, and the analysis of the tangent space projection
in section 3 requires different arguments. Once the projection estimates from section
3 are available, some of the results in section 4 then have essentially the same proof
as the corresponding results in the matrix case (Theorems 4.1 and 4.2), whereas the
proof of Theorem 4.3 follows a different line.
2. Differential equations for dynamical tensor approximation.
2.1. Prerequisites. We use the tensor notation of the review article [8], to
which we refer for further details, results, and numerous references.
Norm and inner product of tensors. The norm of a tensor Y ∈ RI1 ×···×IN is
the Euclidean norm of the vector y that carries the entries yi1 ,...,iN of Y. The inner
product X, Y of two tensors X, Y ∈ RI1 ×···×IN is the Euclidean inner product of the
two corresponding vectors x and y.
Unfolding and reconstruction. The nth unfolding of a tensor Y ∈ RI1 ×···×IN is
the matrix Y (n) ∈ RIn ×In+1 ...IN I1 ...In−1 that aligns all entries yi1 ,...,iN with fixed in in
the in th row of the matrix, ordered lexicographically. We denote
Y (n) = [Y](n) ,
and clearly the tensor Y can be reshaped from its unfolding Y(n) : we write

Y = [Y (n) ](n) .
The n-mode product. For a tensor Y ∈ RI1 ×···×IN and a matrix V ∈ RJn ×In , the
n-mode product
Y ×n V ∈ RI1 ×···×In−1 ×Jn ×In+1 ×···×IN
is defined by the relation
(2.1) [Y ×n V](n) = VY (n) .

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 2362 OTHMAR KOCH AND CHRISTIAN LUBICH

In particular, we then have, for another matrix W of appropriate dimension,

(2.2) (Y ×n V) ×n W = Y ×n (WV) .
Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

2.2. Manifold of rank-(r1 , . . . , rN ) tensors and Tucker decomposition.

The n-rank of a tensor Y ∈ RI1 ×···×IN is
rn = rank (Y(n) ) ,
and the vector r = (r1 , . . . , rN ) is known as the Tucker rank of the tensor. For given
n-rank rn ≤ In (and typically rn  In ), the set
 
Mr = Y ∈ RI1 ×···×IN : Y has n-rank rn for n = 1, . . . , N
is a manifold that will serve as an approximation manifold for general tensors A ∈
RI1 ×···×IN . As is known (see [9]), every tensor Y ∈ Mr can be written as a Tucker
decomposition
N
(2.3) Y = S ×1 U1 · · · ×N UN =: S X Un ,
n=1

where the core tensor S ∈ Rr1 ×···×rN is of full Tucker rank r = (r1 , . . . , rN ), and the
(n)
matrices Un ∈ RIn ×rn have orthonormal columns ujn (jn = 1, . . . , rn ). In terms
of the entries of S = (sj1 ,...,jN ), the above expression can be rewritten as a linear
combination of rank-1 tensors that are formed as the outer products of the column
vectors:
 (1) (N )
Y= sj1 ,...,jN uj1 ◦ · · · ◦ ujN .
j1 ,...,jN

In terms of the nth unfolding, we have the useful matrix formula


(2.4) Y(n) = Un S(n) (UN ⊗ · · · ⊗ Un+1 ⊗ Un−1 ⊗ · · · ⊗ U1 )T =: Un S(n) UTk ,
k=n

where ⊗ denotes the Kronecker matrix product.

The representation (2.3) is not unique: replacing Un by Ũn = Un Qn with orthog-
onal matrices Qn and S by S̃ = S Xn=1 QTn yields the same tensor Y = S Xn=1 Un =
N N

S̃ Xn=1 Ũn .
N

2.3. Tangent tensors. As a substitute for the nonuniqueness in the decompo-

sition (2.3), we will use a unique decomposition in the tangent space. Let VI,r denote
the Stiefel manifold of real I × r matrices with orthonormal columns. The tangent
space at U ∈ VI,r is
T
TU VI,r = {U̇ ∈ RI×r : U̇ U + UT U̇ = 0} = {U̇ ∈ RI×r : UT U̇ ∈ so(r)},
where so(r) denotes the space of skew-symmetric real r × r matrices. Consider the
extended tangent map of (S, U1 , . . . , UN ) → Y = S Xn=1 Un ,
N


N 
N
Rr1 ×···×rN × TUn VIn ,rn → TY Mr × so(rn ),
n=1 n=1
 
N 
N
(Ṡ, U̇1 , . . . , U̇N ) → Ṡ X Un + S ×n U̇n X Uk , UT1 U̇1 , . . . , UTN U̇N .
n=1 k=n
n=1

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 DYNAMICAL TENSOR APPROXIMATION 2363

This linear map turns out to be an isomorphism, since its inverse can be constructed
explicitly by an argument similar to the proof of (2.7) below. Hence, every tangent
tensor Ẏ ∈ TY Mr is of the form
Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

N 
N
(2.5) Ẏ = Ṡ X Un + S ×n U̇n X Uk ,
n=1 k=n
n=1

where Ṡ ∈ Rr1 ×···×rN and U̇n ∈ TUn VIn ,rn . Moreover, it follows that Ṡ and U̇n
are uniquely determined by Ẏ and the chosen S and Un in (2.3) if we impose the
orthogonality constraints

(2.6) UTn U̇n = 0 .

We now show that then Ṡ and U̇n are actually given by the following formulas:
N
Ṡ = Ẏ X UTn ,
(2.7) n=1

U̇n = P⊥
n Ẏ X UTk S†(n)
k=n (n)

† −1
with the projection P⊥
n = I − Un Un and the pseudoinverse S(n) = S(n) S(n) S(n)
T T T
.
The rn × rn matrix S(n) ST(n) is invertible since S(n) is of rank rn by assumption.
Proof of (2.7). When we multiply (2.5) by Xn=1 UTn , then we obtain by (2.2)
N

and the orthogonality relations UTn Un = I and (2.6) the first equation of (2.7). For
the second equation, we multiply (2.5) by Xk=n UTk , and again using (2.2) and the
orthogonality relations we obtain

(2.8) Ẏ X UTk = Ṡ ×n Un + S ×n U̇n .

k=n

By the equation for Ṡ and once again (2.2), the first expression on the right-hand side
becomes

Ṡ ×n Un = Ẏ ×n (Un UTn ) X UTk ,

k=n

or in its n-mode unfolding,

Un Ṡ(n) = Un UTn Ẏ X UTk .

k=n (n)

Taking the n-mode unfolding on both sides of (2.8) and rearranging terms then gives

U̇n S(n) = I − Un UTn Ẏ X UTk .

k=n (n)

Multiplying this relation from the right with S†(n) yields the second equation of (2.7).
2.4. Dynamical tensor approximation. We now turn to the approxima-
tion (1.2) of time-dependent tensors A(t). The minimization condition (1.2) on
the tangent space is equivalent to an orthogonal projection: find Ẏ ∈ TY Mr

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 2364 OTHMAR KOCH AND CHRISTIAN LUBICH

(we omit the argument t) satisfying

(2.9) Ẏ − Ȧ, V = 0 for all V ∈ TY Mr .

Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

This is a Galerkin condition on the tangent space TY Mr . With this formulation we

derive differential equations for the factors in the Tucker decomposition (2.3).
Theorem 2.1. For a tensor Y = S Xn=1 Un ∈ Mr with n-rank rn for n =
N

r1 ×···×rN
1, . . . , N , with core tensor S ∈ R and n-mode factors Un ∈ RIn ×rn having
orthonormal columns, condition (1.2) or (2.9) is equivalent to

N 
N
(2.10) Ẏ = Ṡ X Un + S ×n U̇n X Uk ,
n=1 k=n
n=1

where the factors in the decomposition satisfy the system of differential equations
N
Ṡ = Ȧ X UTn ,
(2.11) n=1

U̇n = ⊥
Pn Ȧ X UTk S†(n)
k=n (n)

with the projection P⊥n = I − Un Un onto the orthogonal complement of the range of
T
† −1
Un and with the pseudoinverse S(n) = ST(n) S(n) ST(n) of the n-mode unfolding S(n)
of S.
Equations (2.11) are formally like (2.7), with Ẏ replaced by Ȧ.
The differential equations (2.11) are solved numerically, starting from an approx-
imation to the tensor A(0) at the initial time given in the low-rank Tucker format.
Notice that a time step with these differential equations requires no decompositions of
large matrices (the matrices S(n) ST(n) are of small dimension rn × rn ). The main com-
putational work is in the computation of the dimension-contracting n-mode products
appearing on the right-hand side of (2.11). These require inner products of length
In and can exploit possible sparsity in Ȧ(t) so that only the actually time-varying
entries of A(t) are addressed.
Proof. We know from section 2.3 that Ẏ can be written in the form (2.10) with
U̇n satisfying the orthogonality condition (2.6). Choosing V = T Xn=1 Un ∈ TY Mr
N

r1 ×···×rN
with an arbitrary T ∈ R , we have
N
 N

Ȧ, T X Un = Ȧ X UTn , T ,
n=1 n=1
N

Ẏ, T X Un = Ṡ, T .
n=1

Since this holds for every T ∈ Rr1 ×···×rN , (2.9) yields the first equation of (2.11).
We now choose V = S ×n Vn Xk=n Uk , which is in TY Mr if Vn satisfies the
orthogonality relation UTn Vn = 0. We then have
   
Ȧ, S ×n Vn X Uk = Ȧ X UTk , S ×n Vn
k=n k=n
   
= Ȧ X UTk , Vn S(n) = Ȧ X UTk ST(n) , Vn ,
k=n (n) k=n (n)

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 DYNAMICAL TENSOR APPROXIMATION 2365

where the matrix inner product in the second line is the Frobenius inner product. On
the other hand we have, using (2.10) and the orthogonality relations,
Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

 
Ẏ, S ×n Vn X Uk = Ṡ ×n Un + S ×n U̇n , S ×n Vn 
k=n

= Un Ṡ(n) + U̇n S(n) , Vn S(n) 

= Ṡ(n) , UTn Vn S(n)  + U̇n S(n) ST(n) , Vn 
= U̇n S(n) ST(n) , Vn  .

By (2.9), we thus have

 
U̇n S(n) ST(n) − Ȧ X UTk ST(n) , Vn = 0
k=n (n)

for all Vn with UTn Vn = 0, that is, for all matrices Vn = (I − Un UTn )Wn with an
arbitrary matrix Wn ∈ RIn ×rn . We therefore obtain

(I − Un UTn )U̇n S(n) ST(n) = (I − Un UTn ) Ȧ X UTk ST(n) .

k=n (n)

Since UTn U̇n = 0 by condition (2.6), this yields the second equation of (2.11).
Remark 2.1. Related differential equations were derived earlier in the chemical
physics literature in [2] for the multiconfiguration time-dependent Hartree method
of quantum mechanics, where an approximation to the multivariate wave function is
sought for in the form of a linear combination of products of univariate functions.
Remark 2.2. In this paper we deal only with the Tucker format of tensor approx-
imation. Another familiar format is the canonical format (for tensors of canonical
rank K)


K
(1) (N )
A(t) ≈ vk (t) ◦ · · · ◦ vk (t),
k=1

(n)
where no orthogonality relation is required between the vectors vk (t). While this
approach appears attractive in that it has no built-in exponential scaling with the
order N of the tensor, it leads to a number of theoretical difficulties even in the time-
independent approximation problem [15] and does not appear to lend itself to the
dynamical tensor approximation approach. The set of all tensors of canonical rank
K does not have a manifold structure that allows us to give differential equations for
(n)
the vectors vk (t) such that the analogue of (1.2) is fulfilled.
Remark 2.3. The dynamical tensor approximation approach for the Tucker format
can be extended to the hierarchical tensor format described in [2, 4, 12, 18], which has
only cubic scaling with N . This is, however, outside the scope of the present paper.
3. Tangent space projection and curvature bounds. In operator notation,
(2.9) can be written as

(3.1) Ẏ = P(Y)Ȧ ,

where P(Y) is the orthogonal projection onto the tangent space TY Mr . For a the-
oretical understanding of the dynamical low-rank approximation, an analysis of the

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 2366 OTHMAR KOCH AND CHRISTIAN LUBICH

properties of this projection operator is therefore essential. We begin by giving an

explicit formula.
Lemma 3.1. Let Y = S Xn=1 Un ∈ Mr with n-mode factors Un having orthonor-
N
Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

mal columns. With the orthogonal projections

Pn = Un UTn , P⊥n = I − Pn ,
 †

Qn = Uk S(n) S(n) UTk ,
k=n k=n

the orthogonal projection onto the tangent space TY Mr is given by

N 
N (n)
(3.2) P(Y)B = B X Pn + P⊥
n B(n) Qn
n=1
n=1

for B ∈ RI1 ×···×IN , with B(n) = [B](n) the nth unfolding of B.

Proof. Theorem 2.1 gives an expression for Ẏ determined by (2.9), that is, for
Ẏ = P(Y)Ȧ. With B in place of Ȧ, this reads

N N 
N    
P(Y)B = B X UTn X Un + S ×n P⊥
n B X UTk S†
(n) (n)
X Uk
n=1 n=1 k=n k=n
n=1
N N 
   (n)
=B X Pn + P⊥
n B(n) Uk S†(n) S(n) UTk ,
n=1
n=1 k=n k=n

where we have used (2.2) and (2.4). This gives us (3.2).

The next lemma is the key tool for the approximation results of the following
section.
Lemma 3.2. There are constants c and C (depending only on the order N and
satisfying Cc ≤ 12 ) such that the following holds true. Let the rank-(r1 , . . . , rN ) tensor
X ∈ Mr be such that the smallest nonzero singular value of the nth unfolding satisfies
σrn (X(n) ) ≥ ρ > 0 for n = 1, . . . , N , and let X̃ ∈ Mr with X̃ − X ≤ cρ. Then, for
all B ∈ RI1 ×···×IN ,

(3.3)  P(X̃) − P(X) B ≤ Cρ−1 X̃ − X · B,

(3.4) P⊥ (X)(X̃ − X) ≤ Cρ−1 X̃ − X2 ,

where P⊥ (X) = I − P(X) is the projection onto the orthogonal complement of the
tangent space TX Mr .
Proof. (a) Writing X ∈ Mr in Tucker form as X = S Xn=1 Un , we note that the
N

nonzero singular values of X(n) are those of S(n) , and hence

1 1 1
S†(n)  = = ≤ .
σrn (S(n) ) σrn (X(n) ) ρ

Since we have, by [6, p. 448],

|σrn (X̃(n) ) − σrn (X(n) )| ≤ X̃(n) − X(n) 2 ≤ X̃(n) − X(n) F = X̃ − X,

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 DYNAMICAL TENSOR APPROXIMATION 2367

we obtain for X̃ − X ≤ 12 ρ that

σrn (X̃(n) ) ≥ σrn (X(n) ) − |σrn (X̃(n) ) − σrn (X(n) )| ≥ 12 ρ ,

Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

†
and hence X̃ = S̃ Xn=1 Ũn with S̃(n)  ≤ 2ρ−1 .
N

(b) We decompose the tensors on the straight line connecting X and X̃ as

X + τ (X̃ − X) = Y(τ ) + Z(τ ) with Y(τ ) ∈ Mr , Z(τ ) ⊥ TX Mr .

A smooth such decomposition exists at least for small τ , but the arguments below
show that it exists in fact for 0 ≤ τ ≤ 1. We denote

Δ = P(X)(X̃ − X) ∈ TX Mr with Δ ≤ δ := X̃ − X.

We then have P(X)(Y(τ ) − X) = τ Δ, which yields

P(X)Ẏ(τ ) = Δ.

Since P(Y)Ẏ = Ẏ, we have

Ẏ(τ ) = Δ + P(Y(τ )) − P(X) Ẏ(τ ) .

(c) As long as the operator norm of P(Y(τ )) − P(X) is bounded by 12 , we thus

have

Ẏ ≤ 2δ .

This yields Y(τ ) − X ≤ 2δ for 0 ≤ τ ≤ 1, and for the Tucker factors of Y(τ ) =
S(τ ) Xn=1 Un (τ ) with UTn U̇n = 0, we therefore have by (2.7)
N

Ṡ ≤ 2δ , U̇n S(n)  ≤ 2δ .

If 2δ ≤ 12 ρ, then the argument of (a) applied to Y(τ ) instead of X̃ shows that

S†(n) (τ ) ≤ 2ρ−1 , and hence

U̇n  ≤ U̇n S(n)  · S†(n)  ≤ 4δρ−1 =: γ .

With this estimate we further obtain for Pn = Un UTn that

Ṗn  ≤ 2γ .

Using the product rule for dτ d

ST(n) (S(n) ST(n) )−1 S(n) and the estimates for the norms
of Ṡ(n) and S†(n) by 2δ and 2ρ−1 , respectively, we find that the norm of this derivative
is bounded by 4γ, and hence we have for the projection Qn of Lemma 3.1

Q̇n  ≤ 2(N − 1)γ + 4γ = 2(N + 1)γ .

τ
Expressing Un (τ ) − Un (0) = 0 U̇n (s) ds and similarly for the increments of S, Pn ,
and Qn , we obtain from the formula of Lemma 3.1 and the above estimates

 P(Y(τ )) − P(X) B ≤ N · 2γ + N · (2γ + 2(N + 1)γ) τ B

= 2N (N + 3)γ · τ B = 8N (N + 3) δρ−1 τ B .

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 2368 OTHMAR KOCH AND CHRISTIAN LUBICH

1
The operator norm of P(Y(τ )) − P(X) thus does not exceed 2 for 0 ≤ τ ≤ 1 if

1
Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

δ ≤ cρ with c = ,
16 N (N + 3)

and at τ = 1 we then obtain the bound (3.3) with C = 8N (N + 3).

(d) We note
 1  1
P⊥ (X)(X̃ − X) = P⊥ (X) Ẏ(τ ) dτ = P(Y(τ )) − P(X) Ẏ(τ ) dτ .
0 0

By the above estimates, this is bounded by

P⊥ (X)(X̃ − X) ≤ 8N (N + 3) ρ−1 δ 2 ,

which yields (3.4).

4. Approximation properties. We give approximation results that are exten-
sions to tensors of results for the matrix case in [7].
4.1. Local quasi-optimality. If the low-rank approximation problem (1.1) has
a continuously differentiable best approximation X(t) ∈ Mr , then the error of (1.2)
can be bounded in terms of the best-approximation error X(t) − A(t). The result
requires a bound on Ȧ(t):

(4.1) Ȧ(t) ≤ μ for 0 ≤ t ≤ t.

Theorem 4.1. Suppose that a continuously differentiable best approximation

X(t) ∈ Mr to A(t) exists for 0 ≤ t ≤ t. Let ρ > 0 be such that the smallest nonzero
singular value of the nth unfolding of X(t) satisfies σrn (X(n) (t)) ≥ ρ for n = 1, . . . , N ,
and assume that the best-approximation error is bounded by X(t) − A(t) ≤ cρ for
0 ≤ t ≤ t, with c of Lemma 3.2. Then, the approximation error of the dynamical
low-rank approximation (1.2) with initial value Y(0) = X(0) is bounded by
 t
Y(t) − X(t) ≤ 2β e βt
X(s) − A(s) ds with β = Cμρ−1
0

for t ≤ t and as long as the right-hand side remains bounded by cρ. Here, C and c
are the constants of Lemma 3.2.
Proof. With Lemma 3.2 at hand, the proof is essentially the same as that of
Theorem 5.1 in [7]. The error of the best approximation X − A is orthogonal to
the tangent space TX Mr , or, equivalently, P(X)(X − A) = 0. We differentiate this
relation with respect to t and denote (P (X) · B)Ẋ = dt
d
P(X(t))B to obtain

P(X)(Ẋ − Ȧ) + P (X) · (X − A) Ẋ = 0.

Since Ẋ ∈ TX Mr , we have P(X)Ẋ = Ẋ, and the equation becomes

(4.2) I − P (X) · (X − A) Ẋ = P(X)Ȧ.

Lemma 3.2 and the condition d := X − A ≤ cρ yield

P (X) · (X − A) ≤ Cρ−1 d ≤ Cc ≤ 12 ,

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 DYNAMICAL TENSOR APPROXIMATION 2369

and hence (4.2) can be solved for Ẋ to yield

Ẋ = P(X)Ȧ + D with D ≤ 2Cρ−1 dμ = 2βd.

Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

We subtract this equation from (3.1), viz., Ẏ = P(Y)Ȧ, and integrate from 0 to t. As
long as e(t) := Y(t) − X(t) ≤ cρ, Lemma 3.2 yields

 P(Y) − P(X) Ȧ ≤ Cρ−1 eμ = βe,

and hence we obtain

 t  t
e(t) ≤ β e(s) ds + 2β d(s) ds.
0 0

The result now follows with the Gronwall inequality.

4.2. A farther-reaching error bound. Smaller errors over longer time inter-
vals are obtained if X − A and its derivative are small. We assume that A(t) is of
the form

(4.3) A(t) = X(t) + E(t), 0 ≤ t ≤ t,

where X(t) ∈ Mr (now this need not necessarily be the best approximation) with

(4.4) Ẋ(t) ≤ μ

and the derivative of the remainder term is bounded by

(4.5) Ė(t) ≤ ε

with a small ε > 0.

Theorem 4.2. In addition to the above assumptions, suppose that the smallest
singular values of the unfoldings of X(t) are bounded from below by ρ > 0. Then, the
approximation error of (1.2) with initial value Y (0) = X(0) is bounded by
ρ
Y(t) − X(t) ≤ 2tε for t≤ √ ,
C με

provided that t ≤ cρ
2ε and t ≤ t. The constants C and c are those of Lemma 3.2.
Proof. The proof follows that of Theorem 5.2 in [7]. We note Ẋ = P(X)Ẋ, rewrite
(3.1) as Ẏ = P(Y)Ẋ + P(Y)Ė, and subtract the two equations. We observe

P(Y) − P(X) Ẋ = − P⊥ (Y) − P⊥ (X) Ẋ = −P⊥ (Y)Ẋ = −P⊥ (Y)2 Ẋ.

We take the inner product with Y − X to obtain

Y − X, P(Y) − P(X) Ẋ = −Y − X, P⊥ (Y)Ẋ = −P⊥ (Y)(Y − X), P⊥ (Y)Ẋ

= P⊥ (Y)(Y − X), P(Y) − P(X) Ẋ.

With Lemma 3.2 and (4.4), (4.5) this equation yields

Y − X, Ẏ − Ẋ = P⊥ (Y)(Y − X), P(Y) − P(X) Ẋ + Y − X, P(Y)Ė

≤ C 2 ρ−2 μ Y − X3 + Y − X · ε,

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 2370 OTHMAR KOCH AND CHRISTIAN LUBICH

and on the other hand we have

1 d d
Y − X, Ẏ − Ẋ = Y − X2 = Y − X Y − X.
Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

2 dt dt
Taken together, we obtain for e(t) = Y(t) − X(t) the differential inequality
ė ≤ γe2 + ε, e(0) = 0,
with γ = C 2 ρ−2 μ. Hence, e(t) is majorized by the solution of
ẏ = γy 2 + ε, y(0) = 0,
 √ √
which equals y(t) = ε/γ tan(t γε) and is bounded by 2tε for t γε ≤ 1. Lemma 3.2
remains applicable as long as 2tε ≤ cρ.
4.3. Systems without gaps between the singular values. The results of
the preceding subsections only give satisfactory error bounds when there is a gap in
the distribution of the singular values of the unfoldings so that essential and inessential
singular values are widely separated. We now consider a situation where such a gap
need not exist. We make the assumptions of Theorem 4.2 and further that X(t) ∈ Mr
with all singular values of all unfoldings greater than ρ > 0 has a Tucker decomposition
N
(4.6) X(t) = S(t) X Un (t)
n=1

such that the following bounds hold for n = 1, . . . , N and 0 ≤ t ≤ t :

 
 † ˙  ˙ (t) ≤ ν .
(4.7) S(n) S(n)  ≤ κ , U n

Under these conditions we can show an O(ε) error over times O(1) even with ρ ∼ ε.
Theorem 4.3. Under the conditions of Theorem 4.2 and with (4.6)– (4.7), the
approximation error of (1.2) with initial value Y(0) = X(0) is bounded by

1 ρ
Y(t) − X(t) ≤ 2tε for t ≤ √
c1 κ + c2 ν ε
with constants c1 and c2 that depend only on the order N . This holds as long as the
right-hand side remains bounded by cρ with c of Lemma 3.2.
Proof. From the proof of Theorem 4.2 we have the equation
(4.8) Y − X, Ẏ − Ẋ = −P⊥ (Y)(Y − X), P⊥ (Y)Ẋ + Y − X, P(Y)Ė.
For e(t) = Y(t)− X(t) ≤ cρ, the proof of Lemma 3.2 shows that there is a homotopy
N
Y(t, τ ) = S(t, τ ) X Un (t, τ )
n=1

with
Y(t, 1) = Y(t), Y(t, 0) = X(t) ,
for which the Tucker factors are bounded by
 ∂S 
 
 (t, τ ) ≤ 2e(t),
∂τ
 ∂U 
 n 
(4.9)  (t, τ ) S(n) (t, τ ) ≤ 2e(t),
∂τ
 ∂U 
 
(t, τ ) ≤ 4e(t)ρ−1 ≤ 4c .
n

∂τ

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 DYNAMICAL TENSOR APPROXIMATION 2371

We write

N 
N
Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

⊥ ⊥
P (Y(t))Ẋ(t) = P (Y(t)) Ṡ(t, 0) X Un (t, 0) + S(t, 0) ×n U̇n (t, 0) X Uk (t, 0)
n=1 k=n
n=1

N 
N
− Ṡ(t, 0) X Un (t, 1) + S(t, 1) ×n U̇n (t, 0) X Uk (t, 1) ,
n=1 k=n
n=1

noting that P⊥ (Y(t)) annihilates the terms in the second line, since they are in TY Mr .
We then have
 
 N N 
 
Ṡ(t, 0) X Un (t, 1) − Ṡ(t, 0) X Un (t, 0)
 n=1 n=1 
   
 1 ∂ N 
 
= Ṡ(t, 0) X Un (t, τ ) ∂τ 
 0 ∂τ n=1 
 N   
 1 
 ∂Un 
= Ṡ(t, 0) ×n (t, τ ) X Uk (t, τ ) dτ 
 0 ∂τ k=n 
n=1
N  1
 
  
 ∂Un 
≤  (t, τ ) Ṡ(n) (t, 0) Uk (t, τ )T  dτ
0  ∂τ 
n=1 k=n
   
N  1   
 ∂Un   † 
≤  (t, τ ) S(n) (t, 0) · S(n) (t, 0) Ṡ(n) (t, 0) dτ .
0  ∂τ   
n=1

Writing
 τ
∂Un ∂Un ∂Un ∂S(n)
(t, τ ) S(n) (t, 0) = (t, τ ) S(n) (t, τ ) − (t, τ ) (t, σ) dσ
∂τ ∂τ ∂τ 0 ∂σ
and using the bounds (4.7) and (4.9) thus yields
 
 N N 
 
Ṡ(t, 0) X Un (t, 1) − Ṡ(t, 0) X Un (t, 0) N ≤ 2e(t) + 8ce(t) κ = c̃1 κe(t)
 n=1 n=1 

with c̃1 = N · (2 + 8c). In the same way, the remaining terms in the above expression
for P⊥ (Y(t))Ẋ(t) are bounded by c̃2 νe(t), and hence we have

P⊥ (Y(t))Ẋ(t) ≤ (c̃1 κ + c̃2 ν)e(t) .

Using this inequality in (4.8) and also the second estimate of Lemma 3.2, we obtain
the differential inequality, as long as e(t) ≤ cρ,

(4.10) ė ≤ γe2 + ε

with γ = ρ−1 C(c̃1 κ + c̃2 ν). With c1 = Cc̃1 and c2 = Cc̃2 , the result now follows as at
the end of the proof of Theorem 4.2.
4.4. Low-rank approximation of tensor differential equations. For the
low-rank approximation to a solution of the tensor differential equation

(4.11) Ȧ = F (A),

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 2372 OTHMAR KOCH AND CHRISTIAN LUBICH

condition (1.2) is replaced, at every time t, by

(4.12) Ẏ ∈ TY Mr such that Ẏ − F (Y) = min!.

Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

Equivalently, condition (2.9) is replaced by the Galerkin condition

(4.13) Ẏ − F (Y), V = 0 for all V ∈ TY Mr ,

and correspondingly, the expression Ȧ is replaced by F (Y) for Y = S Xn=1 UN ∈ Mr

N

in the differential equations (2.11).

Theorems 4.1–4.3 extend to the low-rank approximation of matrix differential
equations (4.11). We assume that F has a moderate bound along the approximations,

(4.14) F (X(t)) ≤ μ, F (Y(t)) ≤ μ for 0 ≤ t ≤ t,

and satisfies a one-sided Lipschitz condition: there is a real λ (positive or negative or

zero) such that

(4.15) F (Y) − F (X), Y − X ≤ λ Y − X2

for all tensors X, Y ∈ Mr . We further assume that for the best approximation X(t),

(4.16) F (X(t)) − F (A(t)) ≤ L X(t) − A(t) for 0 ≤ t ≤ t,

which is, in particular, satisfied if F is Lipschitz continuous with Lipschitz constant L.

Furthermore, assume that the best-approximation error is bounded by X(t)−A(t) ≤
cρ for 0 ≤ t ≤ t, with c of Lemma 3.2.
We then have the following extension of the quasi-optimality result of Theo-
rem 4.1.
Theorem 4.4. Suppose that a continuously differentiable best approximation
X(t) ∈ Mr to a solution A(t) of (4.11) exists for 0 ≤ t ≤ t, and assume the bounds
(4.14)–(4.16). Let ρ > 0 be such that the smallest nonzero singular value of the nth
unfolding of X(t) satisfies σrn (X(n) (t)) ≥ ρ for n = 1, . . . , N , and assume that the
best-approximation error is bounded by X(t) − A(t) ≤ cρ with c of Lemma 3.2 for
0 ≤ t ≤ t. Then, the approximation error of (4.13) with initial value Y(0) = X(0) is
bounded by
 t
(2β+λ)t
Y(t) − X(t) ≤ (2β + L) e X(s) − A(s) ds with β = Cμρ−1
0

for t ≤ t and as long as the right-hand side is bounded by cρ.

Proof. Equation (4.13) rewritten as in (3.1) reads

(4.17) Ẏ = P(Y)F (Y).

As in the proof of Theorem 4.1, we have the equation

Ẋ = P(X)F (A) + D with D ≤ 2βd

for d = X − A. We subtract the two equations, write

P(Y)F (Y) − P(X)F (A) − D = (P(Y) − P(X))F (X) + P(X)(F (X) − F (A))
+ (F (Y) − F (X)) − P⊥ (Y)(F (Y) − F (X)) − D,

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 DYNAMICAL TENSOR APPROXIMATION 2373

and take the inner product with Y − X. With Lemma 3.2 we obtain

Ẏ − Ẋ, Y − X ≤ β Y − X2 + Ld Y − X

Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

+ λ Y − X2 + β Y − X2 + 2βd Y − X.

For e = Y − X this gives the differential inequality

(4.18) ė ≤ (2β + λ)e + (2β + L)d, e(0) = 0,

which yields the result.

We refer to [11, Theorem 4.1] for a related quasi-optimality result in a situation
of a linear differential equation with an unbounded operator.
In the differential equation analogue of Theorem 4.2 with the splitting (4.3), we
start from the equations Ẏ − Ẋ = P(Y)F (Y) − P(X)Ẋ and Ẋ = F (A) − Ė, yielding

Ẏ − Ẋ = (P(Y) − P(X))Ẋ − P⊥ (Y)(F (Y) − F (X))

+ (F (Y) − F (X)) + P(Y)(F (X) − F (A)) + P(Y)Ė,

where we now take the inner product with Y − X. If it is additionally assumed that
F has Lipschitz constant L, then this leads to the differential inequality

(4.19) ė ≤ Cρ−1 (β + L)e2 + λe + Ld + ε, e(0) = 0.

With γ = Cρ−1 (β + L) and ε = ε + L max0≤t≤t d(t), and with ϕ(x) = (ex − 1)/x, this
yields the error bound

(4.20) Y(t) − X(t) ≤ 2t ϕ(λt) ε γ ε)−1/2

for tϕ(λt) ≤ 12 (

as long as t ≤ t and 2tϕ(λt)ε ≤ cρ.

Theorem 4.3 is similarly extended to tensor differential equations.
5. Numerical experiment. To illustrate the theoretical error analysis, we
compute the dynamical low-rank approximation of a time-dependent tensor
A(t) ∈ R15×15×15×15 , which is constructed to give control over the order of magnitude
of the approximation error: with U = (U1 , . . . , U4 ) a set of matrices in R15×10 with or-
thonormal columns, B a random core tensor in R10×10×10×10 , and C ∈ R15×15×15×15
a random perturbation, the time-dependent data tensor is constructed as

(5.1) A(t) = exp(t)B ×1 U1 ×2 · · · ×4 U4 + ε(t + 1 + sin(3t))C, t ∈ [0, 1].

We compute the dynamical tensor approximation (2.9), where Y is defined with a

core tensor of dimension 10 × 10 × 10 × 10. We numerically solve the differential
equations (2.11), taking as initial value an approximate best approximation computed
by alternating least squares [1]. As time integrator, we use the implicit midpoint rule
(yn+1 = yn + hf 12 (yn + yn+1 ) for a differential equation y  = f (y)) with step size
h = 10−4 , where the nonlinear equations in each step are solved by fixed point iteration
with a stopping criterion in a norm combining the tensor norm for the core tensor
and the Frobenius norms of the factor matrices. Figure 5.1 shows the approximation
errors Y(t)−A(t) as functions of t for (5.1) with ε ∈ {10−5, 10−4 , 10−3 , 10−2 , 10−1 }.
We observe that indeed the error of the approximation is proportional to ε and grows
only moderately as a function of t (note the logarithmic scaling of the vertical axis).
To demonstrate the potential computational advantage of dynamical tensor ap-
proximation versus pointwise approximation, we compare the computation time for

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 2374 OTHMAR KOCH AND CHRISTIAN LUBICH

100

eps=1e−1
Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

eps=1e−2
10 eps=1e−3
eps=1e−4
eps=1e−5

1
Error

1e−1

1e−2

1e−3

0 0.2 0.4 0.6 0.8 1

t

Fig. 5.1. Error of the approximation of rank (10, 10, 10, 10) for A from (5.1) with ε ∈
{10−5 , 10−4 , 10−3 , 10−2 , 10−1 }. Interval t ∈ [0, 1], step size h = 10−4 .

1800

1600 DLRTen
ALS
1400

1200
computation time

1000

800

600

400

200

0
0 2000 4000 6000 8000 10000 12000
time step

Fig. 5.2. Comparison of the computation times for dynamical low-rank tensor approximation
(DLRTen) and pointwise approximation by alternating least squares (ALS). Approximations of rank
(10, 10, 10, 10) for A from (5.1) with ε = 10−2 . Interval t ∈ [0, 1], step size h = 10−4 .

our dynamical tensor approximation to pointwise approximation by alternating least

squares as implemented in the MATLAB routines [1]. In the same setting as above,
the tensor from (5.1) (with ε = 10−2 ) is approximated by a tensor with a core tensor
of dimension 10 × 10 × 10 × 10. Figure 5.2 shows that the computation time is about
halved by adopting our approach. In other examples that we computed, this gain is
larger when the size of the core tensor is smaller as compared to the original tensor
size or when Ȧ is significantly sparser than A, as has similarly been observed in the
dynamical low-rank approximation of large data matrices [14].
For further numerical examples and comparisons, we refer to [14], where in par-
ticular the dynamical low-rank approximation of the 3-tensors arising from the spatial
discretization of a 3-dimensional reaction-diffusion PDE is considered.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 DYNAMICAL TENSOR APPROXIMATION 2375

Acknowledgment. We thank Tamara Kolda for helpful comments on a prelim-

inary version of this paper and for useful hints to the literature.
Downloaded 03/19/24 to [Link] . Redistribution subject to SIAM license or copyright; see [Link]

REFERENCES

[1] B. W. Bader and T. G. Kolda, MATLAB Tensor Toolbox Version 2.3, 2009, available online
at [Link]
[2] M. H. Beck, A. Jäckle, G. A. Worth, and H.-D. Meyer, The multiconfiguration time-
dependent Hartree (MCTDH) method: A highly efficient algorithm for propagating
wavepackets, Phys. Rep., 324 (2000), pp. 1–105.
[3] L. Grasedyck, Existence and computation of low Kronecker-rank approximations for large
linear systems of tensor product structure, Computing, 72 (2004), pp. 247–265.
[4] L. Grasedyck, Hierarchical Singular Value Decomposition of Tensors, preprint 20,
DFG SPP 1324, 2009, available online from [Link]
preprints/[Link].
[5] W. Hackbusch and B. N. Khoromskij, Tensor-product approximation to operators and func-
tions in high dimensions, J. Complexity, 23 (2007), pp. 697–714.
[6] R. A. Horn and C. R. Johnson, Matrix Analysis, Cambridge University Press, Cambridge,
UK, 1985.
[7] O. Koch and C. Lubich, Dynamical low-rank approximation, SIAM J. Matrix Anal. Appl., 29
(2007), pp. 434–454.
[8] T. G. Kolda and B. W. Bader, Tensor Decompositions and Applications, Sandia Report
SAND2007-6702, Sandia National Laboratories, Albuquerque, NM and Livermore, CA,
2007.
[9] L. de Lathauwer, B. de Moor, and J. Vandewalle, A multilinear singular value decompo-
sition, SIAM J. Matrix Anal. Appl., 21 (2000), pp. 1253–1278.
[10] L. de Lathauwer, B. de Moor, and J. Vandewalle, On the best rank-1 and rank-
(R1 , . . . , Rn ) approximation of higher-order tensors, SIAM J. Matrix Anal. Appl., 21
(2000), pp. 1324–1342.
[11] C. Lubich, On variational approximations in quantum molecular dynamics, Math. Comp., 74
(2005), pp. 765–779.
[12] C. Lubich, From Quantum to Classical Molecular Dynamics: Reduced Models and Numerical
Analysis, Zur. Lect. Adv. Math., European Mathematical Society, Zürich, 2008.
[13] H.-D. Meyer, F. Gatti, and G. A. Worth, eds., Multidimensional Quantum Dynamics:
MCTDH Theory and Applications, Wiley-VCH, Weinheim, Berlin, 2009.
[14] A. Nonnenmacher and C. Lubich, Dynamical low rank approximation: Applications and
numerical experiments, Math. Comput. Simulation, 79 (2009), pp. 1346–1357.
[15] V. de Silva and L.-H. Lim, Tensor rank and the ill-posedness of the best low-rank approxima-
tion problem, SIAM J. Matrix Anal. Appl., 30 (2008), pp. 1084–1127.
[16] J. Sun, S. Papadimitriou, and P. S. Yu, Window-based tensor analysis on high-dimensional
and multi-aspect streams, in ICDM ’06: Proceedings of the Sixth International Conference
on Data Mining, IEEE Computer Society, Washington, DC, 2006, pp. 1076–1080.
[17] J. Sun, D. Tao, and C. Faloutsos, Beyond streams and graphs: Dynamic tensor analysis, in
KDD ’06: Proceedings of the 12th ACM SIGKDD International Conference on Knowledge
Discovery and Data Mining, ACM, New York, 2006, pp. 374–383.
[18] H. Wang and M. Thoss, Multilayer formulation of the multiconfiguration time-dependent
Hartree theory, J. Chem. Phys., 119 (2003), pp. 1289–1299.
[19] T. Zhang and G. H. Golub, Rank-one approximation to high order tensors, SIAM J. Matrix
Anal. Appl., 23 (2001), pp. 534–550.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.