{'fraction_non_alphanumeric': 0.039582310678745146, 'fraction_numerical': 0.010933107233700745, 'mean_word_length': 4.045024763619991, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient evaluation of AGP reduced density matrices

We propose and implement an algorithm to calculate the norm and reduced density matrices of the antisymmetrized geminal power (AGP) of any rank with polynomial cost. Our method scales quadratically per element of the reduced density matrices. Numerical tests indicate that our method is very fast and capable of treating systems with a few thousand orbitals and hundreds of electrons reliably in double-precision. In addition, we present reconstruction formulae that allows one to decompose higher order reduced density matrices in terms of linear combinations of lower order ones, thereby reducing the computational cost significantly.


I. INTRODUCTION
In electronic structure theory, a geminal is a wavefunction for two electrons. Geminals are central to the concept of bonding and have a long history in quantum chemistry. 1 The geminal creation operator can, in general, be written as where η pq is antisymmetric, 2M is the total number of spin-orbitals and c † p is the creation operator of a fermion in spin-orbital p.
The simplest geminal wave function is perhaps the antisymmetrized geminal power (AGP) where all geminals are the same 2 Here |− is the physical vacuum containing no electrons, N is the number of pairs (2N electrons), and the factor 1/N ! is introduced for convenience. Without loss of generality, we choose to work in the natural orbital basis of the geminal. This is accomplished by applying a unitary transformation that brings η pq into a block diagonal form, wherein the one-body density matrix is diagonal and all spin-orbitals are paired. In this basis, it is mathematically simpler to work with hardcore boson operators, so that where such that c † p is the fermion creation operator in orbital p, andp is the "paired" companion of p. The pair creation and annihilation operators, P † p , and P p along with N p are the generators of a su(2) algebra N p , P † q = 2δ pq P † q . Notice that, written in its natural orbital basis, AGP preserves the so-called seniority symmetry. That is, orbitals with bars and no bars are either both occupied or empty. While in a general Hamiltonian the senioritydefining pairing scheme is arbitrary but can be optimized, in Hamiltonians for which seniority is a symmetry the pairs are naturally defined.
In the nuclear structure and condensed matter physics communities, AGP is known as the number projected Bardeen-Cooper-Schrieffer (PBCS) wavefunction. 3,4 The claim to fame of AGP is its ability to describe off-diagonal long range order, a criterion for superconductivity, 5 without breaking number symmetry as in BCS. 6 AGP is not a great wavefunction per se in quantum chemistry or condensed matter physics, because in most situations electron pairs are very different from each other. However, AGP is potentially an excellent starting point for geminal correlation models because it encompasses a combinatorial number of Slater determinants. Indeed, it is easy to see that This is a superposition of all possible seniority zero (paired) determinants involving N pairs in M orbitals but the coefficients are factorized by η pi rather than being a tensor. In the latter case, when coefficients are general rather than factorized, the wavefunction is known as doubly occupied configuration interaction (DOCI). 7-10 DOCI is obtainable by exact diagonalization over the space of seniority zero determinants with combinatorial cost as a function of M and N . However, AGP (PBCS) can be solved variationally and optimized using symmetry breaking and restoration techniques with mean-field O M 3 cost. 11,12 Other classes of geminal theories like Bethe ansatz (BA) for solving the integrable Richardson-Gaudin Hamiltonians 13,14 or APIG (antisymmetrized product of interacting geminals) wherein all geminals are different 15 are very interesting models but have combinatorial cost when applied to general (non-integrable) Hamiltonians. AGP is qualitatively correct for attractive pairing interactions at all correlation regimes where many other models have serious difficulties. [16][17][18][19] This makes AGP an attractive starting point for more accurate correlated geminal theories. Recent work in our research group 20 has shown that an AGP-based configuration interaction (CI) is a promising step in this direction. For AGP-based correlated theories to be computationally affordable for large systems, it is crucial that the AGP reduced density matrices (RDMs), loosely defined here as the expectation value of strings of ordered generators P † p ...N q ...P r ... over AGP, be obtainable with low cost. This is a necessary and crucial ingredient for developing successful geminal theories based on AGP.
In order to meet this goal, we introduce two techniques in this paper. First, we develop an efficient way of calculating the individual elements of RDMs to all ranks. For this, we formulate all RDMs in terms of the elementary symmetric polynomials and use the sumESP algorithm 21-23 to compute the sum efficiently. We argue that, with appropriate normalization, this is a reliable, fast, and stable method capable of treating large systems. Secondly, we rigorously prove that all AGP RDMs are expressible in terms of linear combinations of lower rank RDMs, using what we call reconstruction formulae. This is the most significant and novel contribution of this paper and has important theoretical and numerical implications. The significance of this relies in the fact that all correlated theories require high rank RDMs-often as high as 5-and 6-body. Therefore the ability to breakdown high-rank RDMs makes AGP a good starting point for correlated methods. Indeed, our result is reminiscent of Hartree-Fock theory wherein all high rank RDMs can be written as products of 1-body RDMs. From a computational perspective, the reconstruction formulae presented here reduce the cost and scaling of correlated AGP calculations significantly as demonstrated in Sec. V.

II. BASIC EXPRESSIONS
The analytic expressions for the norm and RDMs of AGP already exist in the literature and have been expressed in many different forms. 2,[24][25][26] For the sake of completeness, and to familiarize the reader with our notation, we present our own version of these derivations. In so doing, we introduce a Lie algebraic approach to understanding how P † p , P p , and N p act on the manifold of AGP and its excitations. And, we purposefully formulate all the matrix elements in such a way that they can be directly computed by the elementary symmetric polynomials.

A. Norm of the AGP wavefunction
The norm of the AGP wavefunction corresponding to N pairs and M orbitals can be obtained by calculating the contractions explicitly. From the commutation relations in Eq. (6) and the fact that (P qi ) 2 = 0, it follows that where S M N is the elementary symmetric polynomial (ESP) of degree N with M variables associated with the vector {η 2 1 , ..., η 2 M }. It is often convenient and numerically better posed to work with normalized AGP, i.e. AGP |AGP = 1. One can easily verify that the following choice does the job:

B. Differential Representation
Before we derive the expressions for the RDMs, we need to understand how P † p , P p , and N p act on AGP. The AGP state and its excitations describe a Riemann manifold with a positive semidefinite metric. In Lie algebra terms 27 this implies that generators acting on AGP may be represented as differential operators. By direct calculation, one can show that from which we get where |N is the shorthand notation for |AGP with N pairs. On the other hand, using explicit derivatives with respect to η p on Γ † N , one obtains This results in A similar derivation for P p yields and therefore Eq. (13), and Eq. (15) are the differential representation of operators N p and P p over AGP. Noting that P † p is the Hermitian conjugate of P p and acts to the left, we have the differential representation for all the generators as desired. As a corollary to these equations and the commutations relations of Eq. (6), it is easy to show that (N p ) 2 = 2N p and N p = 2P † p P p over AGP. We frequently use these properties in this paper.

C. AGP reduced density matrices
Consider a many-body system of fermions. To evaluate the energy or other observables thereof, one needs to compute the many-body RDMs. Since AGP preserves seniority and the total number of pairs is fixed, only two kinds of contractions are nonzero Clearly, all terms like that of Eq. (16b) can be written as P † p and P q by definition; and all like that of Eq. (16a) can be replaced by N p /2 as all electrons come in pairs. Therefore, all electron density matrices can be written as linear combinations of terms like P † p ...N q ...P r ... . Recall that N p = 2P † p P p over AGP; therefore, we can ultimately normal order everything such that all P † are on the left and all P are on the right, i.e. P † p ...P q ... . We refer to these density matrices as pair RDMs.
We set out by getting the expression for the n = 1 pair RDM Notice that and from this it follows where Generalization of this to higher rank RDMs follows the same reasoning. By induction on n, one can show that the matrix elements of the n-pair RDM, γ (n) , is γ p1...pn where d is the number of unique indices among {p i }∪{q i }.
Note that γ p1...pn q1...qn is symmetric with respect to permutations of p i 's and q i 's. If p i = p j , or q i = q j for some i, j, then the corresponding matrix element is instead zero. Therefore, in practice, one needs to compute only the unique elements, e.g. γ p1<...<pn q1<...<qn , which requires at most M n 2 elements-as opposed to M 2n . There is an ex- We make use of this property later in Sec. IV.
There is a special case of γ p1...pn q1...qn that we need later in this paper; that is when p i = q i for all i. Again, by N p = 2P † p P p , this can be written in terms of N p1 N p2 ...N pn . We refer to this as a number RDM. Formally, we define the number RDM of rank n, ν (n) , as where all the indices are assumed to be differentotherwise it would be equal a to lower rank number RDM by N 2 p = 2N p . From Eq. (20) this can be computed by III. NUMERICAL ALGORITHM From Sec. II it is clear that computing ESP efficiently is imperative to using AGP for any realistic system. This is because there are M N summands in every S M N and a straightforward summation could grow combinatorially with system size. In this section, we introduce the sumESP algorithm 23 that calculates ESP with polynomial cost. First, we briefly review some important remarks about the error analysis that were extensively studied in Ref. 22,23 In Sec. III B and III C, we study the time scales of various computations involving the norm of AGP and the corresponding RDMs as a function of N and M . The computer environments used for the runtime measurements can be found in Appendix C.

A. sumESP and error analysis
The algorithm we use to calculate the ESP is one initially proposed by Fischer. 21,28 This is the same summation algorithm used today in MATLAB's poly function. Following the notation of Ref. 23 we call it sumESP. The sumESP algorithm takes advantage of the following property of ESP Intuitively, Eq. (23) says that we can split any ESP into two sums such that one contains some arbitrary term η 2 p in all of its summands and one that does not. In the context of quantum chemistry, Eq. (23) is mentioned in Ref. 29 However, there has been substantial progress in the computational and applied mathematics community to better understand and craft the algorithm. In particular, Ref. 22 performed an error analysis and proved the stability of the algorithm. Ref. 23 made a slight improvement to the roundoff error bound and introduced new ways of calculating sumESP with enhanced precision. From these analyses it follows that sumESP is highly accurate and stable for positive summands 22 -as is the case for the overlaps of AGP. Quantitatively, one can say that, given {η 2 1 , ..., η 2 M } as a set of floating-point numbers, the "worst case" error due to sumESP is bounded above as follows 23 where S * M N is the exact value of ESP and S M N is computed using sumESP in floating-point arithmetic; denotes machine epsilon (unit roundoff). The error bound assumes that there are no numerical overflow and underflow occurring anywhere in the calculation. 22 In practice, however, poor choice of normalization of η p could lead to overflow and/or underflow. For example, when η p > 1, the sum could overflow even for small values of M . Generally, sumESP is better conditioned when 0 < η p < 1; this is because the underlying summations in the algorithm prevents the summands from getting too small. Nevertheless, the sum is prone to underflow when N is close to M and N is large. (Fortunately, this is not a concern in physical problems since the number of orbitals is always much greater than the number of electrons.) And it is prone to overflow when 1 N < M/2 and M is large. To make this analysis more quantitative, for 0 < η p < 1, we have plotted the magnitude of sumESP (the same as the overlap of AGP) as a function of M in Fig. 1. In Fig. 2 (the right axes), we fix M = 10 2 , 10 3 and vary N . It follows from the plots that, in very large systems, one must take extra caution when performing calculations away from N M or half-filling. But in moderate or small systems overflow and underflow should not be a concern in physical problems. It is noteworthy that, often times, the overflow and underflow is not due to sumESP per se as the exact values of the ESP itself is not expressible in double-precision. Issues of this kind can generally be avoided by scaling η's appropriately for a certain regime.

B. Runtime cost of an individual matrix element
The cost of calculating a single matrix element of npair RDM is bounded above by the cost of the norm of AGP, i.e. AGP |AGP . This is because all higher order RDMs require evaluation of lower degree ESPs by Eq. (20) and the cost of the prefactor is negligible. As such, to get the upper bound of the cost, we only report the elapsed time for evaluating the norm of AGP, which we refer to as the overlap. we remind the reader that it is sufficient to merely compute and store γ p1<...<pn q1<...<qn . Since the calculations of the matrix elements are independent from each other, this is highly parallelizable. (See Appendix A for the pseudocode used to calculate the matrix elements.) The theoretical cost of constructing an n-pair RDM on a single core is  (N (M − N )).
In Fig. 4 we report the average elapsed time of n-pair RDM for n = 1, 2, and 3, with M ≤ 100. The plots show the sample average of 100 observations for n = 1, 2 and  10 observations for n = 3 in which η p ∼ Unif(0,1).

IV. RECONSTRUCTION FORMULAE
In Sec. III C we argued that the asymptotic cost of constructing an n-pair RDM is O (N M 2n+1 ). Here, we show a way of cutting down the cost by expressing higher order RDMs as a linear combination of lower order ones. This is an important ingredient for developing efficient correlated theories based on AGP, as most techniques require evaluation of high order RDMs. For example, in AGP based configuration interaction (AGP-CI) calculations done in our group 20 up to to 5-body density matrices were needed (the Hamiltonian and the correlators contain 2-body operators each, and using the killers one can introduce a commutator to reduce the rank by 1). Therefore, a systematic way of reducing the scaling and the cost of calculating many-body RDMs is of great interest. Rosina mathematically anticipated 26 that the 2-RDMs of AGP determine all of its higher order RDMs. As it turns out, we here show that 1-RDM occupation numbers RDMs (ν (1) ) and η's are sufficient to determine all RDMs of AGP.
Note that the decomposition method presented here differs from the decomposition of PBCS in which the density matrix is a weighted integral over the projection grid of a factorized transition density matrix. 12

A. Direct decomposition
Our goal, here, is to express the nonzero elements of any n-pair RDM, γ (n) , as a linear combination of ν (n ) where n ≤ n ≤ 2n; then we show that we can further decompose each ν (n ) into a sum of ν (1) , hence γ (n) − →  ν (1) . Consider an n-pair RDM, γ p1...pn q1...qn . In general, it can be that p i = q j for some i, j. Since N p = 2P † p P p , we can write ...(P † p n−k P q n−k + P † q n−k P p n−k )|AGP .
Now, by manipulating the killer of AGP, i.e. K pq |AGP = 0 reported in Ref. 20 we can write By plugging Eq. (28) into Eq. (27) we arrive at an expression for γ (n) that is written purely as a linear combination of ν (n ) as desired. Now, we decompose ν (n ) into a sum of ν (1) and factors of η p . For this, we introduce a closed form expression which we prove in Appendix B To make this discussion concrete, consider γ (2) as an example, whose non-zero elements can be written as For each of the three cases, we obtain an expression entirely in terms of the number RDMs and η's using Eq. : Then each of the N p N q , N p N q N r , etc., can be inserted in Eq. (29) to produce an expression in terms of a linear combination of N p and η's only, as desired. This result is profound as it implies that we can obtain all higher rank pair RDMs by merely computing ν (1) whose cost grows asymptotically as O(M 2 N ) and a prefactor of n(n − 1) due to the direct decomposition. (For all practical purposes n N, M ; thus the cost of the prefactor is manageable.) However, we pay the price of introducing 1/(η 2 p − η 2 q ) factors that can be numerically ill-posed when η 2 p ≈ η 2 q or when η p , η q − → 0. In the regime that these factors are not problematic, the decomposition γ (n) − → ν (1) is a major improvement over computing all the matrix elements of an n-pair RDM. We must note that in our own implementation of these equations for practical problems (e.g. the attractive pairing Hamiltonian in Ref. 20 ) we have not observed these potential numerical issues.

B. Stepwise decomposition
In practice, there are situations in which it is more advantageous to break down high rank RDMs in terms of "slightly" lower rank ones. This also makes the issue of having too many 1/(η 2 p − η 2 q ) factors less severe. To this end, we need a new notation for P † ...N...P... . Define where k is the number of N p operators in the middle, and m is the total number of indices such that k ≤ m. And define t = (m − k)/2. For example, Notice that the subscript indices on each Z (k,m) are all different; otherwise they are either zero or reducible to some other Z (k ,m ) by N p = 2P † p P p and N 2 p = 2N p .
Obviously all γ (n) can be mapped to Z (k,2n−k) for some k and n, and vice versa. Our goal, here, is to show that As we prove in Appendix B, this can be accomplished by using the following formula at every step here, i, j ∈ {1, 2, ...n}. The advantage of breaking down the density matrices like this is that, at every step, we reduce the dimension by one, thereby reducing the asymptotic scaling of computing it by a factor of M 2 and a negligible prefactor. Moreover, we can stop at any step of our choice based on the cost that we are willing to tolerate.

V. RUNTIME OF ENERGY AND AGP-CI
In this section, we benchmark the speed gain when using the reconstruction formulae. Here, we report the runtime measurements of two calculations: (1)  calculations as reported in Ref. 20 . See Appendix C for the computer environments used in these benchmarks. We start with the pairing Hamiltonian. Recall that the attractive pairing Hamiltonian can be written as follows: 18 The expected value of the energy over AGP in terms of the pair RDMs is And using the reconstruction formulae, we can get Fig. 5 shows the elapsed time differences between implementing Eq. (37) and Eq. (38) as a function of number of levels at half-filling. For each case, the corresponding RDM is stored and then called in the calculation. By inspection, it is easy to see that the reconstruction formula speed up the calculations often by an order of magnitude. Obviously, the improvement becomes more noticeable as the number of levels increases.
At last, we report the elapsed time of performing AGP-CI calculations with and without the reconstruction formulae in Fig. 6. For these calculations, we stored up to 4-indexed RDMs in memory and calculated all the rank 5 RDMs using the direct decomposition formula Eq. (29). Here, η's are optimized beforehand for the pairing Hamiltonian (G = 2) for various values of M . Similar to the energy calculations, we observe that the reconstruction formulae lead to a substantial improvement.

VI. CONCLUSIONS
Analytic expressions of the norm and reduced density matrices of AGP wavefunction are proportional to ESP. We used the sumESP algorithm 22,23 to efficiently calculate ESPs and argued that, with appropriate normalization, it is well suited for physical problems wherein N M . We have shown that our method can reliably calculate the norm and elements of RDMs to all ranks for systems as large as a few thousand orbitals and hundreds of electrons. Our runtime measurements indicate that the asymptotic cost per element of an n-pair RDMs is at most quadratic, and the cost of building an n-pair RDM grows asymptotically as O(N M 2n+1 ).
However, to reduce the cost of computing the RDMs even further, we have derived reconstruction formulae that allow decomposition of any RDM into a linear combinations of lower rank ones. We introduced two methods: (1) Direct decomposition-breaks down a high rank pair RDM in terms of linear combination of rank-1 occupation RDMs; (2) Stepwise decomposition-reduces the dimension of a RDM by one at every step thereby reducing the cost of computing it by a factor of M 2 .
We demonstrate the advantage of using our reconstruction formulae by benchmarking it against the energy of the pairing Hamiltonian and AGP-CI calculation without the reconstruction formulae. The numerical results indicate that, indeed, the reconstruction formulae lead to a substantial speed-up, especially in systems with large number of orbitals. One can store any RDM as a 2 dimensional array by linear indexing p i 's and q i 's such that γ p1<...<pn q1<...<qn → γ µ ν . The reader may find the following relation handy in the parallel implementation µ = p 1 + n i=2 p i − 1 i ; 1 ≤ p 1 < ... < p n ≤ M.