O(dlog?N)-Quantics Approximation of N-d Tensors in High-Dimensional Numerical Modeling

In the present paper, we discuss the novel concept of super-compressed tensor-structured data formats in high-dimensional applications. We describe the multifolding or quantics-based tensor approximation method of O(dlog N)-complexity (logarithmic scaling in the volume size), applied to the discrete functions over the product index set {1,…,N}^⊗d , or briefly N-d tensors of size N ^d , and to the respective discretized differential-integral operators in ℝ ^d . As the basic approximation result, we prove that a complex exponential sampled on an equispaced grid has quantics rank 1. Moreover, a Chebyshev polynomial, sampled over a Chebyshev Gauss–Lobatto grid, has separation rank 2 in the quantics tensor format, while for the polynomial of degree m over a Chebyshev grid the respective quantics rank is at most 2m+1. For N-d tensors generated by certain analytic functions, we give a constructive proof of the O(dlog Nlog ε ⁻¹)-complexity bound for their approximation by low-rank 2-(dlog N) quantics tensors up to the accuracy ε>0. In the case ε=O(N ^−α ), α>0, our approach leads to the quantics tensor numerical method in dimension d, with the nearly optimal asymptotic complexity O(d/αlog ² ε ⁻¹). From numerical examples presented here, we observe that the quantics tensor method has proved its value in application to various function related tensors/matrices arising in computational quantum chemistry and in the traditional finite element method/boundary element method (FEM/BEM). The tool apparently works.     

Blended kernel approximation in the ℋ︁‐matrix techniques

Several types of ??‐matrices were shown to provide a data‐sparse approximation of non‐local (integral) operators in FEM and BEM applications. The general construction is applied to the operators with asymptotically smooth kernel function provided that the Galerkin ansatz space has a hierarchical structure. The new class of ??‐matrices is based on the so‐called blended FE and polynomial approximations of the kernel function and leads to matrix blocks with a tensor‐product of block‐Toeplitz (block‐circulant) and rank‐k matrices. This requires the translation (rotation) invariance of the kernel combined with the corresponding tensor‐product grids. The approach allows for the fast evaluation of volume/boundary integral operators with possibly non‐smooth kernels defined on canonical domains/manifolds in the FEM/BEM applications. (Here and in the following, we call domains canonical if they are obtained by translation or rotation of one of their parts, e.g. parallelepiped, cylinder, sphere, etc.) In particular, we provide the error and complexity analysis for blended expansions to the Helmholtz kernel. Copyright © 2002 John Wiley & Sons, Ltd.     

Data-sparse approximation to a class of operator-valued functions

In earlier papers we developed a method for the data-sparse approximation of the solution operators for elliptic, parabolic, and hyperbolic PDEs based on the Dunford-Cauchy representation to the operator-valued functions of interest combined with the hierarchical matrix approximation of the operator resolvents. In the present paper, we discuss how these techniques can be applied to approximate a hierarchy of the operator-valued functions generated by an elliptic operator .

     

Tensor method for optimal control problems constrained by fractional three-dimensional elliptic operator with variable coefficients

We introduce the tensor numerical method for solving optimal control problems that are constrained by fractional two- (2D) and three-dimensional (3D) elliptic operators with variable coefficients. We solve the governing equation for the control function which includes a sum of the fractional operator and its inverse, both discretized over large 3D $n\times n\times n$ spacial grids. Using the diagonalization of the arising matrix-valued functions in the eigenbasis of the one-dimensional Sturm–Liouville operators, we construct the rank-structured tensor approximation with controllable precision for the discretized fractional elliptic operators and the respective preconditioner. The right-hand side in the constraining equation (the optimal design function) is supposed to be represented in a form of a low-rank canonical tensor. Then the equation for the control function is solved in a tensor structured format by using preconditioned CG iteration with the adaptive rank truncation procedure that also ensures the accuracy of calculations, given an $\epsilon$-threshold. This method reduces the numerical cost for solving the control problem to $O(n\log n)$ (plus the quadratic term $O(n^2)$ with a small weight), which outperforms traditional approaches with $O(n^3\log n)$ complexity in the 3D case. The storage for the representation of all 3D nonlocal operators and functions involved is also estimated by $O(n\log n)$. This essentially outperforms the traditional methods operating with fully populated $n^3\times n^3$ matrices and vectors in ${?}^{n^3}$. Numerical tests for 2D/3D control problems indicate the almost linear complexity scaling of the rank truncated preconditioned conjugate gradient iteration in the univariate grid size n.     

Approximate iterations for structured matrices

Important matrix-valued functions f (A) are, e.g., the inverse A ⁻¹, the square root and the sign function. Their evaluation for large matrices arising from pdes is not an easy task and needs techniques exploiting appropriate structures of the matrices A and f (A) (often f (A) possesses this structure only approximately). However, intermediate matrices arising during the evaluation may lose the structure of the initial matrix. This would make the computations inefficient and even infeasible. However, the main result of this paper is that an iterative fixed-point like process for the evaluation of f (A) can be transformed, under certain general assumptions, into another process which preserves the convergence rate and benefits from the underlying structure. It is shown how this result applies to matrices in a tensor format with a bounded tensor rank and to the structure of the hierarchical matrix technique. We demonstrate our results by verifying all requirements in the case of the iterative computation of A ⁻¹ and . This work was performed during the stay of the third author at the Max-Planck-Institute for Mathematics in the Sciences (Leipzig) and also supported by the Russian Fund of Basic Research (grants 05-01-00721, 04-07-90336) and a Priority Research Grant of the Department of Mathematical Sciences of the Russian Academy of Sciences.     

Fast tensor method for summation of long‐range potentials on 3D lattices with defects

In this paper, we present a method for fast summation of long‐range potentials on 3D lattices with multiple defects and having non‐rectangular geometries, based on rank‐structured tensor representations. This is a significant generalization of our recent technique for the grid‐based summation of electrostatic potentials on the rectangular L × L × L lattices by using the canonical tensor decompositions and yielding the O(L) computational complexity instead of O(L³) by traditional approaches. The resulting lattice sum is calculated as a Tucker or canonical representation whose directional vectors are assembled by the 1D summation of the generating vectors for the shifted reference tensor, once precomputed on large N × N × N representation grid in a 3D bounding box. The tensor numerical treatment of defects is performed in an algebraic way by simple summation of tensors in the canonical or Tucker formats. To diminish the considerable increase in the tensor rank of the resulting potential sum, the ?‐rank reduction procedure is applied based on the generalized reduced higher‐order SVD scheme. For the reduced higher‐order SVD approximation to a sum of canonical/Tucker tensors, we prove the stable error bounds in the relative norm in terms of discarded singular values of the side matrices. The required storage scales linearly in the 1D grid‐size, O(N), while the numerical cost is estimated by O(NL). The approach applies to a general class of kernel functions including those for the Newton, Slater, Yukawa, Lennard‐Jones, and dipole‐dipole interactions. Numerical tests confirm the efficiency of the presented tensor summation method; we demonstrate that a sum of millions of Newton kernels on a 3D lattice with defects/impurities can be computed in seconds in Matlab implementation. The tensor approach is advantageous in further functional calculus with the lattice potential sums represented on a 3D grid, like integration or differentiation, using tensor arithmetics of 1D complexity. Copyright © 2015 John Wiley & Sons, Ltd.     

Robust Schur complement method for strongly anisotropic elliptic equations

We present robust and asymptotically optimal iterative methods for solving 2D anisotropic elliptic equations with strongly jumping coefficients, where the direction of anisotropy may change sharply between adjacent subdomains. The idea of a stable preconditioning for the Schur complement matrix is based on the use of an exotic non‐conformal coarse mesh space and on a special clustering of the edge space components according to the anisotropy behavior. Our method extends the well known BPS interface preconditioner [2] to the case of anisotropic equations. The technique proposed also provides robust solvers for isotropic equations in the presence of degenerate geometries, in particular, in domains composed of thin substructures. Numerical experiments confirm efficiency and robustness of the algorithms for the complicated problems with strongly varying diffusion and anisotropy coefficients as well as for the isotropic diffusion equations in the ‘brick and mortar’ structures involving subdomains with high aspect ratios. Copyright © 1999 John Wiley & Sons, Ltd.     

Fast and accurate tensor approximation of a multivariate convolution with linear scaling in dimension

In the present paper we present the tensor-product approximation of a multidimensional convolution transform discretized via a collocation-projection scheme on uniform or composite refined grids. Examples of convolving kernels are provided by the classical Newton, Slater (exponential) and Yukawa potentials, 1/‖x‖, and with x∈R^d. For piecewise constant elements on the uniform grid of size n^d, we prove quadratic convergence O(h²) in the mesh parameter h=1/n, and then justify the Richardson extrapolation method on a sequence of grids that improves the order of approximation up to O(h³). A fast algorithm of complexity O(dR₁R₂nlogn) is described for tensor-product convolution on uniform/composite grids of size n^d, where R₁,R₂ are tensor ranks of convolving functions. We also present the tensor-product convolution scheme in the two-level Tucker canonical format and discuss the consequent rank reduction strategy. Finally, we give numerical illustrations confirming: (a) the approximation theory for convolution schemes of order O(h²) and O(h³); (b) linear-logarithmic scaling of 1D discrete convolution on composite grids; (c) linear-logarithmic scaling in n of our tensor-product convolution method on an n×n×n grid in the range n≤16384.     

Data-sparse approximation to the operator-valued functions of elliptic operator

In previous papers the arithmetic of hierarchical matrices has been described, which allows us to compute the inverse, for instance, of finite element stiffness matrices discretising an elliptic operator The required computing time is up to logarithmic factors linear in the dimension of the matrix. In particular, this technique can be used for the computation of the discrete analogue of a resolvent

In the present paper, we consider various operator functions, the operator exponential negative fractional powers , the cosine operator function and, finally, the solution operator of the Lyapunov equation. Using the Dunford-Cauchy representation, we get integrals which can be discretised by a quadrature formula which involves the resolvents mentioned above. We give error estimates which are partly exponentially, partly polynomially decreasing.