In computational schemes including bilinear tensor-tensor or matrix-tensor operations the increase of tensor ranks leads to the critical loss of efficiency. Moreover, in many applications, for example in electronic structure calculations, the canonical tensors with large rank parameters arise as the result of polynomial type or convolution transforms of some function related tensors (say, electron density, the Hartree potential, etc.) In what follows, we present the new look on the direct method of rank reduction for the canonical tensors with large initial rank, the reduced HOSVD, first introduced and analyzed in [46].

In what follows, we consider the vector space of d-fold real-valued data arrays ℝn1×⋯×nd endorsed by the Euclidean scalar product 〈·, ·〉 with the related norm ||u|| = 〈u, u〉^1/2. We denote by Tr,n the class of tensors A ∈ ℝ^n⊗d parametrized in the rank-r, r = (r₁, …, r_d) orthogonal Tucker format,

with the orthogonal side-matrices V(ℓ)=[v1(ℓ)…vrℓ(ℓ)]∈ℝn×rℓ and with the core coefficient tensor β∈ℝr1×...×rd. Here and thereafter × _ℓ denotes the contracted tensor-matrix product in the dimension ℓ, and ℝ^n⊗d denotes the Euclidean vector space of n₁ × ⋯ × n_d-tensors with equal mode size n_ℓ = n, ℓ = 1, …, d.

Likewise, CR,n denotes the class of rank-R canonical tensors. For given A∈CR,n in the rank-R canonical format,

with normalized canonical vectors, i.e., ||uν(ℓ)||=1 for ℓ = 1, …, d, ν = 1, …, R.

The standard algorithm for the Tucker tensor decomposition [7] is based on HOSVD applied to full tensors of size n^d which exhibits O(n^d+1) computational complexity. The question is how to simplify the HOSVD Tucker approximation in the case of canonical input tensor in the form Equation (1) without use of the full format representation of A, and in the situation when the CP rank parameter R and the mode sizes n of the input can be sufficiently large.

First, let us use the equivalent (nonorthogonal) rank-r = (R, …, R) Tucker representation of the tensor Equation (1),

via contraction of the diagonal tensor ξ=diag{ξ1,…,ξR}∈ℝR⊗d with ℓ-mode side matrices U(ℓ)=[u1(ℓ),…,uR(ℓ)]∈ℝn×R (see Figure 1). By definition the tensor diag{ξ1,…,ξR}∈ℝR⊗d is called diagonal if it has all zero entries except the diagonal elements given by ξ(i_ℓ, …, i_ℓ) = ξ_iℓ, ℓ = 1, …, R.Then the problem of canonical to Tucker approximation can be solved by the method of reduced HOSVD (RHOSVD) introduced in [46]. The basic idea of the reduced HOSVD is that for large (function related) tensors given in the canonical format their HOSVD does not require the construction of a tensor in the full format and SVD based computation of its matrix unfolding. Instead, it is sufficient to compute the SVD of the directional matrices U^(ℓ) in Equation (2) composed by only the vectors of the canonical tensor in every dimension separately, as shown in Figure 1. This will provide the initial guess for the Tucker orthogonal basis in the given dimension. For the practical applicability, the results of the approximation theory on the low-rank approximation to the multi-variate functions, exhibiting exponential error decay in the Tucker rank, are of the principal significance [51].

In the following, we suppose that n ≤ R and denote the SVD of the side-matrix U^(ℓ) by

with the orthogonal matrices Z(ℓ)=[z1(ℓ),…,zn(ℓ)]∈ℝn×n, and V(ℓ)=[v1(ℓ),…,vn(ℓ)]∈ℝR×n, ℓ = 1, …, d. We use the following notations for the vector entries, vk(ℓ)(ν)=vk,ν(ℓ) (ν = 1, …, R).

To fix the idea, we introduce the vector of rank parameters, r = (r₁, …, r_d), and let

be the rank-r_ℓ truncated SVD of the side-matrix U^(ℓ) (ℓ = 1, …, d). Here, the matrix D_{ℓ, 0} = diag{σ_{ℓ, 1}, σ_{ℓ, 2}, …, σ_ℓ,rℓ} is the submatrix of D_ℓ in Equation (3) and

represent the respective dominating (n × r_ℓ)-submatrices of the left and right factors in the complete SVD decomposition in Equation (3).

Definition 2.1. (Reduced HOSVD, [46]). Given the canonical tensor A∈CR,n, the truncation rank parameter r, (r_ℓ ≤ R), and rank-r_ℓ truncated SVD of U^(ℓ), see Equation (4), then the RHOSVD approximation of A is defined by the rank-r orthogonal Tucker tensor

obtained by the projection of canonical side matrices U^(ℓ) onto the left orthogonal singular matrices Z0(ℓ), defined in Equation (4).

Notice that the general error bound for the RHOSVD approximation will be presented by Theorem 2.3, see also the discussion afterwards. Corollary 2.4 provides the conditions which guarantee the stability of RHOSVD.

The sub-optimal Tucker approximand Equation (5) is simple to compute and it provides accurate approximation to the initial canonical tensor even with rather small Tucker rank. Moreover, this provides the good initial guess to calculate the best rank-r Tucker approximation by using the ALS iteration. In our numerical practice, usually, only one or two ALS iterations are required for convergence. For example, in case d = 3, algorithmically, the one step of the canonical-to-Tucker ALS algorithm reduces to the following operations. Substituting the orthogonal matrices Z0(1) and Z0(3) from Equation (5) into Equation (2), we perform the initial step of the first ALS iteration

where A₂ is given by the contraction

as illustrated in Figure 2. Then we optimize the orthogonal subspace in the second variable by calculating the best rank-r₂ approximation to the r₁r₃ × n₂ matrix unfolding of the tensor A₂. The similar contracted product representation can be used when d > 3, as well as for the construction of the TT representation for the canonical input.

Here, we notice that the core tensor in the RHOSVD decomposition can be represented in the CP data-sparse format.

Proposition 2.2. The core tensor

in the orthogonal Tucker representation Equation (5), A(r)0=β0×1Z0(1)×2⋯×dZ0(d)∈Tr,n, can be recognized as the rank-R canonical tensor of size r₁ × ⋯ × r_d with the storage request R(∑ℓ=1drℓ), which can be calculated entry-wise in O(Rr₁ ⋯ r_d) operations.

Indeed, introducing the matrices U0(ℓ)T=Dℓ,0V0(ℓ)T∈ℝrℓ×R, for ℓ = 1, …d, we conclude that the canonical core tensor β₀ is determined by the ℓ-mode side matrices U0(ℓ)T. In the other words, the tensor A(r)0 is represented in the mixed Tucker-canonical format getting rid of the “curse of dimensionality” (see also Section 2.2 below).

The accuracy of the RHOSVD approximation can be controlled by the given ε-threshold in truncated SVD of side matrices U^(ℓ). The following theorem proves the absolute error bound for the RHOSVD approximation.

Theorem 2.3. (RHOSVD error bound, [46]). For given A∈CR,n in Equation (1), let σ_ℓ,1 ≥ σ_ℓ,2 … ≥ σ_{ℓ,min(n, R)} be the singular values of ℓ-mode side matrices U^(ℓ) ∈ ℝ^n×R (ℓ = 1, …, d) with normalized skeleton vectors. Then the error of RHOSVD approximation, A(r)0, is bounded by

The complete proof can be found in Section 8 (see Appendix).

The accuracy of the RHOSVD can be controlled in terms of the ε-criteria. To that end, given ε > 0, chose the Tucker ranks such that ∑ℓ=1d(∑k=rℓ+1min(n,R)σℓ,k2)1/2≤ε is satisfied, then Theorem 2.3 provided the error bound adapted to the ε-threshold.

The error estimate in Theorem 2.3 differs from the case of complete HOSVD by the extra factor ||ξ||, which is the payoff for the lack of orthogonality in the canonical input tensor. Hence, Theorem 2.3 does not provide, in general, the stable control of relative error since for the general canonical tensors there is no uniform upper bound on the constant C in the estimate

The problem is that Equation (8) applies to the general non-orthogonal canonical decomposition.

The stable RHOSVD approximation can be proven in the case of the so-called partially orthogonal or monotone decompositions. With partially orthogonal decomposition we mean that for each pair of indexes ν, μ in Equation (1) there holds prodℓ=1d〈uν(ℓ),uμ(ℓ)〉=0. For monotone decompositions we assume that all coefficients and skeleton vectors in Equation (1) have non-negative values.

Corollary 2.4. (Stability of RHOSVD) Assume the conditions of Theorem 2.3 are satisfied. (A) Suppose that at least one of the side matrices U^(ℓ), ℓ = 1, ⋯ , d, in Equation (2), is orthogonal or the decomposition Equation (1) is partially orthogonal. Then the RHOSVD error can be bounded by

(B) Let decomposition Equation (1) be monotone. Then (9) holds.

Proof. (A) The partial orthogonality assumption combined with normalization constraints for the canonical skeleton vectors imply

The above relation also holds in the case of orthogonality of the side matrix U^(ℓ) for some fixed ℓ. Then the result follows by (7).

(B) In case of monotone decomposition we conclude that the pairwise scalar product of all summands in Equation (1) is non-negative, while the norm of each ν-term is equal to ξ_ν. Then the upper bound

holds for vectors uν=ξνuν(1)⊗…⊗uν(d), ν = 1, ⋯ , R, with non-negative entries applied to the case of R summands, thus implying ||ξ||² ≤ ||A||². Now, the result follows.

Clearly, the orthogonality assumption may lead to slightly higher separation rank, however, this constructive decomposition stabilizes the RHOSVD approximation method applied to the canonical format tensor (i.e., it allows the stable control of relative error). The case of monotone canonical sums typically arises in the sinc-based canonical approximation to radially symmetric Green's kernels by a sum of Gaussians. On the other hand, in long term computational practice the numerical instability of RHOSVD approximation was not observed in case of physically relevant data.

In the procedure for the canonical tensor rank reduction the goal is to have a result in a canonical tensor format with a smaller rank. By converting the core tensor to CP format, one can use the mixed two-level Tucker data format [12, 27], or canonical CP format. Figure 3 illustrates the computational scheme of the two-level Tucker approximation.

We define by n¯ℓ the single-hole product of dimension-modes,

The same definition applies to the quantity r¯ℓ.

Next lemma describes the approximation of the Tucker tensor by using canonical representation [12, 27].

Lemma 2.5. (Mixed Tucker-to-canonical approximation, [27]).

(A) Let the target tensor A have the form A=β×1V(1)×2…×dV(d)∈Tr,n, with the orthogonal side-matrices V(ℓ)=[v1(ℓ)…vrℓ(ℓ)]∈ℝn×rℓ and β∈ℝr1×...×rd. Then, for a given R≤min1≤ℓ≤dr¯ℓ,

(B) Assume that there exists the best rank-R approximation A(R)∈CR,n of A, then there is the best rank-R approximation β(R)∈CR,r of β, such that

The complete proof can be found in Section 8 (see Appendix). Notice that condition R≤min1≤ℓ≤dr¯ℓ simply means that the canonical rank does not exceed the maximal CP rank of the Tucker core tensor.

Combination of Theorem 2.3 and Lemma 2.5 paves the way to the rank optimization of canonical tensors with the large mode-size arising, for example, in the grid-based numerical methods for multi-dimensional PDEs with non-regular (singular) solutions. In such applications the univariate grid-size (i.e., the mode-size) may be about n = 10⁴ and even larger.

Notice that the Tucker (for moderate d) and canonical formats allow to perform basic multi-linear algebra using one-dimensional operations, thus reducing the exponential scaling in d. Rank-truncated transforms between different formats can be applied in multi-linear algebra on mixed tensor representations as well, see Lemma 2.5. The particular application to tensor convolution in many dimensions was discussed, for example, in [1, 2].

We summarize that the direct methods of tensor approximation can be classified by:

Direct analytic approximation methods by sinc quadrature/interpolation are of principal importance. Basic examples are given by the tensor representation of Green's kernels, the elliptic operator inverse and analytic matrix-valued functions. In all cases, the algebraic methods for rank reduction by the ALS-type iterative Tucker/canonical approximation can be applied.

Further improvement and enhancement of algebraic tensor approximation methods can be based on the combination of advanced nonlinear iteration, multigrid tensor methods, greedy algorithms, hybrid tensor representations, and the use of new problem adapted tensor formats.

In the rank reduction scheme for the canonical rank-R tensors, we use successively the canonical-to-Tucker (C2T) transform and then the Tucker-to-canonical (T2C) tensor approximation.

First, we notice that the canonical rank of a tensor A ∈ V_n has the upper bound (see [27, 46]),

where nℓ¯ is given by Equation (10). Rank bound (13) applied to the Tucker core tensor of the size r × r × r, indicates that the ultimate canonical rank of a large-size tensor in V_n has the upper bound r². Notice that for function related tensors the Tucker rank scales logarithmically in both approximation accuracy and the discretization grid size (see the proof for some classes of function in [51]).

The following remark shows that the maximal canonical rank of the Tucker core of 3rd order tensor can be easily reduced to the value less than r² by the SVD-based procedure applied to the matrix slices of the Tucker core tensor β. Though, being not practically attractive for arbitrary high order tensors, the simple algorithm described in Remark 2.6 below is proved to be useful for the treatment of small size 3rd order Tucker core tensors within the rank reduction algorithms described in the previous sections.

Remark 2.6. Let d = 3 for the sake of clarity [27, 46]. There is a simple procedure based on SVD to reduce the canonical rank of the core tensor β, within the accuracy ε > 0. Denote by Bm∈ℝr×r, m = 1, …, r the two-dimensional slices of β in each fixed mode and represent

where z_m(m) = 1, z_m(j) = 0 for j = 1, …, r, j ≠ m (there are exactly d possible decompositions). Let p_m be the minimal integer, such that the singular values of B_m satisfy σk(m)≤εr3/2 for k = p_m + 1, …, r (if σr(m)>εr3/2, then set p_m = r). Then, denoting by

the corresponding rank-p_m approximation to B_m (by truncation of σpm+1(m),…,σr(m)), we arrive at the rank-R canonical approximation to β,

providing the error estimate

Representation (15) is a sum of rank-p_m terms so that the total rank is bounded by R≤p1+...+pr≤r2. The approach can be extended to arbitrary d ≥ 3 with the bound R ≤ r^d−1.

Figure 4 illustrates the canonical decomposition of the core tensor by using the SVD of slices B_m of the core tensor β, yielding matrices Um={ukm}k=1pm, Vm={vkm}k=1pm and a diagonal matrix of small size p_m × p_m containing the truncated singular values. It also shows the vector z_m = [0, …, 0, 1, 0, …, 0], containing all entries equal to 0 except 1 at the mth position.

It is worse to note that the rank reduction for the rank-R core tensor of small size r₁ × ⋯ × r_d, can be also performed by using the cascading ALS algorithms in CP format applied to the canonical input tensor, as it was applied in [50]. Moreover, a number of numerical examples presented in the present paper and in the included literature (applied to function generated tensors) demonstrate the substantial reduction of the initial canonical rank R.

First, we recall the grid-based method for the low-rank canonical representation of a spherically symmetric kernel functions p(||x||), x ∈ ℝ^d for d = 2, 3, …, by its projection onto the finite set of basis functions defined on tensor grid. The approximation theory by a sum of Gaussians for the class of analytic potentials p(||x||) was presented in [1, 11, 51, 52]. The particular numerical schemes for rank-structured representation of the Newton and Slater kernels

discretized on a fine 3D Cartesian grid in the form of low-rank canonical tensor was described in [11, 51].

In what follows, for the ease of exposition, we confine ourselves to the case d = 3. In the computational domain Ω = [−b, b]³, let us introduce the uniform n × n × n rectangular Cartesian grid Ω_n with mesh size h = 2b/n (n even). Let {ψi=∏ℓ=13ψiℓ(ℓ)(xℓ)} be a set of tensor-product piecewise constant basis functions, labeled by the 3-tuple index i = (i₁, i₂, i₃), i_ℓ ∈ I_ℓ = {1, …, n}, ℓ = 1, 2, 3. The generating kernel p(||x||) is discretized by its projection onto the basis set {ψ_i} in the form of a third order tensor of size n × n × n, defined entry-wise as

The low-rank canonical decomposition of the 3rd order tensor P is based on using exponentially fast convergent sinc-quadratures for approximating the Laplace-Gauss transform to the analytic function p(z), z ∈ ℂ, specified by a certain weight p^(t)>0,

with the proper choice of the quadrature points t_k and weights p_k. The sinc-quadrature based approximation to generating function by using the short-term Gaussian sums in Equation (23) are applicable to the class of analytic functions in certain strip |z| ≤ D in the complex plane, such that on the real axis these functions decay polynomially or exponentially. We refer to basic results in [11, 52, 53], where the exponential convergence of the sinc-approximation in the number of terms (i.e., the canonical rank) was analyzed for certain classes of analytic integrands.

Now, for any fixed x=(x1,x2,x3)∈ℝ3, such that ||x|| > a > 0, we apply the sinc-quadrature approximation Equation (23) to obtain the separable expansion

providing an exponential convergence rate in M,

In the case of Newton kernel, we have p(z) = 1/z, p^(t)=2π, so that the Laplace-Gauss transform representation reads

which can be approximated by the sinc quadrature Equation (24) with the particular choice of quadrature points t_k, providing the exponential convergence rate as in Equation (25) [11, 51].

In the case of Yukawa potential the Laplace Gauss transform reads

The analysis of the sinc quadrature approximation error for this case can be found, in particular, in [1, 51], section 2.4.7.

Combining (22) and (24), and taking into account the separability of the Gaussian basis functions, we arrive at the low-rank approximation to each entry of the tensor P = [p_i],

Define the vector (recall that p_k > 0)

then the 3rd order tensor P can be approximated by the R-term (R = 2M + 1) canonical representation

Given a threshold ε > 0, in view of Equation (25), we can choose M = O(log²ε) such that in the max-norm

In the case of continuous radial function p(||x||), say the Slater potential, we use the collocation type discretization at the grid points including the origin, x = 0, so that the univariate mode size becomes n → n₁ = n + 1. In what follows, we use the same notation P_R in the case of collocation type tensors (for example, the Slater potential) so that the particular meaning becomes clear from the context.

The range separated (RS) tensor format was introduced in [32] for efficient representation of the collective free-space electrostatic potential of large biomolecules. This rank-structured tensor representation of the collective electrostatic potential of many-particle systems of general type allows to reduce essentially computation of their interaction energy, and it provides convenient form for performing other algebraic transforms. The RS format proved to be useful for range-separated tensor representation of the Dirac delta [34] in ℝ^d and based on that, for regularization of the Poisson-Boltzmann equation (PBE) by decomposition of the solution into short- and long-range parts, where the short-range part of the solution is evaluated by simple tensor operations without solving the PDE. The smooth long-range part is calculated by solving the PBE with the modified right-hand side by using the RS decomposition of the Dirac delta, so that now it does not contain singularities. We refer to papers [33, 35] describing the approach in details.

First, we recall the definition of the range separated (RS) tensor format, see [32], for representation of d-tensors A∈ℝn1×⋯×nd. The RS format is served for the hybrid tensor approximation of discretized functions with multiple cusps or singularities. This allows the splitting of the target tensor onto the highly localized components approximating the singularity and the component with global support that allows the low-rank tensor approximation. Such functions typically arise in computational quantum chemistry, in many-particle modeling and in the interpolation of multi-dimensional data measured at certain set of spatial points in ℝ^{n × n × n}.

In the following definition of RS-canonical tensor format, we use the notion of localized canonical tensor U₀, which is characterized by the small support whose diameter has a size of a few grid points. This tensor will be used as the reference one for presentation of the short-range part in the RS tensor. To that end we use the operation Replica_xν(U₀) which replicates U₀ into some given grid point x_ν. In this construction, we assume that the chosen grid points x_ν are well separated, i.e., the distance between each pair of points is not less then some given threshold n_δ > 0.

Definition 4.1. (RS-canonical tensors, [32]). Given the rank-R_s reference localized CP tensor U₀. The RS-canonical tensor format defines the class of d-tensors A∈ℝn1×⋯×nd, represented as a sum of a rank-R_l CP tensor Ulong=∑k=1Rlξkuk(1)⊗⋯⊗uk(d), and a cumulated CP tensor Ushort=∑ν=1N0cνUν, such that

where U_short is generated by the localized reference CP tensor U₀, i.e., U_ν = Replica_xν(U₀), with rank(U_ν) = rank(U₀) ≤ R_s, where, given the threshold n_δ > 0, the effective support of U_ν is bounded by diam(suppU_ν) ≤ 2n_δ in the index size.

Each RS-canonical tensor is, therefore, uniquely defined by the following parametrization: rank-R_l canonical tensor U_long, the rank-R_s reference canonical tensor U₀ with the small mode size bounded by 2n_δ, list J of the coordinates and weights of N₀ particles in ℝ^d. The storage size is linear in both the dimension and the univariate grid size,

The main benefit of the RS-canonical tensor decomposition is the almost uniform bound on the CP/Tucker rank of the long-range part Ulong=∑k=1Rlξkuk(1)⊗⋯⊗uk(d), in the multi-particle potential discretized on fine n × n × n spatial grid. It was proven in [32] that the canonical rank R scales logarithmically in both the number of particles N₀ and the approximation precision, see also Lemma 4.5.

Given the rank-R CP decomposition Equation (29) based on the sinc-quadrature approximation Equation (24) of the discretized radial function p(||x||), we define the two subsets of indices, K_l: = {k:t_k ≤ 1} and K_s: = {k:t_k > 1}, and then introduce the RS-representation of this tensor as follows,

where

This representation allows to reduce the calculation of the multi-particle interaction energy of the many-particle system. Recall that the electrostatic interaction energy of N charged particles is represented in the form

and it can be computed by direct summation in O(N²) operations. The following statement is the modification of Lemma 4.2 in [32] (see [54] for more details).

Lemma 4.2. [54] Let the effective support of the short-range components in the reference potential P_R for the Newton kernel does not exceed the minimal distance between particles, σ > 0. Then the interaction energy E_N of the N-particle system can be calculated by using only the long range part in the tensor P representing on the grid the total potential sum,

in O(R_lN) operations, where R_l is the canonical rank of the long-range component in P, P_l.

Here, z ∈ ℝ^N is a vector composed of all charges of the multi-particle systems, and pl∈ℝN is the vector of samples of the collective electrostatic long-range potential P_l in the nodes corresponding to particle locations. Thus, the term 12〈z,pl〉 denotes the “non–calibrated” interaction energy associated with the long-range tensor component P_l, while P_Rl denotes the long-range part in the tensor representing the single reference Newton kernel, and P_Rl(0) is its value at the origin.

Lemma 4.2 indicates that the interaction energy does not depend on the short-range part in the collective potential, and this is the key point for the construction of energy preserving regularized numerical schemes for solving the basic equations in bio-molecular modeling by using low-rank tensor decompositions.

In what follows, we address the important question on how the basic multi-linear operations can be implemented in the RS tensor format by using the RHOSVD rank compression. The point is that various tensor operations arise in the course of commonly used numerical schemes and iterative algorithms which usually include many sums and products of functions as well as the actions of differential/integral operators, always making the tensor structure of input data much more complicated requiring the robust rank reduction schemes.

The other important aspect is related to the use of large (fine resolution) discretization grids which is limited by the restriction on the size of the full input tensors, O(n^d) (curse of dimensionality), representing the discrete functions and operators to be approximated in low rank tensor format. Remarkably, that tensor decomposition for special class of functions, which allow the sinc-quadrature approximation, can be performed on practically indefinitely large grids because the storage and numerical costs of such numerical schemes scale linearly in the univariate grid size, O(dn). Hence, having constructed such low rank approximations for certain set of “reproducing” radial functions, makes it possible to construct the low rank RS representation at linear complexity, O(dn), for the wide class of functions and operators by using the rank truncated multi-linear operations. The examples of such “reproducing” radial functions are commonly used in our computational practice.

First, consider the Hadamard product of two tensors P_R and Q_R1 corresponding to the pointwise product of two generating multi-variate functions centered at the same point. The RS representation of the product tensor is based on the observation that the long-range part of the Hadamard product of two tensors in RS-format is basically determined by the product of their long-range parts.

Lemma 4.3. Suppose that the RS representation Equation (31) of tensors P_R and Q_R1 is constructed based on the sinc-quadrature CP approximation Equation (29). Then the long-range part of the Hadamard product of these RS-tensors,

can be represented by the product of their long-range parts, Z_l = P_l ⊙ Q_l, with the subsequent rank reduction. Moreover, we have rank(Z_l) ≤ R_lQ_l.

Proof. We consider the case of collocation tensors and suppose that each skeleton vector in CP tensors P_R and Q_R1 is given by the restriction of certain Gaussians to the set of grid points. Chose the arbitrary short-range components in P_R and some component in Q_R1, generated by Gaussians e-tkxℓ2 and e-tmxℓ2, respectively. Then the effective support of the product of these two terms becomes smaller than that for each of the factors in view of the identity e-tkxℓ2e-tmxℓ2=e-(tk+tm)xℓ2 considered for arbitrary t_k, t_m > 0. This means that each term that includes the short-range multiple remains to be in the short range. Then the long range part in Z takes a form Z_l = P_l ⊙ Q_l with the subsequent rank reduction.

The sums of several tensors in RS format can be easily split into short- and long-range parts by grouping the respective components in the summands. The other important operation is the operator-function product in RS tensor format (see the example in [34] related to the action of Laplacian with the singular Newton kernel resulting in the RS decomposition of the Dirac delta). This topic will be considered in detail elsewhere.

In what follows, we consider the RS-canonical tensor format for the rank-structured representation of the Slater function

which has the principal significance in electronic structure calculations (say, based on the Hartree-Fock equation) since it represents the cusp behavior of electron density in the local vicinity of nuclei. This function (or its approximation) is considered as the best candidate to be used as the localized basis function for atomic orbitals basis sets. Another direction is related to the construction of the accurate low-rank global interpolant for big scattered data to be considered in the next section. In this way, we calculate the data adaptive basis set living on the fine Cartesian grid in the region of target data. The main challenge, however, is due to the presence of point singularities which are hard to approximate in the problem independent polynomial or trigonometric basis sets.

The construction of low-rank RS approximation to the Slater function is based on the generalized Laplace transform representation for the Slater function written in the form G(ρ)=e-2αρ, ρ(x)=x12+...+xd2, reads

which corresponds to the choice G^(τ)=απτ-3/2e-α/τ in the canonical form of the Laplace transform representation for G(ρ),

Denote by G_R the rank-R canonical approximation to the function G(ρ) discretized on the n × n × n Cartesian grid.

Lemma 4.4. ([51]) For given threshold ε > 0 let ρ ∈ [1, A]. Then the (2M + 1)-term sinc-quadrature approximation of the integral in (34) with

ensures the max-error of the order of O(ε) for the corresponding rank-(2M + 1) CP approximation G_R to the tensor G.

Figure 5 illustrates the RS splitting for the tensor G_R = G_Rl + G_Rs representing the Slater potential G(x) = e^−λ||x||, λ = 1, discretized on the n × n × n grid with n = 1024. The rank parameters are chosen by R = 24, R_l = 6 and R_s = 18. Notice that for this radial function the long-range part (Figure 5, left) includes much less canonical vectors comparing with the case of Newton kernel. This anticipates the smaller total canonical rank for the long-range part in the large sum of Slater-like potentials arising, for example, in the representation of molecular orbitals and the electron density in electronic structure calculations. For instance, the wave function for the Hydrogen atom is given by the Slater function e^−μ||x||. In the following section, we consider the application of RS tensor format to interpolation of scattered data in ℝ^d.

In scattered data modeling the problem is in a low parametric approximation of multi-variate functions f:ℝ^d → ℝ by sampling at a finite set X={x1,…,xN}⊂ℝd of piecewise distinct points. Here, the function f might be the surface of a solid body, the solution of a PDE, many-body potential field, multi-parametric characteristics of physical systems, or some other multi-dimensional data, etc.

Traditional ways of recovering f from a sampling vector f_|X = (f(x₁), …, f(x_N)) is the constructing a functional interpolant PN:ℝd→ℝ such that PN|X=f|X=:f∈ℝN, i.e.,

Using radial basis (RB) functions one can find interpolants P_N in the form

where p = p(r):[0, ∞) → ℝ is a fixed RB function, and r = ||·|| is the Euclidean norm on ℝ^d. In further discussion, we set Q(x) = 0. For example, the following RB functions are commonly used