Skip to content
Licensed Unlicensed Requires Authentication Published by De Gruyter May 31, 2017

Block Circulant and Toeplitz Structures in the Linearized Hartree–Fock Equation on Finite Lattices: Tensor Approach

Venera Khoromskaia and Boris N. Khoromskij


This paper introduces and analyzes the new grid-based tensor approach to approximate solutions of the elliptic eigenvalue problem for the 3D lattice-structured systems. We consider the linearized Hartree–Fock equation over a spatial L1×L2×L3 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 L1×L2×L3 lattices due to reduced numerical cost for 3D problems. The numerical simulations for both box-type and periodic L×1×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.

A Appendix

A.1 Overview on Block Circulant Matrices

We recall that a one-level block circulant matrix A𝒞(L,m0) is defined by


where Akm0×m0 for k=0,1,,L-1 are matrices of general structure (see [14]). The equivalent Kronecker product representation is defined by the associated matrix polynomial,


where π=πLL×L is the periodic downward shift (cycling permutation) matrix,


and denotes the Kronecker product of matrices.

In the case m0=1, a matrix A𝒞(L,1) defines a circulant matrix generated by its first column vector a^=(a0,,aL-1)T. The associated scalar polynomial then reads


so that (A.2) simplifies to


Let ω=ωL=exp(-2πiL) and denote by

FL={fk}L×Lwith fk=1LωL(k-1)(-1),k,=1,,L,

the unitary matrix of Fourier transform. Since the shift matrix πL is diagonalizable in the Fourier basis,


the same holds for any circulant matrix,




Conventionally, we denote by diag{x} a diagonal matrix generated by a vector x. Let X be an Lm0×m0 matrix obtained by concatenation of m0×m0 matrices Xk, k=0,,L-1,


For example, the first block column in (A.1) has the form conc(A0,,AL-1). We denote by bdiag{X} the Lm0×Lm0 block-diagonal matrix of block size L generated by m0×m0 blocks Xk.

It is known that similarly to the case of circulant matrices (A.5), a block circulant matrix in 𝒞(L,m0) is unitary equivalent to the block diagonal one by means of Fourier transform via representation (A.2), see [14]. In the following, we describe the block-diagonal representation of a matrix A𝒞(L,m0) 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 𝒯L that maps an Lm0×m0 matrix X (i.e., the first block column in A) to the L×m0×m0 tensor B=𝒯LX by plugging the ith m0×m0 block in X into a slice B(i,:,:). The respective unfolding returns the initial matrix X=𝒯LB. We denote by A^Lm0×m0 the first block column of a matrix A𝒞(L,m0), with a shorthand notation


so that the L×m0×m0 tensor 𝒯LA^ represents slice-wise all generating m0×m0 matrix blocks.

Proposition A.1.

For ABC(L,m0) we have




can be recognized as the j-th m0×m0 matrix block in block column TL(FL(TLA^)), such that


A set of eigenvalues λ of A is then given by


The eigenvectors corresponding to the spectral sets


can be represented in the form

Uj,m=(FLIm)U¯j,m,where U¯j,m=E[j]vec[u0,m,u1,m,,uL-1,m],

with E[j]=diag{ej}Im0RLm0×Lm0, and ejRL being the jth Euclidean basis vector.


We combine representations (A.2) and (A.4) to obtain


where the final step follows by the definition of FT matrix and by the construction of 𝒯L. The structure of eigenvalues and eigenfunctions then follows by straightforward calculations with block-diagonal matrices. This concludes the proof. ∎

The next statement describes the block-diagonal form for a class of symmetric BC matrices, 𝒞s(L,m0). It is a simple corollary of [14] and Proposition A.1. In this case we have A0=A0T, and AkT=AL-k, k=1,,L-1.

Corollary A.2.

Let ABCs(L,m0) be symmetric. Then A is unitary similar to a Hermitian block-diagonal matrix, i.e., A is of the form


where Im0 is the m0×m0 identity matrix. The matrices A~jCm0×m0, j=0,1,,L-1, are defined for even n2 as


Corollary A.2 combined with Proposition A.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 block A~jm0×m0, j=0,1,,L-1, in (A.6).

A.2 Multilevel Block Circulant/Toeplitz Matrices

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 𝐋=(L1,L2,L3), we denote |𝐋|=L1L2L3. A matrix class 𝒞(d,𝐋,m0) (d=1,2,3) of d-level block circulant matrices can be introduced by the following recursion.

Definition A.3.

For d=1, define a class of one-level block circulant matrices by 𝒞(1,𝐋,m)𝒞(L1,m) (see Appendix A.1), where 𝐋=(L1,1,1). For d=2, we say that a matrix A|𝐋|m0×|𝐋|m0 belongs to a class 𝒞(d,𝐋,m0) if

A=bcirc(A1,,AL1)with Aj𝒞(d-1,𝐋[1],m0),j=1,,L1,

where 𝐋[1]=(L2,L3)d-1. A 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𝒞(d,𝐋,m0), d=1,2,3, of size |𝐋|m0×|𝐋|m0 is completely defined (parametrized) by a dth order matrix-valued tensor 𝐀=[Ak1kd] of size L1××Ld (k=1,,L, =1,,d) with m0×m0 matrix entries Ak1kd, obtained by folding the generating first column vector in A.

Recall that a symmetric block Toeplitz matrix A𝒯s(L,m0) is defined by


where Akm0×m0 for k=0,1,,L-1 is a matrix of a general structure (see [14]).

A matrix class 𝒯s(d,𝐋,m0) of symmetric d-level block Toeplitz matrices can be introduced by the following recursion, similarly to Definition A.3.

Definition A.4.

For d=1, 𝒯s(1,𝐋,m0)𝒯s(L1,m0) is the class of one-level symmetric block circulant matrices with 𝐋=(L1,1,1). For d=2 we say that a matrix A|𝐋|m×|𝐋|m0 belongs to a class 𝒯s(d,𝐋,m0) if

A=btoepls(A1,,AL1)with Aj𝒯s(d-1,𝐋[𝟏],m0),j=1,,L1.

A similar recursion applies to the case d=3.

A.3 Rank-Structured Tensor Formats

We consider a tensor of order d, as a multidimensional array numbered by a d-tuple index set, 𝐀=[ai1,,id]n1××nd. A tensor is an element of a linear vector space equipped with the Euclidean scalar product. In particular, a tensor with equal sizes n=n, =1,,d, is called an nd tensor. The required storage for entry-wise representation of tensors scales exponentially in the dimension, nd (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 𝐀=𝐮(1)𝐮(d)n1××nd with entries ai1,,id=ai1(1)aid(d) requires only dn numbers to store it. A tensor in the R-term canonical format (CP tensors) is defined by the parametrization


where uk() are normalized vectors, and R is called the canonical rank of a tensor. The storage size is bounded by dnRnd.

Given the rank parameter 𝐫=(r1,,rd), a tensor in the rank-𝐫 Tucker format is defined by the parametrization


completely specified by a set of orthonormal vectors 𝐯ν()n, and the Tucker core tensor 𝜷=[βν1,,νd]. The storage demand is bounded by |𝐫|+(r1++rd)n.

The remarkable approximating properties of the Tucker and canonical tensor decomposition applied to the wide class of function-related tensors were revealed in [39, 27, 42], promoting using tensor tools for the numerical treatment of the multidimensional PDEs. It was proved for some classes of function related tensors that the rank-𝐫 Tucker approximation provides the exponentially small error of the order of e-αr with r=minr, where r=O(logn) (see [39]).

In the case of many spacial dimensions, the product type tensor formats provide the stable rank-structured approximation. The matrix-product states (MPS) decomposition has been for a long time 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 [51, 50]. The quantics-TT (QTT) tensor approximation method for functional n-vectors was introduced in [40] and shown to provide the logarithmic complexity, O(dlogn), 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 [41, 15].


[1] P. Benner, S. Dolgov, V. Khoromskaia and B. N. Khoromskij, Fast iterative solution of the Bethe-Salpeter eigenvalue problem using low-rank and QTT tensor approximation, J. Comput. Phys. 334 (2017), 221–239. 10.1016/ in Google Scholar

[2] P. Benner, H. Faßbender and C. Yang, Some remarks on the complex J-symmetric eigenproblem, preprint (2015), 10.1016/j.laa.2018.01.014Search in Google Scholar

[3] P. Benner, V. Khoromskaia and B. N. Khoromskij, A reduced basis approach for calculation of the Bethe–Salpeter excitation energies using low-rank tensor factorizations, Mol. Phys. 114 (2016), no. 7–8, 1148–1161. 10.1080/00268976.2016.1149241Search in Google Scholar

[4] P. Benner, V. Mehrmann and H. Xu, A new method for computing the stable invariant subspace of a real Hamiltonian matrix, J. Comput. Appl. Math. 86 (1997), 17–43. 10.1016/S0377-0427(97)00146-5Search in Google Scholar

[5] C. Bertoglio and B. N. Khoromskij, Low-rank quadrature-based tensor approximation of the Galerkin projected Newton/Yukawa kernels, Comput. Phys. Commun. 183 (2012), no. 4, 904–912. 10.1016/j.cpc.2011.12.016Search in Google Scholar

[6] A. Bloch, Les théorèmes de M. Valiron sur les fonctions entières et la théorie de l’uniformisation, Ann. Fac. Sci. Toulouse Math. 17 (1925), no. 3, 1–22. 10.5802/afst.335Search in Google Scholar

[7] D. Braess, Asymptotics for the approximation of wave functions by exponential-sums, J. Approx. Theory 83 (1995), 93–103. 10.1006/jath.1995.1110Search in Google Scholar

[8] A. Bunse-Gerstner, R. Byers and V. Mehrmann, A chart of numerical methods for structured eigenvalue problems, SIAM J. Matrix Anal. Appl. 13 (1992), 419–453. 10.1007/978-3-642-75536-1_39Search in Google Scholar

[9] A. Bunse-Gerstner and H. Faßbender, Breaking Van Loan’s curse: A quest for structure-preserving algorithms for dense structured eigenvalue problems, Numerical Algebra, Matrix Theory, Differential-Algebraic Equations and Control Theory, Springer, Cham (2015), 3–23. 10.1007/978-3-319-15260-8_1Search in Google Scholar

[10] E. Cancés, A. Deleurence and M. Lewin, A new approach to the modeling of local defects in crystals: The reduced Hartree–Fock case, Comm. Math. Phys. 281 (2008), 129–177. 10.1007/s00220-008-0481-xSearch in Google Scholar

[11] E. Cancés, V. Ehrlacher and Y. Maday, Periodic Schrödinger operator with local defects and spectral pollution, SIAM J. Numer. Anal. 50 (2012), no. 6, 3016–3035. 10.1137/110855545Search in Google Scholar

[12] A. Cichocki, N. Lee, I. Oseledets, A. H. Phan, Q. Zhao and D. P. Mandic, Tensor networks for dimensionality reduction and large-scale optimization: Part 1 low-rank tensor decompositions, Found. Trends Mach. Learn. 9 (2016), no. 4–5, 249–429. 10.1561/2200000059Search in Google Scholar

[13] T. Darten, D. York and L. Pedersen, Particle mesh Ewald: An O(NlogN) method for Ewald sums in large systems, J. Chem. Phys. 98 (1993), 10089–10092. 10.1063/1.464397Search in Google Scholar

[14] J. P. Davis, Circulant Matrices, John Wiley & Sons, New York, 1979. Search in Google Scholar

[15] S. Dolgov and B. N. Khoromskij, Two-level QTT-Tucker format for optimized tensor calculus, SIAM J. Matrix Anal. Appl. 34 (2013), no. 2, 593–623. 10.1137/120882597Search in Google Scholar

[16] S. Dolgov, B. N. Khoromskij, D. Savostyanov and I. Oseledets, Computation of extreme eigenvalues in higher dimensions using block tensor train format, Comput. Phys. Commun. 185 (2014), no. 4, 1207–1216. 10.1016/j.cpc.2013.12.017Search in Google Scholar

[17] R. Dovesi, R. Orlando, C. Roetti, C. Pisani and V. R. Sauders, The periodic Hartree–Fock method and its implementation in the CRYSTAL code, Phys. Stat. Sol. (b) 217 (2000), 63–88. 10.1002/3527603107.ch4Search in Google Scholar

[18] T. H. Dunning, Jr., Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen, J. Chem. Phys. 90 (1989), 1007–1023. 10.1063/1.456153Search in Google Scholar

[19] V. Ehrlacher, C. Ortner and A. V. Shapeev, Analysis of boundary conditions for crystal defect atomistic simulations, Arch. Ration. Mech. Anal. 222 (2016), no. 3, 1217–1268. 10.1007/s00205-016-1019-6Search in Google Scholar

[20] P. P. Ewald, Die Berechnung optische und elektrostatischer Gitterpotentiale, Ann. Phys. 369 (1921), no. 3, 253–287. 10.1002/andp.19213690304Search in Google Scholar

[21] H. Faßbender and D. Kressner, Structured eigenvalue problem, GAMM-Mitt. 29 (2006), no. 2, 297–318. 10.1002/gamm.201490035Search in Google Scholar

[22] L. Frediani, E. Fossgaard, T. Flå and K. Ruud, Fully adaptive algorithms for multivariate integral equations using the non-standard form and multiwavelets with applications to the Poisson and bound-state Helmholtz kernels in three dimensions, Mol. Phys. 111 (2013), 9–11. 10.1080/00268976.2013.810793Search in Google Scholar

[23] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci and G. A. Petersson, Gaussian Development Version Revision H1, Gaussian Inc., Wallingford, 2009. Search in Google Scholar

[24] I. P. Gavrilyuk, W. Hackbusch and B. N. Khoromskij, Hierarchical tensor-product approximation to the inverse and related operators in high-dimensional elliptic problems, Computing 74 (2005), 131–157. 10.1007/s00607-004-0086-ySearch in Google Scholar

[25] I. V. Gavrilyuk and B. N. Khoromskij, Quantized-TT-Cayley transform to compute dynamics and spectrum of high-dimensional Hamiltonians, Comput. Methods Appl. Math. 11 (2011), no. 3, 273–290. 10.2478/cmam-2011-0015Search in Google Scholar

[26] L. Greengard and V. Rochlin, A fast algorithm for particle simulations, J. Comput. Phys. 73 (1987), 325–348. 10.1016/0021-9991(87)90140-9Search in Google Scholar

[27] W. Hackbusch and B. N. Khoromskij, Low-rank Kronecker product approximation to multi-dimensional nonlocal operators. Part I. Separable approximation of multi-variate functions, Computing 76 (2006), 177–202. 10.1007/s00607-005-0144-0Search in Google Scholar

[28] W. Hackbusch, B. N. Khoromskij, S. Sauter and E. Tyrtyshnikov, Use of tensor formats in elliptic eigenvalue problems, Numer. Linear Algebra Appl. 19 (2012), no. 1, 133–151. 10.1002/nla.793Search in Google Scholar

[29] R. J. Harrison, G. I. Fann, T. Yanai, Z. Gan and G. Beylkin, Multiresolution quantum chemistry: Basic theory and initial applications, J. Chem. Phys. 121 (2004), no. 23, 11587–11598. 10.1063/1.1791051Search in Google Scholar PubMed

[30] D. R. Hartree, The Calculation of Atomic Structure, Wiley, New York, 1957. Search in Google Scholar

[31] T. Helgaker, P. Jørgensen and J. Olsen, Molecular Electronic-Structure Theory, Wiley, New York, 1999. 10.1002/9781119019572Search in Google Scholar

[32] T. Kailath and A. Sayed, Fast Reliable Algorithms for Matrices with Structure, SIAM, Philadelphia, 1999. 10.1137/1.9781611971354Search in Google Scholar

[33] V. Khoromskaia, Black-box Hartree–Fock solver by tensor numerical methods, Comput. Methods Appl. Math. 14 (2014), no. 1, 89–111. 10.1515/cmam-2013-0023Search in Google Scholar

[34] V. Khoromskaia, D. Andrae and B. N. Khoromskij, Fast and accurate 3D tensor calculation of the Fock operator in a general basis, Comput. Phys. Commun. 183 (2012), 2392–2404. 10.1016/j.cpc.2012.06.007Search in Google Scholar

[35] V. Khoromskaia and B. N. Khoromskij, Grid-based lattice summation of electrostatic potentials by assembled rank-structured tensor approximation, Comput. Phys. Commun. 185 (2014), 3162–3174. 10.1016/j.cpc.2014.08.015Search in Google Scholar

[36] V. Khoromskaia and B. N. Khoromskij, Tensor approach to linearized Hartree–Fock equation for Lattice-type and periodic systems, preprint (2014), Search in Google Scholar

[37] V. Khoromskaia and B. N. Khoromskij, Tensor numerical methods in quantum chemistry: From Hartree–Fock to excitation energies, Phys. Chem. Chem. Phys. 17 (2015), 31491–31509. 10.1039/C5CP01215ESearch in Google Scholar PubMed

[38] V. Khoromskaia and B. N. Khoromskij, Fast tensor method for summation of long-range potentials on 3D lattices with defects, Numer. Linear Algebra Appl. 23 (2016), 249–271. 10.1002/nla.2023Search in Google Scholar

[39] B. N. Khoromskij, Structured rank-(r1,,rd) decomposition of function-related operators in Rd, Comput. Methods Appl. Math. 6 (2006), no. 2, 194–220. 10.2478/cmam-2006-0010Search in Google Scholar

[40] B. N. Khoromskij, O(dlogN)-quantics approximation of N-d tensors in high-dimensional numerical modeling, Constr. Approx. 34 (2011), no. 2, 257–289. 10.1007/s00365-011-9131-1Search in Google Scholar

[41] B. N. Khoromskij, Tensors-structured numerical methods in scientific computing: Survey on recent advances, Chemometr. Intell. Lab. Syst. 110 (2012), 1–19. 10.1016/j.chemolab.2011.09.001Search in Google Scholar

[42] B. N. Khoromskij and V. Khoromskaia, Multigrid tensor approximation of function related multi-dimensional arrays, SIAM J. Sci. Comput. 31 (2009), no. 4, 3002–3026. 10.1137/080730408Search in Google Scholar

[43] B. N. Khoromskij and S. Repin, A fast iteration method for solving elliptic problems with quasi-periodic coefficients, Russian J. Numer. Anal. Math. Modelling 30 (2015), no. 6, 329–344. 10.1515/rnam-2015-0030Search in Google Scholar

[44] B. N. Khoromskij and S. Repin, Rank structured approximation method for quasi-periodic elliptic problems, preprint (2016), 10.1515/cmam-2017-0014Search in Google Scholar

[45] T. G. Kolda and B. W. Bader, Tensor decompositions and applications, SIAM Rev. 51 (2009), no. 3, 455–500. 10.1137/07070111XSearch in Google Scholar

[46] L. Lin, C. Yang, J. C. Meza, J. Lu, L. Ying and E. Weinan, SelInv–An Algorithm for selected inversion of a sparse symmetric matrix, ACM Trans. Math. Software 37 (2011), no. 4, Aricle No. 40. 10.21236/ADA522688Search in Google Scholar

[47] S. A. Losilla, D. Sundholm and J. Juselius, The direct approach to gravitation and electrostatics method for periodic systems, J. Chem. Phys. 132 (2010), no. 2, Article ID 024102. 10.1063/1.3291027Search in Google Scholar PubMed

[48] M. Luskin, C. Ortner and B. Van Koten, Formulation and optimization of the energy-based blended quasicontinuum method, Comput. Methods Appl. Mech. Engrg. 253 (2013), 160–168. 10.1016/j.cma.2012.09.007Search in Google Scholar

[49] D. S. Mackey, N. Mackey and F. Tisseur, Structured tools for structured matrices, Electron. J. Linear Algebra 10 (2003), 106–145. 10.13001/1081-3810.1101Search in Google Scholar

[50] I. V. Oseledets, Approximation of 2d×2d matrices using tensor decomposition, SIAM J. Matrix Anal. Appl. 31 (2010), no. 4, 2130–2145. 10.1137/090757861Search in Google Scholar

[51] I. V. Oseledets and E. E. Tyrtyshnikov, Breaking the curse of dimensionality, or how to use SVD in many dimensions, SIAM J. Sci. Comput. 31 (2009), no. 5, 3744–3759. 10.1137/090748330Search in Google Scholar

[52] P. Parkkinen, S. A. Losilla, E. Solala, E. A. Toivanen, W. Xu and D. Sundholm, A generalized grid-based fast multipole method for integrating Helmholtz kernels, J. Chem. Theory Comput. 13 (2017), 10.1021/acs.jctc.6b01207. 10.1021/acs.jctc.6b01207Search in Google Scholar PubMed

[53] C. Pisani, M. Schütz, S. Casassa, D. Usvyat, L. Maschio, M. Lorenz and A. Erba, CRYSCOR: A program for the post-Hartree–Fock treatment of periodic systems, Phys. Chem. Chem. Phys. 14 (2012), 7615–7628. 10.1039/c2cp23927bSearch in Google Scholar PubMed

[54] M. V. Rakhuba and I. V. Oseledets, Calculating vibrational spectra of molecules using tensor train decomposition, J. Chem. Phys. 145 (2016), no. 12, Article ID 124101. 10.1063/1.4962420Search in Google Scholar PubMed

[55] M. V. Rakhuba and I. V. Oseledets, Grid-based electronic structure calculations: The tensor decomposition approach, J. Comput. Phys. 312 (2016), 19–30. 10.1016/ in Google Scholar

[56] Y. Saad, J. R. Chelikowsky and S. M. Shontz, Numerical methods for electronic structure calculations of materials, SIAM Rev. 52 (2010), no. 1, 3–54. 10.1137/060651653Search in Google Scholar

[57] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Ann. Phys. 51 (2011), no. 326, 96–192. 10.1016/j.aop.2010.09.012Search in Google Scholar

[58] F. Stenger, Numerical Methods Based on Sinc and Analytic Functions, Springer, New York, 1993. 10.1007/978-1-4612-2706-9Search in Google Scholar

[59] A. Szabo and N. Ostlund, Modern Quantum Chemistry, Dover Publication, New York, 1996. Search in Google Scholar

[60] H.-J. Werner and P. J. Knowles, Molpro version 2010.1, a package of Ab-Initio programs for electronic structure calculations. Search in Google Scholar

Received: 2017-1-30
Revised: 2017-3-25
Accepted: 2017-3-27
Published Online: 2017-5-31
Published in Print: 2017-7-1

© 2017 Walter de Gruyter GmbH, Berlin/Boston