Boris N. Khoromskij

Sergey Dolgov, Boris N. Khoromskij, Alexander Litvinenko et al.

We apply the Tensor Train (TT) decomposition to construct the tensor product Polynomial Chaos Expansion (PCE) of a random field, to solve the stochastic elliptic diffusion PDE with the stochastic Galerkin discretization, and to compute some quantities of interest (mean, variance, exceedance probabilities). We assume that the random diffusion coefficient is given as a smooth transformation of a Gaussian random field. In this case, the PCE is delivered by a complicated formula, which lacks an analytic TT representation. To construct its TT approximation numerically, we develop the new block TT cross algorithm, a method that computes the whole TT decomposition from a few evaluations of the PCE formula. The new method is conceptually similar to the adaptive cross approximation in the TT format, but is more efficient when several tensors must be stored in the same TT representation, which is the case for the PCE. Besides, we demonstrate how to assemble the stochastic Galerkin matrix and to compute the solution of the elliptic equation and its post-processing, staying in the TT format. We compare our technique with the traditional sparse polynomial chaos and the Monte Carlo approaches. In the tensor product polynomial chaos, the polynomial degree is bounded for each random variable independently. This provides higher accuracy than the sparse polynomial set or the Monte Carlo method, but the cardinality of the tensor product set grows exponentially with the number of random variables. However, when the PCE coefficients are implicitly approximated in the TT format, the computations with the full tensor product polynomial set become possible. In the numerical experiments, we confirm that the new methodology is competitive in a wide range of parameters, especially where high accuracy and high polynomial degrees are required.

Peter Benner, Sergey Dolgov, Venera Khoromskaia et al.

In this paper, we study and implement the structural iterative eigensolvers for the large-scale eigenvalue problem in the Bethe-Salpeter equation (BSE) based on the reduced basis approach via low-rank factorizations in generating matrices, introduced in the previous paper. The approach reduces numerical costs down to $\mathcal{O}(N_b^2)$ in the size of atomic orbitals basis set, $N_b$, instead of practically intractable $\mathcal{O}(N_b^6)$ complexity scaling for the direct diagonalization of the BSE matrix. As an alternative to rank approximation of the static screen interaction part of the BSE matrix, we propose to restrict it to a small active sub-block, with a size balancing the storage for rank-structured representations of other matrix blocks. We demonstrate that the enhanced reduced-block approximation exhibits higher precision within the controlled numerical cost, providing as well a distinct two-sided error estimate for the BSE eigenvalues. It is shown that further reduction of the asymptotic computational cost is possible due to ALS-type iteration in block tensor train (TT) format applied to the quantized-TT (QTT) tensor representation of both long eigenvectors and rank-structured matrix blocks. The QTT-rank of these entities possesses almost the same magnitude as the number of occupied orbitals in the molecular systems, $N_o$, hence the overall asymptotic complexity for solving the BSE problem can be estimated by $\mathcal{O}(\log(N_o) N_o^{2})$. We confirm numerically a considerable decrease in computational time for the presented iterative approach applied to various compact and chain-type molecules, while supporting sufficient accuracy.

Boris N. Khoromskij

We present a brief survey on the modern tensor numerical methods for multidimensional stationary and time-dependent partial differential equations (PDEs). The guiding principle of the tensor approach is the rank-structured separable approximation of multivariate functions and operators represented on a grid. Recently, the traditional Tucker, canonical, and matrix product states (tensor train) tensor models have been applied to the grid-based electronic structure calculations, to parametric PDEs, and to dynamical equations arising in scientific computing. The essential progress is based on the quantics tensor approximation method proved to be capable to represent (approximate) function related $d$-dimensional data arrays of size $N^d$ with log-volume complexity, $O(d \log N)$. Combined with the traditional numerical schemes, these novel tools establish a new promising approach for solving multidimensional integral and differential equations using low-parametric rank-structured tensor formats. As the main example, we describe the grid-based tensor numerical approach for solving the 3D nonlinear Hartree-Fock eigenvalue problem, that was the starting point for the developments of tensor-structured numerical methods for large-scale computations in solving real-life multidimensional problems. We also address new results on tensor approximation of the dynamical Fokker-Planck and master equations in many dimensions up to $d=20$. Numerical tests demonstrate the benefits of the rank-structured tensor approximation on the aforementioned examples of multidimensional PDEs. In particular, the use of grid-based tensor representations in the reduced basis of atomics orbitals yields an accurate solution of the Hartree-Fock equation on large $N\times N \times N$ grids with a grid size of up to $N= 10^{5}$.

Venera Khoromskaia, Boris N. Khoromskij, Felix Otto

We describe the numerical scheme for the discretization and solution of 2D elliptic equations with strongly varying piecewise constant coefficients arising in the stochastic homogenization of multiscale composite materials. An efficient stiffness matrix generation scheme based on assembling the local Kronecker product matrices is introduced. The resulting large linear systems of equations are solved by the preconditioned CG iteration with a convergence rate that is independent of the grid size and the variation in jumping coefficients (contrast). Using this solver we numerically investigate the convergence of the Representative Volume Element (RVE) method in stochastic homogenization that extracts the effective behavior of the random coefficient field. Our numerical experiments confirm the asymptotic convergence rate of systematic error and standard deviation in the size of RVE rigorously established in [6]. The asymptotic behavior of covariances of the homogenized matrix in the form of a quartic tensor is also studied numerically. Our approach allows laptop computation of sufficiently large number of stochastic realizations even for large sizes of the RVE.

Alexander Litvinenko, David Keyes, Venera Khoromskaia et al.

In this work, we describe advanced numerical tools for working with multivariate functions and for the analysis of large data sets. These tools will drastically reduce the required computing time and the storage cost, and, therefore, will allow us to consider much larger data sets or finer meshes. Covariance matrices are crucial in spatio-temporal statistical tasks, but are often very expensive to compute and store, especially in 3D. Therefore, we approximate covariance functions by cheap surrogates in a low-rank tensor format. We apply the Tucker and canonical tensor decompositions to a family of Matern- and Slater-type functions with varying parameters and demonstrate numerically that their approximations exhibit exponentially fast convergence. We prove the exponential convergence of the Tucker and canonical approximations in tensor rank parameters. Several statistical operations are performed in this low-rank tensor format, including evaluating the conditional covariance matrix, spatially averaged estimation variance, computing a quadratic form, determinant, trace, loglikelihood, inverse, and Cholesky decomposition of a large covariance matrix. Low-rank tensor approximations reduce the computing and storage costs essentially. For example, the storage cost is reduced from an exponential $\mathcal{O}(n^d)$ to a linear scaling $\mathcal{O}(drn)$, where $d$ is the spatial dimension, $n$ is the number of mesh points in one direction, and $r$ is the tensor rank. Prerequisites for applicability of the proposed techniques are the assumptions that the data, locations, and measurements lie on a tensor (axes-parallel) grid and that the covariance function depends on a distance, $\Vert x-y \Vert$.

Peter Benner, Venera Khoromskaia, Boris N. Khoromskij

We introduce and analyze the new range-separated (RS) canonical/Tucker tensor format which aims for numerical modeling of the 3D long-range interaction potentials in multi-particle systems. The main idea of the RS tensor format is the independent grid-based low-rank representation of the localized and global parts in the target tensor which allows the efficient numerical approximation of $N$-particle interaction potentials. The single-particle reference potential like $1/\|x\|$ is split into a sum of localized and long-range low-rank canonical tensors represented on a fine 3D $n\times n\times n$ Cartesian grid. The smoothed long-range contribution to the total potential sum is represented on the 3D grid in $O(n)$ storage via the low-rank canonical/Tucker tensor. We prove that the Tucker rank parameters depend only logarithmically on the number of particles $N$ and the grid-size $n$. Agglomeration of the short range part in the sum is reduced to an independent treatment of $N$ localized terms with almost disjoint effective supports, calculated in $O(N)$ operations. Thus, the cumulated sum of short range clusters is parametrized by a single low-rank canonical reference tensor with a local support, accomplished by a list of particle coordinates and their charges. The RS canonical/Tucker tensor representations reduce the cost of multi-linear algebraic operations on the 3D potential sums arising in modeling of multi-dimensional data by radial basis functions, say, in computation of the electrostatic potential of a protein, in 3D integration and convolution transforms, computation of gradients, forces and the interaction energy of a many-particle systems, and in low parametric fitting of multi-dimensional scattered data by reducing all of them to 1D calculations.

Peter Benner, Venera Khoromskaia, Boris N. Khoromskij et al.

In this paper, we introduce a new interpolation scheme to approximate the density of states (DOS) for a class of rank-structured matrices with application to the Tamm-Dancoff approximation (TDA) of the Bethe-Salpeter equation (BSE). The presented approach for approximating the DOS is based on two main techniques. First, we propose an economical method for calculating the traces of parametric matrix resolvents at interpolation points by taking advantage of the block-diagonal plus low-rank matrix structure described in [6, 3] for the BSE/TDA problem. Second, we show that a regularized or smoothed DOS discretized on a fine grid of size $N$ can be accurately represented by a low rank quantized tensor train (QTT) tensor that can be determined through a least squares fitting procedure. The latter provides good approximation properties for strictly oscillating DOS functions with multiple gaps, and requires asymptotically much fewer ($O(\log N)$) functional calls compared with the full grid size $N$. This approach allows us to overcome the computational difficulties of the traditional schemes by avoiding both the need of stochastic sampling and interpolation by problem independent functions like polynomials etc. Numerical tests indicate that the QTT approach yields accurate recovery of DOS associated with problems that contain relatively large spectral gaps. The QTT tensor rank only weakly depends on the size of a molecular system which paves the way for treating large-scale spectral problems.

Peter Benner, Boris N. Khoromskij, Venera Khoromskaia et al.

We introduce and analyze a mesh-free two-level hybrid Chebyshev-Tucker tensor representation for approximating multivariate functions, which combines tensor-product Chebyshev interpolation with the low-rank Tucker decomposition of the tensor of Chebyshev coefficients. This construction allows to avoid the expensive rank-structured grid-based approximation of function-related tensors on large spatial grids, while benefiting from the Tucker decomposition of the moderate-sized core tensor of Chebyshev coefficients. Thus, we can compute the nearly optimal Tucker decomposition of the 3D function with controllable accuracy $\varepsilon >0$ without discretizing the function on a full fine grid in the domain, but only using its values at a small set of Chebyshev nodes computed either from the explicit analytic expression of the target function or from its data-sparse representation in a rank-structured tensor format with moderate rank parameter. Finally, we can represent the function in the algebraic Tucker format with optimal $\varepsilon$-rank on an arbitrarily large 3D tensor grid in the computational domain by discretizing the Chebyshev polynomials on that grid. The rank parameters of the nonlinear Tucker-ALS decomposition of the coefficient tensor can be much smaller than the polynomial degrees of the initial Chebyshev linear interpolation in the function independent polynomial basis set. It is shown that our techniques can be gainfully applied to the long-range part of the singular electrostatic potential of multi-particle systems represented on a fine grid in the range-separated (RS) tensor format. We provide error and complexity estimates and demonstrate the computational efficiency of the proposed techniques on challenging examples, including the collective electrostatic potential for large bio-molecular systems and lattice-type compounds.

Boris N. Khoromskij

In this paper, we introduce the operator dependent range-separated tensor approximation of the discretized Dirac delta in $\mathbb{R}^d$. It is constructed by application of the discrete elliptic operator to the range-separated decomposition of the associated Green kernel discretized on the Cartesian grid in $\mathbb{R}^d$. The presented operator dependent local-global splitting of the Dirac delta can be applied for solving the potential equations in non-homogeneous media when the density in the right-hand side is given by the large sum of pointwise singular charges. We show how the idea of the operator dependent RS splitting of the Dirac delta can be extended to the closely related problem on the range separated tensor representation of the elliptic resolvent. The numerical tests confirm the expected localization properties of the obtained operator dependent approximation of the Dirac delta represented on a tensor grid. As an example of application, we consider the regularization scheme for solving the Poisson-Boltzmann equation for modeling the electrostatics in bio-molecules.

Venera Khoromskaia, Boris N. Khoromskij

This paper introduces and analyses the new grid-based tensor approach for approximate solution of the eigenvalue problem for linearized Hartree-Fock equation applied to the 3D lattice-structured and periodic systems. The set of localized basis functions over spatial $(L_1,L_2,L_3)$ lattice in a bounding box (or supercell) is assembled by multiple replicas of those from the unit cell. All basis functions and operators are discretized on a global 3D tensor grid in the bounding box which enables rather general basis sets. In the periodic case, the Galerkin Fock matrix is shown to have the three-level block circulant structure, that allows the FFT-based diagonalization. The proposed tensor techniques manifest the twofold benefits: (a) the entries of the Fock matrix are computed by 1D operations using low-rank tensors represented on a 3D grid, (b) the low-rank tensor structure in the diagonal blocks of the Fock matrix in the Fourier space reduces the conventional 3D FFT to the product of 1D FFTs. We describe fast numerical algorithms for the block circulant representation of the core Hamiltonian in the periodic setting based on low-rank tensor representation of arising multidimensional functions. Lattice type systems in a box with open boundary conditions are treated by our previous tensor solver for single molecules, which makes possible calculations on large $(L_1,L_2,L_3)$ lattices due to reduced numerical cost for 3D problems. The numerical simulations for box/periodic $(L,1,1)$ lattice systems in a 3D rectangular "tube" with $L$ up to several hundred confirm the theoretical complexity bounds for the tensor-structured eigenvalue solvers in the limit of large $L$.

Boris N. Khoromskij, Kishore K. Naraparaju, Jan Schneider

In this article we consider the iterative schemes to compute the canonical (CP) approximation of quantized data generated by a function discretized on a large uniform grid in an interval on the real line. This paper continues the research on the QTT method [16] developed for the tensor train (TT) approximation of the quantized images of function related data. In the QTT approach the target vector of length $2^{L}$ is reshaped to a $L^{th}$ order tensor with two entries in each mode (Quantized representation) and then approximated by the QTT tenor including $2r^2 L$ parameters, where $r$ is the maximal TT rank. In what follows, we consider the Alternating Least-Squares (ALS) iterative scheme to compute the rank-$r$ CP approximation of the quantized vectors, which requires only $2 r L\ll 2^L$ parameters for storage. In the earlier papers [17] such a representation was called Q$_{Can}$ format, while in this paper we abbreviate it as the QCP representation. We test the ALS algorithm to calculate the QCP approximation on various functions, and in all cases we observed the exponential error decay in the QCP rank. The main idea for recovering a discretized function in the rank-$r$ QCP format using the reduced number the functional samples, calculated only at $O(2rL)$ grid points, is presented. The special version of ALS scheme for solving the arising minimization problem is described. This approach can be viewed as the sparse QCP-interpolation method that allows to recover all $2r L$ representation parameters of the rank-$r$ QCP tensor. Numerical examples show the efficiency of the QCP-ALS type iteration and indicate the exponential convergence rate in $r$.

Boris N. Khoromskij, Sergey I. Repin

The paper suggests a preconditioning type method for fast solving of elliptic equations with oscillating quasiperiodic coefficients $A_ε$ specified by the small parameter $ε>0$. We use an iteration method generated by an elliptic operator, associated with a certain simplified (e.g., homogenized) problem. On each step of this procedure it is required to solve an auxiliary elliptic boundary value problem with non--oscillating coefficients $A_0$. All the information related to complicated coefficients of the original differential problem is encompasses in the linear functional, which forms the right hand side of the auxiliary problem. Therefore, explicit inversion of the original operator associated with oscillating coefficients is avoided. The only operation used instead is multiplication of the operator by a vector (vector function), which can be efficiently performed due to the low-rank QTT tensor operations with the rank parameter controlled by the given precision $δ>0$ independent on the parameter $ε$. We deduce two--sided a posteriori error estimates that do not use $A^{-1}_ε$ and provide guaranteed two sided bounds of the distance to the exact solution of the original problem for any step of the iteration process. The second part is concerned with realisations of the iteration method. For a wide class of oscillating coefficients, we obtain sharp QTT rank estimates for the stiffness matrix in tensor representation. In practice, this leads to the logarithmic complexity scaling of the approximation and solution process in both the FEM grid-size, and $O(\vert\logε\vert)$ cost in terms of $ε$. Numerical tests in 1D confirm the logarithmic complexity scaling of our method applied to a class of complicated quasiperiodic coefficients.

Peter Benner, Venera Khoromskaia, Boris N. Khoromskij

The Bethe-Salpeter equation (BSE) is a reliable model for estimating the absorption spectra in molecules and solids on the basis of accurate calculation of the excited states from first principles. This challenging task includes calculation of the BSE operator in terms of two-electron integrals tensor represented in molecular orbital basis, and introduces a complicated algebraic task of solving the arising large matrix eigenvalue problem. The direct diagonalization of the BSE matrix is practically intractable due to $O(N^6)$ complexity scaling in the size of the atomic orbitals basis set, $N$. In this paper, we present a new approach to the computation of Bethe-Salpeter excitation energies which can lead to relaxation of the numerical costs up to $O(N^3)$. The idea is twofold: first, the diagonal plus low-rank tensor approximations to the fully populated blocks in the BSE matrix is constructed, enabling easier partial eigenvalue solver for a large auxiliary system relying only on matrix-vector multiplications with rank-structured matrices. And second, a small subset of eigenfunctions from the auxiliary eigenvalue problem is selected to build the Galerkin projection of the exact BSE system onto the reduced basis set. We present numerical tests on BSE calculations for a number of molecules confirming the $\varepsilon$-rank bounds for the blocks of BSE matrix. The numerics indicates that the reduced BSE eigenvalue problem with small matrices enables calculation of the lowest part of the excitation spectrum with sufficient accuracy.

Venera Khoromskaia, Boris N. Khoromskij

We resume the recent successes of the grid-based tensor numerical methods and discuss their prospects in real-space electronic structure calculations. These methods, based on the low-rank representation of the multidimensional functions and integral operators, led to entirely grid-based tensor-structured 3D Hartree-Fock eigenvalue solver. It benefits from tensor calculation of the core Hamiltonian and two-electron integrals (TEI) in $O(n\log n)$ complexity using the rank-structured approximation of basis functions, electron densities and convolution integral operators all represented on 3D $n\times n\times n $ Cartesian grids. The algorithm for calculating TEI tensor in a form of the Cholesky decomposition is based on multiple factorizations using algebraic 1D ``density fitting`` scheme. The basis functions are not restricted to separable Gaussians, since the analytical integration is substituted by high-precision tensor-structured numerical quadratures. The tensor approaches to post-Hartree-Fock calculations for the MP2 energy correction and for the Bethe-Salpeter excited states, based on using low-rank factorizations and the reduced basis method, were recently introduced. Another direction is related to the recent attempts to develop a tensor-based Hartree-Fock numerical scheme for finite lattice-structured systems, where one of the numerical challenges is the summation of electrostatic potentials of a large number of nuclei. The 3D grid-based tensor method for calculation of a potential sum on a $L\times L\times L$ lattice manifests the linear in $L$ computational work, $O(L)$, instead of the usual $O(L^3 \log L)$ scaling by the Ewald-type approaches.

Venera Khoromskaia, Boris N. Khoromskij

We introduce and study a novel tensor approach for fast and accurate assembled summation of a large number of lattice-allocated potentials represented on 3D $N\times N \times N$ grid with the computational requirements only \emph{weakly dependent} on the number of summed potentials. It is based on the assembled low-rank canonical tensor representations of the collected potentials using pointwise sums of shifted canonical vectors representing the single generating function, say the Newton kernel. For a sum of electrostatic potentials over $L\times L \times L$ lattice embedded in a box the required storage scales linearly in the 1D grid-size, $O(N )$, while the numerical cost is estimated by $O(N L)$. For periodic boundary conditions, the storage demand remains proportional to the 1D grid-size of a unit cell, $n=N/L$, while the numerical cost reduces to $O(N)$, that outperforms the FFT-based Ewald-type summation algorithms of complexity $O(N^3 \log N)$. The complexity in the grid parameter $N$ can be reduced even to the logarithmic scale $O(\log N)$ by using data-sparse representation of canonical $N$-vectors via the quantics tensor approximation. For justification, we prove an upper bound on the quantics ranks for the canonical vectors in the overall lattice sum. The presented approach is beneficial in applications which require further functional calculus with the lattice potential, say, scalar product with a function, integration or differentiation, which can be performed easily in tensor arithmetics on large 3D grids with 1D cost. Numerical tests illustrate the performance of the tensor summation method and confirm the estimated bounds on the tensor ranks.

Venera Khoromskaia, Boris N. Khoromskij

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\times L \times L$ lattices by using the canonical tensor decompositions and yielding the $O(L)$ computational complexity instead of $O(L^3)$ 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\times N \times 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 $\varepsilon$-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(N L)$. 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.