Increasing the Representation Accuracy of Quantum Simulations of Chemistry without Extra Quantum Resources

Proposals for experiments in quantum chemistry on quantum computers leverage the ability to target a subset of degrees of freedom containing the essential quantum behavior, sometimes called the active space. This approximation allows one to treat more difficult problems using fewer qubits and lower gate depths than would otherwise be possible. However, while this approximation captures many important qualitative features, it may leave the results wanting in terms of absolute accuracy (basis error) of the representation. In traditional approaches, increasing this accuracy requires increasing the number of qubits and an appropriate increase in circuit depth as well. Here we explore two techniques requiring no additional qubits or circuit depth that are able to remove much of this approximation in favor of additional measurements. The techniques are constructed and analyzed theoretically, and some numerical proof-of-concept calculations are shown. As an example, we show how to achieve the accuracy of a 20-qubit representation using only four qubits and a modest number of additional measurements for a hydrogen molecule. We close with an outlook on the impact such techniques may have on both near-term and fault-tolerant quantum simulations.


I. INTRODUCTION
Quantum computers promise to make a dramatic impact on a number of fields including optimization and materials simulation.Since the initial proposal by Feynman to simulate quantum systems via quantum systems [1], there has been a wealth of developments studying these applications both theoretically and experimentally.One particular application of note is the simulation of chemical and material systems by quantum computers, as it represents a natural application of this idea in practice.The combination of the practical potential for quantum chemistry as well as its low overhead have made it a target for the first near-term quantum computers.
The progress in quantum chemistry on quantum computers has been extremely rapid over the last few years.Starting from the original proposal to use quantum phase estimation for chemical problems [2], the precise costs of these algorithms and methods to reduce these costs by orders of magnitude have now been developed in detail [3][4][5][6][7][8][9][10][11][12][13][14][15].Precise gate counts are now known for hard chemical systems in implementations on the earliest fault-tolerant computers under a realistic model of error correction using the surface code [16][17][18][19][20]. Methods based on quantum phase estimation, however, are believed by some to require quantum error correction, which places their experimental implementation some time beyond the current noisy intermediate-scale quantum (NISQ) devices.
As NISQ devices progress to a level of development capable of doing computations a classical computer cannot, or quantum supremacy [21], the question arises if one can do a practical application such as chemistry on such a near-term quantum computer.Prime candidates for this possibility have been variational algorithms, such as the variational quantum eigensolver (VQE) [22,23] or the quantum approximation optimization algorithm [24].These methods exhibit a natural adaptation to device parameters as well as intrinsic robustness to systematic errors that make them attractive candidates.Since their inception, they have been extended to treat excited states [25][26][27] and different problem areas [28,29] and have been demonstrated on numerous experimental architectures [22,26,[30][31][32][33].
While hardware and theoretical developments have been rapid, quantum resources are expected to remain costly for some time.As a result, proposals for doing quantum computing for chemical problems have focused on isolating the essential strongly correlated component of a physical system for simulation on a quantum computer.A simple division of the system into this select subsystem is often called an active space within chemistry, and more detailed treatments may incorporate it as an impurity model as well.Despite the ability of these methods to treat qualitative phenomena that would not otherwise be accessible, this division introduces quantitative approximations that can be unacceptable when quantitative accuracy is demanded.Previously, when one wanted to lift this approximation, the consequence was both an increased number of qubits and gate complexity, ruling out compatibility with a NISQ device, and unnecessary qubit overhead in fault-tolerant applications.
Here, we discuss two methods that use no additional qubits or gate complexity to lift the active-space approximation in a systematic way.These methods are expected to impact both near-term applications as well as fault-tolerant applications where the number of qubits may remain a limitation in reaching the desired accuracy.One of these procedures is a novel way to leverage the quantum subspace expansions [25], also known to mitigate errors [34] and provide excited states in experiments [31], that relies only on additional measurements.The other involves orbital relaxations, a common and well-known procedure in classical electronic structure, to remove the active-space approximation and reduce the circuit depth for some classes of VQE circuits.We cost out these methods, provide theoretical justification, and show how these methods work in practice with simple numerical simulations.A framework to develop cost-effective approximations is discussed that makes use of manipulation of marginals [35], and we conclude with an outlook on the impact for near-term quantum simulations.

II. BACKGROUND AND DEFINITIONS
We begin by setting up the problem, establishing notation, and reviewing the quantum subspace expansion method [25] as it was originally devised.The problem of electronic structure in quantum chemistry is typically cast as the problem of determining the electronic ground state in the field of fixed nuclei, or the Born-Oppenheimer approximation [36].From here, one discretizes space, where the accuracy of this discretization determines the ultimate accuracy one can achieve.The abstract blocks one uses to divide space are called basis functions, and a number of canonical choices are known such as linear combinations of atomic orbitals and plane waves.Once the basis is chosen, the electronic structure Hamiltonian may be written in its canonical form as where each index corresponds to one of these basis functions, the h ij and h ijkl are standard integrals over the involved basis functions, and the ladder operators satisfy the canonical fermionic anticommutation relations fa † i ; a j g ¼ δ ij , fa i ; a j g ¼ fa † i ; a † j g ¼ 0. As we discuss above, the number of basis functions used is tied to the ultimate accuracy one can achieve for the problem; however, using too many basis functions may make the problem impractical or be a waste of resources when more clever treatments may be used.A method that was first widely used in traditional quantum chemistry and that has been widely adopted by the quantum-computing community is the active-space approximation.
The physical intuition behind the active-space approximation is that the space may be divided into a portion which exhibits strong correlations or entanglement, the active space, and a portion that while important, exhibits only low-rank contributions that are extremely well treated perturbatively, the virtual space.This methodology has been proven out numerous times in the classical literature by methods such as multireference configuration interaction and perturbation theory [37,38].However, the size of the essential quantum component or active space remains limited on a classical computer, and to date the contributions of virtuals have remained absent on a quantum computer.We aim to show how one can regain virtuals on top of quantum active spaces without the need for additional qubits or gate complexity.
As depicted in Fig. 1, in typical chemistry calculations on quantum computers, the set of basis functions is divided FIG. 1. Schematic virtual quantum subspace expansion to increase accuracy without additional qubits.The method separates the orbitals of a fermionic system into their core, active, and virtual components.The quantum device solves the active-space problem, and additional quantum measurements are taken within this space.These additional measurements are combined with classical postprocessing from data on the virtual space to correct the solutions using no additional qubits or circuit depth.The resulting corrected wave function(s) are stored in a mixed quantum-classical representation that may be used to derive any desired properties of these wave functions.into three sets, which we denote C, A, and V for core, active, and virtual.The core orbitals are assumed to be doubly occupied, and their contributions are integrated out to an effective field felt by the active space and virtual space.The virtual orbitals are ignored, and the problem is solved exactly within the dressed active space, A.
Our approach utilizes some of the machinery of an approach designed originally to provide excited states and error mitigation within the context of VQE, which utilizes quantum subspace expansions (QSEs) [25].These techniques are also applicable in fault-tolerant approaches and leverage the ability to expand and manipulate representations of operators within a subspace using measurements and classical computation without knowing the details of the states themselves.In this framework, one assumes the ability to prepare a wave function within the active space that we denote as jΨ ref i.We then choose a set of expansion operators fO i g which act on this reference in combination with another operator, such as the Hamiltonian, to form a representation of that operator in the basis given by fO i jΨ ref ig.As this basis may be nonorthogonal, we also measure the overlap or metric matrix within this basis in order to ensure the problem is well defined.The matrices may be formed through additional measurements only and have matrix elements given by With these matrices, one then uses canonical diagonalization to remove the approximately 0 eigenvalues of S and solve the generalized eigenvalue problem in the wellconditioned subspace given by where C is the matrix of eigenvectors in the basis fO i jΨ ref ig, and E is the diagonal matrix of eigenvalues.This approach can both improve the accuracy of the ground state and provide approximations to excited states.In the original work [25], it was suggested to use O i approximating fermionic excitations of the form which when considered with the Hamiltonian composed of only up to two-particle operators, means that the matrix elements can be evaluated as sums over subsets of the fourelectron reduced density matrix (RDM) where the sum weights are determined by integrals in the Hamiltonian.The four-electron density matrix in the active space is Each of these elements can be evaluated on a quantum computer by repeated preparation of jΨ ref i and measurement of the Pauli operators corresponding to the transformation of this matrix element by a Jordan-Wigner or equivalent transformation.To perform the procedure exactly, the number of terms one must measure scales as N 8 A where N A is the number of active-space orbitals.

III. VIRTUAL QUANTUM SUBSPACE EXPANSION
The original formulation of the QSE remained in the active space and did not include contributions from the virtual orbitals.Here we show how simple classical postprocessing and measurements can be used to reintroduce contributions from the virtuals in a systematic way to approach quantitative accuracy for the chemical problem.
Consider a reference function which is constructed only on the original active space A. We now introduce a set of expansion operators that act, in principle, on additional qubits that would define the virtual space.However, by definition, the reference wave function has no components in the virtual space, and hence, the contraction over the virtual space can be done using only efficient classical computation on the Hamiltonian in addition to the information from the active-space RDM.The only information that is required is the integrals over the full space, which are classical inputs to the problem, as well as the appropriate density matrix elements.We note that the quantum advantage over classical multireference configuration interaction (MRCI) methods here is that the reference wave function may contain an exponential number of determinants, which gives classical MRCI methods exponential scaling in the size of the active space.In contrast, our method scales only polynomially in the size of the active space, and if one measures the appropriate RDM of the active space, the quantumcomputing resources required are independent of the size of the virtual space.For both classical MRCI and the method proposed here, the classical computation required scales polynomially in the size of the virtual space.Thus, our method allows for very large basis sets to be treated on a small quantum computer.
In this set of operators, we restrict ourselves to single and double excitations from the references within the active space and virtuals.The single and double excitations have empirically been shown to be the dominant contributions in classical multireference methods.Our method is general enough to go to higher-order approximations; however, with each higher order of excitation, a larger reduced density matrix must be measured.If we denote single and doubles as excitation level 2, at excitation level k, one must measure the (2 þ k) RDM, which is expected to have a cost that scales as N 2ð2þkÞ A in the number of terms that must be measured.When the excitation level k matches the number of electrons, the method is formally exact.
To illustrate how we eliminate the virtual space, we take the highest-order matrix element as an example.The relevant quantity to measure for the case where both expansion operators in Eq. ( 2) are double excitations is Using Wick's theorem where greek indices denote virtual orbitals, we contract the operators on the virtual space, where Δ ξη;μν ¼ δ ξν δ ημ − δ ξμ δ ην , x (ȳ) means x (y) changing to the other value in the sum (e.g., ī ¼ j for x), and hÁ Á Ái denotes the expectation values with respect to the reference state jΨ ref i in the active space.The reduction in Eq. ( 9) shows that the expectation value M can be derived by knowing the fourelectron RDM of the reference state jΨ ref i in the active space.Similar but simpler results hold by replacing one or both double-excitation operators by single-excitation operators.
Hence, the inclusion of virtual orbitals amounts to simple classical postprocessing, and no additional qubits after the appropriate measurements are performed.

IV. CUMULANT AND RESTRICTED ACTIVE-SPACE APPROXIMATIONS
One may introduce a number of approximations to make the construction more efficient.The simplest approximation is the division of the active space into a part which is excited into the virtuals A v and a part of the active space which will be treated as correlated core orbitals A c .This partitioning reduces the scaling in the number of measurements to the cost of originally estimating the energy plus the size of the virtuals active space A v , N 8  Av , which for small sizes can dramatically reduce the cost.
The second class of approximations involves estimating the matrix elements of the reduced density matrix via cumulant approximations.For example, one can form a series of approximations to the four-electron RDM using products of lower RDMs and perturbative corrections.The 4-RDM may be expressed in terms of cumulants as may be denoted by where ðlÞ Δ is the l-particle cumulant expressing irreducible l-body correlations in the density matrix, and ∧ is the Grassmann wedge product.The simplest form of approximation is given by setting ðlÞ Δ ¼ 0 for l > 2 which allows one to express the 4-RDM with only measurements from the original 2-RDM.This truncation greatly reduces the number of terms to measure back to N 4 A , but it introduces some considerable approximation to the energetic values.Much work has been done on improving these approximations as well.An alternative scheme we do not exploit here may stochastically sample the elements to measure with hopefully increasing degrees of accuracy as time proceeds within a calculation.

V. NUMERICAL EXPERIMENTS
In this section, we show on a prototype system the accuracy gains one may expect from using this method on a real system.We assess the performance of the procedure on a simple molecule and examine the ground-state energy as a function of the number of additional virtual orbitals considered in the system.The first system we look at is H 2 in a variety of different basis sets.These numerics are enabled by the OpenFermion software package for doing quantum chemistry on a quantum computer [39] in conjunction with the Psi4 [40] open source package for electronic structure.The orbitals are obtained from a Hartree-Fock calculation, and the reference function jΨ ref i is the exact solution within the active space consisting of FIG. 2. Ground-state energy as a function of the internuclear separation for the H 2 molecule at different levels of theory.The notation ðn ¼ n q Þ after each line indicates that n q qubits are needed for that representation.The use of the virtual quantum subspace expansion (VQSE) technique using only four qubits attains the same accuracy in this case as a representation using 20 qubits, with an error much smaller than chemical accuracy (approximately 10 −5 E h ).The minimal basis here denotes a STO-3G calculation and 6-31G is an intermediate level of accuracy between the minimal basis and a double-zeta basis cc-pVDZ (here, VDZ). the single occupied orbital and the lowest-energy virtual orbital.Only single and double excitations from the active space to the additional virtual orbitals are included in these calculations.A correlation-consistent basis set of doublezeta quality (cc-pVDZ) [41] is used, as well as the 6-31G basis set [42] for reference.We compare with the exact results with both basis sets.
The results of the numerical simulations are shown in Fig. 2, and the energy differences with respect to exact/cc-pVDZ are shown in Fig. 3.We see that the addition of virtuals to an active space at a higher level of theory quickly surpasses what is possible at a lower level of theory, using no additional qubits.In fact, using only four qubits, which is the same as a typical minimal basis for H 2 , the VQSE procedure attains an accuracy commensurate with the exact solution on a basis that would require 20 qubits.The curves smoothly improve as a function of the number of virtuals included and show excellent accuracy across the range of calculations.

VI. FULL-SPACE ORBITAL RELAXATION
A common classical electronic structure method to improve active-space calculations is to allow for orbitals to relax in the presence of the newly optimized active-space wave function.Variants of this idea fall into the family of methods known as the multiconfigurational self-consistentfield (MCSCF) method.These algorithms work by iteratively solving the active-space Schrödinger equation and then finding a single-particle rotation on the full space that minimizes the energy.Integrating a quantum resource as an active-space solver in the MCSCF framework has been previously suggested for use with the phase estimation approach to quantum simulation [16] and also motivated other embedding methods in the quantum-computing context [43].Here we demonstrate how a single step of a MCSCF orbital relaxation can be used to further improve energies with a procedure that is compatible with the limitations of NISQ simulations and also readily applied to fault-tolerant implementations where qubit limitations are a concern.This technique requires measuring only the two-electron reduced density matrix with no additional quantum resources.
Given a 2-RDM from the ground state of the dressed active-space Hamiltonian, we seek to minimize the energy of the full-space Hamiltonian Eq. ( 1) through one-body rotations U on the full space.Because of Thouless's theorem, the one-body rotations can be efficiently implemented as a rotation of the underlying basis [44].This can be formulated as the following nonlinear optimization problem: X ¼ X p;q t p;q a † p a q ; t p;q ¼ −t Ã q;p ; where ha † i a † j a k a l i A and ha † i a j i A are the 2-and 1-RDM of the ground-state active-space wave function.
In the classical literature, it is common to use the secondorder approximation of u, to minimize energy in Newton-Raphson style [45][46][47][48].
A well-known alternative to parametrizing the unitary as an exponentiated anti-Hermitian matrix is the use of Givens rotations [49,50].This parametrization uses a set of angles fθg associated with the set of nonredundant orbital rotation generators.For 2-RDMs obtained from an exact diagonalization of the active-space Hamiltonian, the only nonredundant parameters are one-body generators associated with pairs of orbitals involving rotations from the active space to the virtual space and active space to the core space.Therefore, the unitary in Eq. ( 11) can be expressed as a product of Givens matrices FIG. 3. Ground-state energy error with respect to the exact solution in a cc-pVDZ (abbreviated here as VDZ) as a function of internuclear separation for the H 2 molecule for different levels of theory.The notation ðn ¼ n q Þ after each line indicates that n q qubits are needed for that representation.We see that the VQSE technique using only four qubits attains the same accuracy as an exact solution in a 20-qubit basis, up to an error of 10 −5 E h which is far below chemical accuracy.
The optimal rotation of a single angle with respect to the input 2-RDM, one-, and two-electron integrals has been derived [50] along with a sweep procedure to find an energy minimizing U.This algorithm was used for developing 2-RDM MCSCF methods [51] and the first orbitaloptimized coupled-cluster doubles methods [52].
One consideration in configuration-interaction-type approaches is the lack of a property known as size extensivity, which is related to the proper scaling of energy as a function of system size [36].While coupled-cluster approaches are designed to avoid this flaw, their extra complexity is often cumbersome in multireference schemes, and for modest system sizes, the effects may not be observable.Indeed, many modern approaches used in traditional computation, such as density matrix renormalization group, tensor network methods, or select configuration interaction approaches are not strictly size extensive, but so long as the accuracy attained is sufficient, this is not a concern.We note that for our orbital-relaxation technique, if the active-space solution is exact, size extensivity is satisfied exactly due to the exponential form of the ansatz correction.A related consideration worthy of further investigation is that of the size consistency of the methods used [36] when applied on quantum devices; however, we do not consider this in more detail here.While the configuration interaction form of VQSE is not strictly size extensive, for the types of systems likely to be investigated on quantum computers, this is not expected to be a concern.
One can iterate between solving the active-space Schrödinger equation and full-space one-body rotations, similar to a general MCSCF routine, or perform a singleorbital-relaxation procedure as postprocessing.Because of the variational principle, relaxing the orbitals in a single step as postprocessing is guaranteed to reduced the energy.Furthermore, if the ground state of the active space is not achieved due to an approximate wave-function ansatz, additional one-body orbital rotations between orbital pairs inside the active space can be included in the relaxation.Including these rotations would correspond to an additional linear depth circuit that would have been executed perfectly on the quantum computer [53].Any circuit that contains one-body rotations at the end can be made shorter with classical postprocessing by instead rotating the operator to be measured.
To demonstrate the utility of the orbital-relaxation step, we examine the energy improvement upon orbital relaxation for molecular hydrogen computed in the 6-31G basis set.We select the lowest-energy Hartree-Fock orbitals as the active space and perform orbital optimization with the SciPy COBYLA optimizer.Though this is not guaranteed to find the lowest-energy orbital rotations, it is sufficient for validating the performance one can achieve through orbital relaxation as postprocessing.In Fig. 4, we plot the energy of molecular hydrogen and observe a significant recovery of energy when performing a four-qubit active-space calculation with orbital relaxation on the full eight-qubit space.The six-qubit active-space calculation is included to demonstrate a smooth interpolation of this method to exact diagonalization in the full eight-qubit space.Figure 5 shows the energy difference relative to an eight-qubit exact diagonalization, revealing we reach within 5 × 10 −5 E h of the exact solution with half the number of qubits.
As an additional demonstration of our technique, we consider the removal of a hydrogen atom from a water molecule in a cc-pVTZ basis in Fig. 6.This molecule has been considered in the context of quantum computation before but only in minimal basis sets.We begin with same the number of qubits as the minimal basis (n ¼ 14) as an active space and attain a solution that is improved beyond the minimal basis set by over a hartree, which is several orders of magnitude larger than chemical accuracy.Using the original basis set without this methodology requires a staggering 72 qubits, which is a demonstration of why the active-space method is so important.Using the OO approach, we go beyond the minimal basis accuracy without requiring additional quantum resources and match the accuracy of a solution with six additional qubits in the active space in the bonding region.Moreover, the inclusion of orbital relaxation improves the description of the process considerably in the dissociation region, which has a critical influence on the estimate of the dissociation energy.Hence, our method improves the accuracy while still requiring FIG. 4. Ground-state energy as a function of internuclear separation for the H 2 molecule at different levels of theory.The notation ðn ¼ n q Þ after each line indicates that n q qubits are needed for that representation.The orbital optimization postprocessing is labeled by "OO."We compare active-space calculations on the full space with the 6-31G basis set (requiring eight qubits) to active-space calculations on the same model involving only four and six qubits.fewer qubits.We note that this curve is not yet fully converged with respect to the basis set limit, but our technique can be applied at larger sizes as well.

VII. DISCUSSION
As devices progress, there will continue to be a push to demonstrate their power for interesting chemical problems.As the number of qubits and coherence time are likely to remain limited, one expects that approaches such as activespace divisions of the orbitals will continue to play a dominant role.While powerful, these techniques introduce a level of approximation that may prevent us from reaching the desired accuracy for problems of interest without additional qubits or gate depth.
Here we introduce a method for going beyond the activespace approximation with no additional qubits or gate complexity.We show how this method may be constructed and demonstrate its power for simple test systems.While the number of measurements is increased, there are a number of promising approximation schemes that may allow one to avoid this additional overhead in an intelligent way.Moreover, these methods maintain an exponential advantage over their classical counterparts in the treatment of the active space and active-space reference.We believe methods such as these will be crucial for obtaining accurate solutions on both NISQ and fault-tolerant devices.Additionally, in the long term, these techniques should be useful for the competitiveness of fault-tolerant approaches to simulating quantum chemistry.Studies such as Refs.[16,18] have focused on performing error-corrected quantum computations of chemistry within an active space, and the methods that we develop here should help significantly to capture dynamic correlation in those simulations; in fact, a MCSCF procedure was even suggested in Ref. [16].However, more research would be needed to make such schemes optimal within error correction.This is because both schemes that we present here require a large (although polynomial scaling) number of measurements, and error-corrected logical gates are many orders of magnitude slower than the physical gates used in NISQ devices.FIG. 5. Ground-state energy error with respect to the exact solution in a 6-31G basis as a function of internuclear separation for H 2 .The notation ðn ¼ n q Þ after each line indicates that n q qubits are needed for that representation.The orbital optimization postprocessing is labeled by OO.We see that orbital optimization greatly reduces the need for perturbative corrections in the stretched bonding region.FIG. 6. Ground-state energy as a function of the OH-H bond distance for the water molecule.Using an active space within the large cc-pVTZ (VTZ) basis, we show that the OO step improves upon the active-space solution with (n ¼ 14) qubits (CAS) achieving the accuracy of a larger active space (n ¼ 20) in the bonding region, and superior accuracy in the dissociation region due to better treatment of dynamical correlation.The solution in a minimal basis exact solution, which requires the same number of qubits, is not shown here, due to being in error by over a hartree.