Stability of Low-Rank Tensor Representations and Structured Multilevel Preconditioning for Elliptic PDEs

Foundations of Computational Mathematics https://doi.org/10.1007/s10208-020-09446-z Stability of Low-Rank Tensor Representations and Structured Multilevel Preconditioning for Elliptic PDEs Markus Bachmayr1 · Vladimir Kazeev2 Received: 5 March 2018 / Revised: 18 June 2019 / Accepted: 5 November 2019 © The Author(s) 2020 Abstract Folding grid value vectors of size 2 L into Lth-order tensors of mode size 2 × · · · × 2, combined with low-rank representation in the tensor train format, has been shown to result in highly efficient approximations for various classes of functions. These include solutions of elliptic PDEs on nonsmooth domains or with oscillatory data. This tensor-structured approach is attractive because it leads to highly compressed, adaptive approximations based on simple discretizations. Standard choices of the underlying bases, such as piecewise multilinear finite elements on uniform tensor product grids, entail the well-known matrix ill-conditioning of discrete operators. We demonstrate that, for low-rank representations, the use of tensor structure itself additionally introduces representation ill-conditioning, a new effect specific to computations in tensor networks. We analyze the tensor structure of a BPX preconditioner for a second-order linear elliptic operator and construct an explicit tensor-structured representation of the preconditioner, with ranks independent of the number L of discretization levels. The straightforward application of the preconditioner yields discrete operators whose matrix conditioning is uniform with respect to the discretization parameter, but in decompositions that suffer from representation ill-conditioning. By additionally eliminating certain redundancies in the representations of the preconditioned discrete operators, we obtain reduced-rank decompositions that are free of both matrix and representation ill-conditioning. For an iterative solver based on soft thresholding of low-rank tensors, we obtain convergence and complexity estimates and demonstrate its reliability and efficiency for discretizations with up to 250 nodes in each dimension. Keywords Elliptic boundary value problems · Multilevel preconditioning · Tensor decompositions · Representation condition number · Solver complexity Communicated by Endre Süli. M.B. acknowledges support by the Hausdorff Center of Mathematics, University of Bonn. M.B. was partially supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Projektnummer 211504053 - SFB 1060. Extended author information available on the last page of the article 123 Foundations of Computational Mathematics Mathematics Subject Classification 15A69 · 35J25 · 65N12 · 65N30 · 65N55 · 65F08 · 65F35 · 65Y20 1 Introduction The direct textbook treatment of elliptic PDEs by low-order discretizations on uniform grids becomes unaffordable for many important problem classes. The high computational costs are due to the prohibitively large number of degrees of freedom required to resolve specific features of solutions, such as singularities and high-frequency oscillations, that arise in problems with nonsmooth or oscillatory data. More efficient discretizations can be obtained with basis functions that are adapted to the given problem and require fewer degrees of freedom. However, the construction and analysis of such methods (for instance, of hp-adaptive solvers) generally depends on specific features of the considered problem classes and accordingly specialized analytical tools. By the approach considered in this work, efficiency is achieved in a different way: extremely large arrays of coefficients parametrizing simple, uniformly refined loworder discretizations are themselves parametrized as nonlinear functions of relatively few effective degrees of freedom. The latter parametrization is based on representing the coefficient arrays, reshaped into high-order tensors, in the tensor train decomposition with low ranks. This representation exploits low-rank structure with respect to a hierarchy of dyadic scales, providing, at each scale, a problem-adapted basis that can be computed using standard techniques of numerical linear algebra. In other words, for the identification of suitable degrees of freedom, this approach avoids relying on problem-specific a priori information; instead, suitable degrees of freedom are found by the low-rank tensor compression of generic, conceptually straightforward discretizations. In numerical solvers for PDE problems that operate on such highly compressed, nonlinear representations of basis coefficients, new difficulties arise compared to a standard entrywise representation. As we demonstrate in this contribution, specific types of ill-conditioning in such tensor representations can dramatically affect the numerical stability of solvers. We show how a special low-rank representation of a BPX preconditioner allows to overcome these difficulties and obtain estimates for the total computational complexity of computing solutions with low-rank tensor train structure. 1.1 Low-Rank Tensor Approximations The development of low-rank tensor representations [18,25,45,47,50], such as the tensor train format, has originally been motivated by applications to high-dimensional PDEs. As observed in [19,37,43,44], the artificial treatment of coefficient vectors in lower-dimensional problems as high-dimensional quantities, known in the literature as quantized tensor train (QTT) decomposition or tensorization, leads to highly efficient approximations in many problems of interest. See [38] for a general overview and, for instance, [29,36] for further applications. 123 Foundations of Computational Mathematics To briefly illustratethis concept, let us suppose that a function u has an accurate N approximation u ≈ j=1 u j φ j in terms of the basis functions {φ j } j=1,...,N with the coefficient vector u = (u j ) j=1,...,N ∈ R N . The basic idea is to reinterpret u L as a higher-order tensor of mode sizes n 1 × · · · × n L with =1 n  = N via the identification j ↔ (i 1 , . . . , i L ) ∈ {0, . . . , n 1 − 1} × · · · × {0, . . . , n L − 1} provided by the unique decomposition j −1= L  =1 L  i n k with i  ∈ {0, . . . , n  − 1} for all  = 1, . . . , L . k=+1 We assume a simple choice of basis functions, such as low-order splines, combined with a compressed, nonlinearly parametrized approximation of the corresponding coefficients u in the tensor train format, ui1 ,...,i L ≈ r1  α1 =1  r L−1 ··· U1 (1, i 1 , α1 ) U2 (α1 , i 2 , α2 ) · · · U L (α L−1 , i L , 1). (1) α L−1 =1 The actual degrees of freedom are now the entries of the third-order tensors U ∈ Rr−1 ×n  ×r with  ∈ {1, . . . , L}, which are referred to as cores (where r0 = r L = 1 for notational convenience). In the case of n  = n ∈ N for all , which we consider defining this approximation equals  L in this work, the total number2 of parameters 2 =1 n  r−1 r  (log N ) max{r1 , . (...truncated)