
==== Front
J Chem Theory Comput
J Chem Theory Comput
ct
jctcce
Journal of Chemical Theory and Computation
1549-9618
1549-9626
American Chemical Society

38788209
10.1021/acs.jctc.4c00352
Article
Reducing the Runtime of Fault-Tolerant Quantum Simulations in Chemistry through Symmetry-Compressed Double Factorization
https://orcid.org/0000-0003-2122-6933
Rocca Dario *†
Cortes Cristian L. †
https://orcid.org/0000-0002-2933-4085
Gonthier Jérôme F. †
Ollitrault Pauline J. †
Parrish Robert M. †
Anselmetti Gian-Luca ‡
https://orcid.org/0000-0002-8850-7708
Degroote Matthias *‡
Moll Nikolaj ‡
Santagati Raffaele ‡
Streif Michael ‡
† QC Ware Corporation, Palo Alto, California 94306, United States
‡ Quantum Lab, Boehringer Ingelheim, 55218 Ingelheim am Rhein, Germany
* Email: dario.rocca@qcware.com.
* Email: matthias.degroote@boehringer-ingelheim.com.
24 05 2024
11 06 2024
20 11 46394653
20 03 2024
13 05 2024
08 05 2024
© 2024 American Chemical Society
2024
American Chemical Society
https://creativecommons.org/licenses/by-nc-nd/4.0/ Permits non-commercial access and re-use, provided that author attribution and integrity are maintained; but does not permit creation of adaptations or other derivative works (https://creativecommons.org/licenses/by-nc-nd/4.0/).

Quantum phase estimation based on qubitization is the state-of-the-art fault-tolerant quantum algorithm for computing ground-state energies in chemical applications. In this context, the 1-norm of the Hamiltonian plays a fundamental role in determining the total number of required iterations and also the overall computational cost. In this work, we introduce the symmetry-compressed double factorization (SCDF) approach, which combines a CDF of the Hamiltonian with the symmetry shift technique, significantly reducing the 1-norm value. The effectiveness of this approach is demonstrated numerically by considering various benchmark systems, including the FeMoco molecule, cytochrome P450, and hydrogen chains of different sizes. To compare the efficiency of SCDF to other methods in absolute terms, we estimate Toffoli gate requirements, which dominate the execution time on fault-tolerant quantum computers. For the systems considered here, SCDF leads to a sizable reduction of the Toffoli gate count in comparison to other variants of DF or even tensor hypercontraction, which is usually regarded as the most efficient approach for qubitization.

document-id-old-9ct4c00352
document-id-new-14ct4c00352
ccc-price
==== Body
pmc1 Introduction

Quantum chemistry simulations hold a significant potential to advance many industry-relevant applications, including the development of new drugs,1 catalysts,2 and materials.3 The simulation of chemical systems from first-principles requires the solution of the Schrödinger equation, a task particularly challenging for classical approaches because of the exponential growth of the computational cost with system size. Quantum computing provides a promising solution to address this scalability issue, with significant ongoing efforts focused on developing resource-efficient algorithms.4 Much of this work has been dedicated to approaches tailored for early stage noisy quantum devices, such as the variational quantum eigensolver (VQE).5,6 Besides the challenges of working with noisy hardware, optimizing the parameters in the VQE ansatz is nontrivial, and the number of required measurements grows rapidly with the system size.

The significant challenges to achieving a quantum advantage in near-term noisy devices motivate current efforts to transition toward fault-tolerant quantum computing (FTQC). Early hardware demonstrations of error-corrected logical qubits have already been achieved,7,8 and many companies, including IBM9 and Google,10 have announced roadmaps to build FTQ computers in the next few years. At the same time, a parallel effort is underway to develop quantum algorithms that can efficiently exploit error-corrected qubits.

Quantum phase estimation (QPE) can be considered as the prototypical algorithm for chemistry simulations within the FTQC framework.11,12 Within the standard formulation of QPE, the calculation of the ground state energy of a given chemical system relies on the implementation of the Hamiltonian evolution operator for some duration τ; this operator can be approximated in practice using the Trotter–Suzuki13 formula or Taylor series expansion.14 More recently, an alternative approach for QPE has been proposed based on the quantum walk operator .15−17 Instead of directly returning the ground state energy, the algorithm outputs the arccosine of the ground state energy. The advantage of this procedure is that the quantum circuit corresponding to the quantum walk operator can be implemented exactly using qubitization.18 The parameter λ in the definition of corresponds to the 1-norm of the Hamiltonian, and its value is influenced by the specific representation of Ĥ and the strategy used to block encode it. This parameter plays a crucial role in the QPE efficiency, and its optimization is one of the main topics of this work.

Within the qubitization-based QPE approach, the total number of Toffoli gates scales as . Here, ϵ represents the accuracy required for the ground state energy (typically, this should be within the chemical accuracy threshold of 1.6 mHa), and the ratio λ/ϵ determines the total number of iterations; is the Toffoli gate cost per iteration and depends on the specific approach used for implementing . Implementing Toffoli gates or, similarly, T gates on quantum hardware requires a procedure known as magic state distillation.17,19,20 This process demands a considerably large number of qubits and takes significantly more time than other operations in the computation. For this reason, a reduction in the overall Toffoli gate count for a quantum algorithm is expected to lead to an equivalent reduction in the overall runtime. While Toffoli gates indeed dominate the computational cost and our discussion will primarily revolve around these gates, it is important to recognize the impact of the 1-norm in determining the total number of required QPE iterations and circuit depths. Accordingly, a method that decreases the 1-norm, such as the Hamiltonian factorization proposed in this study, not only reduces the demand for Toffoli gates but also for all other gates necessary for circuit compilation.

While the 1-norm plays a fundamental role in determining the total number of iterations, also has a significant contribution to the overall computational cost. Specifically, the computational complexity of realizing depends on Γ, the amount of information needed to specify the Hamiltonian, and the specific approach employed for the quantum implementation. As discussed in ref (21), the combination of tensor factorizations with techniques such as unary iteration17 and optimized Quantum Read-Only Memory (QROM) assisted by ancillae22,23 leads to a Toffoli gate and logical qubit scaling. Further details on the origin of this square root dependence will be provided in Section 3.1. A summary of the computational complexity of state-of-the-art approaches for qubitization-based QPE is presented in Table 1. Beyond a different cost in the implementation, these approaches also involve different definitions of the 1-norm λ, whose scaling varies between and , in which N is the number of spatial orbitals the Hamiltonian is expressed in.

Table 1 Asymptotic Scaling of the Computational Resources Required by Different Approaches Used in Qubitization-Based Quantum Phase Estimationa

approach	logical qubits	Toffoli gates	
sparse method23	

	

	
single factorization23	

	

	
XDF24	

	

	
tensor hypercontraction21	

	

	
a The definition of the 1-norm λ depends on the specific implementation. The ratio of λ to the required precision ϵ in the final result determines the total number of iterations. N denotes the number of orbitals, S is the sparsity of the Hamiltonian, and Ξ is the average rank of the second factorization.

A straightforward implementation based on the electronic Hamiltonian in second quantization involves terms. To improve over this complexity, a sparse method was introduced that truncates the components of the two-electron tensor according to a chosen threshold.23 The main limitation of this approach is that the number of remaining terms S in the Hamiltonian cannot be systematically predicted and, in some cases, still behaves as . The single factorization (SF) approach applies an eigendecomposition to the two-electron tensor (see eq 17 below) and this effectively decreases the number of terms in the Hamiltonian to .23 The explicit double factorization (XDF) approach introduces a second factorization on top of the SF (see eq 18 below).23,25−27 This reduces the Hamiltonian to pieces of information, where Ξ is the average rank of the second tensor factorization. The numerical experiments for hydrogen chains considered in Section 3.4 show that Ξ itself is characterized by a behavior. As discussed in Section 2.2, an alternative approach known as compressed DF (CDF) builds the tensors in the factorization by optimizing a suitable cost function.28−30 The new methodology presented in this paper will be based on a variant of the CDF approach. The recent work of von Burg et al.24 has introduced an efficient quantum algorithm to implement the double-factorized Hamiltonian in the QPE framework by employing Givens rotations and qubitization.24 Compared to a straightforward qubitization of the double-factorized Hamiltonian, this formulation also benefits from significantly reducing the 1-norm.

The tensor hypercontraction (THC) approach decomposes the two-electron tensor in the Hamiltonian as1

where χpμ and ζμν denote the components of the tensors used for this decomposition, NTHC is the THC rank, and p, q, r, and s are indices identifying the spatial orbitals.31,32 The tensors are obtained by minimizing a cost function that determines the deviation of the decomposition from the exact two-electron tensor. The application of this approach in the context of QPE was first proposed in ref (21). To effectively decrease the 1-norm, the χ tensors were used as basis set rotations, applying them to redefine the representation of the corresponding creation and annihilation operators in the second-quantized Hamiltonian. In practice, this amounts to reformulating the Hamiltonian in a larger nonorthogonal basis set and new techniques were developed to block encode and qubitize it.21 Beyond decreasing the 1-norm, the THC approach provides a very compact representation of the Hamiltonian, with and, correspondingly, an improved asymptotic computational complexity. Since this method has systemically provided the most favorable resource estimations for many examples of electronic Hamiltonians,21,33 it will serve as the main benchmark for the methodological developments proposed in this work.

This work is largely focused on the reduction of the 1-norm that has a strong impact on the number of iterations and, accordingly, on the total runtime of the QPE algorithm. Different approaches have been proposed in the literature to optimize the 1-norm. The XDF in the implementation of von Burg et al.24 and the THC21 benefit themselves from formulations that significantly reduce the 1-norm as compared to a straightforward transformation of the electronic Hamiltonian into Pauli words. The optimization of the 1-norm for quantum simulations has been considered in previous work. Orbital transformations were proven to improve the 1-norm values significantly.34 In ref (35), several different approaches (including orbital transformation) were compared by considering small molecules in the minimal STO-3G basis set; it was shown that DF coupled with a symmetry shift provides the best results in terms of 1-norm reduction and scaling with the system size. This symmetry shift approach, described in detail in Section 2.3, effectively decreases the 1-norm by subtracting a function of the number operator of electrons from the electronic Hamiltonian. Since the number operator of electrons commutes with the Hamiltonian, the eigenvectors of the Hamiltonian are not affected by this shift, and the correct ground state energy can be obtained by applying a simple a posteriori correction.

The new symmetry-CDF (SCDF) approach introduced here exploits the symmetry shift idea but additionally optimizes the DF tensor decomposition to further decrease the 1-norm. Numerical demonstrations of this method include active space models of the FeMoco molecule and cytochrome P450, and hydrogen chains with up to 80 atoms. For all of these systems, SCDF, to the best of our knowledge, provides the smallest values of the 1-norm reported in the literature. This leads to a Toffoli gate count and runtime that are sizeably reduced with respect to THC (for example, by one-half for FeMoco and P450). The SCDF has the same structure of XDF and can be implemented using the techniques proposed by von Burg et al.24 Accordingly, SCDF inherits an analogous computational complexity both in terms of Toffoli gate and logical qubit requirements (see Table 1) but with a 1-norm that scales more favorably with the number of orbitals compared to XDF (see Section 3.4). Despite a slightly worse asymptotic behavior than THC, the numerical applications considered in this work show that SCDF provides more systematic accuracy for ground state energies owing to a simpler numerical optimization scheme. This feature is crucial to address systems of large size and to obtain reliable properties in the thermodynamic limit.

2 Methodological Approach

To introduce our new SCDF approach (Section 2.4), it is important to first review some of the main ideas at the base of this factorization. Section 2.1 introduces the general DF framework and discusses some of its properties and implications for the quantum implementation. Section 2.2 describes some numerical techniques used to compute the DF in practice, one based on eigenvalue decomposition (XDF) and one on the optimization of a cost function (CDF). Section 2.3 discusses the symmetry shift approach, a general technique recently developed to reduce the 1-norm of the electronic Hamiltonian. This approach is applied in practice by shifting the Hamiltonian by a function of the number operator of electrons. Section 2.4 finally introduces our original contribution, the SCDF, that combines ideas from CDF and symmetry shift. By optimizing a regularized cost function, the SCDF tensor decomposition is built to effectively compensate the symmetry shift operator while preserving a good level of accuracy for two-electron integrals.

2.1 General Double Factorization Framework

Within the second quantization formalism, the electronic Hamiltonian is expressed as2

where p, q, r, and s are indices identifying the N spatial orbitals, and the singlet spin-summed one-particle substitution operators are defined as . In this definition, and â denote creation and annihilation operators, respectively, and the bar on top of the orbital indexes indicates a ↓ spin orbital.

In 2 the constant term Enuc corresponds to the nuclear repulsion energy3

is the two-electron tensor, and the modified one-electron tensor is defined in terms of the one-electron integrals4

that include the kinetic and electron–nucleus interaction energies; ZI here denotes the atomic number of the I-th atom and RI denotes the position of its nucleus. Without loss of generality for molecular systems, the spatial orbitals ϕ(r) have been chosen to be real.

The second quantized Hamiltonian can then be expressed in a quantum computing amenable form by expanding it in terms of Pauli words using, for example, the Jordan–Wigner or Bravyi–Kitaev transformations.36−38 These approaches lead to number of terms in the Hamiltonian. Despite the polynomial growth, this number of terms poses practical challenges for noisy near-term and fault-tolerant algorithms, and, in this context, the DF of the Hamiltonian can provide several advantages.

The main idea of DF consists in decomposing the two-electron tensor in the following way23,25−28,39,405

where the Ut tensors are orthonormal, namely6

and the “core” tensors7

are symmetric for all t’s. The sum over t in eq 5 runs up to a maximum value NDF that depends on the specific approach used to build the tensor factorization and crucially determines the trade-off between the accuracy and efficiency of the method. Depending on the specific DF implementation, the sums over k and l can also be limited to values Ξ(t) ≤ N. The truncation of the tensors within the DF approach will be further discussed below.

By inserting the factorized tensors of eq 5 into eq 2, it is possible to reformulate the second quantized Hamiltonian in terms of the operators and , that create and annihilate electrons, respectively, in a new set of rotated orbitals. In order to apply these operators, it is convenient to introduce the operators, which rotate the quantum state in the new orbital basis and, using the Thouless theorem,41 can be expressed as8

These rotations can be formulated in terms of Givens rotation networks that can be efficiently implemented on quantum hardware.26 Within the double-factorized formalism, the Hamiltonian in second quantization can then be expressed as9

since the operators and are diagonal in the Jordan–Wigner representation, in practice, the two-body part of the Hamiltonian has been decomposed in NDF terms that are separately diagonal in the basis set representation obtained by applying . We can now apply explicitly a Fermion-to-qubit mapping based on the Jordan–Wigner transformation.36 Since within this framework and10

the Hamiltonian can be finally expressed as11

In this equation, the one-electron tensor has been redefined as and a “single factorized” (eigenvalue) decomposition has been applied to obtain12

The constant term contains Enuc and additional constant terms originating from the one and two-body operators.

The block encoding of the double-factorized Hamiltonian can be obtained by straightforward application of the linear combination of unitaries (LCU) approach to the Hamiltonian in the form of eq 11.18,42 In this case, the 1-norm takes the following value13

An alternative method to block encode the double-factorized Hamiltonian has been introduced by von Burg et al.24 To introduce this approach, we assume that the Vt tensors have rank one and are positive definite for every value of t. As discussed in the next section, not all the approaches to building the DF satisfy these properties, and this has important repercussions on the efficiency of a specific method. A positive-definite rank-one Vt can always be decomposed as14

by replacing this factorization in eq 9 and reapplying the Jordan–Wigner transformation of the Hamiltonian takes the form15

where the sums over k have been truncated at Ξ(t) ≤ N by eliminating the elements of the Wt tensors below a certain threshold δDF; Ξ, which denotes the average of the Ξ(t) values, has been used in the Introduction to discuss the computational complexity of DF (see Table 1). For this formulation of the electronic Hamiltonian, we have , since in the two-body term, we have a sum over NDF (which is itself a sum over Ξ, and N additional degrees of freedom in the basis rotations. Indeed, as shown in ref (24), an efficient quantum implementation of and can be obtained by decomposing these operators in single term rotations, requiring N angles each to be defined. Within this formulation, the 1-norm of the Hamiltonian is given by16

This approach has two main advantages: (1) the Hamiltonian in eq 15 can be efficiently implemented using qubitization;18 (2) the 1-norm in eq 16 is typically significantly smaller with respect to the LCU 1-norm in eq 13. As discussed in the following sections of the paper, the implementation of von Burg et al. benefits significantly from the low rank of the Vt tensor, and this will be an important feature included in our new SCDF methodology.

2.2 Explicit and Compressed Double Factorization

In the previous section, we introduced the general DF formalism and presented the benefits of this approach. Here, we explain the main practical schemes that can be used to build the tensor factorization in eq 5. These approaches fall into two main categories: (1) XDF, which builds the Vt and Ut tensors using a two-step eigenvalue or Cholesky decomposition;23,25−27 (2) CDF and its variants that build those tensors, optimizing a cost function.28,29

Within the framework of XDF, the two-electron tensor is first decomposed in terms of eigenvalues and eigenvectors17

where we introduced the definition . A second factorization can then be obtained from the eigendecomposition of the Lt tensors18

whose rank has been truncated to a t-dependent value Ξ(t), which is at most equal to N. The combination of these two equations provides a tensor decomposition in the form of eq 5 by defining19

which is clearly consistent with the general definition in eq 14. According to 19 the Vt tensor has rank 1 and, as already mentioned in the previous section, this is an important feature to implement efficiently the double-factorized Hamiltonian using the approach of von Burg et al.24 This XDF procedure is, in principle, exact if no truncation is applied, namely, NDF = N2 and Ξ(t) = N. In practice, it is well-known from widely used techniques such as density fitting or Cholesky decomposition43−46 that the rank of the two-electron integrals to achieve reasonable accuracy is much smaller than N2 and, in practice, NDF scales as . Concerning the truncation of the second factorization, while the average rank Ξ can be smaller than N its complexity still behaves as .

To further improve the scaling or at least the numerical complexity prefactor of the DF, it is alternatively possible to build the Hamiltonian from Vt and Ut tensors obtained by minimizing a suitable cost function. While Ut and Vt are constrained to still be orthogonal and symmetric, respectively, this type of approach takes advantage of the full rank of these tensors. This is the main idea at the base of CDF, which determines the DF tensors by minimizing the cost function2820

where denotes the Frobenius norm. With respect to XDF, a smaller value of NDF can typically achieve the same level of accuracy in the final result. Minimization of the cost function can be achieved more efficiently by alternating the optimization of the Vt and Ut tensors. The optimization with respect to Vt can be recast in the form of a linear system and solved with standard linear algebra libraries. The orthogonality of the Ut tensors should be constrained during the optimization. As explained in Appendix A, this problem can be reformulated in an unconstrained form by introducing antisymmetric orbital rotation generators. The derivatives of the cost function with respect to the components of the generators can be evaluated analytically and, since, the minimization problem is nonlinear in this case, a numerical unconstrained continuous optimizer such as limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm (L-BFGS) is used.

The CDF approach provides a more compact representation in terms of the NDF rank required for a given level of accuracy but usually converges to tensors Vt with large components, which in turn leads to 1-norms comparable to or even larger with respect to XDF. For this purpose, a regularized CDF (RCDF) has been recently introduced that is based on the following cost function29,3321

where the components of the tensor ρtkl are usually fixed to a constant value and γ takes the values 1 or 2 for L1 and L2 regularization, respectively. The last term in the cost function can be considered a penalty function that prevents the elements of the Vt tensor from becoming too large and, accordingly, limits the growth of the 1-norm. In ref (29) all numerical applications were based on the L2 regularization, and this approach showed a sizable decrease of the λDFBurg norm with respect to XDF and THC.

It is important to notice that in both the CDF and RCDF approaches, the Vt tensors are not necessarily positive definite, and their rank is unconstrained during the optimization. As discussed in ref (29), to generalize the approach of von Burg et al. to the (R)CDF case, it is possible to introduce the factorization22

However, this approach has two major disadvantages. First, the components of Wt are complex, requiring some modifications of the original implementation of DF based on qubitization. Second, with respect to eq 14, the tensor decomposition in eq 22 involves an additional sum over N terms. This implies that the amount of information Γ defining the Hamiltonian has a complexity , which is of the same order as the unfactorized Hamiltonian and strongly affects the computational resource requirements. Using the RCDF tensor factorization provided in ref (29), in Section 3.3, we will show how Toffoli gate and logical qubit requirements are affected in practice for the case of cytochrome P450. Since our new SCDF method is based on a cost function analogous to those used in CDF and RCDF, this issue is overcome by constraining the Vt tensor to be rank 1 during the optimization.

2.3 Symmetry Shift Approach

The main idea of symmetry shift involves replacing the Hamiltonian Ĥ with , where Ŝ(a) is a generic symmetry operator that satisfies .35,47 The set of parameters a = {a1, a2, ...} is chosen to minimize the 1-norm. Since the symmetry operator commutes with the full electronic Hamiltonian, has the same eigenstates of Ĥ and can be directly used in QPE, taking advantage of the smaller 1-norm. The generic operator Ŝ(a) can be built as a function of one or multiple reciprocally commuting basic symmetry operators, of which we know the eigenvalue for the ground state ahead of time. The possible choices of these symmetry operators include the number operator of electrons , the z-projection of the spin, the total spin, and the molecular point group symmetries. The original work of Loaiza and Izmaylov35,47 focused on the shift, which was shown to be effective to significantly reduce the 1-norm. The use of to at most the power 2 preserves the structure of the Hamiltonian in second quantization (eq 2). It adds only spin-independent terms, where the electronic tensors run exclusively on spatial orbitals, and introduces no terms that include more than four orbitals (two-body). Accordingly, the quantum implementation of the shifted Hamiltonian can be based on the same methodologies used in the unshifted case. The use of other symmetry operators would instead modify the structure of the Hamiltonian and possibly impact the cost of the quantum implementation.

The use of other symmetry shift operators in Hamiltonian could involve the In practice, this method is based on the shifted Hamiltonian23

where the two terms corresponding to and are intended to decrease the 1-norm of the one and two-body components of the Hamiltonian, respectively. The strategy to optimize the parameters a is detailed in ref (35). For the one-body term of the Hamiltonian , we have24

where we have used the fact that commutes with the orbital rotation operator; the a1′ = -2a1 definition is not strictly necessary but has been introduced to simplify the equation. From the definition of the 1-norm in eq 16, it is clear that the optimal a1′ has to be chosen to minimize . The optimal value can be simply obtained from the median of the fkø coefficients.

Starting from eq 2, the symmetry shift for the two-body term can be written as25

where a2′ = 2a2. Similarly to the one-body case, the optimal value can be found by minimizing . Once the two-electron integrals have been redefined, including the symmetry shift, the double-factorized two-body Hamiltonian is obtained as in the regular XDF case, but the 1-norm is typically significantly reduced.

For several small molecules in the minimal STO-3G basis set, the XDF approach coupled with symmetry shift was shown to be a very promising method both in terms of the 1-norm values and overall scaling of the 1-norm as a function of the system size.35 Indeed, this approach outperformed many others, such as orbital optimization, anticommuting Pauli product grouping, and greedy Cartan subalgebra decomposition both with and without symmetry shift.

Our new methodology, which we will introduce in the next section, is exclusively focused on the 1-norm reduction of the two-body part of the Hamiltonian. The approach of eq 24 without modifications will be used for the one-body term.

2.4 Symmetry-Compressed Double Factorization

In this section, we introduce our new approach, which will be denoted as SCDF. This method combines some of the advantages of RCDF29 and symmetry shift35,47 to significantly decrease the 1-norms of the Hamiltonian, which results in lower Toffoli gate counts. As a first step to introducing the SCDF approach, we focus on the two-body part of the Hamiltonian and consider the identity26

where the Kronecker deltas in the third line were resolved using the orthogonality of the Ut tensors (see eq 6). The coefficient in front of has been decomposed as , where, in the ideal case, the differences should be as small as possible to effectively decrease the 1-norm. There are two main differences between the symmetry shift technique proposed by Loaiza et al. (eq 25) and the approach that we are proposing in eq 26: in eq 25, the symmetry shift is applied before the DF (namely, gpqrs – a2′δpqδrs is factorized), while in our approach it is applied after (namely, gpqrs is factorized and then the shift is applied); in eq 26, the “global” symmetry shift is decomposed into different t-dependent contributions. We empirically observed from numerical experiments that eq 26 is less effective than eq 25 in reducing the 1-norm, when applied within the XDF framework. However, the Vt tensors themselves can be optimized to decrease the fluctuations of and this is the main idea of the SCDF approach. In practice, this can be achieved by minimizing the cost function27

with respect to the Vt and Ut tensors and the prefactors αt of the symmetry shift. The tuning of the regularization coefficient ρ determines the trade-off between the accuracy in the tensor decomposition of the two-electron integrals and the decrease of .

Similarly to RCDF, the formulation of the cost function in eq 27 leads to full rank Vt tensors that increase the number of terms in the double factorized Hamiltonian by a factor N with respect to the XDF approach. To constrain the rank of Vt to be equal to 1, we factorize the tensor as in eq 14 and modify the SCDF cost function to be28

which is now optimized with respect to Wt rather than Vt. The downside is that the easy single-step update of Vt is now replaced with a higher order dependence on Wt. Because of the constraint on the rank of Vt, the cost function of SCDF has less variational freedom compared to other CDF methods. However, the numerical applications of Section 3 show that minima with very small 1-norms can be achieved at the price of a large number of iterations in the optimization of the cost function.

It is important to notice that while Vklt = WktWlt is rank 1, its shifted counterpart WktWlt – αt is rank 2, unless αt is 0 (see Appendix B for a detailed discussion). From a numerical standpoint, this is effectively equivalent to doubling NDF, with potential negative consequences on quantum computing resources. In practice, only a limited number Nα of αt’s have an optimized value different from 0 and, accordingly, the cost of the implementation does not significantly change with respect to the basic XDF. In the resource estimation presented in Section 3, the influence of these additional terms is taken into account. In the same section, we will provide additional quantitative details on the number of required αt’s and their influence on the 1-norm.

An additional important observation is that, while the Vt tensors are positive definite by construction, Vt – αt can have negative eigenvalues. Accordingly, implementing the double-factorized Hamiltonian based on eq 15 requires some modifications. This point is discussed in Appendix B.

The cost in the evaluation of the SCDF cost function (eq 28) or other CDF variants is dominated by the reconstruction of the approximate two-electron tensor from the Ut and Vt tensors (eq 5). Since NDF is itself , the computational complexity of this operation is . By assuming that the number of iterations required in the optimization does not depend on N, this numerical cost is expected to dominate the scaling of the SCDF tensor decomposition. In a similar way, the evaluation of eq 1 dominates the computational cost of the THC factorization, leading to a complexity. Alternative algorithms to compute the THC factorization have also been proposed in the literature (refs (31 and 32)) with complexity. However, their application is not straightforward if regularization is required to limit the value of the 1-norm (this is the case, for example, of ref (33)). The polynomial complexity to build the DF and THC factorizations is still significantly smaller than the exponential complexity involved in the classical solution of the Schrödinger equation.

3 Results

3.1 Computational Details

The SCDF numerical results for all the other systems considered in this work are based on a Python implementation that uses the JAX library.48 A description of the optimization procedure and the parameters used are provided in Appendix A. For all the numerical applications considered in this work, a value of 10–5 for the regularization coefficient ρ was found to be a reliable option. Indeed, if, for example, ρ is set to 10–4, the norm can be effectively optimized, but, at least in certain cases, chemical accuracy for ground state energies is not achieved; if instead the value of ρ is decreased to 10–6, the minimization of the cost function is slow, and this option should be considered only if a high level of accuracy is required.

The Toffoli gate and logical qubit requirements for the quantum implementation of XDF and SCDF are estimated using the OpenFermion library.49 For the FeMoco molecule and all the hydrogen chains independently of the size, 10 bits for state preparation and 16 bits for rotations were used;21 for cytochrome P450, we considered 10 bits for state preparation and 20 bits for rotations.33 As discussed in Appendix B, the quantum implementation of SCDF is analogous to XDF, and the same approach for resource estimation can be used.

It is important to consider that OpenFermion’s resource estimation for double-factorized Hamiltonians is based on a QROM algorithm that optimizes the number of Toffoli gates at the expense of a higher number of logical qubits.21−23 This subroutine is used in different steps of the quantum implementation and tends to dominate the estimate of the total computational cost (this is especially the case for the rotation angle lookup). By using auxiliary ancillae, the data lookup can be implemented using L/k + b(k – 1) Toffoli gates and b(k – 1) + log(L/k) ancillae, where L is the number of entries to load and b is the number of bits used to represent each entry. The parameter k, which must be a power of 2, controls the trade-off between the number of the required Toffoli gates and of logical qubits. If k = 1, the conventional implementation is recovered, requiring L Toffoli gates and log(L) auxiliary qubits. To minimize the Toffoli gate count, the k is chosen as close as possible to the optimal value ; with this choice both the Toffoli gate and logical qubit requirements behave as . The use of QROM for the rotation angle lookup tends to dominate the computational cost and overall scaling. The specific cost for this task is (NDFΞ)/kr + Nβ(kr – 1) for the Toffoli gates and Nβ(kr – 1) + log(NDFΞ/kr) for the auxiliary logical qubits, where β is the number of bits used to represent each single angle and kr denotes the specific k parameter used for this task. In this case, the number of Toffoli gates is optimized by , which leads to a cost for both the number of required logical qubits and Toffoli gates. By considering that NDF grows itself as , this discussion explains the computational complexity of DF, as reported in Table 1. Since NDFΞN corresponds to the amount of information Γ contained in the Hamiltonian, the use of this approach to trade off Toffoli gates for logical qubits explains the complexity discussed in the Introduction. While the choice of the optimal value of kr is crucial to reducing the Toffoli gate requirements and scaling, depending on the specific system and the number of logical qubits available, it may be of interest to find a different space-time balance. To this purpose, for the FeMoco and P450 active space models considered below, we also present results for “suboptimal” values of the kr parameter, which sizably decrease the number of required logical qubits while still providing a low number of Toffoli gates.

The rank of the first factorization NDF is chosen to be a multiple of the number of orbitals N, from a minimum of 4N to a maximum of 6N. The rank of the second factorization Ξ is determined by eliminating the components of the Wt tensor below the threshold δDF = 10–4; this choice ensures a significant reduction of the terms in the Hamiltonian while preserving a high level of accuracy. In applying the symmetry shift, only the components αt above the threshold δα = 10–3 are included.

3.2 Active Space Model of the FeMoco Molecule

We begin by applying our approach to simulate the ground state of an active space model of the FeMoco active site of nitrogenase, which plays a crucial role in understanding the mechanism of biological nitrogen fixation.50 This system was identified as a potential killer application of quantum computing because of its strong static correlation, which is challenging to simulate with classical techniques.51 The original active space studied by Reiher et al.51 was later improved by Li et al.52 For the purpose of demonstrating the efficiency of the SCDF approach, we focus here on the Reiher Hamiltonian. In ref (51), the geometry of the FeMoco model was first optimized using the B3LYP density functional53 with the def2-TZVP Ahlrichs triple-ζ basis set plus polarization functions on all atoms54 (see Figure 1). The ground state quantum calculation for the corresponding active space involves 54 electrons with N = 54 spatial complete active space self-consistent field (CASSCF) orbitals, amounting to 108 qubits or spin orbitals. This active space was obtained by first performing a smaller singlet CASSCF calculation with 24 electrons in 16 orbitals; this CASSCF calculation was then used to compute the natural orbitals that were selected to build the larger active space based on occupation number criteria.

Figure 1 Model systems of FeMoco from ref (51) (C, black; O, red; H, white; S, yellow; N, blue; Fe, orange; and Mo, cyan).

All of the methods compared in Table 2 involve truncations of the tensor ranks that control the trade-off between accuracy and computational efficiency. For the methods based on DF, the behavior with respect to the truncation of the first factorization is assessed by comparing the three values, NDF = 4N, 5N, and 6N. The average rank of the second factorization Ξ is also shown in Table 2. Interestingly, for XDF without or with symmetry shift, the application of the threshold δDF = 10–4 to eliminate the components of the Wt tensors has minimal effects on the total number of terms in the Hamiltonian and, accordingly, on Ξ. This is different for SCDF, where Ξ decreases when increasing NDF, leading to a total number of terms in the Hamiltonian that grows slowly with NDF (this implies also a rather steady number of Toffoli gates). This behavior is likely to be related to the characteristics of the SCDF cost function (eq 28), where large components of the Wt tensor are penalized. It is worth mentioning that ref (21) also presents results for XDF, but a different procedure was used to truncate the XDF tensor decomposition. All the N2 terms were initially maintained in the first factorization, and the pruning was exclusively performed for the second factorization, removing the jth components that satisfy (δDF′ serves the same purpose of δDF, but they are not strictly equivalent). Within this procedure, the truncation of the second factorization effectively decreases also the rank of the first factorization, NDF. The XDF resource estimation for FeMoco in ref (21) was performed choosing δDF′ = 0.00125, which leads to NDF = 360 and Ξ = 36. At first sight, the behavior of Ξ could seem radically different with respect to our results reported in Table 2. In practice, for the Wt tensors with the largest contribution to the tensor factorization, the two procedures provide similar values of the Ξ(t) rank. For Wt’s of decreasing importance, the procedure of ref (21) tends to keep more tensors, even with small values of Ξ(t); this explains the larger number of NDF and the smaller average rank Ξ form in ref (21). In practice, for FeMoco, our truncation scheme for XDF provides slightly more accurate ground state energies (0.24 mHa CCSD(T) error for NDF = 4N vs 0.44 mHa in ref (21)) and similar total numbers of terms in the Hamiltonian and resource estimations (11,596 terms for NDF = 4N vs 13,031 in ref (21)).

Table 2 Resource Estimates, Correlation Energy Errors, and 1-Norm Values (λ) for the Active Space Model of the FeMoco Moleculea

approach	NDF or NTHC	Ξ	CCSD(T) error (mHa)	λ (Ha)	Toffoli gates	logical qubits	
XDF	4N	54	0.24	293.9	9.6 × 109	3722	
XDF	5N	54	0.28	295.3	1.0 × 1010	3724	
XDF	6N	54	0.12	296.0	1.1 × 1010	3724	
 	 	 	 	 	 	 	
XDF + sym. shift	4N	54	0.25	182.9	6.0 × 109	3722	
XDF + sym. shift	5N	54	0.28	184.3	6.5 × 109	3724	
XDF + sym. shift	6N	54	0.13	185.0	7.1 × 109	3724	
 	 	 	 	 	 	 	
THCb	350 ≈ 6.5N	 	–0.29	306.3	5.3 × 109	2142	
 	 	 	 	 	 	 	
SCDF	4N	39	0.60	79.9	2.4 × 109	3719	
SCDF	5N	35	0.32	78.0	2.4×109	3722	
SCDF	6N	31	0.36	77.9	2.5 × 109	3722	
 	 	 	 	 	 	 	
SCDF (kr = 2)	5N	35	0.32	78.0	2.6 × 109	1994	
a Different approaches based on DF and THC are compared for different values of the rank used in the factorization. The best-performing methods in terms of Toffoli gates and logical qubits are highlighted in bold (while SCDF with NDF = 4N provides a slightly smaller Toffoli gate count than NDF = 5N, the result with the highest accuracy in the correlation energy is considered). If not explicitly indicated, the SCDF results were obtained with the optimal value kr = 4.

b Results from ref (21).

The errors in the ground state energy and the resource estimations for the different approaches considered here are presented in Table 2. A reliable tensor factorization of the Hamiltonian should preserve a high level of accuracy and, following ref (21), we consider the correlation energy of the coupled cluster with singles, doubles, and perturbative triples (CCSD(T)) as an error metric. As shown in the fourth column of Table 2, the error of all the different approaches is well below the chemical accuracy threshold of 1.6 mHa for all the values of NDF. The resources required by XDF agree with previous findings in the literature.21,24 The application of the symmetry shift significantly decreases the 1-norm to about 62% of the initial XDF value, and this is reflected in a very similar way in the number of required Toffoli gates. However, the symmetry-shifted XDF is still not competitive with the THC approach for both the numbers of required Toffoli gates and logical qubits. For the FeMoco model, the new SCDF approach provides a significant additional reduction of the 1-norm, with values amounting to about one-quarter of those obtained with the XDF and THC methods. Because of this significant 1-norm decrease, the SCDF approach requires less than half the number of Toffoli gates (and, accordingly, runtime) compared to the state-of-the-art THC method. Since the number of logical qubits weakly depends on the 1-norm (log λDF dependence), the SCDF and all the DF-based approaches tend to be equivalent with respect to the qubit requirements. However, if for SCDF, the kr parameter is decreased from its optimal value of 4 to 2, the number of Toffoli gates increases only slightly, while the logical qubit count becomes the smallest among all methods, as listed in Table 2. The same procedure could be applied to the other approaches to further decrease their logical qubit requirements, but this would further increase their Toffoli gate count. According to these observations, for the FeMoco active space model, the SCDF method can achieve a significant speed-up with respect to all the other approaches, and also provide an excellent balance between the number of Toffoli gates and logical qubits.

As discussed in Section 2.4, only a limited number of symmetry shift prefactors αt actually contribute to the 1-norm reduction. Using the δα = 10–3 threshold for FeMoco only 5 αt’s are retained for NDF = 4N = 216, 13 for NDF = 5N = 270, and 25 for NDF = 6N = 324; while this number grows with NDF, it remains limited and, as shown in Table 2, does not have a significant impact on the computational cost. This truncation has a negligible effect on the 1-norm (differences of the order of 10–8 Ha). Still, the symmetry shift plays a fundamental role in the 1-norm reduction: for example, if all the symmetry shift prefactors are set to 0 for NDF = 5N, the 1-norm increases to 145.8 Ha.

As a concluding note of this section, it is interesting to explore the effects of thresholds δα and δDF on the accuracy and computational requirements of SCDF. If for NDF = 5N, these two thresholds are set to 0 and all the terms of the factorization are included, the 1-norm and the correlation energy error are not affected (the error in the correlation energy actually slightly increases to 0.33 mHa). Instead, the computational requirements sizeably increase to 3.8 × 109 Toffoli gates and 7180 logical qubits. This finding demonstrates the importance of the sparsity in the SCDF tensor factorization, beyond the beneficial effects of the 1-norm reduction.

3.3 Active Space Model of Cytochrome P450

To further assess the efficiency and accuracy of the SCDF method, we consider here the cytochrome P450, which has been proposed as a benchmark system for fault-tolerant quantum algorithms in ref (33). The (34↑ + 29↓e, 58o) active space model of the Cpd I species (see Figure 2) was chosen as an example. It was obtained by performing a Restricted Open Shell Hartree–Fock (ROHF) calculation in the cc-pVDZ basis for the high-spin instance to generate the molecular orbitals. After that, the orbitals were split-localized with the Pipek-Mezey55 procedure, and the 58 active orbitals were selected based on their likeness to the atomic orbitals that are important for the correlations in this molecule. Resource estimates for the ground state calculation of the P450 model are reported in Table 3. The behavior of the different methods is similar to what we already discussed for the FeMoco case. In particular, our new SCDF method decreases to less than one-half the Toffoli gate requirements of THC, with an equivalent speed-up expected for the runtime of the quantum algorithm. It is also important to mention that for this application, the accuracy of THC is less systematic. As shown in the Supporting Information of ref (33), the CCSD(T) error tends to oscillate as a function of the THC rank. For example, if the rank is increased to 380 (to be compared to 320, the value chosen as optimal in ref (33)), the CCSD(T) error increases in absolute value to −0.83 mHa. The correlation energy error of SCDF and other DF-based approaches is instead significantly smaller and less-dependent on NDF. Compared to other DF variants, SCDF also decreases the requirements in terms of logical qubits. Since the qubit requirements have only a weak dependence on the 1-norm, this reduction is mainly due to the decreased number of terms in the Hamiltonian due to the lower Ξ rank of the SCDF. Similarly to the FeMoco active space model, by changing the kr parameter from its optimal value of 2 to 1, the number of logical qubits can be further decreased at the price of a higher number of Toffoli gates. In this case, for kr = 1, the number of logical qubits required by SCDF lies between the two THC results, as reported in Table 3.

Figure 2 Model systems of the Cpd I species in cytochrome P450 from ref (33) (C, black; O, red; H, white; S, yellow; N, blue; and Fe, orange).

Table 3 Resource Estimates, Correlation Energy Errors, and 1-Norm Values (λ) for the Active Space Model of Cytochrome P450a

approach	NDF or NTHC	Ξ	CCSD(T) error (mHa)	λ (Ha)	Toffoli gates	logical qubits	
XDF	4N	57	0.12	472.2	1.9 × 1010	4922	
XDF	5N	57	0.069	472.7	2.1 × 1010	4926	
XDF	6N	57	0.060	472.9	2.2 × 1010	4925	
 	 	 	 	 	 	 	
XDF + sym. shift	4N	57	0.12	298.9	1.2 × 1010	4920	
XDF + sym. shift	5N	57	0.066	299.4	1.3 × 1010	4924	
XDF + sym. shift	6N	57	0.057	299.6	1.4 × 1010	4923	
 	 	 	 	 	 	 	
RCDFb	100 ≈ 1.7 N	58	0.019	284.1	4.6 × 1010	18,856	
 	 	 	 	 	 	 	
THCc	320 ≈ 5.5 N	 	0.10	388.9	7.8 × 109	1434	
THCc	380 ≈ 6.6 N	 	–0.83	392.5	8.3 × 109	2158	
 	 	 	 	 	 	 	
SCDF	4N	33	–0.10	112.3	3.9 × 109	2590	
SCDF	5N	26	0.044	111.3	3.8×109	2596	
SCDF	6N	24	0.069	111.0	4.0 × 109	2594	
 	 	 	 	 	 	 	
SCDF (kr = 1)	5N	26	0.044	111.3	4.8 × 109	1706	
a Different approaches based on DF and THC are compared for different values of the rank used in the factorization. The best performing methods in terms of Toffoli gates and logical qubits are highlighted in bold. If not explicitly indicated, the SCDF results were obtained with the optimal value kr = 2.

b The resource estimation was obtained using the tensors provided with ref (29).

c Results from ref (33).

For the P450 application with NDF = 5N = 290, only 7 nonzero αt’s are included using the δα = 10–3 threshold. This truncation changes the 1-norm only by 0.3 Ha, corresponding to 0.3%. Without the symmetry shift contribution, the 1-norm of the SCDF approach would increase from 111.3 to 216.5. Similarly to the case of the FeMoco molecule, if the δα and δDF thresholds are set to zero, the 1-norm and correlation energy error is minimally affected, but the resource requirements sizeably increase to 6.7 × 109 Toffoli gates and 4924 logical qubits.

Since the RCDF tensor decomposition for the same P450 active space model was provided with ref (29), in Table 3, we also include the resource estimation for this methodology. To obtain the Toffoli gate and logical qubit requirements, we assume that the implementation of the complex components of the Wt tensors (see discussion about eq 22) does not involve any overhead costs with respect to the real case. As discussed in Section 2.2, the RCDF involves a factor N more terms in the Hamiltonian as compared to the XDF and, similarly, the SCDF. In order to decrease the number of terms, we used also in this case, the same δDF = 10–4 threshold, but this procedure is not effective in this case. Despite the sizable decrease of the 1-norm with respect to XDF, the cost of implementing is much higher for RCDF, and this explains the high computational requirements, as shown in Table 3.

3.4 Hydrogen Chains

In this section, we discuss the scaling of the SCDF method for systems of growing size by considering hydrogen chain models with up to 80 atoms. To compare with previous results in the literature, we use the same interatomic distance (d = 1.4 Bohr) of ref (21) and the STO-6G basis set. For this specific choice of the basis set, the number of spatial orbitals N is equal to the number of hydrogen atoms NH. For all the different DF-based approaches, NDF is set to 4N and the δDF = 10–4 threshold is used to truncate the components in the second factorization. Our results are compared with the THC results from ref (21), where the tensor decomposition rank was set to 7N.

First of all, for systems of growing size, it is of fundamental importance to establish the accuracy in the prediction of the ground-state energy. Table 4 shows the error in the CCSD(T) correlation energy per atom in Hartree. The SCDF method is characterized by errors of the order of 10–7 Ha per atom and a maximum total (absolute) deviation of 5.06 × 10–5 for H80. The XDF method, without and with symmetry shift, achieves an even higher level of accuracy. In ref (21) the authors aimed at achieving an accuracy in the energy per atom within 50–60 μHa. This choice leads to an energy error per atom for THC that is about one or 2 orders of magnitude larger than what we report here for SCDF. Importantly, these errors present a rather erratic behavior as a function of NH and, for H60, the total deviation in the CCSD(T) correlation energy even achieves the value of 3.1 × 10–3 Ha, which largely exceeds the chemical accuracy threshold (1.6 × 10–3 Ha). While this level of accuracy for THC might not be sufficient to obtain accurate energy differences or properties in the thermodynamic limit, we nevertheless consider this factorization for the purpose of comparing the resource requirements of THC and DF-based methods.

Table 4 CCSD(T) Correlation Energy Error Per Atom (in Ha) for Hydrogen Chains with up to 80 Atomsa

NH	XDF	XDF + sym. shift	THCb	SCDF	
10	7.2 × 10–10	–5.8 × 10–10	4.4 × 10–6	–2.1 × 10–7	
20	8.1 × 10–10	8.3 × 10–10	1.7 × 10–5	8.7 × 10–7	
30	–4.6 × 10–8	–3.4 × 10–8	–2.9 × 10–7	7.7 × 10–7	
40	1.1 × 10–7	–2.9 × 10–9	1.2 × 10–5	–2.8 × 10–7	
50	5.5 × 10–8	8.7 × 10–9	5.4 × 10–6	–2.0 × 10–7	
60	–6.8 × 10–9	–1.8 × 10–8	5.1 × 10–5	1.5 × 10–7	
70	–6.1 × 10–8	–2.1 × 10–8	1.9 × 10–5	–3.2 × 10–7	
80	–1.1 × 10–7	–3.63 × 10–10	4.9 × 10–6	–6.3 × 10–7	
a Four methods are compared: XDF, XDF with symmetry shift, THC, and SCDF.

b Results from ref (21).

The 1-norm of the Hamiltonian plays a fundamental role in determining the number of QPE iterations, and its scaling has strong implications for the efficiency of a specific method. Figure 3 and Table 5 show the empirical behavior of the 1-norm as a function of NH. The 1-norms of the standard XDF approach and its symmetry-shifted counterpart do not significantly differ for the hydrogen chains and are characterized by a asymptotic complexity. This has decreased considerably for SCDF, which has a similar scaling as THC and is expected to significantly decrease the computational requirements as compared to other DF-based methods. To this purpose, the left panel of Figure 4 and Table 5 show the dependence of the number of the Toffoli gates on the number of hydrogen atoms. The SCDF approach requires significantly less Toffoli gates than XDF (with or without symmetry shift) and improves over their computational complexity by approximately a factor N, as expected by the 1-norm behavior. Within this range, SCDF also outperforms THC, but with a tendency to grow more rapidly. Beyond the slightly different growth of λ, this should be expected as the theoretical scaling of the Toffoli gate requirements per iteration for the DF-based methods behaves as while for THC behaves as ; although the hydrogen chains considered here are still too small to reproduce exactly these asymptotic behaviors, the slopes in Table 5 already reflect the trends correctly.

Figure 3 Scaling of the 1-norm λ (in Ha) as a function of the number of hydrogen atoms NH for XDF, XDF with symmetry shift, THC, and SCDF.

Table 5 Slopes of the Linear Fits on a Logarithmic Scale of the 1-Norm λ, the Number of Toffoli Gates, and the Number of Logical Qubits As a Function of the Number of Hydrogen Atoms NH

approach	λ	number Toffolis	number logical qubits	
 	slope	R2	slope	R2	slope	R2	
XDF	1.87	0.9998	3.08	0.9997	1.27	0.9647	
XDF + sym. shift	1.98	0.9999	3.18	0.9998	1.27	0.9647	
THCa	1.11	0.9991	2.08	0.9998	1.01	0.9653	
SCDF	1.24	0.9998	2.38	0.9998	1.21	0.9249	
a Results obtained using the data provided with ref (21).

Figure 4 Left: scaling of the number of Toffoli gates required by the XDF, XDF with symmetry shift, THC, and SCDF approaches as a function of the number of hydrogens NH. Right: scaling of the number of logical qubits.

The behavior of the number of logical qubits as a function of the system size is less smooth, and the fitting on a logarithmic scale in Table 5 and the right panel of Figure 4 has to be considered only indicative of the overall trends. Since the number of logical qubits has only a weak dependence on the 1-norm, in this case, the behaviors of the SCDF do not significantly differ with respect to the XDF, and, as in the previous examples, THC has lower requirements than the DF methods.

Finally, it is interesting to discuss how the SCDF performs for hydrogen chains in a strongly correlated regime with large values of the interatomic distance. This is particularly important for the calculation of accurate binding curves. To this purpose, we considered H10 with d = 2.8 Bohr, which is double the value used for the hydrogen chains considered above. The XDF and SCDF were computed for this system with NDF = 4N. The SCDF approach for the stretched H10 (d = 2.8 Bohr) leads to a smaller value of the 1-norm, 5.72 Ha, to be compared to 12.75 Ha of the original H10 chain (d = 1.4 Bohr). Although still small and well below chemical accuracy, the error in CCSD(T) correlation energy per atom is −4.0 × 10–6 Ha, to be compared to −2.1 × 10–7 Ha in Table 4. Considering that the stretched system is strongly correlated but small enough to be solved exactly, it is interesting to determine the full configuration interaction (fullCI) error; the value obtained, −4.3 × 10–6 Ha, is in excellent agreement with the CCSD(T) estimate. The trend for XDF is similar, with an increase (in absolute value) of the CCSD(T) correlation energy error to −4.75 × 10–9 Ha as compared to 7.2 × 10–10 Ha in Table 4. This behavior suggests that, at least in certain cases, stretched geometries could require larger values of NDF to maintain a consistent level of accuracy with respect to configurations closer to equilibrium. Since the 1-norm is also smaller in the stretched regime, this should not involve a significant increase in the resource requirements.

4 Conclusions

In conclusion, we have introduced the SCDF approach that couples the symmetry shift technique with regularized DF to decrease the 1-norm of the electronic Hamiltonian significantly. As the 1-norm determines the total number of iterations needed in QPE, its reduction is important in decreasing the resource requirements and the runtime of fault-tolerant quantum simulations in chemistry. The effectiveness of the SCDF method in reducing the 1-norm is demonstrated numerically with applications to different chemical systems, including active space models of the FeMoco molecule, cytochrome P450, and hydrogen chains up to 80 atoms. For these systems, the 1-norm values achieved by SCDF are significantly lower than those achieved by other methods, including XDF, XDF with symmetry shift, RCDF, and THC.

Despite the fundamental role of the 1-norm, other factors should also be taken into account to assess the performance of a specific Hamiltonian factorization in the context of qubitization-based QPE.

First, the cost of a single QPE iteration, which depends on the specific implementation and the number of terms in the Hamiltonian, can significantly impact the global computational cost. To perform an unbiased comparison of methods, the Toffoli gate and logical qubit requirements were estimated for the chemical systems considered here. This shows that SCDF still outperforms all the other methods regarding Toffoli gate counts. For example, for the FeMoco molecule and cytochrome P450, SCDF requires fewer than 50% of the Toffoli gates compared to THC, thus far recognized as the best-performing method in the literature. However, the number of terms in the THC-factorized Hamiltonian grows with lower complexity as compared to all other methods, and this approach will tend to become the most efficient for systems of growing size. Concerning the logical qubit count, SCDF inherits the same behavior as the other DF-based approaches and tends to require relatively large numbers of qubits.

Second, given that Hamiltonian factorizations frequently employ tensor truncations to enhance efficiency, it is crucial to determine their impact on the ground-state total energy. We considered CCSD(T) correlation energies as an accuracy metric consistent with prior literature. Despite the potentially disruptive effect of the regularization, SCDF maintains a high level of accuracy comparable to other DF methods. The behavior of THC is more problematic. For example, in the case of cytochrome P450, the error in the correlation energy exhibits an erratic behavior as a function of the THC rank, with a tendency to significantly increase after reaching very small values. The behavior of the correlation energy is also nonsystematic as a function of the size of the hydrogen chains, and, in at least one case, the error significantly exceeds the chemical accuracy threshold.

In this work, the significant 1-norm reduction provided by SCDF has been exclusively exploited in the context of QPE. However, different quantum algorithms could benefit from the 1-norm decrease. This is the case, for example, of the stochastic compilation protocol known as qDRIFT,56 which requires O(λ2) repetitions to approximate the time evolution. The number of shots required to measure the expectation value of the Hamiltonian also grows with the 1-norm and, accordingly, also VQE and other near-term quantum algorithms could benefit from its reduction. Establishing the accuracy and efficiency of SCDF, also for near-term quantum algorithms, will be the subject of future work. In the context of fault-tolerant quantum computing, the ideas developed in this work may also apply within the THC framework. While not straightforward, this extension could lead to a highly efficient methodology with optimal scaling.

Appendix A

Optimization of the SCDF Cost Function

The SCDF cost function in eq 28 has to be minimized with respect to the components of the Wt, Ut, and αt tensors. As previously discussed for CDF,28 the simultaneous optimization of the different tensors is inefficient, and a nested approach with multiple steps should be preferred. For the SCDF cost function, the optimization is performed through the following steps:1. Optimization of the Wt tensor using the unconstrained minimization L-BFGS algorithm;

2. Update of the αt values as the median of the new Vt that is obtained from Wt using eq 14;

3. Optimization of the Ut tensor using the L-BFGS algorithm;

4. Repeat from step 1 until convergence is achieved.

For the optimizations in steps 1 and 3, we have developed two different implementations based on analytical gradients and automatic differentiation using the JAX library.48 While the former implementation is mainly intended to provide a reference and to test the soundness of the numerics, the latter is significantly more efficient and is used in production runs. The convergence criterion in step 4 is naturally defined in terms of thresholds on gradients or on the decrease of the cost function between subsequent iterations. However, since the 1-norm reduction is crucial to enhancing the QPE efficiency and its value is found to decrease monotonically as a function of the iteration count, the optimization procedure was stopped only when the 1-norm value was steadily decreasing well below the 0.05 Ha threshold.

The analytical gradient of the SCDF cost function with respect to the components of the Wt tensor is given by29

where30

The analytical gradient of the SCDF cost function with respect to the components of Ut is given by31

Since the regularization term does not depend on the Ut tensors, this derivative is the same as for the CDF and RCDF approaches. The optimization of Ut is more complex than in the Wt case since this tensor has to be constrained to be orthogonal. In practice, this problem can be formulated as an unconstrained optimization by introducing the antisymmetric orbital rotation generator matrices Xt to define Ut = exp(Xt). The minimization of the cost function is then performed with respect to the components of Xt; the procedure is described in detail in ref (28).

The minimization of the SCDF cost function is susceptible to the presence of local minima. Significantly lower values of the SCDF cost function and 1-norm can be found if very tight values (10–12) are used for the tolerance of the stopping criterion of the L-BFGS algorithm. This tight threshold was found to be effective for all the systems considered in this work at the price of a large number of steps to minimize Wt and Ut at each iteration. The presence of local minima affects the values of the 1-norm rather than the accuracy of the energy. For example, if the L-BFGS threshold is increased to 10–8 for P450 (NDF = 4N case), the 1-norm converges to 176.9 Ha, which is more than 50% higher than the value reported in Table 3. While this value of the 1-norm is still significantly smaller than the values of THC and other DF variants, this increase has a sizable impact on quantum resource requirements. However, the CCSD(T) error with this higher threshold is only −0.13 mHa, not significantly differing from the result in Table 3, −0.10 mHa. In a similar way, by increasing the L-BFGS threshold to 10–8 for FeMoco, the 1-norm converges to a larger value, 135.7 Ha (to be compared to 79.9 Ha in Table 2), but the accuracy in the CCSD(T) correlation energy does not significantly deteriorate, with an error of 0.73 mHa (to be compared to 0.60 mHa in Table 2). To find lower 1-norm minima and optimize the quantum resource requirements, all the results presented in this work were obtained with the 10–12 threshold for the L-BFGS minimization.

Appendix B

Quantum Implementation of the Symmetry-Compressed Double-Factorized Hamiltonian

The implementation of the (rank 1) double-factorized Hamiltonian in the form of eq 15 has been discussed in detail in ref (24). Specifically, this approach is valid for a positive definite/rank one Vt tensor, as in the XDF case. In this appendix, we show that the implementation of SCDF requires minimal modifications with respect to the original work of von Burg et al. Within the SCDF framework, eq 26 redefines the two-body part of the Hamiltonian by introducing a symmetry shift of each term in the sum over the index t, leading to the following form of eq 932

As shown by the numerical applications part in Section 3, the values of αt are actually different from zero only for a small number Nα of indexes t. In SCDF, Vt = Wt ⊗ Wt is rank 1 and the element-wise constant shift αt is equivalent to subtracting the rank 1 tensor – αt1 ⊗ 1 (here ⊗ denotes the outer product and 1 is a vector whose N components are all 1). This implies that Vt – αt has rank two and, by applying an eigenvalue decomposition, can be expressed as33

where Pt and Qt are 1D tensors. To keep the formalism consistent, in this appendix, we will define Pt = Wt and Qt = 0 for the t indexes corresponding to αt = 0 (or, more precisely, corresponding to the αt’s smaller than a given threshold δα). The application of the Jordan–Wigner transformation to the Hamiltonian in eq 32 leads to34

where the primed sum ∑′ indicates that only the Nα nonzero terms are actually included. Similarly to eq 15, a threshold δDF is applied to eliminate the small components of the Pt and Qt tensors; this leads to sums over k that run up to values Ξ(t) and Θ(t) that are less than or equal to N. This equation is analogous to the formulation of von Burg et al. in eq 15, with only Nα additional terms with a negative sign. The implementation of terms with a negative sign requires minimal modifications with respect to the implementation described in the Supporting Information of ref (24). The minus sign can be included when combining the two body terms evaluated from qubitization. Using the same notation of the Supporting Information of ref (24), this can be achieved by replacing with in the circuit in eqs 7 and 8 therein. The notation with an additional line on top is used in ref (24) to distinguish eqs 7 and 8, which correspond to the prepare operator for block encoding without or with sign, respectively. The quantum circuits for and are very similar and require the same number of Toffoli gates. Accordingly, the required resources are not affected by the negative sign, and the resource estimation was carried out using the standard OpenFermion implementation.49

Data Availability Statement

The SCDF tensor factorizations discussed in this work are available on a public Zenodo repository with https://doi.org/10.5281/zenodo.10836681.

The authors declare the following competing financial interest(s): D.R., C.L.C., J.F.G., P.J.O., and R.M.P. own stock/options in QC Ware Corp.

Acknowledgments

We thank Nicholas Rubin for explaining the Openfermion resource estimation code and the details of THC, William Poll and Mark Steudtner for insightful discussions on DF and THC circuits, and Oumarou Oumarou and Christian Gogolin for valuable discussions on the RCDF method.
==== Refs
References

Heifetz A. Quantum Mechanics in Drug Discovery; Springer, 2020.
Nørskov J. K. ; Bligaard T. ; Rossmeisl J. ; Christensen C. H. Towards the computational design of solid catalysts. Nat. Chem. 2009, 1 , 37–46. 10.1038/nchem.121.21378799
Jain A. ; Ong S. P. ; Hautier G. ; Chen W. ; Richards W. D. ; Dacek S. ; Cholia S. ; Gunter D. ; Skinner D. ; Ceder G. ; Persson K. A. Commentary: The Materials Project: A materials genome approach to accelerating materials innovation. APL Mater. 2013, 1 , 011002 10.1063/1.4812323.
Cao Y. ; Romero J. ; Olson J. P. ; Degroote M. ; Johnson P. D. ; Kieferová M. ; Kivlichan I. D. ; Menke T. ; Peropadre B. ; Sawaya N. P. D. ; Sim S. ; Veis L. ; Aspuru-Guzik A. Quantum Chemistry in the Age of Quantum Computing. Chem. Rev. 2019, 119 , 10856–10915. 10.1021/acs.chemrev.8b00803.31469277
Peruzzo A. ; McClean J. ; Shadbolt P. ; Yung M.-H. ; Zhou X.-Q. ; Love P. J. ; Aspuru-Guzik A. ; O’brien J. L. A variational eigenvalue solver on a photonic quantum processor. Nat. Commun. 2014, 5 , 4213 10.1038/ncomms5213.25055053
Tilly J. ; Chen H. ; Cao S. ; Picozzi D. ; Setia K. ; Li Y. ; Grant E. ; Wossnig L. ; Rungger I. ; Booth G. H. ; Tennyson J. The variational quantum eigensolver: a review of methods and best practices. Phys. Rep. 2022, 986 , 1–128. 10.1016/j.physrep.2022.08.003.
Suppressing quantum errors by scaling a surface code logical qubit. Nature 2023, 614 , 676–681. 10.1038/s41586-022-05434-1.36813892
Bluvstein D. ; Evered S. J. ; Geim A. A. ; Li S. H. ; Zhou H. ; Manovitz T. ; Ebadi S. ; Cain M. ; Kalinowski M. ; Hangleiter D. ; Bonilla Ataides J. P. ; et al. Logical quantum processor based on reconfigurable atom arrays. Nature 2023, 626 , 58–65. 10.1038/s41586-023-06927-3.38056497
https://newsroom.ibm.com/2023-12-04-IBM-Debuts-Next-Generation-Quantum-Processor-IBM-Quantum-System-Two,-Extends-Roadmap-to-Advance-Era-of-Quantum-Utility (accessed on May 1, 2024).
https://quantumai.google/qecmilestone (accessed on May 1, 2024).
Abrams D. S. ; Lloyd S. Quantum algorithm providing exponential speed increase for finding eigenvalues and eigenvectors. Phys. Rev. Lett. 1999, 83 , 5162–5165. 10.1103/PhysRevLett.83.5162.
Aspuru-Guzik A. ; Dutoi A. D. ; Love P. J. ; Head-Gordon M. Simulated quantum computation of molecular energies. Science 2005, 309 , 1704–1707. 10.1126/science.1113479.16151006
Suzuki M. Improved Trotter-like formula. Phys. Lett. A 1993, 180 , 232–234. 10.1016/0375-9601(93)90701-Z.
Berry D. W. ; Childs A. M. ; Cleve R. ; Kothari R. ; Somma R. D. Simulating Hamiltonian dynamics with a truncated Taylor series. Phys. Rev. Lett. 2015, 114 , 090502 10.1103/PhysRevLett.114.090502.25793789
Berry D. W. ; Kieferová M. ; Scherer A. ; Sanders Y. R. ; Low G. H. ; Wiebe N. ; Gidney C. ; Babbush R. Improved techniques for preparing eigenstates of fermionic Hamiltonians. Npj Quantum Inf. 2018, 4 , 22 10.1038/s41534-018-0071-5.
Poulin D. ; Kitaev A. ; Steiger D. S. ; Hastings M. B. ; Troyer M. Quantum algorithm for spectral measurement with a lower gate count. Phys. Rev. Lett. 2018, 121 , 010501 10.1103/PhysRevLett.121.010501.30028152
Babbush R. ; Gidney C. ; Berry D. W. ; Wiebe N. ; McClean J. ; Paler A. ; Fowler A. ; Neven H. Encoding electronic spectra in quantum circuits with linear T complexity. Phys. Rev. X 2018, 8 , 041015 10.1103/physrevx.8.041015.
Low G. H. ; Chuang I. L. Hamiltonian simulation by qubitization. Quantum 2019, 3 , 163 10.22331/q-2019-07-12-163.
Bravyi S. ; Kitaev A. Universal quantum computation with ideal Clifford gates and noisy ancillas. Phys. Rev. A 2005, 71 , 022316 10.1103/PhysRevA.71.022316.
Reichardt B. W. Quantum universality from magic states distillation applied to CSS codes. Quantum Inf. Process. 2005, 4 , 251–264. 10.1007/s11128-005-7654-8.
Lee J. ; Berry D. W. ; Gidney C. ; Huggins W. J. ; McClean J. R. ; Wiebe N. ; Babbush R. Even more efficient quantum computations of chemistry through tensor hypercontraction. PRX Quantum 2021, 2 , 030305 10.1103/PRXQuantum.2.030305.
Low G. H. ; Kliuchnikov V. ; Schaeffer L. Trading T-gates for dirty qubits in state preparation and unitary synthesis. arXiv 2018, arXiv:1812.00954 10.48550/arXiv.1812.00954.
Berry D. W. ; Gidney C. ; Motta M. ; McClean J. R. ; Babbush R. Qubitization of arbitrary basis quantum chemistry leveraging sparsity and low rank factorization. Quantum 2019, 3 , 208 10.22331/q-2019-12-02-208.
von Burg V. ; Low G. H. ; Häner T. ; Steiger D. S. ; Reiher M. ; Roetteler M. ; Troyer M. Quantum computing enhanced computational catalysis. Phys. Rev. Res. 2021, 3 , 033055 10.1103/PhysRevResearch.3.033055.
Motta M. ; Ye E. ; McClean J. R. ; Li Z. ; Minnich A. J. ; Babbush R. ; Chan G. K.-L. Low rank representations for quantum simulation of electronic structure. Npj Quantum Inf. 2021, 7 , 83 10.1038/s41534-021-00416-z.
Kivlichan I. D. ; McClean J. ; Wiebe N. ; Gidney C. ; Aspuru-Guzik A. ; Chan G. K.-L. ; Babbush R. Quantum simulation of electronic structure with linear depth and connectivity. Phys. Rev. Lett. 2018, 120 , 110501 10.1103/PhysRevLett.120.110501.29601758
Huggins W. J. ; McClean J. R. ; Rubin N. C. ; Jiang Z. ; Wiebe N. ; Whaley K. B. ; Babbush R. Efficient and noise resilient measurements for quantum chemistry on near-term quantum computers. Npj Quantum Inf. 2021, 7 , 23 10.1038/s41534-020-00341-7.
Cohn J. ; Motta M. ; Parrish R. M. Quantum filter diagonalization with compressed double-factorized hamiltonians. PRX Quantum 2021, 2 , 040352 10.1103/PRXQuantum.2.040352.
Oumarou O. ; Scheurer M. ; Parrish R. M. ; Hohenstein E. G. ; Gogolin C. Accelerating Quantum Computations of Chemistry Through Regularized Compressed Double Factorization. arXiv 2022, arXiv:2212.07957 10.48550/arXiv.2212.07957.
Rubin N. C. ; Lee J. ; Babbush R. Compressing many-body fermion operators under unitary constraints. J. Chem. Theory Comput. 2022, 18 , 1480–1488. 10.1021/acs.jctc.1c00912.35166529
Hohenstein E. G. ; Parrish R. M. ; Martínez T. J. Tensor hypercontraction density fitting. I. Quartic scaling second-and third-order Møller-Plesset perturbation theory. J. Chem. Phys. 2012, 137 , 044103 10.1063/1.4732310.22852593
Parrish R. M. ; Hohenstein E. G. ; Martínez T. J. ; Sherrill C. D. Tensor hypercontraction. II. Least-squares renormalization. J. Chem. Phys. 2012, 137 , 224106 10.1063/1.4768233.23248986
Goings J. J. ; White A. ; Lee J. ; Tautermann C. S. ; Degroote M. ; Gidney C. ; Shiozaki T. ; Babbush R. ; Rubin N. C. Reliably assessing the electronic structure of cytochrome p450 on today’s classical computers and tomorrow’s quantum computers. Proc. Natl. Acad. Sci. U.S.A. 2022, 119 , e2203533119 10.1073/pnas.2203533119.36095200
Koridon E. ; Yalouz S. ; Senjean B. ; Buda F. ; O’Brien T. E. ; Visscher L. Orbital transformations to reduce the 1-norm of the electronic structure Hamiltonian for quantum computing applications. Phys. Rev. Res. 2021, 3 , 033127 10.1103/PhysRevResearch.3.033127.
Loaiza I. ; Khah A. M. ; Wiebe N. ; Izmaylov A. F. Reducing molecular electronic hamiltonian simulation cost for linear combination of unitaries approaches. Quantum Sci. Technol. 2023, 8 , 035019 10.1088/2058-9565/acd577.
Jordan P. ; Wigner E. On Paul’s prohibition on equivalence. Z. Phys. 1928, 47 , 631–651. 10.1007/BF01331938.
Bravyi S. B. ; Kitaev A. Y. Fermionic quantum computation. Ann. Phys. 2002, 298 , 210–226. 10.1006/aphy.2002.6254.
Seeley J. T. ; Richard M. J. ; Love P. J. The Bravyi-Kitaev transformation for quantum computation of electronic structure. J. Chem. Phys. 2012, 137 , 224109 10.1063/1.4768229.23248989
Poulin D. ; Hastings M. B. ; Wecker D. ; Wiebe N. ; Doherty A. C. ; Troyer M. The Trotter step size required for accurate quantum simulation of quantum chemistry. arXiv 2014, arXiv:1406.4920 10.48550/arXiv.1406.4920.
Matsuzawa Y. ; Kurashige Y. Jastrow-type decomposition in quantum chemistry for low-depth quantum circuits. J. Chem. Theory Comput. 2020, 16 , 944–952. 10.1021/acs.jctc.9b00963.31939668
Thouless D. J. Stability conditions and nuclear rotations in the Hartree-Fock theory. Nucl. Phys. 1960, 21 , 225–232. 10.1016/0029-5582(60)90048-1.
Childs A. M. ; Wiebe N. Hamiltonian simulation using linear combinations of unitary operations. Quant.Inf.Comput. 2012, 12 , 901–924. 10.26421/qic12.11-12-1.
Dunlap B. I. ; Connolly J. ; Sabin J. On some approximations in applications of X α theory. J. Chem. Phys. 1979, 71 , 3396–3402. 10.1063/1.438728.
Werner H.-J. ; Manby F. R. ; Knowles P. J. Fast linear scaling second-order Møller-Plesset perturbation theory (MP2) using local and density fitting approximations. J. Chem. Phys. 2003, 118 , 8149–8160. 10.1063/1.1564816.
Pedersen T. B. ; Aquilante F. ; Lindh R. Density fitting with auxiliary basis sets from Cholesky decompositions. Theor. Chem. Acc. 2009, 124 , 1–10. 10.1007/s00214-009-0608-y.
Hohenstein E. G. ; Sherrill C. D. Density fitting of intramonomer correlation effects in symmetry-adapted perturbation theory. J. Chem. Phys. 2010, 133 , 014101 10.1063/1.3451077.20614953
Loaiza I. ; Izmaylov A. F. Block-Invariant Symmetry Shift: Preprocessing Technique for Second-Quantized Hamiltonians to Improve Their Decompositions to Linear Combination of Unitaries. J. Chem. Theory Comput. 2023, 19 , 8201–8209. 10.1021/acs.jctc.3c00912.37939198
Bradbury J. ; Frostig R. ; Hawkins P. ; Johnson M. J. ; Leary C. ; Maclaurin D. ; Necula G. ; Paszke A. ; VanderPlas J. ; Wanderman-Milne S. ; Zhang Q. JAX: composable transformations of Python+NumPy programs Version 0.3.13, 2018. http://github.com/google/jax (accessed on May 1, 2024).
McClean J. R. ; Rubin N. C. ; Sung K. J. ; Kivlichan I. D. ; Bonet-Monroig X. ; Cao Y. ; Dai C. ; Fried E. S. ; Gidney C. ; Gimby B. ; Gokhale P. ; et al. OpenFermion: the electronic structure package for quantum computers. Quantum Sci. Technol. 2020, 5 , 034014 10.1088/2058-9565/ab8ebc.
Beinert H. ; Holm R. H. ; Münck E. Iron-sulfur clusters: nature’s modular, multipurpose structures. Science 1997, 277 , 653–659. 10.1126/science.277.5326.653.9235882
Reiher M. ; Wiebe N. ; Svore K. M. ; Wecker D. ; Troyer M. Elucidating reaction mechanisms on quantum computers. Proc. Natl. Acad. Sci. U.S.A. 2017, 114 , 7555–7560. 10.1073/pnas.1619152114.28674011
Li Z. ; Li J. ; Dattani N. S. ; Umrigar C. ; Chan G. K. The electronic complexity of the ground-state of the FeMo cofactor of nitrogenase as relevant to quantum simulations. J. Chem. Phys. 2019, 150 , 024302 10.1063/1.5063376.30646701
Stephens P. J. ; Devlin F. J. ; Chabalowski C. F. ; Frisch M. J. Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields. J. Phys. Chem. 1994, 98 , 11623–11627. 10.1021/j100096a001.
Weigend F. ; Ahlrichs R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7 , 3297–3305. 10.1039/b508541a.16240044
Pipek J. ; Mezey P. G. A fast intrinsic localization procedure applicable for ab initio and semiempirical linear combination of atomic orbital wave functions. J. Chem. Phys. 1989, 90 , 4916–4926. 10.1063/1.456588.
Campbell E. Random compiler for fast Hamiltonian simulation. Phys. Rev. Lett. 2019, 123 , 070503 10.1103/PhysRevLett.123.070503.31491106
