Tensor Approximation of Slater-Type Orbital Basis Functions

This paper deals with a tensor representation of the Slater-type orbital basis functions. Localized basis systems are primarily used in electronic structure calculations. A choice of the system is usually limited to Gaussian-type orbitals due to the impossibility of the analytical evaluation of necessary integrals using other basis types. Unfortunately, it is not possible to use direct discretization techniques due to the dimensionality of the problem so the numerical integration is problematic. Tensor Numerical Methods overcome this problem by using special data representations, which are discussed. Finally, it is demonstrated how to use these methods to construct a tensor approximation of Slatertype orbital basis functions including an error estimate and its numerical verification.


Introduction
In electronic structure calculations we encounter several mathematical challenges. Most of the computations (Hartree-Fock, correlated, or those based on the Density Functional Theory) have the following structure. They usually start with a choice of a basis system followed by precomputing phase where integrals that are necessary for further phases are calculated. Finally, a non-linear eigenvalue problem or an optimization one is solved to obtain the ground state energy of the system. At this point the computation can be stopped or followed by other computations to include dynamical electronic correlation (e.g. post Hartree-Fock methods). At each part of the computation we can choose from a variety of mathematical tools.
In case of the electronic structure of molecules most of the standard software packages use the Gaussian-Type basis (GTO) functions in the form: with integer exponents n, l, m and a normalization constant N . A main advantage of this type of basis is the analytical integrability of the two-electron integrals where φ * are basis functions, µ, ν, κ and λ are their indices in a given basis set, x = (x 1 , x 2 , x 3 ) and y = (y 1 , y 2 , y 3 ) are three-dimensional spatial vectors. These integrals represent the basic ingredient of the precomputing phase mentioned above. In recent paper [1], a numerical evaluation of two-electron integrals and a subsequent solution of the HF equations is demonstrated using the GTO basis functions. Unfortunately, GTO basis functions behave poorly in the area of the atomic core, so for the required accuracy of the electronic structure calculation we have to enlarge the number of basis functions. In this paper we discuss a way how to represent physically more appropriate Slater-Type Orbital (STO) basis functions that can be used in a similar way as in [1]. The STO functions can be written as where N is a normalization constant and Y lm denotes a real spherical harmonics in form with a normalization constant C lm and polynomial P lm . However, these basis sets are not very popular due to the impossibility of the analytic integration. Moreover, standard grid-based discretization and numerical integration is too expensive from the perspective of the time and storage demands. In practice, the STOs are substituted by linear combination of GTOs with fixed coefficients. This representation is known as a contracted basis. Similarly to the contracted basis approximating the STOs, the Tensor Numerical Methods (TNM) use several Gaussians to approximate an STO function. The difference is in the way how these functions are generated. Computing of contracted basis coefficients requires a solution of the optimization problem of fitting in a given Slater function in the sense of least squares for each number n of representing Gaussians. These coefficients can be found in many tables (e.g. [2]), but usually with n not exceeding 6. On the other hand, the TNM approach [3] uses a simple quadrature formula, which produces less optimal coefficients, but finding an approximation for arbitrary n is very simple and straightforward with a guaranteed error convergence rate. Another advantage is that within this approach we can estimate an error of ground state energy acquired by computation with accurate STOs using extrapolation techniques.
The paper is divided into three main sections. In Sec. 2. we describe basic tensor formats that are used to represent multidimensional data discretized on very fine grids (of the order of magnitude 10 4 −10 6 discretization points in one dimension). The standard full representation of data suffers from the so called "Curse of Dimensionality", which means that the memory demands grow exponentially with the dimension of the problem. TNM overcome this problem by representing the data in an economical form where the space complexity is linear with respect to the dimension. The section also includes the Sinc approximation method [3] and [4], which is used to construct a tensor-structured representation of multidimensional data representing discretized multivariate functions. Section 3. focuses on STO functions and their tensor approximation by methods described in Sec. 2. The error of the approximation is discussed as well. Section 4. includes numerical experiments that show the accuracy and confirm advantages of the tensor representation. Finally, conclusive remarks are given in Sec. 5.

Tensor Representation of Data
In recent years, several important papers describing the basic principles of TNM have been published (e.g. [1], [5], [6] and [7]). In this section we recall some of these principles and concepts for better readability. Consider a real matrix A ∈ R m×n . The matrix can be represented by its Singular Value Decomposition (SVD) where R is the rank of matrix A, U ∈ R m×R and V ∈ R n×R are orthonormal matrices, S ∈ R R×R is a diagonal matrix which contains singular values of A.
Storage costs for the matrix A are mn whereas the SVD needs to store R · (m + n + 1) values. If the rank R is sufficiently small, the SVD representation is more economical.
Similar ideas can be applied to tensors. A tensor of order d (N -d tensor) is a multidimensional array over a d-tuple index set [1]: A so-called rank-R tensor A ∈ R n1×...×n d (for simplicity let's assume, that ∀ ∈ {1, . . . , d} : n = n) can be written as a tensor product where c k is an expansion coefficient, u ( ) k ∈ R n is -th mode vector (see Fig. 1). This representation is known as canonical R-term representation [1]. Unlike the factors of the SVD the vectors u ( ) k , k = 1 . . . R do not necessarily have to be orthogonal (for a given ). A tensor in canonical format can be stored with d·R·n numbers. In case of low rank R this quantity is much lower than n d . Generally, for an arbitrary tensor R may be high. However, in many cases R can be reduced by omitting rank-1 contributions whose expansion coefficients have absolute values lower than a specified accuracy (analogously to the truncated SVD). This way we can get a good memory-saving approximation of the original tensor. Advantages of canonical representation consist not only of low storage demands. It offers effective evaluation of multilinear algebra operations, too. For example, given two tensors with canonical ranks R 1 , R 2 , the euclidean inner product, or the Hadamard (entrywise) product takes d · R 1 · R 2 · n operations instead of n d [5].
Although using the low rank canonical representation is very effective, no universal methods to calculate it for a general arbitrary full sized tensor are known. An alternative to the canonical representation can be the Tucker format [6]. Given the rank parameter r = (r 1 , . . . , r d ) we define the orthogonal Tucker format as: where k = (k 1 , . . . , k d ) and for any ∈ {1, . . . , d} the set of vectors v ( ) k , 1 ≤ k ≤ r is orthonormal and coefficients b k1,...,k d form a full sized tensor B of order d (see Fig. 2). Unlike the canonical representation there exist robust algorithms for calculating a good Tucker approximation of an arbitrary full sized tensor. However, for some types of function related tensors effective algorithms for finding a canonical approximation are known, e.g. the Sinc approximation.

Sinc Approximation
The Sinc approximation is a robust method suitable for finding a good canonical approximation of functionrelated tensors. Consider a multivariate function u : R d → R that is spherically symmetrical, i.e. it depends on x where x ∈ R d . Assume that u can be written as an integral: with Ω ∈ R, R + , I , where I is closed real interval and G, F : R → R. By applying a suitable quadrature we obtain a separable approximation: with quadrature weights w k . To get a discrete canonical approximation of u on a domain [−a, a] d we just have to represent each one-dimensional discretized interval [−a, a] as a vector y ∈ R n . Then the canonically represented tensor approximating discretized function u can be written as: where and The Sinc-quadrature is based on the interpolation by the Sinc function A Sinc interpolant is defined as: for functions f from the Hardy space ∈ H 1 (D δ ) [6]. The Hardy space H 1 (D δ ) is defined as a set of functions that are holomorphic on the strip with a given δ ∈ 0, π 2 and For f ∈ H 1 (D δ ) the integral can be well approximated by an integral of Sinc interpolant, i.e.
For practical use the sum has to be truncated, i.e.
The estimate of error of such approximation is given by the following theorem [4].
then there exists a constant C 1 depending on f, δ and α such that

Canonical Representation of STO Functions
Let's have a look at the approximation of STO functions. For simplicity the STO function Eq. (3) can be rewritten as: whereỸ lm does not contain the denominator Our goal is to find satisfactory approximations of the exponential function e −α √ x 2 1 +x 2 2 +x 2 3 , the radial function and the modified spherical harmonicsỸ lm (x 1 , x 2 , x 3 ).

Approximation of the Exponential Function
For the exponential function e −α x there exists an integral representation (based on the Laplace transform), Let us denote g(t) the restriction of the integrand to x = 0, i.e.
For this function, it is not straightforward to find a suitable quadrature. The problem is that behavior of g(t) strongly depends on parameter α. Especially for larger values of α the function converges very slowly to 0 as t goes to infinity (see Fig. 3) so the choice of α-independent quadrature is not possible. This problem can be solved by substitution so we can rewrite integral Eq. (26) as:

du. (29)
This representation is more appropriate than Eq. (26). First, if we look at the restriction of integrand to x = 0, i.e.
we can see that it does not depend on α. Moreover, it can be easily shown that g(u) reaches its maximum at u = 0. These properties facilitate a good choice of quadrature points.
To make an error estimation we use Thm. 1. First, we have to show that integrand of Eq. (29) as a function of u belongs to the Hardy space H 1 (D δ ) independently of x and α. Let's define complex function (now, we understand x and α as parameters). Function f is evidently holomorphic on C, so it is holomorphic on D δ for arbitrary δ. Now we will check condition Eq. (18). Let's find an upper bound for integral Eq. (18). Function f satisfies: For the right hand side of Eq. (32) we can easily find an upper bound By choosing δ = π 3 we can estimate integral Eq. (18) by: We have proven that f ∈ H 1 D π 3 , so now, using Thm. 1, we can make an error estimate of the Sinc approximation. It can be easily shown that If we use the quadrature step the error of the Sinc-quadrature can be estimated as: where C 1 > 0 is a real constant that does not depend on x or α. This estimate will be checked against numerical experiments in Sec. 4.

Approximation of the Radial Function
The radial function can be written as: where a = n − l − 1. First, assume that a = 2. Then, the radial function can be written as: Finding a canonical representation of discretization of r 2 is straightforward. It can be simply written as: where i ( ) is vector of ones and z ( ) i = y i 2 . The canonical rank of the tensor R 2 is equal to 3. Now for an arbitrary even exponent a = 2k, k ∈ N, we can simply find a canonical representation by computing the k-th power of R 2 (in the Hadamard sense), i.e.
so the canonical rank is equal to 3 k . In case of odd exponent a = 2k − 1, the situation is a little more complicated. We have to represent our function as: For the first term of the r.h.s. of Eq. (42) (so called Newton potential) 1 x we can find a good approximation using the Sinc-quadrature applied to the following formula [6]: Similarly to the exponential function, there can be shown an exponential convergence of the error. A resulting canonical rank of R 2k−1 is equal to 3 k R N , where R N is the canonical rank of the Newton potential approximation.

Canonical Representation of Modified Spherical Harmonics
The canonical representation ofỸ lm is straightforward. It can be made analogously to the radial function with α = 2. If we look at the table of spherical harmonics [8], we can observe that their canonical rank does not exceed 6 for l ≤ 4 (i.e. including s,p,d,f,g orbitals, without g orbitals the canonical rank does not exceed 3).

Reduction of Rank
Although the canonical rank of each part is quite low, the total canonical rank is equal to the product of those ranks. Too high ranks can negatively affect the computational effectiveness. There are two ways how to overcome this problem. One, physically less correct way is to simplify the STO functions. Instead of radial function and spherical harmonics we can consider (similarly to GTOs) powers of spatial coordinates only, i.e.
In this case, the total canonical rank is equal to the rank of the exponential function approximation. Unfortunately, in this form we have to consider more basis functions due the lower physical correctness. Another option is to use the original Slater functions, but expressed in terms of a high rank canonical representation and subsequently simplified via a rank-reduction algorithm. One of such algorithms is called The Best Tucker Approximation (C_BTA). Roughly described, the algorithm consists of the following steps: • conversion of canonically-represented tensor to the Tucker format with the core tensor represented in the canonical format, • optimization of the Tucker approximation using the least square algorithm, • return to the canonical representation with the reduced rank.
For more details about algorithm see e.g. [7].

Numerical Results
In this section, we will show the results of numerical experiments. The error estimate of the Sinc approximation of the exponential function is verified. An independence of the error on α is also discussed. Finally, an error dependence of the tensor approximation of a selected STO function on the canonical rank is shown.

Error of the Sinc Approximation of the Exponential Function
First, we are going to verify the error estimate given by Eq. (37). The quadrature step h is generated by formula Eq. (36). Figure 4 shows the dependence of the error on the canonical rank of the approximating tensor. The error is measured as the maximum of the difference between the analytical formula Eq. (25) with α = 1, n = 1, l = m = 0 and the Sinc approximation. Furthermore, it shows the theoretical error estimate given by formula Eq. (37) where C 1 has been adjusted to 1.2.
We can see that the error estimate is correct. Table 1 shows experimentally determined values of C 1 for given values of α. We can see that the estimate C 1 = 1.2 is independent on α which corresponds to the theory. Although we have a quadrature step with a theoretical (and numerically tested) error estimate, we use an alternative way for generating the quadrature because of better convergence properties. It is based on the observation that function g (u) given by formula Eq. (30) rapidly and monotonously converges to 0 as u goes to ±∞. It means that an integration over a fixed interval [u min , u max ] well approximates the integral over R. For example by setting u min = −3.1938 and u max = 26.0010 we get umax umin g (u) du So our strategy is to distribute the quadrature points equidistantly over the interval [u min , u max ]. Figure 5 shows the dependence of the error on the canonical rank of the approximating tensor compared to the theoretical error estimate related to the original quadrature.
We can see that for a small values of the canonical rank the error is larger than in case of the original quadrature. However, the decrease of the error is significantly faster. Furthermore, we can notice that the error decreases up to the canonical rank value 40-50, then it stagnates. This is not surprising, because the integration over the finite interval has a limited accuracy. The essential fact is that using the original quadrature an adequate order of error can be achieved at canonical rank values greater than 100.
Using the quadrature step given by formula Eq. (36) the error estimate does not depend on α. Unfortunately, for an alternative quadrature we have no proof of independence on α. However, numerical tests showed that in this case the error is practically independent on α.

Error of the Sinc Approximation of the STO
To approximate the STO function we have to combine the Sinc approximation of the exponential function with the Sinc approximation of the radial function and with the spherical harmonics. The combination is realized via the Hadamard product so the resulting canonical rank is equal to the product of ranks of individual contributions. The test has been performed on STO function with parameters α = 0.1, n = 1, Canonical ranks of both Sinc approximations were generated from a parameter r ∈ N as for preserving a balance of errors of both contributions. The absolute error dependence of the tensor approximation on the parameter r is shown on Fig. 6. The error is measured in the same way as in Subsec. 4.1.
Similarly to the exponential function we can see the exponential convergence rate. It should be noted that the total canonical rank is equal to 3r 2 , which might be too high for practical use. However, the rank can be reduced by discarding rank-1 updates with coefficients c k lower than some accuracy limit and/or by using the reduction algorithm mentioned in Subsec. 3.4.

Conclusion
In this paper, we have discussed a tensor representation of STO functions using the Sinc approximation method. In case of the exponential function the dependence of the error on the canonical rank has been shown. An exponential convergence rate of the error provides a good usability of the approximation of the exponential function at low values of the canonical rank. Especially, it can be used for the canonical representation of the STO basis functions. These basis functions are included in a currently developed Hartree-Fock module. Currently we are working on extensions of our module. A tensor approximation of STO functions can be understood as a first step towards more complicated bases such as hydrogen-like orbitals or even non-localized system independent tensor product bases constructed from 1D basis functions (e.g. polynomials, trigonometric functions or wavelets). Operations performed on tensor-structured data require using of high performance computing infrastructures. Therefore, they need to be effectively parallelized. The ultimate goal of the present research is to extend the use of tensor methods beyond the limits of the HF approximation to the realm of correlated and DFT calculations.