A reduced basis approach for calculation of the Bethe-Salpeter excitation energies using low-rank tensor factorizations

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


Introduction
In modern material science, there is a growing interest to ab initio computation of the absorption spectra for molecules or surfaces of solids. This computational problem can be treated either by using the time-dependent density functional theory (TD-DFT) [1][2][3][4][5][6] or by solving the Bethe-Salpeter equation (BSE) [7,8] based on the Green's function formalism and many-body perturbation theory [9][10][11][12][13]. A specific choice of the approximate computational model may depend on many physical and CONTACT Boris N. Khoromskij bokh@mis.mpg.de * Dedicated to Prof. Andreas Savin on occasion of his th birthday.
implementation aspects, see [9] for the detailed discussion. In particular, the BSE approach leads to the challenging numerical task concerning the solution of large eigenvalue problem for a dense matrix that, in general, is non-symmetric.
In the present paper, we consider the computational aspects of the large-scale algebraic BSE spectral problem when using the data-sparse matrix structures. We follow the particular formulation of the BSE problem based on the non-interacting Green's function via representation in terms of the Hartree-Fock (HF) molecular orbitals (MOs) [14,15], where it was applied to H 2 molecule in the minimal basis of two Slater functions. 1 In the framework of this specific BSE formulation, we focus on the algebraic aspects of solving the computationally extensive spectral problems arising in the case of larger molecular systems. It is demonstrated that this scheme becomes practically applicable to moderate-size molecules when using tensor-structured HF calculations [16][17][18][19], accomplished by efficient representation of the two-electron integrals (TEI) in the MO basis in the form of a low-rank Cholesky factorisation [17,20]. In this way, the low-rank representation of the TEI tensor stipulates the beneficial structure of the BSE matrix blocks, thus enabling numerical algorithms of reduced complexity.
It is worth to note that the size of the BSE matrix scales quadratically with the size of the atomic orbital (AO) basis set, O(N 2 b ), used in ab initio HF calculations. The direct diagonalisation is limited by the O(N 6 b ) complexity, making the problem computationally expensive already for moderately sized molecules with a basis size N b 100. Hence, a procedure that relies entirely on multiplication of the governing BSE matrix, or its approximation, with vectors (in the framework of some iterative procedure) is the only viable approach. In turn, fast matrix computations can be based on the use of low-rank representations since such data structures allow efficient storage and fast algebraic operations with linear complexity scaling in the matrix size.
Methods for solving partial eigenvalue problems for matrices with the special structure as in the BSE eigenvalue problem have been intensively studied in the literature. These structures are related to the so-called Hamiltonian matrices, exposing a particular block structure. Papers and books treating Hamiltonian eigenvalue problems include [21][22][23][24][25]; see also the recent survey [26] and the references therein. Special cases of the BSE and other eigenvalue problems related to HF approximations lead to anti-block-diagonal Hamiltonian eigenproblems that can be solved by special techniques based on minimisation principles [27,28]. The algebraic structure of the BSE matrix is not that of a Hamiltonian matrix in the general case, but yields a so-called complex J-symmetric matrix. Theory and numerical solution of such eigenvalue problems are discussed in [29][30][31][32][33], where the particular instance of the BSE matrix is considered in [33]. Other partial eigensolvers tailored for electronic structure calculations are discussed in [34,35]. The reduced basis method for large-scale systems is described in [36].
In this paper, we study a reduced basis approach to the approximate numerical solution of the BSE eigenvalue problem based on model reduction via projection onto a reduced basis, which is constructed by using the eigenvectors of a simplified system matrix obeying a simpler data-sparse structure. The reduced basis method includes two steps. First, the diagonal plus low-rank approximation to the fully populated blocks in the BSE matrix is calculated, enabling an easier partial eigenvalue solver for a large auxiliary system relying only on matrix-vector multiplications with rank-structured matrices. Second, a small subset of eigenvectors from the auxiliary eigenvalue problem is selected to build the projection of the exact BSE system onto this reduced basis set.
The approximation error incurred by the reduced basis approach depending on the rank truncation parameters is investigated. Theoretical and numerical analysis on the existence of the low-rank approximation and the respective rank bounds for matrix blocks in the BSE system matrix are presented. One of the favourable features of the approach is the quadratic convergence rate in approximate excitation energies compared with the accuracy of the reduced basis set thresholded by a rank truncation parameter ε > 0 (see Remark 3.3).
The reduced basis approach applies to the BSE system with matrix blocks of size denote the number of occupied and virtual HF orbitals, respectively, such that , the direct numerical calculation of the matrix elements, based on the precomputed TEI tensor in the HF MO basis, has a storage and numerical cost of the order of O(N 4 b ). The construction of the reduced basis and the preceding low-rank decomposition of matrix blocks in the Bethe-Salpeter kernel are motivated by the use of truncated Cholesky factorisation of the TEI matrix [17,20]. To that end, the BSE matrix blocks are represented in terms of the precomputed Cholesky factors in the HF MO basis. Along with the diagonal energy matrix, this constitutes the structured representations of the dielectric and response functions, as well as the static screened interaction matrix. Taking into account the rank decomposition of TEI, the above quantities tolerate the low-rank approximation up to a chosen threshold. 2 This yields the construction of a so-called simplified BSE matrix with a diagonal plus low-rank structure in the matrix blocks, thus admitting efficient storage and matrix-vector products. The reduced basis is obtained by calculating several of the lowest eigenvectors of the auxiliary eigenvalue problem for the simplified matrix. A projection of the exact BSE matrix onto the reduced basis set and diagonalisation of the arising small-size matrix completes the reduced basis scheme.
Numerical tests for single molecules and finite chains of hydrogen atoms indicate the convergence in the senior excitation energies by the increase of the separation ranks.
The rest of the paper is organised as follows. Section 2 recalls the truncated Cholesky decomposition scheme for low-rank factorisation of the TEI tensor in the HF MO basis, that is the building block in the construction of the BSE matrix. Section 3 describes the algebraic computational scheme for evaluation of the entries in the BSE matrix, analyses the low-rank structure in different matrix blocks and describes the reduced basis approach. Section 4 presents numerical tests for several compact and extended molecules demonstrating the computational features of the reduced basis method applied to the full BSE system as well as to simplified model by the socalled Tamm-Dancoff approximation (TDA).

Cholesky decomposition of the TEI matrix
The numerical treatment of the TEI tensor (also known as the electron repulsion integrals) is one of the bottlenecks in the numerical solution of the HF equation and in density function theory (DFT) calculations for large molecules [37].
Given the AO basis set The matrix B is proven to be symmetric and positive definite ensuring application of the incomplete Cholesky decomposition [37,[39][40][41]. The tensor-based HF solver [16,19] employs the efficient calculation of the Cholesky factors [17,20] in where the adaptively chosen column vectors in B are calculated in an efficient way. This allows the partial decoupling of the index sets {μν} and {λσ }.
Notice that the Cholesky factorisation (2.2) can be written in the following index form: where the second factor corresponds to the transposed matrix L T k . Here L k = L k (μ; ν), k = 1, …, R B , denotes the N b × N b matrix unfolding of the column vector L(:, k) from the Cholesky factor L ∈ R N 2 b ×R B . Numerical experiments indicate that the truncated Cholesky decomposition with the separation rank O(N b ) ensures a satisfactory numerical precision ε > 0 of order 10 −5 to 10 −6 . The refined rank estimate O(N b |log ε|) was observed for every molecular system considered so far [17,20].
In the standard quantum chemical implementations in the Gaussian-type AO basis, the numerically confirmed rank bound, rank(B) ≤ C B N b (C B is of order 1-10), allows to reduce the complexity of building up the Fock matrix F to O(N 3 b ), which is by far dominated by the computational cost for the exchange term K(D) (see Section 2.2 ).

Rank bounds for the TEI matrix V
The 2N o -electron HF equation for N o pairwise L 2orthogonal occupied electronic orbitals, ψ i : Gaussian-type AOs. We consider the complete set of HF MOs {C p ∈ R N b }, i.e. the pth column vectors in the coefficients matrix C ∈ R N b ×N b , and the corresponding energies {ε p }, p = 1, 2, …, N b . The occupied MOs ψ i are represented (approximately) by the coefficient The coefficients matrix C ∈ R N b ×N b solves the nonlinear eigenvalue problem (2.5) where H is the core Hamiltonian matrix, S is the mass matrix and In BSE calculations, the TEI tensor B = [b μνλσ ], corresponding to the AO basis set, is represented in the MO basis: The BSE calculations utilise the two subtensors of V specified by the index sets The first subtensor is defined as in the MP2 calculations [17]: while the second one lives on the extended index set: In the following, {C i } and {C a } denote the sets of occupied and virtual orbitals, respectively. We shall also use the notation (2.9). The straightforward computation of the matrix V by the above representations accounts for the dominating impact on the overall numerical cost of order O(N 5 b ) in the evaluation of the block entries in the BSE matrix. A method of complexity O(N 4 b ) based on the low-rank tensor decomposition of the matrix V on the full index set was described in [17].
It can be shown that the rank R B = O(N b ) approximation to matrix B LL T with the N × R B Cholesky factor, L, allows to introduce the low-rank representation of the tensor V, and then to reduce the asymptotic complexity of calculations to O(N 4 b ) (see [17]). Indeed, let C m be the mth column of the coefficient matrix (2.10) A similar factorisation can be derived in the case of (2.9). The precise formulation is given by the following lemma [17], which will be used in further considerations.
where the columns of L V are given by On the index set (2.9), The numerical cost is determined by the computation complexity and storage size for the factors L V , U V and W V in the above rank-structured factorisations.
Lemma 2.1 provides the upper bounds on rank(V ) in the representation (2.10) which might be reduced by the ε-rank truncation. It can be shown that the ε-rank of the matrix V remains of the same magnitude as that for the TEI matrix B obtained by its ε-rank truncated Cholesky factorisation (see the numerical illustration in Section 3.2).
Numerical tests in [17] indicate that the singular values of the TEI matrix B decay exponentially as where the constant z > 0 depends weakly on the molecule configuration. If we define R B (ε) as the minimal number satisfying the condition then the estimate (2.11) leads to the ε-rank bound R B (ε) ࣘ CN b |log ε|, which will be postulated in the following discussion.
Our goal is to justify that R V (ε) increases only logarithmically in ε, similar to the bound for R B (ε). To that end, we introduce the singular value decomposition (SVD) decomposition of the matrix B, which can be written in the following index form:

Lemma 2.2:
For given ε > 0, there exists a rank-r approximation V r to the matrix V, and a constant C > 0 not depending on ε, such that r ࣘ R B (ε) and Proof. We estimate the R B (ε)-term truncation error by using the representation (2.13), which can be presented in the matrix form V = taking into account that U k = 1, k = 1, …, R B , and the Frobenius norm estimate , then the multiple of ε|log ε| in (2.15) does not depend on ε, which proves our lemma.
The storage cost of these decompositions restricted to the active index set Figure 1 represents the singular values of the matrix V for H 32 chain, N 2 H 4 and C 2 H 5 NO 2 (Glycine amino acid) molecules with the size of the basis set (N b , N o ) equal to (128, 16), (82,9) and (170, 20), respectively. They indicate that R V (ε) is linearly proportional to |log ε|. Moreover, it is of the same order of magnitude as R B (ε) (see [20]).
The calculation of V r is based on the reduced truncated SVD algorithm applied to the initial rank-R B Cholesky decomposition of the matrix V inherited from that for the TEI matrix B (see Lemma 2.1). Complexity of this straightforward computation on the active index set can be estimated by

Tensor factorisation of the BSE matrix blocks
Here we discuss the main ingredients for calculation of blocks in the BSE matrix and their reduced rank approximate representation. We compose the 2N ov × 2N ov BSE matrix by Equations (46a) and (46b) in [14], though the construction of static screened interaction matrix w(ij, ab) in Equation (3.4) may slightly differ.

Tensor representations using TEI matrix in MTO basis
Construction of the BSE matrix includes computation of several auxiliary quantities. First, introduce a fourthorder diagonal 'energy' matrix by that can be represented in the Kronecker product form where I o and I v are the identity matrices on respective index sets. It is worth noting that if the so-called homolumo gap of the system is positive, i.e.
then the matrix ε is invertible. Using the matrix ε and the N ov × N ov TEI matrix V = [v ia, jb ] represented in the MO basis as in (2.7), the dielectric function (N ov × N ov matrix) Z = [z pq, rs ] is defined by z pq,rs := δ pr δ qs − v pq,rs [χ 0 (ω = 0)] rs,rs , with χ 0 (ω) being the matrix form of the so-called Lehmann representation to the response function. In turn, the matrix representation of the inverse of χ 0 (ω) is known to have the following form: Let 1 ∈ R N ov and d ε = diag{ ε −1 } ∈ R N ov be the allones and diagonal vectors of ε −1 , respectively, specifying the rank-1 matrix 1ࣹd ε . In this notation, the matrix Z = [z pq, rs ] takes a compact form where V is calculated by (2.8). Lemma 2.1 suggests the existence of a low-rank factorisation for the matrix W defined above.
which can be calculated by solving linear system with structured data (see Section 3.2).
Furthermore, Equation (46a) in [14] includes matrix entries w ij, ab for a, b ∈ I v , i, j ∈ I o . To this end, the modified matrix W = [ w pq,rs ] is computed by (3.2) on the index set {p, q ∈ I o } × {r, s ∈ I v } by using entries v i j,ab in the matrix unfolding of the tensor V in (2.9) multiplied from the left with the N 2 o × N 2 o submatrix of Z −1 . Now the matrix representation of the BSE in the (ov, vo) subspace reads as the following eigenvalue problem determining the excitation energies ω n : where the matrix blocks are defined in the index notation by (see (46a) and (46b) in [14] for more details): In the matrix form, we obtain where the matrix elements in W = [w ia, jb ] are defined by w ia, jb = w i j,ab , computed by (3.2) and (2.9) as described above. Here the diagonal plus low-rank sparsity structure in ε + V can be recognised in view of Lemma 2.1. For the matrix block B, we have where the matrix V , corresponding to the partly transposed tensor, is defined entrywise by thus coinciding with V in (2.8) due to the symmetry properties. Here W is defined by permutation, In the following, we investigate the ε-rank structure in the matrix blocks A and B resulting from the corresponding factorisations of V. Solutions of Equation (3.3) come in pairs: excitation energies ω n with eigenvectors (x n , y n ), and de-excitation energies −ω n with eigenvectors (x * n , y * n ). The block structure in the matrices A and B is inherited from the symmetry of the TEI matrix V, v ia, jb = v * ai,bj and the matrix W, w ia, jb = w * bj,ai . In particular, it is known from the literature that the matrix A is Hermitian and the matrix B is (complex) symmetric (since v ia, bj = v jb, ai and w ib, aj = w ja, bi ), which we control in the matrix construction (see also [33] for implications on the algebraic properties of the BSE matrix).
In the following, we confine ourselves to the case of real spin orbitals, i.e. the matrices A and B remain real. It is known that for the real spin orbitals and if A + B and A − B are positive definite, the problem can be transformed into a half-size symmetric eigenvalue equation [3]. Indeed, in this case for every eigenpair, we have , depending on how to compute the matrix W. Here, the low-rank structure in the matrix V can be adapted.
The challenging computational tasks arise in the case of lattice-structured compounds, where the number of basis functions increases proportionally to the lattice size L × L × L, i.e. N b ∼ n 0 L 3 , that quickly leads to intractable problems even for small lattices.

The reduced basis approach using low-rank approximations
The large matrix size in Equation (3.3) makes the solution of the full eigenvalue problem computationally intractable even for moderate-size molecules, not to speak of lattice-structured compounds. Hence, in realistic quantum chemical simulations of excitation energies, the calculation of several (tens of) eigenpairs may be sufficient.
In the following, we show that the part ε + V in the matrix block A has diagonal plus low-rank (DPLR) structure, while the submatrix V in the block B exhibits lowrank approximability. Taking into account these structures, we propose a special partial eigenvalue problem solver based on the use of a reduced basis obtained from the eigenvectors of the reduced matrix that picks up only the essential part of the initial BSE matrix with the DPLR structure. The iterative solver is based on fast matrix-vector multiplication and efficient storage of all data involved in the computational scheme. Using the reduced basis approach, we then approximate the initial problem by its projection onto a reduced basis of moderate size.
We begin with the low-rank decomposition of the matrix V, where the rank parameter R V = R V (ε) = O(N b |log ε|) can be optimised depending on the truncation error ε > 0 (see [17] and Section 2.2). First, we represent all matrix blocks and intermediate matrices included in the representation of the BSE matrix by using the above decomposition and diagonal matrices as follows. The properties of the Hadamard product imply that the matrix Z exhibits the representation where the rank of the second summand does not exceed R V . Hence, the linear system solver W = Z −1 V can be implemented by algorithms tailored to the DPLR structure by adapting the Sherman-Morrison formula. The computational cost for setting up the full BSE matrix F in (3.3) can be estimated by O(N 2 ov ), which includes the cost O(N ov R B ) for generation the matrix V and the dominating cost O(N 2 ov ) for setting up W . In the following, we rewrite the spectral problem (3.3) in the following equivalent form: The main idea of the reduced basis approach proposed in this paper is as follows. Instead of solving the partial eigenvalue problem for finding of, say, m 0 eigenpairs satisfying Equation (3.7), we first solve the slightly simplified auxiliary spectral problem with a modified matrix F 0 . The approximation F 0 is obtained from F 1 by using lowrank approximations of the parts W and W of the matrix blocks A and B, respectively, i.e. A and B are replaced by respectively. Here we assume that the matrix V is already represented in the low-rank format inherited from the Cholesky decomposition of the TEI matrix B.
The modified auxiliary problem reads This eigenvalue problem is much simpler than that in (3.3), since now the matrix blocks A 0 and B 0 , defined in (3.8), are composed of diagonal and low-rank matrices.
Having at hand the set of m 0 eigenpairs computed for the modified (reduced model) problem (3.9), {(λ n , ψ n ) =  (λ n , (u n , v n ) T )}, we solve the full eigenvalue problem for the reduced matrix obtained by projection of the initial equation onto the problem adapted small basis set {ψ n } of size m 0 . Define a matrix G 1 = ψ n (:, 1 : m 0 ) ∈ R 2N ov ×m 0 whose columns present the spanning vectors of the reduced basis, compute the stiffness and mass matrices by projection onto the reduced basis specified by the columns in G 1 , and then solve the projected generalised eigenvalue problem of small size m 0 × m 0 , The portion of small eigenvalues γ n , n = 1, …, m 0 , is thought to be very close to the corresponding excitation energies ω n , (n = 1, …, m 0 ) in the initial spectral problem (3.3). Table 1 illustrates that the larger the size m 0 of the reduced basis is, the better is the accuracy of the lowest excitation energy γ 1 , as to be expected.

Remark 3.2:
Notice that the matrix W might have rather large ε-rank for small values of ε, which increases the cost of high-accuracy solutions. Numerical tests show (see Table 3) that the ε-rank approximation to the matrix W with a moderate rank parameter allows for a numerical error in the excitation energies of the order of few percents. For this reason, we study another approximation strategy in which the rank approximation of the matrix W remains fixed, while the matrices V and W are substituted by their adaptive ε-rank approximations (see Table 3).
Matrix blocks in the auxiliary equation (3.9) are obtained by rather rough ε-rank approximation to the initial system matrix. However, we observe much smaller approximations error γ n − ω n for the solution of the projected reduced basis system (3.10) compared with that for auxiliary equation (3.9). Numerical tests indicate that Table . Accuracy (in eV) for the first eigenvalue, |γ  − ω  |, and norms of the differences between the exact and reducedrank matrices, F  − F  , vs. ε-rank for V, W and W . the difference γ n −ω n behaves merely quadratically in the rank truncation parameter ε.

Remark 3.3:
In the case of a symmetric matrix, the above-mentioned effect of 'quadratic' convergence rate can be justified by a well-known property of the quadratic error behaviour in the approximate eigenvalue, computed by the Rayleigh quotient with respect to the perturbed eigenvector (vectors of the reduced basis ψ n in our construction), compared with the perturbation error in the eigenvector, which is of order O(ε). This beneficial property may explain the efficiency of the proposed reduced basis approach.
In the particular BSE formulation based on the HF MO basis, we may have a slight perturbation of the symmetry in the matrix block W , i.e. the above argument does not apply directly. However, we observe the same quadratic error decay in all numerical experiments implemented so far.
It is also worth to note that due to the symmetry features of the eigenproblem, the approximation computed by the reduced basis approach is always an upper bound of the true excitation energies obtained from the full BSE model. Again, this is a simple consequence of the variational properties of the Ritz values being upper bounds on the smaller eigenvalues for symmetric matrices. The 'upper bound' character is also clearly visible in the figures in Section 4.

Numerical tests for the reduced basis method
In this section, we present numerical illustrations of the reduced basis approach applied to the BSE problem for single molecules and finite chains of hydrogen atoms. The TEI tensor and MOs are obtained from ab initio HF calculations using tensor-structured solver [16][17][18][19] implemented in MATLAB R .
Both the core Hamiltonian and two-electron (repulsion) integrals are computed by rank-structured algorithms using the discrete representation of the AO basis functions on fine n × n × n three-dimensional (3D) Cartesian grids [17,20]. In TEI calculations for molecules, the basis functions and convolution kernels involved are represented on fine 3D grids of size up to 131, 726, 3 which guaranties the sufficient accuracy of numerical quadratures. The TEI matrix is precomputed in the form of lowrank Cholesky factorisation by tensor-based algorithm incorporating 1D density fitting [20].

Reduced basis method for the BSE system
The numerical examples below demonstrate that a small reduced basis set, obtained by separable approximation of the BSE matrix blocks with rank parameters of about several tens, allows to reveal several of the lowest excitation energies. Accuracy is controlled by the rank truncation threshold. Examples below utilise the grid representation of the Gaussian basis sets of type cc-pDVZ (see e.g. [42,43]). Table 2 presents the size of GTO basis set, N b , and the number of MOs, N o , in the numerical examples considered. Table 3 shows numerics for H 2 O (360 × 360), N 2 H 4 (1314 × 1314) and C 2 H 5 OH (2860 × 2860), where the numbers in brackets specify the BSE matrix size. It demonstrates the quadratic decay of the error |γ 1 − ω 1 | in the lowest excitation energy with respect to the approximation error |λ 1 − ω 1 | for the modified auxiliary BSE problem (3.9). Errors for eigenvalues are given in eV. The numerical error is controlled by a tolerance ε in the rank truncation procedure applied to the BSE submatrices V, W and W . The resulting ε-ranks for the corresponding matrices are presented. Notice that the rank decomposition of the matrix V can be derived from the respective Cholesky factorisation of the TEI matrix B accomplished by the simple rank reduction. The rank approximation for the symmetric matrices W and W can be calculated by pivoted Cholesky factorisation. Table 2 demonstrates that the approximation error in the reduced basis, |γ 1 − ω 1 |, is at least one order of magnitude smaller than that for auxiliary problem, |λ 1 − ω 1 |, i.e. |γ n − ω n | |λ n − ω n |, which justifies the use of the reduced basis equation (3.10).
This effect can be also seen in Figure 2 for N 2 H 4 molecule demonstrating the convergence γ n → ω n and λ n → ω n with respect to the increasing rank parameter determining the auxiliary problem (the size of the reduced basis set is m 0 = 30). The left and right figures correspond to the rank truncation thresholds, ε = 0.6 and ε = 0.1, respectively. The quantities λ n , γ n and ω n are marked by black, blue and red colours, respectively. Figure 3 represents similar results for amino acid glycine, C 2 H 5 NO 2 , with the BSE matrix size 6000 × 6000. In this case, the truncation thresholds ε = 0.2 leads to the rank parameters R V = 54 , R W = 50, R W = 50 and the error for the minimal eigenvalue, ω 1 = 0.72 eV. For  The lowest values of the BSE excitation energy for H 2 O molecule computed by solving our exact system is 8.72 eV which agrees with the value 8.7 eV for ice water presented in [44]. The reduced basis method using the rank truncation threshold ε = 10 −1 provides the value 8.95 eV. Figure 4, left and right, illustrates the BSE energy spectrum of the H 2 O molecule for the lowest N red = 30 eigenvalues vs. the rank truncation parameters ε = 0.6 and 0.1, where the ranks of V and the BSE matrix block W are equal to 4, 5 and 28, 30, respectively, while the block W remains unchanged. For the choice ε = 0.6 and ε = 0.1, the error in the first (lowest) eigenvalue for the solution of the problem using the reduced basis is about 0.11 and 0.025 eV, respectively.
Next, we present BSE calculations for chains of 16 and 32 hydrogen atoms placed in a 3D bounding box with the size 64 3 bohr. 3 The interatomic interval equals to 1.39 bohr. The HF calculations were performed with 64 and 128 Gaussian-type basis functions using grids of size 32, 768 3 and 16, 384 3 , for computation of the core Hamiltonian and TEI, respectively. Cholesky factors of TEI matrix are of size 4096 × 175 and 16, 384 × 348 for the chains with 16 and 32 hydrogen atoms. Table 4 demonstrates the decay of the error in the lowest eigenvalues |γ 1 − ω 1 | with respect to the tolerance ε in the rank approximation of the BSE matrix for the chains of hydrogen atoms, H 16 (896 × 896) and H 32 (3584 × 3584), where the numbers in brackets specify the BSE matrix size. We observe linear scaling of the corresponding ranks of V,W and W with respect to the size of the system as expected. Since the rank of W decays slowly, we studied the case, with fixed rank(W ) = max {rank(V), rank W } (usually it coincides with rank(V)). The improved accuracy of the resulting spectrum is achieved even for rather large ε = 2 × 10 −1 .

Reduced basis approach to the Tamm-Dancoff model
It is interesting to apply the reduced basis approach described above to the so-called TDA [3], which corresponds to setting the matrix B = 0 in Equation (3.3). It also allows to estimate the difference between the excitation energies from the full BSE model and those obtained by the TDA, which introduces an additional small model error.
The TDA model simplifies Equation (3.3) to a standard Hermitian eigenvalue problem with the factor-two smaller matrix size N ov . The reduced basis approach via low-rank approximation described in Section 3.2 can be applied directly to the TDA equation. Below we present numerical tests indicating that the approximation error introduced by the TDA compared with the initial BSE system (3.3) remains on the level of 0.003 Hartree for several compact molecules (see Table 5). This table indicates a tendency to decrease the TDA model error for larger molecules, say 0.0017 Hartree (0.045 eV) for glycine amino acid.

Conclusions
This paper introduces and analyses the reduced basis method for computation of several lowest eigenvalues in the BSE, based on the solution of an auxiliary, simplified eigenvalue problem via diagonal plus low-rank approximations to the BSE matrix blocks. The reduced spectral problem of small size is derived via projection of the full BSE matrix onto the reduced basis set, composed of several dominant eigenvectors of the simplified problem. The ε-rank bounds for the requested subtensors of the TEI tensor represented in the basis set of HF MOs are proved in Lemmas 2.1 and 2.2. Asymptotic estimates on storage demands are provided. Numerical tests confirm merely quadratic error behaviour in the excitation energies with respect to the accuracy of the rank approximation. The particular construction of the BSE matrix is based on the BSE-GW approximation with noninteracting HF Green's function (we follow the construction considered in [14]).
The main goal of the present paper is the development and numerical verification of the reduced basis method applied to large-scale BSE eigenvalue problems. Potential efficiency of the approach is illustrated numerically on the examples of single molecules and finite hydrogen chains. The numerical studies demonstrate that eigenvalues of the reduced spectral problem provide a sufficient approximation to the lowest excitation energies of the exact BSE system. For all the examples considered so far, the accuracy of the order of 0.15 − 0.3 eV was achieved 3 by the reduced basis method with rather moderate ranks as indicated in Tables 3 and 4. The behaviour of the approximation error vs. rank parameters remains similar for compact molecules and for finite chains of atoms.
Future work will be focused on application and development of efficient linear algebra algorithms for fast and accurate solution of the arising large eigenvalue problems taking into account data-sparse block structures. In particular, the complementary structural representations to the matrix W may enhance the accuracy of the reduced basis method. These issues will be considered in a forthcoming paper.
Another possible research direction is concerned with the quantised tensor train (QTT) approximation [45] of the long vectors and large matrices involved, in order to perform the fast matrix-vector calculations in the QTT tensor arithmetics with ε-rank truncation, see [46]. Notes 1. In [14], it was demonstrated that in the case of small interatomic distance this model remains in a good agreement with the reference data for full configuration-interaction calculations though it does not describe the dissociation case. The BSE model was shown to provide satisfactory results in the latter case when using the exact Green's function. 2. Though one can notice that the static-screened interaction matrix, responsible for the exchange interaction of electrons, does not allow the high-accuracy low-rank approximation, as it is already known about the HF exchange. 3. For example, paper [11] assumes the BSE model errors in the range of 0.1-0.3 eV as the acceptable precision.