Tensor Numerical Methods: Actual Theory and Recent Applications

Numerical treatment of the high-dimensional problems by using the traditional numerical methods suffers from the so-called “curse of dimensionality”. This exponential growth O⁢(nd) of the complexity of numerical algorithms in the dimension parameter d can be only slightly relaxed by parallel computations and high performance computing. The tensor numerical methods, bridging the multilinear algebra and nonlinear approximation theory, allow to reduce solution of multidimensional problems to one-dimensional calculations. They rely on an efficient separable representation of multivariate functions and operators on large n⊗d grids, leading to algorithms of low computational cost that scale polynomially or linearly in the dimension parameter d. Thus, tensor methods may be understood as a discrete analogue of the separation of variables, which may be efficiently maintained at all steps of calculations. The target function can be the solution of some operator equation A⁢u=f, in particular PDE, and it can be represented through the so-called solution operator u=S⁢(A)⁢f. In tensor numerical methods the entities u, f, A and also S⁢(A) are gainfully approximated in the low rank tensor formats.

Traditional methods of separable approximation combine the canonical, Tucker, as well as the matrix product states (MPS) formats, the latter known as the tensor train (TT) decomposition [22, 21]. The recent tensor methods in combination with exponentially accurate sinc-based approximations are proven to provide the 𝒪⁢(d) data-compression on a wide class of functions and operators [5, 6, 9, 7, 4]. The quantized-TT (QTT) tensor approximation of functions [13] makes it possible to solve high-dimensional PDEs in quantized tensor spaces, with the log-volume complexity scaling in the full-grid size, i.e. 𝒪⁢(d⁢log⁡n), instead of 𝒪⁢(nd).

At present, tensor numerical methods and multilinear algebra continue to expand rapidly to a wide range of theoretical and applied fields, see for example [11, 14, 10]. We also refer to the recent research monographs [15, 12], where the tensor numerical methods in scientific computing with the particular focus on multi-dimensional PDEs and electronic structure calculations have been presented. These trends are reflected also in the papers from the present issue of CMAM.

This special issue is a collection of papers which demonstrate that the tensor techniques allow to solve various hard theoretical and computational problems including approximation of multi-dimensional elliptic/parabolic PDEs. This issue includes ten invited contributions on theoretical analysis and applications of tensor-based numerical methods. These papers cover a broad range of topics including construction of computational schemes for steady-state and dynamical problems as well as for stochastic and parametric equations, separation rank estimates for classes of functions and operators, numerical simulations etc. Below we briefly describe the content of the special issue.

The goal of the paper [1] is the efficient numerical solution of stochastic eigenvalue problems. Such problems often lead to prohibitively high-dimensional systems with tensor product structure when discretized with the stochastic Galerkin method. The authors exploit this inherent tensor product structure to develop a globalized low-rank inexact Newton method with which they tackle the stochastic eigenvalue problem. The effectiveness of the solver is illustrated by numerical experiments.

The paper [2] deals with an algorithm for solution of high-dimensional evolutionary equations (ODEs and discretized time-dependent PDEs) in the TT decomposition, assuming that the solution and the right-hand side of the ODE admit such a decomposition with a low rank parameter. A linear ODE, discretized via one-step or Chebyshev differentiation schemes, turns into a large linear system. The tensor decomposition allows to solve this system for several time points simultaneously. In numerical experiments with the transport and the chemical master equations, the author demonstrates that the new method is faster than traditional time stepping and stochastic simulation algorithms.

The paper [3] examines a completely non-intrusive, sample-based method for the computation of functional low-rank solutions of high-dimensional parametric random PDEs which have become an area of intensive research in Uncertainty Quantification. In order to obtain a generalized polynomial chaos representation of the approximate stochastic solution, a novel black-box rank-adapted tensor reconstruction procedure is proposed. The performance of the described approach is illustrated with several numerical examples and compared to Monte Carlo sampling.

The authors of [8] consider the abstract differential equations of the heat and Schrödinger type and discuss various N-parametric approximations on the base of the Cayley transform and of the Laguerre expansion providing a sub-exponential accuracy, i.e. the accuracy of the order 𝒪⁢(e-N⁢log⁡N). They propose a new approximation using the combination of the Gauss–Lobatto–Chebyshev interpolation and the Cayley transform and obtain a purely exponential accuracy of the order 𝒪⁢(e-N). The rank-structured tensor form of this approximation for a d-dimensional spatial operator coefficient results in an algorithm having a linear complexity in d.

The paper [16] study a dynamical low-rank approximation on the manifold of fixed-rank tensor trains, and analyze projection methods for the time integration of such problems. The authors prove error estimates for the explicit Euler method, amended with quasi-optimal projections to the manifold, under suitable approximability assumptions. Then they discuss the possibilities and difficulties with higher order explicit and implicit projected Runge–Kutta methods, in particular, the ways for limiting rank growth in the increments, and robustness with respect to small singular values.

The paper [17] deals with a new algorithm for spectral learning of Hidden Markov Models (HMM). In contrast to standard approach, the parameters of the HMM are not approximated directly, but through an estimate for the joint probability distribution. Using TT-format, the authors get an approximation by minimizing the Frobenius distance between the empirical joint probability distribution and tensors with low TT-ranks with core tensors normalization constraints. An algorithm for the solution of optimization problem that is based on the alternating least squares (ALS) approach is proposed and its fast version for sparse tensors is developed. The authors compare the performance of the proposed algorithm with the existing schemes and found that it is much more robust if the number of hidden states is overestimated.

The paper [18] describes advanced numerical tools for working with multivariate functions and for the analysis of large data sets. In particular, covariance matrices are crucial in spatio-temporal statistical tasks, but are often very expensive to compute and store, especially in 3D. Therefore, one can alternatively use a low-rank tensor formats, which reduce the computing and storage costs essentially. The authors apply the Tucker and canonical tensor decompositions to a family of Matérn-type radial functions with varying parameters and demonstrate theoretically and numerically that their tensor approximations exhibit exponentially fast convergence in the rank parameter, thus providing low computational complexity.

The paper [19] deals with a space-time isogeometric analysis scheme for the discretization of parabolic evolution equations with diffusion coefficients depending on both time and space variables. The problem is considered in a space-time cylinder in ℝd+1, with d=2,3 and is discretized using higher-order and highly-smooth spline spaces. This makes the matrix formation task very challenging from a computational point of view. The authors overcome this problem by introducing a low-rank decoupling of the operator into space and time components. Numerical experiments demonstrate the efficiency of this approach.

In [20] the authors propose an efficient algorithm to compute a low-rank approximation to the solution of so-called “Laplace-like” linear systems. The idea is to transform the problem into the frequency domain, and then to use cross approximation. In this case, we do not need to form explicit approximation to the inverse operator and can approximate the solution directly, which leads to reduced complexity. It is demonstrated that the proposed method is fast and robust by using it as a solver inside Uzawa iterative method for solving the Stokes problem.

The problem of approximately solving a system of univariate polynomials with one or more common roots and its coefficients corrupted by noise is studied in [23]. New Rayleigh quotient methods are proposed and evaluated for estimating the common roots. Using tensor algebra, reasonable starting values for the Rayleigh quotient methods can be computed. The new methods are compared to Gauss–Newton, solving an eigenvalue problem obtained from the generalized Sylvester matrix, and building a cluster among the roots of all polynomials. It is shown in a simulation study that Gauss–Newton and a new Rayleigh quotient method perform best, where the latter is more accurate when other roots than the true common roots are close together.