Block circulant and Toeplitz structures in the linearized Hartree-Fock equation on finite lattices: tensor approach

TL;DR

This paper presents a tensor-based grid approach for efficiently solving the linearized Hartree-Fock equation on 3D lattice systems, leveraging block-circulant and Toeplitz structures for computational speedup.

Contribution

It introduces a novel tensor technique exploiting block-circulant and Toeplitz structures to efficiently solve 3D lattice Hartree-Fock problems with reduced computational complexity.

Findings

Efficient FFT-based diagonalization for periodic systems.

Fast matrix-vector multiplication for finite systems.

Successful numerical validation on large lattice chains.

Abstract

This paper introduces and analyses the new grid-based tensor approach to approximate solution of the elliptic eigenvalue problem for the 3D lattice-structured systems. We consider the linearized Hartree-Fock equation over a spatial $L_1\times L_2\times L_3$ lattice for both periodic and non-periodic problem setting, discretized in the localized Gaussian-type orbitals basis. In the periodic case, the Galerkin system matrix obeys a three-level block-circulant structure that allows the FFT-based diagonalization, while for the finite extended systems in a box (Dirichlet boundary conditions) we arrive at the perturbed block-Toeplitz representation providing fast matrix-vector multiplication and low storage size. The proposed grid-based 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) in the periodic case 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. Lattice type systems in a box with Dirichlet boundary conditions are treated numerically by our previous tensor solver for single molecules, which makes possible calculations on rather large $L_1\times L_2\times L_3$ lattices due to reduced numerical cost for 3D problems. The numerical simulations for both box-type and periodic $L\times 1\times 1$ lattice chain in a 3D rectangular "tube" with $L$ up to several hundred confirm the theoretical complexity bounds for the block-structured eigenvalue solvers in the limit of large $L$.

Tables

Table 1. Table 3.1: CPU times (sec.): full eig-solver vs. FFT-based MBC diagonalization for L × 1 × 1 𝐿 1 1 L\times 1\times 1 lattice system, and with m 0 = 4 subscript 𝑚 0 4 m_{0}=4 , L = 2 p 𝐿 superscript 2 𝑝 L=2^{p} , p = 7 , 8 , … , 15 𝑝 7 8 … 15 p=7,8,...,15 .

Matrix size $N_b=m_{0}L$	$512$	$1024$	$2048$	$4096$	$8192$	$16384$	$32768$	$65536$	$131072$
Full EIG-solver	$0.67$	$5.49$	$48.6$	$497.4$	$--$	$--$	$--$	$--$	$--$
MBC diagonalization	$0.10$	$0.09$	$0.08$	$0.14$	$0.44$	$1.5$	$5.6$	$22.9$	$89.4$

Equations

$$\mathcal{H}\phi (x)\equiv [-\Delta +v(x)]\phi (x)=\lambda \phi (x),x\in \Omega \subset {\mathbb{R}}^d,d=1,2,3,$$

$$\left[-\frac{1}{2}\Delta -v_c(x)+{\int }_{{\mathbb{R}}^3}\frac{\rho (y)}{\Vert x-y\Vert }dy\right]{\phi }_i(x)-{\int }_{{\mathbb{R}}^3}\frac{\tau (x,y)}{\Vert x-y\Vert }{\phi }_i(y)dy={\lambda }_i{\phi }_i(x),x\in {\mathbb{R}}^3,$$

$$\mathcal{H}=-\frac{1}{2}\Delta -v_c,$$

$$v_c(x)=\sum _{\nu =1}^M\frac{Z_{\nu }}{\Vert x-a_{\nu }\Vert },Z_{\nu }>0,x,a_{\nu }\in {\mathbb{R}}^3,$$

$${\psi }_i=\sum _{\mu =1}^{N_b}{C}_{\mu i}{g}_{\mu },i=1,...,N_{orb},$$

$$h_{\mu \nu }=\frac{1}{2}{\int }_{{\mathbb{R}}^3}\nabla g_{\mu }\cdot \nabla g_{\nu }dx-{\int }_{{\mathbb{R}}^3}{v}_c(x)g_{\mu }{g}_{\nu }dx,1\le \mu ,\nu \le N_b,$$

$$HC$$

$$C^{T}SC$$

$$V_c=[v_{\mu \nu }]:v_{\mu \nu }={\int }_{{\mathbb{R}}^3}{v}_c(x)g_{\mu }{g}_{\nu }dx,1\le \mu ,\nu \le N_b,$$

$$\mathbf{P}:=[p_{\mathbf{i}}]\in {\mathbb{R}}^{n\times n\times n},p_{\mathbf{i}}={\int }_{{\mathbb{R}}^3}\frac{{\psi }_{\mathbf{i}}(x)}{\Vert x\Vert }\mathrm{d}x.$$

$$\frac{1}{z}=\frac{2}{\sqrt{\pi }}{\int }_{{\mathbb{R}}_{+}}{e}^{-z^2{t}^2}dt,z>0,$$

$$\mathbf{P}\approx {\mathbf{P}}_R=\sum _{q=1}^R{\mathbf{p}}_q^{(1)}\otimes {\mathbf{p}}_q^{(2)}\otimes {\mathbf{p}}_q^{(3)}\in {\mathbb{R}}^{n\times n\times n}.$$

$${\tilde{\mathbf{P}}}_R=\sum _{q=1}^R{\tilde{\mathbf{p}}}_q^{(1)}\otimes {\tilde{\mathbf{p}}}_q^{(2)}\otimes {\tilde{\mathbf{p}}}_q^{(3)}\in {\mathbb{R}}^{\tilde{n}\times \tilde{n}\times \tilde{n}}.$$

$${\mathbf{P}}_c=\sum _{\nu =1}^{M_0}{Z}_{\nu }{\mathbf{P}}_{c,\nu }\approx {\hat{\mathbf{P}}}_c\in {\mathbb{R}}^{n\times n\times n},$$

$${\mathbf{P}}_{c,\nu }={\mathcal{W}}_{\nu }{\tilde{\mathbf{P}}}_R=\sum _{q=1}^R{\mathcal{W}}_{\nu }^{(1)}{\tilde{\mathbf{p}}}_q^{(1)}\otimes {\mathcal{W}}_{\nu }^{(2)}{\tilde{\mathbf{p}}}_q^{(2)}\otimes {\mathcal{W}}_{\nu }^{(3)}{\tilde{\mathbf{p}}}_q^{(3)}\in {\mathbb{R}}^{n\times n\times n}.$$

$$v_{km}\approx {\int }_{\Omega }{V}_c(x){\overline{g}}_k(x){\overline{g}}_m(x)dx\approx ⟨{\mathbf{G}}_k\odot {\mathbf{G}}_m,{\hat{\mathbf{P}}}_c⟩=:V_{km},1\le k,m\le m_0.$$

$${\Omega }_L=B_1\times B_2\times B_3,\text{mbox}with{B}_{\ell }=\frac{1}{2}[-b_0{L}_{\ell }-b,b_0{L}_{\ell }+b]$$

$$k_{\ell }\in {\mathcal{K}}_{\ell }:=\{0,1,...,L_{\ell }-1\}.$$

$$v_{c_L}(x)=\sum _{\nu =1}^{M_0}\sum _{\mathbf{k}\in {\mathcal{K}}^3}\frac{Z_{\nu }}{\Vert x-a_{\nu }-b\mathbf{k}\Vert },x\in {\Omega }_L.$$

$${\mathbf{P}}_{c_L}=\sum _{\nu =1}^{M_0}{Z}_{\nu }\sum _{q=1}^R(\sum _{k_1\in {\mathcal{K}}_1}{\mathcal{W}}_{\nu (k_1)}{\tilde{\mathbf{p}}}_q^{(1)})\otimes (\sum _{k_2\in {\mathcal{K}}_2}{\mathcal{W}}_{\nu (k_2)}{\tilde{\mathbf{p}}}_q^{(2)})\otimes (\sum _{k_3\in {\mathcal{K}}_3}{\mathcal{W}}_{\nu (k_3)}{\tilde{\mathbf{p}}}_q^{(3)}),$$

$$V_{\mathbf{k}\mathbf{m}}(\mu ,\nu )=⟨{\mathbf{G}}_{\mathbf{k},\mu }\odot {\mathbf{G}}_{\mathbf{m},\nu },{\mathbf{P}}_{c_L}⟩,\mathbf{k},\mathbf{m}\in {\mathcal{K}}^3.$$

$$v_{km}={\int }_{{\mathbb{R}}^3}{v}_c(x){\overline{g}}_k(x){\overline{g}}_m(x)dx\approx ⟨{\mathbf{G}}_k\odot {\mathbf{G}}_m,{\mathbf{P}}_{c_L}⟩=:v_{km},1\le k,m\le N_b,$$

$$\begin{array}{rl}{V}_{\mathbf{k}\mathbf{m}}& =\sum _{\nu =1}^{M_0}{Z}_{\nu }\sum _{q=1}^R⟨{\mathbf{G}}_{\mathbf{k}}\odot {\mathbf{G}}_{\mathbf{m}},(\sum _{k_1\in \mathcal{K}}{\mathcal{W}}_{\nu (k_1)}{\tilde{\mathbf{p}}}_q^{(1)})\otimes (\sum _{k_2\in \mathcal{K}}{\mathcal{W}}_{\nu (k_2)}{\tilde{\mathbf{p}}}_q^{(2)})\otimes (\sum _{k_3\in \mathcal{K}}{\mathcal{W}}_{\nu (k_3)}{\tilde{\mathbf{p}}}_q^{(3)})⟩\\ & =\sum _{\nu =1}^{M_0}{Z}_{\nu }\sum _{q=1}^R\prod _{\ell =1}^3⟨{\mathbf{g}}_{\mathbf{k}}^{(\ell )}\odot {\mathbf{g}}_{\mathbf{m}}^{(\ell )},\sum _{k_{\ell }\in \mathcal{K}}{\mathcal{W}}_{\nu (k_{\ell })}{\tilde{\mathbf{p}}}_q^{(\ell )}⟩.\end{array}$$

$${\mathbf{G}}_{\mathbf{k}}\odot {\mathbf{G}}_{\mathbf{m}}=0\text{mbox}if|k_{\ell }-m_{\ell }|\ge L_0,$$

$$(I_{1}w)(x_{\ell }):=\sum _{i=1}^{\hat{n}}w(x_{i_{\ell }}){\xi }_i(x_{\ell }),x_{\mathbf{i}}\in {\Omega }_{\hat{n}},\ell =1,2,3.$$

$$A^{(\ell )}:=\{⟨\frac{d}{d{x}_{\ell }}{\xi }_i(x_{\ell }),\frac{d}{d{x}_{\ell }}{\xi }_j(x_{\ell })⟩{\}}_{i,j=1}^{\hat{n}}=\frac{1}{h}\text{mbox}tridiag\{-1,2,-1\},$$

$$S^{(\ell )}=\{⟨{\xi }_i,{\xi }_j⟩{\}}_{i,j=1}^{\hat{n}}=\frac{h}{6}\text{mbox}tridiag\{1,4,1\},\ell =1,2,3.$$

$${\otimes }_{(d⋎k)}(S^{(\ell )},A^{(k)})=S^{(1)}\otimes ...\otimes S^{(k-1)}\otimes A^{(k)}\otimes S^{(k+1)}\otimes ...\otimes S^{(d)}.$$

$$A_3:=A^{(1)}\otimes S^{(2)}\otimes S^{(3)}+S^{(1)}\otimes A^{(2)}\otimes S^{(3)}+S^{(1)}\otimes S^{(2)}\otimes A^{(3)}\equiv \sum _{k=1}^d{\otimes }_{(d⋎k)}(S^{(\ell )},A^{(k)}).$$

$$S_3=S^{(1)}\otimes S^{(2)}\otimes S^{(3)}\in {\mathbb{R}}^{{\hat{n}}^3\times {\hat{n}}^3}.$$

$$g_k(x)\approx {\hat{g}}_k(x):=\prod _{\ell =1}^3{\hat{g}}_k^{(\ell )}(x_{\ell })=\prod _{\ell =1}^3\sum _{i=1}^{\hat{n}}{g}_k^{(\ell )}(x_{i_{\ell }}){\xi }_i(x_{\ell }).$$

Full text

AMS Subject Classification:* * 65F30, 65F50, 65N35, 65F10

Key words: Tensor structured numerical methods for PDEs, 3D grid-based tensor approximation, Hartree-Fock equation, linearized Fock operator, periodic systems, lattice sum of potentials, block circulant/Toeplitz matrix, fast Fourier transform.

1 Introduction

Efficient numerical simulation of lattice systems for both periodic and non-periodic settings in application to crystalline, metallic, polymer-type compounds, and nano-structures is one of the challenging tasks in computational quantum chemistry. The reformulation of the nonlinear Hartree-Fock equation for periodic molecular systems based on the Bloch theory [1] has been addressed in the literature for more than forty years ago, and nowadays there are several implementations mostly relying on the analytic treatment of arising integral operators [19, 53, 23]. Mathematical analysis of spectral problems for PDEs with the periodic-type coefficients was an attractive topic in the recent decade, see [10, 20, 48] and the references therein. However, the systematic optimization of the basic numerical algorithms in the ab initio Hartree-Fock calculations for large lattice structured compounds with perturbed periodicity are largely unexplored.

The real space methods for single molecules based on the locally adaptive grids and multiresolution techniques have been discussed in [29, 56, 22, 11]. The grid-based tensor-structured approach for solving the Hartree-Fock nonlinear spectral problem approximated in the basis set of localized Gaussian-type-orbitals (GTO) has been developed and proved to be efficient for moderate size molecular systems [35, 42, 38].

This paper presents the grid-based tensor approach to the solution of elliptic eigenvalue problem for the 3D lattice-structured systems in a bounding box. We focus on the basic application to the linearized Hartree-Fock equation for extended systems composed of atoms or molecules located at nodes of a $L_{1}\times L_{2}\times L_{3}$ finite lattice, for both open boundary conditions and periodic supercell. The latter is useful because the structure of the respective Galerkin matrix (i.e., the Fock matrix) in the presence of defects can be treated as a small (local) perturbation to an ideally periodic system. We consider the 3D model eigenvalue problem for the Fock operator confined to the core Hamiltonian part, composed of the 3D Laplacian and the nuclear potential operator describing the Coulomb interaction of electrons and nuclei, which requires a sum of the total electrostatic potential of nuclei in the considered extended system. This is the typical example of a nontrivial elliptic eigenvalue problem arising in the numerical modeling of electronic structure in large almost periodic molecular systems. We observed in numerical experiments that there is an irreducible difference between the spectral data for periodic and non-periodic settings.

Computation of 3D lattice sums of a large number of long-distance Coulomb interaction potentials is one of the severe difficulties in the Hartree-Fock calculations for lattice-structured periodic or box-restricted systems. Traditionally this problem was treated by the so-called Ewald-type summation techniques [21, 17] combined with the fast multipole expansion or/and FFT methods [26, 52], which scale as $O(L^{3}\log L)$, $L=\max\{L_{1},L_{2},L_{3}\}$, for both periodic and box-type lattice sums. In this paper we apply the new, recently introduced method for summation of long-range potentials on lattices [38, 40] by using the assembled rank-structured tensor decomposition. This approach reduces the cost of summation for $L\times L\times L$ lattices to linear scaling in $L$, i.e. $O(L)$.

In the presented approach the Fock matrix is calculated directly by grid-based tensor numerical operations in the basis set of localized Gaussian-type-orbitals111GTO basis can be viewed as the special reduced basis constructed on the base of physical insight. (GTO) first specified by $m_{0}$ elements in the unit cell and then finitely replicated on 3D extended lattice structure in a box [35, 36]. For numerical integration by using low-rank tensor formats all basis functions are represented on the fine rectangular grid covering the whole computational box, where we introduce either the Dirichlet or periodic boundary conditions.

We show that in the case of finite lattices in a box (the Dirichlet boundary conditions) the core Hamiltonian exhibits the $C(2d+1)$-diagonal block sparsity, see Lemma 2.1. In particular both the discrete Laplacian and the mass matrix reveal the block-Toeplitz structure. The nuclear potential operator can be constructed, in general, as for the large single molecule in a box or by replication from the central unit cell to the whole lattice. In the latter case we arrive at the block-Toeplitz structure which allows the fast FFT-based matrix vector product.

For periodic boundary conditions (periodic supercell) we do not impose explicitly the periodicity-like features of the solution by means of the approximation ansatz that is the common approach in the Bloch formalism. Instead, the periodic properties of the considered system appear implicitly through the Toeplitz or circulant block structures in the Fock matrix. In case of periodic supercell the Fock matrix is proved to inherit the $d$-level symmetric block circulant form, that allows its diagonalization in the Fourier basis [32, 12] at the expense $O(m_{0}^{2}L^{d}\log L)$, $d=1,2,3$, see Lemma 3.4. In the case of $d$-dimensional lattice, the weak overlap between lattice translated basis functions leads to banded block sparsity thus reducing the storage cost. Furthermore, we introduce the low-rank tensor structure to the diagonal blocks of the Fock matrix represented in the Fourier space which allows to reduce the numerical cost to handle the block-circulant Galerkin matrix to linear scaling in $L$, $O(m_{0}^{2}L\log L)$, see Theorem 3.3.

The presented numerical scheme can be further investigated in the framework of the reduced Hartree-Fock model [10], where the similar block-structure in the Coulomb term of the Fock matrix can be observed. The Wannier-type basis functions constructed by the lattice translation of the localized molecular orbitals precomputed on the reference unit cell, can be also adapted to this algebraic framework.

The arising block-structured matrix representing the discretized core Hamiltonian, as well as some auxiliary function-related tensors arising, can be considered for further optimization by imposing the low-rank tensor formats, and in particular, the quantics-TT (QTT) tensor approximation [39] of long vectors and large matrices, which especially benefits in the limiting case of large $L\times L\times L$ perturbed periodic systems. In the QTT approach the algebraic operations on the 3D $n\times n\times n$ representation Cartesian grid can be implemented with logarithmic cost $O(\log n)$. Literature surveys on tensor algebra and rank-structured tensor methods for multi-dimensional PDEs can be found in [43, 42, 36, 13], see also [28, 25, 46, 14, 5, 4] and [36, 54, 55] concerning the low-rank decompositions in eigenvalue and electronic structure calculations. The present paper represents the revised and essentially extended version of the previous preprint [37].

Notice that in the recent years the analysis of eigenvalue problem solvers for large structured matrices has been widely discussed in the linear algebra community [8, 2, 9, 3]. Tensor structured approximation of elliptic equations with quasi–periodic coefficients has been considered in [44, 45].

The rest of the paper is organized as follows. Section 2 includes the main results on the analysis of core Hamiltonian on lattice structured compounds. In particular, §2.1 describes the tensor-structured calculation of the core Hamiltonian for large lattice-type molecular/atomic systems. We recall tensor-structured calculation of the Laplace operator and fast summation of lattice potentials by assembled canonical tensors. The complexity reduction due to low-rank tensor structures in the matrix blocks is discussed, see Remark 3.5). Section 3 discusses in detail the block circulant structure of the core Hamiltonian and presents numerical illustrations for a rectangular 3D “tube” of size $L\times 1\times 1$ with $L$ for large $L$. In particular, §3.2 introduces the new block structures by imposing the low-rank factorizations within multi-indexed blocks of the diagonalized three-level block-circulant matrix. We present a number of numerical experiments illustrating the pollution effect on the spectrum of periodic system compared with the system in a finite box. We also demonstrate the optimal performance for the direct FFT-based solver that implements the one-level block-circulant matrix structure describing the $L\times 1\times 1$ lattice systems for large $L$ (polymer-type compounds). Appendix recalls the classical results on the properties of block circulant/Toeplitz matrices and describes the basic tensor formats.

2 Elliptic operators with lattice-structured potentials

In this section we analyze the matrix structure of the Galerkin discretization for the elliptic eigenvalue problem in the form

[TABLE]

where the potential $v(x)$ is constructed by replication of those in the rectangular unit cell $\Omega_{0}$ over a $d$-dimensional rectangular $L_{1}\times L_{2}\times L_{3}$ lattice in a box, such that $\varphi\in H^{1}_{0}(\Omega)$, or in a rectangular supercell $\Omega$ with periodic boundary conditions. We focus on the important particular case of $v(x)=v_{c}(x)$ corresponding to the core Hamiltonian part in the Fock operator that constitutes the Hartree-Fock spectral problem arising in electronic structure calculations. In this case the electrostatic potential $v_{c}(x)$ is obtained as the large lattice sum of long-range interactions defined by the Newton kernel.

2.1 The Hartree-Fock core Hamiltonian in a GTO basis set

The nonlinear Fock operator ${\cal F}$ in the governing Hartree-Fock eigenvalue problem, describing the ground state energy for $2N_{b}$-electron system, is defined by

[TABLE]

where $i=1,...,N_{orb}$ and $\|\cdot\|$ means the distance function in $\mathbb{R}^{3}$. The linear part in the Fock operator is presented by the core Hamiltonian

[TABLE]

while the nonlinear Hartree potential and exchange operators depend on the unknown eigenfunctions (molecular orbitals) comprising the electron density, $\rho({y})=2\tau(y,y)$, and the density matrix, $\tau(x,y)=\sum\limits^{N_{orb}}_{i=1}\varphi_{i}(x)\varphi_{i}(y),\quad x,y\in\mathbb{R}^{3}.$ The electrostatic potential generated by the core Hamiltonian is defined by a sum

[TABLE]

where $M$ is the total number of nuclei in the system, $a_{\nu}$, $Z_{\nu}$, represent their Cartesian coordinates and the respective charge numbers.

Given a set of localized GTO basis functions $\{g_{\mu}\}$ ($\mu=1,...,N_{b}$), the occupied molecular orbitals $\psi_{i}$ are approximated in the form

[TABLE]

with the unknown coefficients matrix $C=[C_{\mu i}]\in\mathbb{R}^{N_{b}\times N_{orb}}$ obtained as the solution of the discretized Hartree-Fock equation with respect to the Galerkin basis $\{g_{\mu}\}$, and governed by $N_{b}\times N_{b}$ Fock matrix. The stiffness matrix $H=[h_{\mu\nu}]\in\mathbb{R}^{N_{b}\times N_{b}}$ of the core Hamiltonian (2.2) is represented by the single-electron integrals,

[TABLE]

such that the resulting eigenvalue equations governed by the reduced Fock matrix, $H$, read

[TABLE]

where the mass (overlap) matrix $S=[s_{\mu\nu}]_{1\leq\mu,\nu\leq N_{b}}$, is given by $s_{\mu\nu}=\int_{\mathbb{R}^{3}}g_{\mu}g_{\nu}dx.$

In the case of $L\times L\times L$ lattice system in a box the number of basis functions scales cubically in $L$, $N_{b}=m_{0}L^{3}$, hence the evaluation of the Fock matrix and further computations may become prohibitive as $L$ increases ($m_{0}$ is the number of basis functions in the unit cell). Moreover, the numerically extensive part in the matrix evaluation (2.5) is related to the integration with the large sum of lattice translated Newton kernels. Indeed, let $M_{0}$ be the number of nuclei in the unit cell, then the expensive calculations are due to the summation over $M_{0}L^{3}$ Newton kernels, and further spacial integration of this sum with the large set of localized atomic orbitals $\{g_{\mu}\}$, ($\mu=1,...,N_{b}$), where $N_{b}$ is of order $m_{0}L^{3}$.

In what follows, we describe the grid-based tensor approach for the block-structured representation of the core Hamiltonian in the Fock matrix for the lattice system in a box or in a periodic supercell. The main ingredients of the present approach include:

(a) the fast and accurate grid-based tensor method for evaluation of the electrostatic potential $v_{c}$ defined by the lattice sum in (2.3), see [38, 40];

(b) fast computation of the entries in the stiffness matrix $V_{c}$,

[TABLE]

by numerical grid-based integration using the low-rank tensor representation of all functions involved,

(c) block-structured factorized representation of the large and densely populated matrix $V_{c}$ in the form of perturbed multilevel block-circulant matrix, and

(d) block representation of the Galerkin matrix for the Laplacian.

The approach provides the opportunities to reduce computational costs in the case of large $L\times L\times L$ lattice systems. In the next sections, we show that in the periodic setting the resultant stiffness matrix $H=[h_{\mu\nu}]$ of the core Hamiltonian can be parametrized in the form of a symmetric, three-level block circulant matrix that allows further structural improvements by introducing tensor factorizations of the matrix blocks. In the case of lattice system in a box the block structure of $H$ is obtained by a small perturbation of the block Toeplitz matrix. These matrix structures allow the efficient storage and fast matrix-vector multiplication within iterations on a subspace for solving partial eigenvalue problem.

2.2 Nuclear potential operator for a single molecule

In this paragraph, we describe the evaluation of the stiffness matrix $V_{c}$ by tensor operations. It is based on the low-rank separable approximation to the nuclear (core) potential $v_{c}(x)$ representing the Coulomb interaction of the electrons with the nuclei, see (2.3).

In the case $L=1$ we have the single (discrete) molecule embedded into the scaled unit cell $\Omega=[-\frac{b}{2},\frac{b}{2}]^{3}$. In the computational domain $\Omega$, we introduce the uniform $n\times n\times n$ rectangular Cartesian grid $\Omega_{n}$ with the mesh size $h=b/n$, and define the set of tensor-product piecewise constant finite element basis functions $\{\psi_{\textbf{i}}\}$, which are supposed to be separable, i.e., $\psi_{\textbf{i}}(\textbf{x})=\prod_{\ell=1}^{d}\psi_{i_{\ell}}^{(\ell)}(x_{\ell})$ for ${\bf i}=(i_{1},i_{2},i_{3})$, $i_{\ell}\in I=\{1,...,n\}$.

Following [6, 35], the Newton kernel is discretized by the projection/collocation method in the form of a third order tensor $\mathbb{R}^{n\times n\times n}$, defined by

[TABLE]

The low-rank canonical decomposition of the $3$rd order tensor $\mathbf{P}$ is based on using exponentially convergent $\operatorname*{sinc}$-quadratures approximation of the Laplace-Gauss transform, [7, 58, 24, 27],

[TABLE]

which can be adapted to the Newton kernel by substitution $z=\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}$. We denote the resulting $R$-term canonical representation by

[TABLE]

We further suppose that all atomic centers are located strictly within subdomain $\Omega_{0}=[-\frac{b_{0}}{2},\frac{b_{0}}{2}]^{3}\subset\Omega$, $b_{0}<b$, called formation domain, and define the auxiliary (bounding) box $\widetilde{\Omega}\supset\Omega$, associated with the grid parameter $\widetilde{n}=n_{0}+n$ (say, $\widetilde{n}=2n$), see Figure 2.1.

Similar to (2.8), we introduce the auxiliary “reference tensor” $\widetilde{\bf P}_{R}\in\mathbb{R}^{\widetilde{n}\times\widetilde{n}\times\widetilde{n}}$, living on the grid $\Omega_{\widetilde{n}}$ and approximating the Newton kernel in $\widetilde{\Omega}$,

[TABLE]

The core potential $v_{c}(x)$ for a single molecule is approximated by a weighted sum of canonical tensors

[TABLE]

where the rank-$R$ tensor ${\bf P}_{{c},\nu}=={\cal W}_{\nu}\widetilde{\bf P}_{R}$ represents the single reference Coulomb potential in the form (2.9) shifted and restricted to $\Omega_{n}$ via the windowing operator ${\cal W}_{\nu}={\cal W}_{\nu}^{(1)}\otimes{\cal W}_{\nu}^{(2)}\otimes{\cal W}_{\nu}^{(3)}$, [38],

[TABLE]

Here every rank-$R$ canonical tensor ${\cal W}_{\nu}\widetilde{\bf P}_{R}\in\mathbb{R}^{n\times n\times n}$, $\nu=1,...,M_{0}$, is understood as a sub-tensor of the reference tensor obtained by a shift and restriction (windowing) of $\widetilde{\bf P}_{R}$ onto the $n\times n\times n$ grid $\Omega_{n}$ in the unit cell $\Omega$, $\Omega_{n}\subset\Omega_{\widetilde{n}}$. A shift from the origin is specified according to the coordinates of the corresponding nuclei, $a_{\nu}$, counted in the $h$-units.

The initial rank bound $rank({\bf P}_{c})\leq M_{0}R$ for the direct sum of canonical tensors in (2.10) can be improved (see [38], Remark 2.2). In the following, we denote by $\widehat{\bf P}_{c}$ the rank-$R_{c}$ ($R_{c}\leq M_{0}R$) canonical tensor obtained from ${\bf P}_{c}$ by the rank optimization procedure subject to certain threshold (in numerical tests we have $R_{c}\approx R$).

For the tensor representation of the Newton potentials, ${\bf P}_{{c},\nu}$, we make use of the piecewise constant discretization on the equidistant tensor grid, where, in general, the univariate grid size $n$ can be noticeably smaller than that used for the piecewise linear discretization applied to the Laplace operator. Indeed, the Galerkin approximation to the eigenvalue problem is constructed by using the global basis functions (reduced basis set $\{{g}_{k}\}$, $k=1,...,m_{0}$), hence the grid-based representation of these basis functions can be different in the calculation of the kinetic and potential parts in the Fock operator. The grid size $n$ is the only controlled by the approximation error for the integrals in (2.5) and by the numerical efficiency depending on the separation rank parameters.

Given tensor ${\bf P}_{c}$, the entries in the stiffness matrix $V_{c}$ in (2.6) can be evaluated by simple tensor operations. Fixed the GTO-type basis set $\{{g}_{k}\}$, $k=1,...,m_{0}$, i.e. $N_{b}=m_{0}$, defined in the scaled unit cell ${\Omega}$, where for ease of presentation functions ${g}_{k}$ are supposed to be separable. Introduce the corresponding rank-$1$ coefficients tensors ${\bf G}_{k}={\bf g}_{k}^{(1)}\otimes{\bf g}_{k}^{(2)}\otimes{\bf g}_{k}^{(3)}$ representing their piecewise constant approximations $\{\widehat{g}_{k}\}$ on the fine ${n}\times{n}\times{n}$ grid. Then the entries of the respective Galerkin matrix ${V}_{c}=[{v}_{km}]$ in (2.6) approximating the core potential operator $v_{c}$ in (2.3) are represented (approximately) by the following tensor operations,

[TABLE]

The error arising due to the separable $\varepsilon$-approximation of the discretized nuclear potential is controlled by the rank parameter $R_{c}=rank(\widehat{\bf P}_{c})$. Now letting $rank({\bf G}_{m})=1$ implies that each matrix element is to be computed with linear complexity in the univariate grid-size $n$, $O(R_{c}\,n)$. The almost exponential convergence of the tensor approximation in the separation rank $R_{c}$ leads to the asymptotic behavior of the $\varepsilon$-rank, $R_{c}=O(|\log\varepsilon|)$.

2.3 Nuclear potential operator for a lattice system in a box

Here we apply the previous constructions to the lattice structured location of nuclei. Given the potential sum $v_{c}(x)$ defined by (2.3) in the scaled unit cell $\Omega=[-\frac{b}{2},\frac{b}{2}]^{3}$ of size $b\times b\times b$, see Figure 2.1, we specify the smaller subdomain $\Omega_{0}=[-\frac{b_{0}}{2},\frac{b_{0}}{2}]^{3}\subset\Omega$ (called the formation cell) whose interior contains all atomic centers in the unit cell included into the summation in (2.3).

Let us consider an interaction potential in a symmetric computational box (supercell)

[TABLE]

consisting of a union of $L_{1}\times L_{2}\times L_{3}$ unit cells $\Omega_{\bf k}$, obtained by a shift of the reference domain $\Omega$ along the lattice vector ${\bf b_{0}}/2+b_{0}{\bf k}$, where ${\bf k}=(k_{1},k_{2},k_{3})\in\mathbb{Z}^{3}$, such that for $\ell=1,2,3$,

[TABLE]

In this notation the choice $L_{\ell}=1$ corresponds to the 3D one-layer system in the respective variable as illustrated in Figure 2.3. Figure 2.2 represents the $2$D projection of the 3D computational domain $\Omega_{L}$ for the $L_{1}\times 1\times 1$ molecular chain with $L_{1}=5$. Dashed regions correspond to the overlapping parts between shifted unit cells.

Figure 2.3 represents the geometry of the $3$D chain-type computational ”tube” $\Omega_{L}$.

For the discussion of complexity issues, we often consider a cubic lattice of equal sizes $L_{1}=L_{2}=L_{3}=L$. By the construction, we set $b=nh$ and $b_{0}=n_{0}h$, where the mesh-size $h>0$ is chosen the same for all spacial variables.

In the most interesting case of extended system in a box, further called case (B), the potential $v_{c_{L}}(x)$, for $x\in\Omega_{L}$, is obtained by summation over all unit cells $\Omega_{\bf k}$ in $\Omega_{L}$,

[TABLE]

Note that the direct calculation by (2.13) is performed at each of $L^{3}$ unit cells $\Omega_{\bf k}\subset\Omega_{L}$, ${\bf k}\in{\cal K}^{3}$, on a 3D lattice, which presupposes substantial numerical costs at least of the order of $O(L^{3})$ per unit cell.

The fast calculation of (2.13) is implemented by using the tensor summation method introduced in [38, 40] which can be described as follows. Let $\Omega_{N_{L}}$ be the $N_{L}\times N_{L}\times N_{L}$ uniform grid on $\Omega_{L}$ with the same mesh-size $h$ as above, and introduce the corresponding space of piecewise constant basis functions of the dimension $N_{L}^{3}$, where we have $N_{L}=n_{0}L+n-n_{0}$. Given the reference tensor in (2.9), the resultant lattice sum is presented by the canonical tensor ${\bf P}_{c_{L}}$

[TABLE]

whose rank is uniformly bounded $R_{c}\leq M_{0}R$. The numerical cost and storage size are bounded by $O(M_{0}RLN_{L})$, and $O(M_{0}RN_{L})$, respectively (see [38], Theorem 3.1), where $N_{L}=O(n_{0}L)$. The lattice sum in (2.14) converges only conditionally as $L\to\infty$. This aspect will be addressed in Section 3.4 following the approach discussed in [38, 40].

In the case of lattice system in a box, we define the basis set on a supercell $\Omega_{L}$ (and on $\widetilde{\Omega}_{L}$) by translation of the generating basis, defined in $\Omega_{0}$ for the single molecule, by the lattice vector $b{\bf k}$, i.e., $\{g_{\mu}({x})\}\mapsto\{g_{\mu}({x+b{\bf k}})\}$, $\mu=1,...,m_{0}$, where ${\bf k}=(k_{1},k_{2},k_{3})\in{\cal K}^{3}$, assuming zero extension of $\{g_{\mu}({x+b{\bf k}})\}$ beyond each local bounding box $\widetilde{\Omega}_{\bf k}$. The corresponding tensor representation of such functions is denoted by ${\bf G}_{{\bf k},\mu}$. The total number of basis functions for the lattice system is equal to $N_{b}=m_{0}L^{3}$.

In what follows, the matrix block entries of the $N_{b}\times N_{b}$ stiffness matrix $V_{c_{L}}$, corresponding to large basis set on a supercell $\Omega_{L}$, will be numbered by a pair of multi-indices, $V_{c_{L}}=[V_{{\bf k}{\bf m}}]$, where each $m_{0}\times m_{0}$ matrix block $V_{{\bf k}{\bf m}}$ is defined by

[TABLE]

This definition introduces the three-level block structure in the matrix $V_{c_{L}}$, which will be discussed in what follows.

In the practically interesting case of localized atomic orbitals (AO) basis, the matrix $V_{c_{L}}$ exhibits the banded block sparsity pattern since the effective support of localized AO associated with every unit cell $\Omega_{\bf k}\subset\widetilde{\Omega}_{\bf k}$ overlaps only fixed (small) number of neighboring cells. We call the number of overlapping neighboring cells by the overlap constant, $L_{0}$. The constant $L_{0}$ measures the essential overlap between basis functions in each spacial direction. For example, Figure 2.2 corresponds to the choice $L_{0}=2$.

Lemma 2.1

Assume that the overlap constant does not exceed $L_{0}$, then:

(a) The number of non-zero blocks in each block row (column) of the symmetric Galerkin matrix $V_{c_{L}}$ does not exceed $(2L_{0}+1)^{3}$.

(b) The storage size is bounded by $m_{0}^{2}[(L_{0}+1)L]^{3}$.

(c) The cost for evaluation of each $m_{0}\times m_{0}$ matrix block is bounded by $O(m_{0}^{2}M_{0}RN_{L})$.

Proof. In case (B), the matrix elements of $V_{c_{L}}=[v_{km}]\in\mathbb{R}^{N_{b}\times N_{b}}$ represented in (2.12), or in the block form in (2.15), can be expressed by the following tensor operations

[TABLE]

where again $\{\overline{g}_{k}\}$ denotes the piecewise constant representations to the respective Galerkin basis functions. This leads to the block representation (2.15) in terms of univariate vector operations

[TABLE]

Combining the block representation (2.15) and taking into account the overlapping property

[TABLE]

we are able to analyze the block sparsity pattern in the Galerkin matrix $V_{c_{L}}$. Given $3M_{0}R$ canonical vectors $\sum\limits_{k_{\ell}\in{\cal K}}{\cal W}_{\nu({k_{\ell}})}\widetilde{\bf p}^{(\ell)}_{q}\in\mathbb{R}^{N_{L}}$, where $N_{L}$ denotes the total number of grid points in $\Omega_{L}$ in each space variable. Now the numerical cost to compute ${v}_{km}$ for every fixed index $(k,m)$ is estimated by $O(M_{0}RN_{L})$ indicating linear scaling in the large grid parameter $N_{L}$ (but not cubic).

Fixed the row index in $(k,m_{\ast})$, then item (b) follows from the bound on the total number of overlapping cells $\Omega_{\bf k}$ in the effective integration domain in (2.16), that is $(2L_{0}+1)^{3}$, and from the symmetry of $V_{c_{L}}$.

Figure 2.4 illustrates the sparsity pattern of the nuclear potential contribution $V_{c_{L}}$ in the matrix $H$, computed for $L\times 1\times 1$ lattice in a $3$D supercell with $L=48$ and $m_{0}=4$, and the overlapping parameter $L_{0}=3$. In Figure (2.4), right one can observe the nearly-boundary effects due to the non-equalized contributions from the left and from the right (supercell in a box).

Notice that the quantized tensor approximation (QTT) of canonical vectors involved in ${\bf G}_{k}$ and ${\bf P}_{c_{L}}$ reduces this cost to the logarithmic scale, $O(M_{0}R\log N_{L})$, that is important in the case of large $L$ in view of $N_{L}=O(L)$, see the discussion in [38].

The block $L_{0}$-diagonal structure of the matrix $V_{c_{L}}=[V_{{\bf k}{\bf m}}]$, ${\bf k},{\bf m}\in{\cal K}^{3}$ described by Lemma 2.1 allows the essential saving in the storage costs.

However, the polynomial complexity scaling in $L$ leads to severe limitations on the number of unit cells. These limitations can be relax if we look more precisely on the defect between matrix ${V}_{c_{L}}$ and its block-circulant version corresponding to the periodic boundary conditions (see §3.3). This defect can be split into two components corresponding to their local and non-local features:

1. (A)

The non-local effect indicates the asymmetry in the interaction potential sum on the lattice in a box. 2. (B)

The near boundary (local) defect effects only those blocks in $V_{c_{L}}=\{V_{{\bf k}{\bf m}}\}$ lying in the $L_{0}$-width of $\partial\Omega_{L}$.

The defect in item (A) can be diminished by a slight modification of the core potential to the shift invariant Toeplitz-type form $V_{{\bf k}{\bf m}}=V_{|{\bf k}-{\bf m}|}$ by replication of the central unit cell to the whole lattice, as considered in Section 3. In this way the overlap condition (2.17) for the tensor ${\bf G}_{\bf k}$ will impose the $(2L_{0}+1)$ block diagonal sparsity in the block-Toeplitz matrix.

The boundary effect in item (B) becomes relatively small for large number of unit cells so that the block-circulant part of the matrix $V_{c_{L}}$ is getting dominating as $L\to\infty$.

The full diagonalization for the above mentioned matrix $V_{c_{L}}$ can be prohibitively expensive. However, the efficient storage and fast matrix-vector multiplication algorithms can be applied in the framework of structured iteration on subspace for calculation of a small subset of eigenvalues, see [4].

2.4 Discrete Laplacian and the mass matrix

In the case of a single molecule, the Laplace operator in (2.1), (2.2) is posed in the unit cell $\Omega=[-b/2,b/2]^{3}\in\mathbb{R}^{3}$, subject to the homogeneous Dirichlet boundary conditions on $\partial\Omega$. Periodic case corresponds to periodic boundary conditions. Given discretization parameter $\widehat{n}\in\mathbb{N}$, we use the equidistant $\widehat{n}\times\widehat{n}\times\widehat{n}$ tensor grid $\Omega_{\widehat{n}}=\{x_{\bf i}\}$, ${\bf i}\in{\cal I}:=\{1,...,\widehat{n}\}^{3}$, defined by the mesh-size $h=b/(\widehat{n}+1)$. This grid might be different from $\Omega_{n}$ introduced in §2.2 for representation of the interaction potential in the set of piecewise constant basis functions (usually, $n\leq\widehat{n}$).

Define a linear tensor-product interpolation operator ${\bf I}$ via the set of product hat functions, $\{\xi_{\bf i}:=\xi_{i_{1}}(x_{1})\xi_{i_{2}}(x_{2})\xi_{i_{3}}(x_{3})$, ${\bf i}\in{\cal I}\}$, associated with the respective grid-cells in $\Omega_{\widehat{n}}$. Here the linear interpolant ${\bf I}={I}_{1}\times{I}_{1}\times{I}_{1}$ is a product of 1D interpolation operators, where ${I}_{1}:C^{0}([-b,b])\to W_{h}:=span\{\xi_{i}\}_{i=1}^{\widehat{n}}$ is defined over the set of piecewise linear basis functions by

[TABLE]

Define the 1D FEM Galerkin stiffness (for Laplacian) and mass matrices $A^{(\ell)},S^{(\ell)}\in\mathbb{R}^{\widehat{n}\times\widehat{n}}$, respectively, by

[TABLE]

[TABLE]

For fixed dimension $d$ and $k\leq d$, introduce the mixed Kronecker product of matrices $S^{(\ell)}$ and $A^{(\ell)}$

[TABLE]

In the following, we apply the similar notations with respect to the Hadamard product of matrices $\odot$ and the usual multiplication operation, $\prod$.

Following [34], the rank-$3$ Kronecker tensor representation of the standard FEM Galerkin stiffness matrix for the Laplacian, $A_{3}\in\mathbb{R}^{\widehat{n}^{3}\times\widehat{n}^{3}}$, in the separable basis $\{\xi_{i}(x_{1})\xi_{j}(x_{2})\xi_{k}(x_{3})\}$, $i,j,k=1,\ldots\widehat{n}$, reads as

[TABLE]

In turn, the mass matrix takes the separable Kronecker product form

[TABLE]

For given GTO-type basis set $\{g_{k}(x)=g_{k}^{(1)}(x_{1})g_{k}^{(2)}(x_{2})g_{k}^{(3)}(x_{3})\}$ define a set of piecewise linear basis functions $\widehat{g}_{k}^{(\ell)}:={I}_{1}g_{k}^{(\ell)}$, $k=1,...,m_{0}$, and introduce the separable grid-based approximation of the initial basis functions $g_{k}(x)$,

[TABLE]

Here the rank-$1$ coefficients tensor ${\bf G}_{k}={\bf g}_{k}^{(1)}\otimes{\bf g}_{k}^{(2)}\otimes{\bf g}_{k}^{(3)}\in\mathbb{R}^{{\widehat{n}}^{\otimes 3}}$ given by the canonical vectors ${\bf g}_{k}^{(\ell)}=\{g_{k}^{(\ell)}(x_{i_{\ell}})\}$, ($k=1,...,m_{0}$) is associated with the Kronecker product of vectors, ${\bf g}_{k}=vec({\bf G}_{k})\in\mathbb{R}^{\widehat{n}^{3}}$. Let us agglomerate vectors ${\bf g}_{k}$ in a Kronecker product matrix $G={G}^{(1)}\otimes{G}^{(2)}\otimes{G}^{(3)}\in\mathbb{R}^{{\widehat{n}}^{3}\times m_{0}}$, where ${G}^{(\ell)}=[{\bf g}_{1}^{(\ell)},...,{\bf g}_{m_{0}}^{(\ell)}]\in\mathbb{R}^{{\widehat{n}}\times m_{0}}$, ($\ell=1,2,3$), is constructed by concatenation of vectors ${\bf g}_{k}^{(\ell)}$. Then the Galerkin stiffness matrix for the Laplacian and the mass matrix in the GTO basis set $\{{\bf G}_{k}\}$ can be written as

[TABLE]

corresponding to the standard matrix-matrix transform under the change of basis.

Applying the above representations to the $L\times L\times L$ lattice systems as described in §2.3 leads to the symmetric and sparce block-Toeplitz structure of the $N_{b}\times N_{b}$ Galerkin matrices with the block size $m_{0}\times m_{0}$ and with $N_{b}=m_{0}L^{3}$.

Proposition 2.2

Assume that the overlap constant does not exceed $L_{0}$, then:

(A) The cost for evaluation of each $m_{0}\times m_{0}$ matrix block is bounded by $O(m_{0}^{2}\widehat{n})$.

(B) The number of non-zero blocks in each block row (column) of the symmetric Galerkin matrices $A_{G}$ and $S_{G}$ does not exceed $(2L_{0}+1)^{3}$.

(C) Both $A_{G}$ and $S_{G}$ are symmetric $3$-level block-Toeplitz matrices. The storage size is bounded by $m_{0}^{2}(L_{0}+1)^{3}L^{3}$.

Proof. First, notice that the matrix entries in $A_{G}=\{a_{km}\}$ and $S_{G}=\{s_{km}\}$, ($k,m=1,...,m_{0}$) can be represented in the product form. For example, we have

[TABLE]

Combining this representation with (2.18) implies the matrix factorization

[TABLE]

where $\odot$ means the Hadamard product of matrices. The similar $d$-term sum of products representing matrix elements $a_{km}=\langle A_{3}{\bf g}_{k},{\bf g}_{m}\rangle$,

[TABLE]

leads to the $d$-term factorized representation of $A_{G}$ (say, $d=3$),

[TABLE]

This proves the numerical cost for the matrix evaluation. Items (B) and (C) can be justified by similar arguments as in Lemma 2.1.

Remark 2.3

Notice that in the periodic case both matrices, $A_{G}$ and $S_{G}$, possess the three-level block circulant structure as discussed in §3.3 (see Appendix for definitions).

3 Tensor factorization meets FFT block-diagonalization

There are two basic approaches to mathematical modeling of the $L$-periodic molecular systems composed of $L\times L\times L$ elementary unit cells [59]. In the first approach, the system is supposed to contain an infinite set of equivalent atoms that map identically into itself under any translation by $L$ units in each spacial direction. The other model is based on the ring-type periodic structures consisting of $L$ identical units in each spacial direction, where every unit cell of the periodic compound will be mapped to itself by applying a rotational transform from the corresponding rotational group symmetry.

The main difference between these two concepts is in the treatment of the lattice sum of Coulomb interactions, thought, at the limit of $L\to\infty$ both models approach each other. In this paper we mainly follow the first approach with the particular focus on the asymptotic complexity optimization for large lattice parameter $L$. The second concept is useful for understanding the block structure of the Galerkin matrices for the Hartree-Fock operator.

The direct Hartree-Fock calculations for lattice structured systems in the localized GTO-type basis lead to the symmetric block circulant/Toeplitz matrices, where the first-level blocks, $A_{0},...,A_{L-1}$, may have further block structures to be discussed in what follows (see also Appendix). In particular, the Galerkin approximation to the 3D Hartree-Fock core Hamiltonian in periodic setting leads to the symmetric, three-level block circulant matrix, see §5.2 concerning the definition of multilevel block circulant (MBC) matrices.

3.1 Block-diagonal form of the system matrix

In this paragraph, we introduce the new data-sparse block structure by imposing the low-rank tensor factorizations within the diagonalized MBC matrix in the matrix class ${\cal BC}(d,{\bf L},m_{0})$, where ${\bf L}=(L_{1},...,L_{d})$.

The block-diagonal form of a MBC matrix is well known in the literature, see e.g. [12]. A diagonalization of a $d$-level MBC matrix is based on representation via a sequence of cycling permutation matrices $\pi_{L_{1}},...,\pi_{L_{d}}$, $d=1,2,3,...$. Recall that the $d$-dimensional Fourier transform (FT) can be defined via the Kronecker product of the univariate FT matrices (Kronecker rank-$1$ operator),

[TABLE]

Here we prove the diagonal representation in a form that is useful for the description of tensor-based numerical algorithms. To that end we generalize the notations ${\cal T}_{L}$ and $\widehat{A}$ (see Appendix, §5.1) to the class of multilevel matrices. We denote by $\widehat{A}\in\mathbb{R}^{|{\bf L}|m_{0}\times m_{0}}$ the first block column of a matrix $A\in{\cal BC}(d,{\bf L},m_{0})$, with a shorthand notation

[TABLE]

and define a $|{\bf L}|\times m_{0}\times m_{0}$ tensor ${\cal T}_{\bf L}\widehat{A}$, which represents slice-wise all generating $m_{0}\times m_{0}$ matrix blocks in $\widehat{A}$ (reshaping of $\widehat{A}$). Notice that in the case $m_{0}=1$, the matrix $\widehat{A}\in\mathbb{R}^{|{\bf L}|\times 1}$ represents the first column of $A$. Now the Fourier transform $F_{\bf L}$ applies to ${\cal T}_{\bf L}\widehat{A}$ column-wise, while the backward reshaping of the resultant tensor, ${\cal T}_{\bf L}^{\prime}$, returns an $|{\bf L}|m_{0}\times m_{0}$ block matrix column. In the following we use the conventional matrix indexing and assume that the lattice ${\bf k}$-index runs as $k_{\ell}=0,1,...,L_{\ell}-1$.

Lemma 3.1

A matrix $A\in{\cal BC}(d,{\bf L},m_{0})$ can be converted to the block-diagonal form by the Fourier transform $F_{\bf L}$,

[TABLE]

where

[TABLE]

Proof. First, we confine ourself to the case of three-level matrices, i.e. $d=3$. We apply the Kronecker product decomposition (5.2) successively to each level of the block-circulant $A$ to obtain (see (5.3) for the definition of $\pi_{L}$)

[TABLE]

where ${A}_{k_{1}}\in\mathbb{R}^{L_{2}L_{3}m_{0}\times L_{2}L_{3}m_{0}}$, ${A}_{k_{1}k_{2}}\in\mathbb{R}^{L_{3}m_{0}\times L_{3}m_{0}}$ and $A_{k_{1}k_{2}k_{3}}\in\mathbb{R}^{m_{0}\times m_{0}}$.

Diagonalizing the periodic shift matrices $\pi_{L_{1}}^{k_{1}},\pi_{L_{2}}^{k_{2}}$, and $\pi_{L_{3}}^{k_{3}}$ via the 1D Fourier transform (see Appendix), we arrive at the block-diagonal representation

[TABLE]

where the monomials of diagonal matrices $D_{L_{\ell}}^{k_{\ell}}\in\mathbb{R}^{L_{\ell}\times L_{\ell}}$, $\ell=1,2,3$ are defined by (5.4). The generalization to the case $d>3$ can be proven by the similar argument.

Taking into account representation (5.11), the multilevel symmetric block circulant matrix can be described in form (3.1), such that all real-valued diagonal blocks remain symmetric.

The following remark compares the properties of circulant and Toeplitz matrices.

Remark 3.2

A block Toeplitz matrix does not allow explicit diagonalization by FT as it is the case for block circulant matrices. However, it is well known that a block Toeplitz matrix can be extended to the double-size (at each level) block circulant that makes it possible the efficient matrix-vector multiplication, and, in particular, the efficient application of power method for finding its senior eigenvalues.

3.2 Low-rank tensor structure within diagonalized block matrix

In the particular case $d=3$, the general block-diagonal representation (3.2) allows the reduced storage cost for the coefficients tensor $[A_{k_{1}k_{2}k_{3}}]$ to the order of $O(|{\bf L}|m_{0}^{2})$, where $|{\bf L}|=L_{1}L_{2}L_{3}$. Introduce the short notation $D_{\bf L}^{\bf k}=D_{L_{1}}^{k_{1}}\otimes D_{L_{2}}^{k_{2}}\otimes\cdots\otimes D_{L_{d}}^{k_{d}}$, then (3.2) takes a form

[TABLE]

For large $L$ the numerical cost becomes prohibitive. However, the above representation indicates that the further storage and complexity reduction can be possible if the third-order coefficients tensor ${\bf A}=[A_{k_{1}k_{2}k_{3}}]$, $k_{\ell}=0,...,L_{\ell}-1$, with the matrix-valued entries $A_{k_{1}k_{2}k_{3}}\in\mathbb{R}^{m_{0}\times m_{0}}$, allows the low-rank tensor factorization (approximation) in the multiindex ${\bf k}=(k_{1},k_{2},k_{3})$, which can be described by a smaller then $L^{3}$ number of parameters.

To fix the idea, let us assume the existence of rank-$1$ separable tensor factorization,

[TABLE]

Given $\ell\in\{1,...,d\}$ and a matrix $G\in\mathbb{R}^{L_{\ell}\times L_{\ell}}$, define the tensor prolongation (lifting) mapping, ${\cal P}_{\ell}:\mathbb{R}^{L_{\ell}\times L_{\ell}}\to\mathbb{R}^{|{\bf L}|\times|{\bf L}|}$, by

[TABLE]

The following theorem introduces the new multilevel block-circulant tensor-structured matrix format, where the coefficient tensor ${\bf A}$ is represented via the low-rank factorization.

Theorem 3.3

Assume the separability of a tensor $[A_{\bf k}]$ in the ${\bf k}$ space in the form (3.3), then the $3$-level block-circulant matrix $A$ can be represented in the factorized block-diagonal form as follows

[TABLE]

where the block-diagonal matrix $D_{A}$ with the block size $m_{0}\times m_{0}$ is given by

[TABLE]

with tri-tensors ${\bf A}^{(\ell)}=[A_{0}^{(\ell)},...,A_{L_{\ell}-1}^{(\ell)}]^{T}\in\mathbb{R}^{L_{\ell}\times m_{0}\times m_{0}}$ defined by concatenation of $\ell$-factors in (3.3).

Proof. The diagonal blocks in (3.2) can be written in the factorized tensor-product form

[TABLE]

Combining this representation with (3.2) leads to the powerful matrix factorization

[TABLE]

where the tensor prolongation operator ${\cal P}_{\ell}$ is given by (3.4).

The expansion (3.5) includes only 1D Fourier transforms thus reducing the representation cost to

[TABLE]

Moreover, and it is even more important, that the eigenvalue problem for the large matrix $A$ now reduces to only $L_{1}+L_{2}+L_{3}\ll L_{1}L_{2}L_{3}$ independent small $m_{0}\times m_{0}$ matrix eigenvalue problems.

The above block-diagonal representation for $d=3$ generalizes easily to the case of arbitrary dimension $d$. Furthermore, the rank-$1$ decomposition (3.3) was considered for the ease of exposition only. For instance, the above low-rank representations can be easily generalized to the case of canonical (CP) or Tucker formats in ${\bf k}$ space (see Proposition 3.5 below). In fact, both CP and Tucker formats provide the additive structure which can be converted to the respective additive structure of the core coefficient in (3.5).

Notice that in the practically interesting 3D case the use of MPS/TT type factorizations does not take the advantage over the Tucker format since the Tucker and MPS ranks in 3D appear to be close to each other. Indeed, the HOSVD for a tensor of order $3$ leads to the same sharp rank estimates for both the Tucker and TT tensor formats.

3.3 Block circulant structure in the periodic core Hamiltonian

In this section we consider the periodic case, further called case (P), and derive the more refined sparsity pattern of the matrix $V_{c_{L}}$ in (2.16) by using the $d$-level ($d=1,2,3$) tensor structure in this matrix. The matrix block entries are numbered by a pair of multi-indices, $V_{c_{L}}=\{V_{{\bf k}{\bf m}}\}$, ${\bf k}=(k_{1},k_{2},k_{3})$, where the $m_{0}\times m_{0}$ matrix block $V_{{\bf k}{\bf m}}$ is defined by (2.15). Figure 3.1 illustrates an example of 3D lattice-type structure of size $4\times 4\times 2$.

Following [38], we introduce the periodic cells ${\cal R}=\mathbb{Z}^{d}$, $d=1,2,3$ for the $\bf k$ index, and consider a 3D $B$-periodic supercell $\Omega_{L}=B\times B\times B$, with $B=\frac{b}{2}[-L,L]$. The total electrostatic potential in the supercell $\Omega_{L}$ is obtained by, first, the lattice summation of the Coulomb potentials over $\Omega_{L}$ for (rather large) $L$, but restricted to the central unit cell $\Omega_{0}$, and then by replication of the resultant function to the whole supercell $\Omega_{L}$. Hence, in this construction, the total potential sum $v_{c_{L}}(x)$ is designated at each elementary unit-cell in $\Omega_{L}$ by the same value (${\bf k}$-translation invariant). The electrostatic potential in each of $B$-periods can be obtained by copying the respective data from $\Omega_{L}$.

The effect of the conditional convergence of the lattice summation as $L\to\infty$ can be treated by using the extrapolation to the limit (regularization) on a sequence of different lattice parameters $L,2L,3L,\ldots$ as described in [38].

Consider the case $d=3$ in the more detail. Recall that the reference value $v_{c_{L}}(x)$ will be computed at the central cell $\Omega_{0}$, indexed by $(0,0,0)$, by summation over all contributions from $L^{3}$ elementary sub-cells in $\Omega_{L}$. For technical reasons here and in the following we vary the summation index by $k_{\ell}=0,...,L-1$, to obtain

[TABLE]

In the following, we use the same notations as in §2.3. The basis set in $\Omega_{L}$ is constructed by replication from the master unit cell $\Omega_{0}$ to the whole periodic lattice. The tensor representation of the local lattice sum on the $n\times n\times n$ grid associated with $\Omega_{0}$ takes a form

[TABLE]

where the tensor ${\bf P}_{\Omega_{0}}$ of size $n\times n\times n$ allows the low-rank expansion as in (2.14) with the reference tensor $\widetilde{\bf P}_{R}$ defined by (2.9). Here, the $\Omega$-windowing operator, ${\cal W}_{\nu({\bf k})}={\cal W}_{\nu(k_{1})}^{(1)}\otimes{\cal W}_{\nu(k_{2})}^{(2)}\otimes{\cal W}_{\nu(k_{3})}^{(3)},$ restricts onto the $n\times n\times n$ unit cell by shifting via the lattice vector ${\bf k}=(k_{1},k_{2},k_{3})$. This reduces both the computational and storage costs by a factor of $L$.

In the 3D case, we set $q=3$ in the notation for multilevel block-circulant (BC) matrix, see Appendix. Similar to the case of one-level BC matrices, we notice that a matrix $A\in{\cal BC}(3,{\bf L},m)$ of size $|{\bf L}|m\times|{\bf L}|m$ is completely defined by a $3$-rd order coefficients tensor ${\bf A}=[A_{k_{1}k_{2}k_{3}}]$ of size $L_{1}\times L_{2}\times L_{3}$, ($k_{\ell}=0,...,L_{\ell}-1$, $\ell=1,2,3$) with $m\times m$ block-matrix entries, obtained by folding of the generating first column vector in $A$.

Lemma 3.4

Assume that in case (P) the number of overlapping unit cells (in the sense of effective supports of basis functions) in each spatial direction does not exceed $L_{0}$. Then the Galerkin matrix $V_{c_{L}}=[V_{{\bf k}{\bf m}}]$ exhibits the symmetric, three-level block circulant Kronecker tensor-product form, i.e. $V_{c_{L}}\in{\cal BC}(3,{\bf L},m_{0})$,

[TABLE]

where the number of non-zero matrix blocks $A_{k_{1}k_{2}k_{3}}$ does not exceed $(L_{0}+1)^{3}$. Similar properties hold for both the Laplacian and the mass matrix.

The required storage is bounded by $m_{0}^{2}(L_{0}+1)^{3}$ independent of $L$. The set of non-zero generating matrix blocks $\{A_{k_{1}k_{2}k_{3}}\}$ can be calculated in $O(m_{0}^{2}(L_{0}+1)^{3}n)$ operations.

Furthermore, assume that the QTT ranks of the assembled canonical vectors do not exceed $r_{0}$. Then the numerical cost can be reduced to the logarithmic scale, $O(m_{0}^{2}(L_{0}+1)^{3}r_{0}^{2}\log n)$.

Proof. First, we notice that the shift invariance property in the matrix $V_{c_{L}}=[V_{{\bf k}{\bf m}}]$ is a consequence of the translation invariance in the canonical tensor ${\bf P}_{c_{L}}$ (periodic case), and in the tensor ${\bf G}_{\bf k}$ representing basis functions (by construction),

[TABLE]

so that we have

[TABLE]

This ensures the perfect three-level block-Toeplitz structure of $V_{c_{L}}$ (compare with the case of a bounded box). Now the block circulant pattern characterizing the class ${\cal BC}(3,{\bf L},m_{0})$ is imposed by the periodicity of a lattice-structured basis set.

To prove the complexity bounds we observe that a matrix $V_{c_{L}}\in{\cal BC}(3,{\bf L},m_{0})$ can be represented in the Kronecker tensor product form (3.7), obtained by an easy generalization of (5.2). In fact, we apply (5.2) by successive slice-wise and fiber-wise unfolding to obtain

[TABLE]

where ${\bf A}_{k_{1}}\in\mathbb{R}^{L_{2}\times L_{3}\times m_{0}\times m_{0}}$, ${\bf A}_{k_{1}k_{2}}\in\mathbb{R}^{L_{3}\times m_{0}\times m_{0}}$, and $A_{k_{1}k_{2}k_{3}}\in\mathbb{R}^{m_{0}\times m_{0}}$. Now the overlapping assumption ensures that the number of non-zero matrix blocks $A_{k_{1}k_{2}k_{3}}$ does exceed $(L_{0}+1)^{3}$.

Furthermore, the symmetric mass matrix, $S_{c_{L}}=\{{s}_{\mu\nu}\}\in\mathbb{R}^{N_{b}\times N_{b}}$, providing the Galerkin representation of the identity operator reads as follows,

[TABLE]

where $N_{b}=m_{0}L^{3}$. It can be seen that in the periodic case the block structure in the ”basis-tensor” ${\bf G}_{\bf k}$ imposes the three-level block circulant structure in the mass matrix $S_{c_{L}}$,

[TABLE]

By the previous arguments we conclude that $S_{k_{1}k_{2}k_{3}}=S^{(1)}_{k_{1}}S^{(2)}_{k_{2}}S^{(3)}_{k_{3}}$ implying the rank-$1$ separable representation in (3.10).

Likewise, it is easy to see that the stiffness matrix representing the (local) Laplace operator in the periodic setting has the similar block circulant structure,

[TABLE]

where the number of non-zero matrix blocks $B_{k_{1}k_{2}k_{3}}$ does not exceed $(L_{0}+1)^{3}$. In this case the matrix block $B_{k_{1}k_{2}k_{3}}$ admits a rank-$3$ product factorization inheriting the tri-term representation of the Laplacian.

In the Hartree-Fock calculations for lattice structured systems we deal with the multilevel, symmetric block circulant/Toeplitz matrices, where the first-level blocks, $A_{0},...,A_{L_{1}-1}$, may have further block structures. In particular, Lemma 3.4 shows that the Galerkin approximation of the 3D Hartree-Fock core Hamiltonian $H$ in periodic setting leads to the symmetric, three-level block circulant matrix structures.

Figure 3.2 shows the difference between matrices $V_{c_{L}}$ in periodic and non-periodic cases.

Figure 3.3 represents the block-sparsity in the core Hamiltonian matrix of a $L\times 1\times 1$ Hydrogen chain in a box with $L=32$ (right), and the matrix profile (left).

In the next section we discuss computational details of the FFT-based eigenvalue solver on the example of 3D linear chain of molecules.

3.4 Spectral problems in different settings: complexity analysis

Combining the block circulant representations (3.7), (3.11) and (3.10), we are able to represent the eigenvalue problem for the Fock matrix in the Fourier space as follows

[TABLE]

where $U=(F_{\bf L}\otimes I_{m})C$ and

[TABLE]

with the diagonal matrices $D_{L_{\ell}}^{k_{\ell}}\in\mathbb{R}^{L_{\ell}\times L_{\ell}}$, $\ell=1,2,3$ and the compact notation

[TABLE]

The equivalent block-diagonal form reads

[TABLE]

The block structure specified by Lemma 3.4 allows to apply the efficient eigenvalue solvers via FFT based diagonalization in the framework of Hartree-Fock calculations with the numerical cost $O(m_{0}^{2}L^{d}\log L)$.

Figure 3.4 visualizes molecular orbitals on fine spatial grid with $n=2^{14}$: the $4$th orbital (left), the $8$th orbital (right). The eigenvectors are computed in GTO basis for $L\times 1\times 1$ system with $L=128$ and $m_{0}=4$.

Remark 3.5

The low-rank structure in the coefficients tensor mentioned above (see Section 3.2) allows to reduce the factor $L^{d}\log L$ to $L\log L$ for $d=2,3$. It was already observed in the proof of Lemma 3.4 that the respective coefficients in the overlap and Laplacian Galerkin matrices can be treated as the rank-$1$ and rank-$3$ tensors, respectively. Clearly, the factorization rank for the nuclear part of the Hamiltonian does not exceed $R$. Hence, Theorem 3.3 can be applied in the generalized form.

Table 3.1 compares CPU times in sec. (Matlab) for the full eigenvalue solver on a 3D $L\times 1\times 1$ lattice in a box, and for the FTT-based MBC diagonalization in the periodic supercell, all computed for $m_{0}=4$, $L=2^{p}$ ($p=7,8,...,15$). The number of basis functions (problem size) is given by $N_{b}=m_{0}L$.

Figure 3.5 represents the spectrum of the core Hamiltonian in a box in comparison with that in a periodic supercell. We consider different number of cells $L=128,256$, where $m_{0}=4$. The systematic difference between the eigenvalues in both cases can be observed even at the limit of large $L$. This kind of spectral pollution effects have been discussed and theoretically analyzed in [10].

Figure 3.6 demonstrates the relaxation (for increasing $L$) of the average energy per unit cell with $m_{0}=4$, for a $L\times 1\times 1$ lattice structure in a 3D rectangular ”tube” up to $L=512$, for both periodic and open boundary conditions.

4 Conclusions

We introduced and analyzed the grid-based tensor approach to discretization and solution of the linearized Hartree-Fock equation in ab initio modeling of the lattice-structured molecular systems. We describe methods and algorithms for banded (finite box) and block-circulant (periodic setting) structured representation of the Fock matrix in GTO basis set (for the core Hamiltonian) and provide the numerical illustrations for both cases by implementing the algorithms in Matlab.

The sparse block structured representation to the Fock matrix combined with tensor techniques manifest several benefits: (a) the entries of the banded block structured Fock matrix are computed by 1D operations using low-rank CP tensors; (b) the storage size in the case of $L\times L\times L$ lattice is reduced to $O(L^{d})\ll L^{2d}$; (c) the $3$-level block-circulant Fock matrix in the periodic setting admits the low-rank tensor structure in the coefficients tensor, thus reducing the matrix diagonalization via conventional 3D FFT to the product of 1D Fourier transforms.

The main contributions include:

• •

Fast computation of the Fock matrix by 1D matrix-vector operations by using low-rank tensors associated with a 3D spacial grid.

• •

Analysis and numerical implementation of the multilevel banded/Toeplitz structure in the Fock matrix for the system in a box, and the block circulant structure in periodic setting.

• •

Establishing the low-rank tensor structure in the diagonal blocks of the Fock matrix represented in the Fourier space, that allows to reduce the storage size and diagonalization costs to $O(m_{0}^{3}\,d\,L\log L)$.

• •

Numerical tests illustrating the computational efficiency of the tensor-structured methods applied to the linearized Hartree-Fock equation for finite lattices and in periodic suprcell. Numerical experiments verify the theoretical results on the asymptotic complexity estimates of the proposed algorithms.

The rigorous numerical study of the nonlinear reduced Hartree-Fock eigenvalue problem for periodic and lattice-structured systems in a box is a subject of the future research.

5 Appendix

5.1 Overview on block circulant matrices

We recall that a one-level block circulant matrix $A\in{\cal BC}(L,m_{0})$ is defined by [12],

[TABLE]

where $A_{k}\in\mathbb{R}^{m_{0}\times m_{0}}$ for $k=0,1,\ldots,L-1$, are matrices of general structure. The equivalent Kronecker product representation is defined by the associated matrix polynomial,

[TABLE]

where $\pi=\pi_{L}\in\mathbb{R}^{L\times L}$ is the periodic downward shift (cycling permutation) matrix,

[TABLE]

and $\otimes$ denotes the Kronecker product of matrices.

In the case $m_{0}=1$ a matrix $A\in{\cal BC}(L,1)$ defines a circulant matrix generated by its first column vector $\widehat{a}=(a_{0},...,a_{L-1})^{T}$. The associated scalar polynomial then reads

[TABLE]

so that (5.2) simplifies to

[TABLE]

Let $\omega=\omega_{L}=\exp(-\frac{2\pi i}{L})$, we denote by

[TABLE]

the unitary matrix of Fourier transform. Since the shift matrix $\pi_{L}$ is diagonalizable in the Fourier basis,

[TABLE]

the same holds for any circulant matrix,

[TABLE]

where

[TABLE]

Conventionally, we denote by $\mbox{diag}\{x\}$ a diagonal matrix generated by a vector $x$. Let $X$ be an $Lm_{0}\times m_{0}$ matrix obtained by concatenation of $m_{0}\times m_{0}$ matrices $X_{k}$, $k=0,...,L-1$, $X=\operatorname{conc}(X_{0},...,X_{L-1})=[X_{0},...,X_{L-1}]^{T}$. For example, the first block column in (5.1) has the form $\operatorname{conc}(A_{0},...,A_{L-1})$. We denote by $\mbox{bdiag}\{X\}$ the $Lm_{0}\times Lm_{0}$ block-diagonal matrix of block size $L$ generated by $m_{0}\times m_{0}$ blocks $X_{k}$.

It is known that similarly to the case of circulant matrices (5.5), block circulant matrix in ${\cal BC}(L,m_{0})$ is unitary equivalent to the block diagonal one by means of Fourier transform via representation (5.2), see [12]. In the following, we describe the block-diagonal representation of a matrix $A\in{\cal BC}(L,m_{0})$ in the form that is convenient for generalization to the multi-level block circulant matrices as well as for the description of FFT based implementation schemes. To that end, let us introduce the reshaping (folding) transform ${\cal T}_{L}$ that maps a $Lm_{0}\times m_{0}$ matrix $X$ (i.e., the first block column in $A$) to $L\times m_{0}\times m_{0}$ tensor $B={\cal T}_{L}X$ by plugging the $i$th $m_{0}\times m_{0}$ block in $X$ into a slice $B(i,:,:)$. The respective unfolding returns the initial matrix $X={\cal T}_{L}^{\prime}B$. We denote by $\widehat{A}\in\mathbb{R}^{Lm_{0}\times m_{0}}$ the first block column of a matrix $A\in{\cal BC}(L,m_{0})$, with a shorthand notation

[TABLE]

so that the $L\times m_{0}\times m_{0}$ tensor ${\cal T}_{L}\widehat{A}$ represents slice-wise all generating $m_{0}\times m_{0}$ matrix blocks.

Proposition 5.1

For $A\in{\cal BC}(L,m_{0})$ we have

[TABLE]

where

[TABLE]

can be recognized as the $j$-th $m_{0}\times m_{0}$ matrix block in block column ${\cal T}_{L}^{\prime}(F_{L}({\cal T}_{L}\widehat{A}))$, such that

[TABLE]

A set of eigenvalues $\lambda$ of $A$ is then given by

[TABLE]

The eigenvectors corresponding to the spectral sets

[TABLE]

can be represented in the form

[TABLE]

with $E_{[j]}=\operatorname{diag}\{e_{j}\}\otimes I_{m_{0}}\in\mathbb{R}^{Lm_{0}\times Lm_{0}}$, and $e_{j}\in\mathbb{R}^{L}$ being the $j$th Euclidean basis vector.

Proof. We combine representations (5.2) and (5.4) to obtain

[TABLE]

where the final step follows by the definition of FT matrix and by the construction of ${\cal T}_{L}$. The structure of eigenvalues and eigenfunctions then follows by straightforward calculations with block-diagonal matrices.

The next statement describes the block-diagonal form for a class of symmetric BC matrices, ${\cal BC}_{s}(L,m_{0})$, that is a simple corollary of [12], Proposition 5.1. In this case we have $A_{0}=A_{0}^{T}$, and $A_{k}^{T}=A_{L-k}$, $k=1,...,L-1$.

Corollary 5.2

Let $A\in{\cal BC}_{s}(L,m_{0})$ be symmetric, then $A$ is unitary similar to a Hermitian block-diagonal matrix, i.e., $A$ is of the form

[TABLE]

where $I_{m_{0}}$ is the $m_{0}\times m_{0}$ identity matrix. The matrices $\tilde{A}_{j}\in\mathbb{C}^{m_{0}\times m_{0}}$, $j=0,1,\ldots,L-1$, are defined for even $n\geq 2$ as

[TABLE]

Corollary 5.2 combined with Proposition 5.1 describes a simplified structure of spectral data in the symmetric case. Notice that the above representation imposes the symmetry of each real-valued diagonal blocks $\tilde{A}_{j}\in\mathbb{R}^{m_{0}\times m_{0}},\;j=0,1,\ldots,L-1$, in (5.10).

5.2 Multilevel block circulant/Toeplitz matrices

Furthermore, we describe the extension of (one-level) block circulant matrices to multilevel structure. First, we recall the main notions of multilevel block circulant (MBC) matrices with the particular focus on the three-level case. Given a multi-index ${\bf L}=(L_{1},L_{2},L_{3})$, we denote $|{\bf L}|=L_{1}\,L_{2}\,L_{3}$. A matrix class ${\cal BC}(d,{\bf L},m_{0})$ ($d=1,2,3$) of $d$-level block circulant matrices can be introduced by the following recursion.

Definition 5.3

For $d=1$, define a class of one-level block circulant matrices by ${\cal BC}(1,{\bf L},m)\equiv{\cal BC}(L_{1},m)$ (see §5.1), where ${\bf L}=(L_{1},1,1)$. For $d=2$, we say that a matrix $A\in\mathbb{R}^{|{\bf L}|m_{0}\times|{\bf L}|m_{0}}$ belongs to a class ${\cal BC}(d,{\bf L},m_{0})$ if

[TABLE]

where ${\bf L}_{[1]}=(L_{2},L_{3})\in\mathbb{N}^{d-1}$. Similar recursion applies to the case $d=3$.

Likewise to the case of one-level BC matrices, it can be seen that a matrix $A\in{\cal BC}(d,{\bf L},m_{0})$, $d=1,2,3$, of size $|{\bf L}|m_{0}\times|{\bf L}|m_{0}$ is completely defined (parametrized) by a $d$th order matrix-valued tensor ${\bf A}=[A_{k_{1}...k_{d}}]$ of size $L_{1}\times...\times L_{d}$, ($k_{\ell}=1,...,L_{\ell}$, $\ell=1,...,d$), with $m_{0}\times m_{0}$ matrix entries $A_{k_{1}...k_{d}}$, obtained by folding of the generating first column vector in $A$.

Recall that a symmetric block Toeplitz matrix $A\in{\cal BT}_{s}(L,m_{0})$ is defined by [12],

[TABLE]

where $A_{k}\in\mathbb{R}^{m_{0}\times m_{0}}$ for $k=0,1,\ldots,L-1$, is a matrix of a general structure.

Similar to Definition 5.3, a matrix class ${\cal BT}_{s}(d,{\bf L},m_{0})$ of symmetric $d$-level block Toeplitz matrices can be introduced by the following recursion.

Definition 5.4

For $d=1$, ${\cal BT}_{s}(1,{\bf L},m_{0})\equiv{\cal BT}_{s}(L_{1},m_{0})$ is the class of one-level symmetric block circulant matrices with ${\bf L}=(L_{1},1,1)$. For $d=2$ we say that a matrix $A\in\mathbb{R}^{|{\bf L}|m\times|{\bf L}|m_{0}}$ belongs to a class ${\cal BT}_{s}(d,{\bf L},m_{0})$ if

[TABLE]

Similar recursion applies to the case $d=3$.

5.3 Rank-structured tensor formats

We consider a tensor of order $d$, as a multidimensional array numbered by a $d$-tuple index set, ${\bf A}=[a_{i_{1},...,i_{d}}]\;\in\mathbb{R}^{n_{1}\times\ldots\times n_{d}}$. A tensor is an element of a linear vector space equipped with the Euclidean scalar product. In particular, a tensor with equal sizes $n_{\ell}=n$, $\ell=1,\ldots d$, is denoted as $n^{\otimes d}$ tensor. The required storage for entry-wise representation of tensors scales exponentially in the dimension, $n^{d}$, (the so-called ”curse of dimensionality“). To get rid of exponential scaling in the dimension, one can apply the rank-structured separable representations of multidimensional tensors.

The rank-$1$ canonical tensor, ${\bf A}={\bf u}^{(1)}\otimes...\otimes{\bf u}^{(d)}\in\mathbb{R}^{n_{1}\times\ldots\times n_{d}}$, with entries $a_{i_{1},\ldots i_{d}}=a^{(1)}_{i_{1}}\cdot\cdot\cdot a^{(d)}_{i_{d}}$ requires only $dn$ numbers to store it. A tensor in the $R$-term canonical format (CP tensors) is defined by the parametrization

[TABLE]

where $u_{k}^{(\ell)}$ are normalized vectors, and $R$ is called the canonical rank of a tensor. The storage size is bounded by $dnR\ll n^{d}$.

Given the rank parameter ${\bf r}=(r_{1},...,r_{d})$, a tensor in the rank-$\bf r$ Tucker format is defined by the parametrization

[TABLE]

completely specified by a set of orthonormal vectors ${\bf v}^{(\ell)}_{\nu_{\ell}}\in\mathbb{R}^{n_{\ell}}$, and the Tucker core tensor $\boldsymbol{\beta}=[\beta_{\nu_{1},...,\nu_{d}}]$. The storage demand is bounded by $|{\bf r}|+(r_{1}+\ldots+r_{d})n$.

The remarkable approximating properties of the Tucker and canonical tensor decomposition applied to the wide class of function related tensors were revealed in [41, 27, 33], promoting using tensor tools for the numerical treatment of the multidimensional PDEs. It was proved for the CP/Tucker decomposition of some classes of function related tensors that the rank-$\bf r$ Tucker approximation with $r=O(\log n)$, i.e., $r_{\ell}\ll n$, provides the exponentially small error of order $e^{-\alpha r}$, [41] (here we discuss for simplicity the equal rank decompositions, $r_{\ell}=r$).

In the case of many spacial dimensions the product type tensor formats provide the stable rank-structured approximation. The matrix-product states (MPS) decomposition was long since used in quantum chemistry and quantum information theory, see the survey paper [57]. The particular case of MPS representation is called a tensor train (TT) format [50]. The quantics-TT (QTT) tensor approximation method for functional $n$-vectors was introduced in [39] and shown to provide the logarithmic complexity, $O(d\log n)$, on the wide class of generating functions. Furthermore, a combination of different tensor formats proved to be successful in the numerical solution of the multidimensional PDEs [42, 15].