A literature survey of low-rank tensor approximation techniques

Introduction

This survey is concerned with tensors in the sense of multidimensional arrays. A general tensor of order $d$ and size $n_{1}\times n_{2}\times\cdots\times n_{d}$ for integers $n_{1},n_{2},\ldots,n_{d}$ will be denoted by

An entry of $\mathcal{X}$ is denoted by $\mathcal{X}_{i_{1},\ldots,i_{d}}$ where each index $i_{\mu}\in\{1,\ldots,n_{\mu}\}$ refers to the $\mu$th mode of the tensor for $\mu=1,\ldots,d$. For simplicity, we will assume that $\mathcal{X}$ has real entries, but it is of course possible to define complex tensors or, more generally, tensors over arbitrary fields.

A wide variety of applications lead to problems where the data or the desired solution can be represented by a tensor. In this survey, we will focus on tensors that are induced by the discretization of a multivariate function; we refer to the survey and to the books for the treatment of tensors containing observed data. The simplest way a given multivariate function $f(x_{1},x_{2},\ldots,x_{d})$ on a tensor product domain $\Omega=^{d}$ leads to a tensor is by sampling $f$ on a tensor grid. In this case, each entry of the tensor contains the function value at the corresponding position in the grid. The function $f$ itself may, for example, represent the solution to a high-dimensional partial differential equation (PDE).

As the order $d$ increases, the number of entries in $\mathcal{X}$ increases exponentially for constant $n=n_{1}=\cdots=n_{d}$. This so called curse of dimensionality prevents the explicit storage of the entries except for very small values of $d$. Even for $n=2$, storing a tensor of order $d=50$ would require 9 petabyte! It is therefore essential to approximate tensors of higher order in a compressed scheme, for example, a low-rank tensor decomposition. Various such decompositions have been developed, see Section 2. An important difference to tensors containing observed data, a tensor $\mathcal{X}$ induced by a function is usually not given directly but only as the solution of some algebraic equation, e.g., a linear system or eigenvalue problem. This requires the development of solvers for such equations working within the compressed storage scheme. Such algorithms are discussed in Section 3.

The range of applications of low-rank tensor techniques is quickly expanding. For example, they have been used for addressing:

the approximation of parameter-dependent integrals , multi-dimensional integrals , and multi-dimensional convolution ;

computational tasks in electronic structure calculations, e.g., based on Hartree-Fock or DFT models ;

the solution of stochastic and parametric PDEs ;

approximation of Green’s functions of high-dimensional PDEs ;

the solution of the Boltzmann equation , and the chemical master / Fokker-Planck equations

the solution of high-dimensional Schrödinger equations ;

the approximation of stationary states of stochastic automata networks ;

multivariate regression and machine learning .

Note that the list above mainly focuses on techniques that involve tensors; low-rank matrix techniques, such as POD and reduced basis methods, have been applied in an even broader setting.

A word of caution is appropriate. Even though the field of low-rank tensor approximation is relatively young, it already appears a daunting task to give proper credit to all developments in this area. This survey is biased towards work related to the TT and hierarchical Tucker decompositions. Surveys with a similar scope are the lecture notes by Grasedyck , Khoromskij , and Schneider , as well as the monograph by Hackbusch . Less detailed attention may be given to other developments, although we have made an effort to at least touch on all important directions in the area of function-related tensors.

Low-rank tensor decompositions

As mentioned in the introduction, it will rarely be possible to store all entries of a higher-order tensor explicitly. Various compression schemes have been developed to reduce storage requirements. For $d=2$ all these schemes boil down to the well known reduced singular value decomposition (SVD) of matrices ; however, they differ significantly for tensors of order $d\geq 3$.

The entries of a rank-one tensor $\mathcal{X}$ can be written as

By defining the vectors u^{(\mu)}:=\big{(}u_{1}^{(\mu)},\ldots,u_{n_{\mu}}^{(\mu)}\big{)}^{T}, a more compact from of this relation is given by

where $\otimes$ denotes the usual Kronecker product and $\operatorname{vec}$ stacks the entries of a tensor into a long column vector, such that the indices are in reverse lexicographical order. (Using the vector outer product $\circ$, this relation takes the form $\mathcal{X}=u^{(1)}\circ u^{(2)}\circ\cdots\circ u^{(d)}$.) When $\mathcal{X}$ represents the discretization of a separable function $f(x_{1},x_{2},\ldots,x_{d})=f_{1}(x_{1})f_{2}(x_{2})\ldots f_{d}(x_{d})$ then $\mathcal{X}$ is a rank-one tensor with each vector $u^{(\mu)}$ corresponding to a discretization of $f_{\mu}$.

The CP (CANDECOMP/PARAFAC) decomposition is a sum of rank-one tensors:

The tensor rank of $\mathcal{X}$ is defined as the minimal $R$ such that $\mathcal{X}$ has a CP decomposition with $R$ terms. Note that, in contrast to matrices, the set CP($R$) of tensors of rank at most $R$ is in general not closed, which renders the problem of finding a best low-rank approximation ill-posed . For more properties of the CP decomposition, we refer to the survey paper .

The CP decomposition requires the storage of $(n_{1}+n_{2}+\cdots+n_{d})R$ entries, which becomes very attractive for small $R$. To be able to use the CP decomposition for the approximation of function-related tensors, robust and efficient compression techniques are essential. In particular, the truncation of a rank-$R$ tensor to lower tensor rank is frequently needed. Nearly all existing algorithms are based on carefully adapting existing optimization algorithms, see, once again, for an overview of the literature until around 2009. More recent developments for general tensors include work on increasing the efficiency and robustness of gradient-based and Newton-like methods , modifying and improving ALS (alternating least squares) , studying the convergence of ALS and reducing the cost of the unfolding operations required during the approximation .

One can impose additional structure on the coefficients of the CP decomposition, such as nonnegativity. As this is of primary interest in data analysis applications, a comprehensive discussion is beyond the scope of this survey, see .

2. Tucker decomposition

A Tucker decomposition of a tensor $\mathcal{X}$ takes the form

Like CP, the Tucker decomposition has a long history and we refer to the survey for a more detailed account. In the following, we briefly summarize the basic techniques, which are needed to motivate the TT and HT decompositions discussed below.

The Tucker decomposition is closely related to the matricizations of $\mathcal{X}$. The $\mu$th matricization $X^{(\mu)}$ is an n_{\mu}\times\big{(}n_{1}\cdots n_{\mu-1}n_{\mu+1}\cdots n_{d}\big{)} matrix formed in a specific way from the entries of $\mathcal{X}$:

In contrast to the tensor rank related to the CP decomposition, the set T($r_{1},\ldots,r_{d}$) of tensors of $\mu$-rank at most $r_{\mu}$ is closed.

Another consequence of the relation (4) is the higher-order SVD (HOSVD) introduced in for approximating a tensor by a Tucker decomposition (3) of lower multilinear rank. In HOSVD, the columns of each factor matrix $U_{\mu}$ are computed as the $k_{\mu}$ dominant left singular vectors of $X^{(\mu)}$. The core tensor is then obtained by forming \text{vec}(\mathcal{C}):=\big{(}U_{d}\otimes\cdots\otimes U_{1}\big{)}^{T}\text{vec}(\mathcal{X}). Eventually, this yields

In contrast to the matrix case, where the SVD yields a best low-rank approximation for all unitarily invariant norms [126, Sec. 7.4.9], the truncated tensor $\widetilde{\mathcal{X}}$ resulting from the HOSVD is usually not optimal. However, we have

This quasi-optimality condition is usually sufficient for the purpose of obtaining an accurate approximation to a function-related tensor.

Various alternatives to improve on the approximation provided by the HOSVD have been developed, see and the references therein. Recent developments include Newton-type methods on manifolds , a Jacobi algorithm for symmetric tensors , generalizations of Krylov subspace methods , and modifications of the HOSVD .

3. Tensor train decomposition

The need for storing the $r_{1}\times\cdots\times r_{d}$ core tensor $\mathcal{C}$ renders the Tucker decomposition increasingly unattractive as $d$ gets larger. This has motivated the search for decompositions which potentially avoid these exponentially growing memory requirements, while still featuring the two most important advantages of the Tucker decomposition: closedness and SVD-based compression.

One well established candidate for such a decomposition is the so called TT (tensor train) decomposition, which takes the form

where $r_{0}=r_{d}=1$. For every mode $\mu$ and every index $i_{\mu}$ the coefficients $G_{\mu}(i_{\mu})$ are matrices. In the context of numerical analysis, a decomposition of the form (5) was first proposed in . However, such a decomposition has been proposed earlier in the density-matrix renormalization group method (DMRG) for simulating quantum systems . In this area, the term matrix product state (MPS) representation for the decomposition (5) has been established . Suitable conditions that imply a unique MPS representation can be found in . The connection between TT and MPS has been explained in .

Similar to the Tucker decomposition, the TT decomposition is closely related to certain matricizations of $\mathcal{X}$. Let $X^{(1,\ldots,\mu)}$ denote the matrix obtained by reshaping the entries of $\mathcal{X}$ into an $(n_{1}\cdots n_{\mu})\times(n_{\mu+1}\cdots n_{d})$ array, such that (5) implies \text{rank}\big{(}X^{(1,\ldots,\mu)}\big{)}\leq r_{\mu} for $\mu=1,\ldots,d$. Consequently, the tuple containing the ranks of these matricizations is called the TT-rank of $\mathcal{X}$. As explained, e.g., in a quasi-best approximation in a TT decomposition for a given TT-rank can be obtained from the SVDs of $X^{(1,\ldots,\mu)}$, similarly to the HOSVD. It is important to avoid the explicit construction of these matrices and the SVDs when truncating a tensor in TT decomposition to lower TT-rank. Such truncation algorithms are described in . On the theoretical side, it turns out that the set $TT(r_{1},\ldots,r_{d-1})$ of tensors with TT-ranks bounded by $r_{\mu}$ is closed, and under a full rank condition it actually forms a smooth manifold . The Kähler manifold structure for complex MPS with open and periodic boundary conditions has been studied in .

Tensor network diagrams, which have been attributed to Penrose , are helpful in visualizing tensor decompositions and their manipulation. Figure 1 gives a few basic examples, see, e.g., for more details.

In particular, Figure 1 (v) gives an illustration of the contraction (5) representing a TT decomposition. In view of this diagram, the TT decomposition is also sometimes called linear tensor network .

In applications related to quantum spin systems, the tensor $\mathcal{X}$ often exhibits symmetries inherited from underlying physical properties. There are variants of MPS/TT that reflect such symmetries in the low-rank decomposition, see and the references therein.

4. Hierarchical Tucker decomposition

An alternative way to reduce the complexity of the Tucker decomposition is given by the hierarchical Tucker (HT) decomposition (also called hierarchical tensor representation). This decomposition is based on the idea of recursively splitting the modes of the tensor, which results in a binary tree $\mathcal{T}$ containing a subset $t\subset\{1,\ldots,d\}$ at each node. An example of such a binary tree is given in the left plot of Figure 2. The matricization $X^{(t)}$ of a tensor $\mathcal{X}$ corresponding to such a subset $t$ merges all modes contained in $t$ into row indices of the matrix, and the other modes into column indices. We then consider a hierarchy of matrices $U_{t}$ whose columns span the image of $X^{(t)}$ for each $t\in\mathcal{T}$. Hence, $U_{t}$ has exactly $r_{t}=\text{rank}(X^{(t)})$ columns. The rank tuple $(r_{t})_{t\in\mathcal{T}}$ is called the HT-rank of $\mathcal{X}$.

The following nestedness property allows for the implicit storage of $(U_{t})_{t\in\mathcal{T}}$, and thus of the tensor $\mathcal{X}$: For $t=t_{l}\cup t_{r}$, $t_{l}\cap t_{r}=\emptyset$, there exists a matrix $B_{t}$ such that

For simplicity, we have assumed that the ordering of the modes in the tree $\mathcal{T}$ is such that all modes contained in $t_{l}$ are smaller than the modes contained in $t_{r}$. The relation (6) implies that it suffices to store the basis matrices $U_{t}$ only for the leaf nodes $t=\{1\},\{2\},\ldots,\{d\}$, and $B_{t}$ for all other nodes in $\mathcal{T}$. The resulting storage requirements are $O(dnr+dr^{3})$, when assuming $r\equiv r_{t}$ and $n\equiv n_{\mu}$.

Similarly to the Tucker and TT decompositions, a quasi-best approximation in the HT decomposition for a given HT-rank can be obtained from the SVDs of $X^{(t)}$. Algorithms that avoid the explicit computation of these SVDs when truncating a tensor that is already in HT decomposition are discussed in . As for the TT decomposition, the set of tensors having fixed HT-rank forms a smooth manifold .

The tensor network corresponding to the HT decomposition is always a binary tree, see also the right plot of Figure 2.

Such tensor tree networks had already been discussed in (without the basis matrices at the leafs). Moreover, the so called multilayer multi-configuration time-dependent Hartree method (ML-MCTDH) introduced in makes use of a decomposition based on general trees instead of binary trees. When allowing for general trees, tensor tree networks include the Tucker decomposition from Section 2.2 as a (quite particular) special case. In the case of a degenerate tree, where at each level, one mode is split from the remaining modes, the HT decomposition becomes equivalent to a variant of the TT decomposition discussed in . In contrast to the TT decomposition defined in (5), this variant features additional basis matrices, which may reduce the storage cost for large $n_{\mu}$. A discussion on the difference between the ranks for the HT and TT decompositions can be found in .

5. More general tensor network formats

6. Hybrid formats

Adding to the diversity of the formats discussed above, it is possible and sometimes useful to combine different low-rank formats. One popular combination is the Tucker format combined with the CP format for the approximation of the core tensor , see [169, Sec. 5.7] for other variations of Tucker and CP. In , combinations of low-rank tensor formats with hierarchical matrices are investigated.

7. A priori approximation results

As an essential prerequisite for the success of tensor-based computations, it is important to decide whether a tensor generated by a certain multivariate function $f(x_{1},\ldots,x_{d})$ can be well approximated by a low-rank tensor decomposition. As discussed in Section 2.1, the tensor rank is closely linked to approximating $f$ by a sum of separable functions. Only in exceptional cases, it will be possible to represent $f$ exactly by such a sum, see .

In general, one is therefore interested in an approximation of the form

For a function of the form $f(x_{1},\ldots,x_{d})=g(x_{1}+\cdots+x_{d})$, such an approximation can be immediately obtained from approximating $g$ by a sum of exponentials. For this purpose, various approaches have been discussed in . Other techniques include applying numerical quadrature to an integral representation of the function, see, e.g., , Taylor series expansion , and polynomial interpolation . Based on results by Temlyakov on bilinear approximation rates, singular value estimates for the hierarchical Tucker decomposition of functions in mixed Sobolev spaces have been obtained in , see also . In fact, a sparse grid approximation to $f$ can be turned quite effectively into a low-rank tensor decomposition . General nonlinear best $R$-term approximation schemes represent another important technique, which we cannot cover in detail.

Even for smooth $f$, it may not always be possible to attain sufficiently low ranks, especially when the variation of $f$ is too strong across its entire domain of definition. In this case, it can be advantageous to subdivide the domain and approximate $f$ on each subdomain separately with a low-rank tensor decomposition, see for examples. As first discussed in , an approximation of the form (7) can also be used to approximate linear operators on tensors, see also Section 2.8 below.

In the other direction, SVD-based approximations of a function-related tensor yield an approximation of the underlying function, where the $L^{2}$-norm approximation error can be directly controlled by the truncated singular values. For smooth functions, best approximations in tensor formats are known to inherit the regularity of the approximated function . This also holds for SVD based quasi-best approximations, for which even the smoothness of the error can be controlled . This can be used, e.g., for deriving $L^{\infty}$ error estimates or approximation results for the basis matrices $U_{\mu}$ .

For quantum many-body systems, the approximability of the ground state by a low-rank TT decomposition is closely linked to the concepts of entropy and entanglement, see for an introduction.

8. Low-rank decomposition of linear operators

The matrix representation of a linear operator

can be viewed as an $m_{1}n_{1}\times m_{2}n_{2}\times\cdots\times m_{d}n_{d}$ tensor after pairing up row and column indices:

This view allows to apply any of the low-rank tensor decompositions discussed above to $\mathcal{A}$. Such a low-rank decomposition of $\mathcal{A}$ is useful, e.g., for performing the matrix-vector product $\mathcal{A}(\mathcal{X})$ efficiently when $\mathcal{X}$ itself admits a low-rank decomposition. This idea appears to be ubiquitous in the literature on low-rank tensor decomposition, see for an early reference. For example, in the study of strongly correlated quantum systems, matrix product operators (MPO) were introduced in , which corresponds to representing $\mathcal{A}$ in the TT decomposition. As pointed out in , having $\mathcal{A}$ in a low-rank decomposition also allows to compute an approximate inverse by combining the Newton-Hotelling-Schulz algorithm with truncation, see Section 3.1.

9. Tensorization

QTT has been applied to the solution of PDEs and eigenvalue problems , evaluation of boundary integrals in BEM , convolution and the FFT . A connection between QTT and the wavelet transform is discussed in .

One important ingredient of QTT is that the involved matrices can be represented in a way that conforms to the format . For the following matrices, QTT representations have been discussed: Toeplitz matrices , (inverse) Laplace operators , linear diffusion operators .

10. Software

There are several Matlab toolboxes available for dealing with tensors in CP and Tucker decomposition, including the Tensor Toolbox , the $N$-way toolbox , the PLS_Toolbox , and the Tensorlab . The TT-Toolbox provides Matlab classes covering tensors in TT and QTT decomposition, as well as linear operators. There is also a Python implementation of the TT-toolbox called ttpy . The htucker toolbox provides a Matlab class representing a tensor in HT decomposition.

The TensorCalculus library is a mathematically oriented C++ library, allowing for computations with general tensor networks. The Heidelberg MCTDH Package is a set of Fortran programs for multi-dimensional quantum dynamics. ALPS provides C++ libraries for simulating strongly correlated quantum mechanical systems, including DMRG. Block is a C++ implementation of the DMRG algorithms discussed in . The tensor contraction engine automatically generates near-optimal code for tensor contractions in many-body electronic structure methods.

Algorithms

In applications for function-related tensors, $\mathcal{X}$ is often given implicitly, e.g., as the solution to a linear system or eigenvalue problem. There are mainly two different types of approaches to obtain an approximation to $\mathcal{X}$. A first class of methods is based on combining classical iterative algorithm with repeated low-rank truncation. A second class is based on reformulating the problem at hand as an optimization problem, constraining the admissible set to low-rank tensors, and applying various optimization techniques.

In principle, any vector iteration for solving a linear algebra problem involving a tensor can be combined with truncation in any of the low-rank decompositions discussed above. To illustrate the basic principle, let us consider the (preconditioned) Richardson iteration for solving a linear system $\mathcal{A}(\mathcal{X})=\mathcal{B}$:

where $\mathcal{P}$ is a preconditioner and $\omega$ is a suitably chosen scalar. Letting $\mathcal{T}$ denote truncation in any of the low-rank tensor decompositions discussed above, we obtain the truncated Richardson iteration

which has been proposed in combination with the CP decomposition.

Other examples of combining iterative methods with low-rank truncation include:

the (shift-and-invert) power method combined with CP ;

the (shift-and-invert) power and Lanczos methods combined with CP and Tucker ;

a restarted Lanczos method combined with TT ;

conjugate gradient type-methods for symmetric eigenvalue problems combined with truncation for low-rank matrices , as well as for TT , HT and QTT .

the Richardson method combined with QTT and low-rank matrix decompositions ;

the conjugate gradient method, BiCGStab and other Krylov subspace methods combined with HT and low-rank matrix decompositions ;

When applying an iterative method in combination with a low-rank decompositions, the ranks will inevitably grow quite quickly in the course of the iteration. For example, the sum of $k$ tensors can multiply the ranks by $k$, while the pointwise (Hadamard) product between two tensors can even square the ranks. To gain efficiency, it is therefore advisable to not let this rank growth happen explicitly and combine such an operation directly with truncation. This has been discussed for sums in and for the Hadamard product (and other bilinear operations) in , see also .

Preconditioners not only help accelerate convergence but numerical evidence suggests that an effective preconditioner leads to a limited intermediate rank growth. At the same time, it is mandatory that the preconditioner can be applied efficiently to a low-rank tensor decomposition. In particular, this is the case when the preconditioner itself admits a low-rank decomposition in the sense of Section 2.8. There are various techniques to construct such preconditioners. An early technique is based on best approximation in the Frobenius norm by a Kronecker product , by a short sum of Kronecker products , or by a more general low-rank tensor decomposition. Other techniques include the use of approximate inverses for high-dimensional Laplace operators , low-rank manipulation of the PDE coefficients , low-rank tensor approximation of multilevel preconditioners , and low-rank tensor diagonal preconditioners for wavelet discretizations .

2. Optimization-based algorithms

In many cases, a linear algebra problems involving a tensor can be posed as an optimization problem. For example, it is well known that a symmetric positive definite linear system $\mathcal{A}(\mathcal{X})=\mathcal{B}$ can be turned into

where $\langle\cdot,\cdot\rangle$ corresponds to the standard inner product for the vectorization of the tensors. For nonsymmetric linear systems, an optimization problem can be obtained by minimizing the norm of the residual. For symmetric eigenvalue problems, the Rayleigh-quotient minimization or, more generally, the trace minimization principle can be used.

Once the optimization problem is set up, the set of admissible tensors $\mathcal{X}$ is then constrained to a low-rank decomposition, for example, to all tensors with fixed tensor rank or fixed multilinear rank. Even when the original optimization problem, such as (8), is convex and in principle simple, the resulting constrained optimization problem is highly nonlinear and non-convex in general. A number of heuristic approaches to the solution of such constrained optimization problems are available, including ALS (alternating linear scheme). The basic principle of ALS is to optimize every factor of the low-rank decomposition separately and to sweep over all factors repeatedly. It is probably most natural to combine ALS with the CP decomposition, see for an application of this idea to linear systems. The combination of ALS with the TT decomposition has been considered in . The convergence of ALS for the TT decomposition has been studied in .

The ALS scheme can be improved in various ways. One quite successful improvement for decompositions described by tensor networks is to join two neighboring factors, optimize the resulting supernode, and split the result into separate factors by a low-rank factorization. Originally, this so called DMRG method had been developed for the simulation of strongly correlated quantum lattice systems, see for an overview. Later on, the ideas of DMRG have been picked up and extended to other applications in the numerical analysis community in a series of papers .

There is a growing interest in applying so called Riemannian optimization techniques to (8). Examples include nonlinear conjugate gradient or Newton-like methods on manifolds of low-rank matrices or low-rank tensors , see also Section 3.4. For tensors in CP decomposition, the approximate solution of (8) by gradient techniques has been discussed in .

3. Successive rank-1 approximation

A tempting and surprisingly successful approach to the solution of high-dimensional problems is to build up a low-rank approximation from successive rank-1 approximations. This idea has been suggested in the context of various applications, including Fokker-Planck equations and stochastic partial differential equations .

such that $\mathcal{X}_{r+1}=\mathcal{X}_{r}+\mathcal{W}$ is an improved approximation, that is,

Analogous to Section 3.2, the unknown vectors $w^{(1)},\ldots,w^{(d)}$ can be determined by turning (10) into a nonlinear (optimization) problem and applying standard methods, such as the alternating direction method. This procedure is repeated until the residual $\mathcal{A}(\mathcal{X}_{r})-\mathcal{B}$ is sufficiently small. Of course, such a greedy approach will not yield the best rank-$R$ approximation after $R$ steps . However, it is important to remember that these methods aim at a more moderate goal, to obtain a reasonable approximation after $R$ steps, with $R$ not too large. Convergence results in this direction can be found in .

A number of improvements to the simple scheme outlined above have been proposed to increase its convergence speed, see, e.g., . A connection between best rank-$1$ or, more generally, best rank-$m$ approximations to a nonlinear eigenvalue problem is explained in . This connection also motivates the use of the term generalized spectral decomposition. Many further developments, improvements, and extensions of successive low-rank approximation techniques have taken place during the last years; we refer to for an overview.

4. Low-rank methods for dynamical problems

for which a typical application is the spatial discretization of a time-dependent $d$-dimensional PDE. Dynamical low-rank methods aim to determine an approximation $\mathcal{Y}(t)$ in a manifold $\mathcal{M}$ of low-rank tensors by restricting the dynamics of (11) to the tangent space $T_{\mathcal{Y}(t)}\mathcal{M}$:

As explained in , this approximation is closely related to the can Dirac-Frenkel-McLachlan variational principle in quantum molecular dynamics.

Initially proposed for low-rank matrix manifolds in , dynamical low-rank methods have been extended to low-rank tensors in Tucker , TT/MPS , and HT decomposition. The efficient and robust numerical integration of (12) is crucial to the success of dynamical low-rank methods; apart from the references above, this aspect has been discussed in .

A more immediate approach to (11) is to combine a standard time stepping method, such as the explicit and implicit Euler methods, with low-rank truncation in every time step . An alternative, which allows to control the error global-in-time, is to apply iterative solvers to a space-time formulation .

5. Black box approximation

Suppose that a matrix $A$ or a tensor $\mathcal{X}$ is defined through a function that returns entries at arbitrary positions. Then the goal of black box approximation is to find a good low-rank approximation based only on relatively few entries. It is important to emphasize that the selection of the entries can be controlled by the user. In this respect, this situation is quite different from the growing area of tensor completion, see, e.g., , where the selection of the entries is usually prescribed by the application.

For an $m\times n$ matrix $A$, the so called cross approximation method produces an approximation of the form

where Matlab notation is used to denote the submatrices of $A$ corresponding to the index sets $I\subset\{1,\ldots,m\}$ and $J\subset\{1,\ldots,n\}$ . In the $p$th step of cross approximation as described in , the entry of largest magnitude in the column $j_{p}$ of $A-A(:,J)A(I,J)^{-1}A(I,:)$ is calculated and its position is denoted by $i_{p}$. Then the entry of largest magnitude in the row $i_{p}$ of that matrix is calculated and its position is denoted by $j_{p+1}$. Moreover, both index sets are updated: $I\leftarrow I\cup\{i_{p}\}$ and $J\leftarrow J\cup\{j_{p}\}$. Volume maximization is an alternative entry selection strategy that attempts to maximize \big{|}\det\big{(}A(I,J)\big{)}\big{|}, see .

A first extension of cross approximation to tensors was proposed in for approximating third-order tensors by a Tucker decomposition. Essentially, this extension consists of applying the algorithm above to an arbitrary matricization of the tensor. However, the rows of this matricization (corresponding to slices of the tensor) are further approximated by a low rank matrix, again using cross approximation. In , this method has been combined with multilevel ideas and applied to quantum chemistry. The adaptive cross approximation for the Tucker decomposition has been analyzed in . Another variant for third and fourth order tensors, focussing on interpolation properties, is discussed in . A cross approximation method for the TT decomposition has been proposed in .

Based on fiber crosses, a quite different extension for tensors of arbitrary order can be found in . In this method a set of multi-indices $I^{1},\ldots,I^{p}\in[1,n_{1}]\times\cdots[1,n_{d}]$ is computed successively. Subspaces $\mathcal{U}_{\mu}$ are constructed approximately containing all $\mu$-mode fibers passing through at least one of the multi-indices. The tensor is then approximated by a CP decomposition (2) under the constraint $u_{j}^{(\mu)}\in\mathcal{U}_{\mu}$, using general optimization methods. The multi-indices are selected by performing an alternating direction search along the fibers of the error tensor. Also based on fiber crosses, a black box approximation for a tensor in the HT decomposition is given in .

Randomized algorithms represent an alternative way to quickly extract a low-rank approximation from partial information on the entries of the matrix, see and the references therein. These ideas have been extended to low-rank tensor decompositions, for the Tucker decomposition in and for the TT decomposition in .

6. Other algorithms

For linear systems and eigenvalue problems with a very particular structure, it is sometime possible to design specialized algorithms that can be more efficient and easier to analyze. This applies, in particular, to discretizations of the multi-dimensional Poisson equation .

Acknowledgments

We thank Jonas Ballani, Wolfgang Hackbusch, Thomas Huckle, Boris Khoromskij, Anthony Nouy, Ivan Oseledets, André Uschmajew, and Bart Vandereycken for many helpful suggestions on a preliminary draft of this survey.

References