Increasing the representation accuracy of quantum simulations of chemistry without extra quantum resources
Tyler Takeshita Nicholas C. Rubin Zhang Jiang Eunseok Lee Ryan Babbush Jarrod R. McClean






tyler.takeshita@daimler.com
Mercedes-Benz Research and Development North America, Sunnyvale, CA 94085

nickrubin@google.com
Google Inc., Venice, CA 90291

Google Inc., Venice, CA 90291

Mercedes-Benz Research and Development North America, Sunnyvale, CA 94085

Google Inc., Venice, CA 90291

jmcclean@google.com
Google Inc., Venice, CA 90291


Proposals for near-term 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 4 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. 


Increasing the representation accuracy of quantum simulations of chemistry
 without extra quantum resources
    Jarrod R. McClean
    December 30, 2023
===========================================================================================================





§ 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 <cit.>, 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 <cit.>, the precise costs of these algorithms and methods to reduce these costs by orders of magnitude have now been developed in detail <cit.>.  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 <cit.>.  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 <cit.>, the question arises of 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) <cit.> or the quantum approximation optimization algorithm <cit.>.  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 <cit.> and different problem areas <cit.>, and have been demonstrated on numerous experimental architectures <cit.>.

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.  

Here, we discuss two methods that use no additional qubits or gate complexity to lift the active space approximation in a systematic way.  One of these procedures is a novel way to leverage the quantum subspace expansions <cit.>, also known to mitigate errors <cit.> and provide excited states in experiments <cit.>, 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 reduces 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 <cit.>, and we conclude with an outlook on the impact for near-term quantum simulations.





§ BACKGROUND AND DEFINITIONS

We begin by setting up the problem, establishing notation, and reviewing the quantum subspace expansion method <cit.> 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 <cit.>.  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

    H = ∑_ij h_ij a_i^† a_j + 1/2∑_ijkl h_ijkl a_i^† a_j^† a_k a_l

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 {a^†_i, a_j} = δ_ij, {a_i, a_j} = {a^†_i, a^†_j}= 0.  As discussed 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 multi-reference configuration interaction and perturbation theory <cit.>.  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. <ref>, in typical chemistry calculations on quantum computers, the set of basis functions is divided into 3 sets, which we denote 𝒞, 𝒜, and 𝒱 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, 𝒜.  

Our approach will utilize 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 (QSE) <cit.>.  These techniques 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 wavefunction within the active space that we denote as |Ψ_ref⟩.  We then choose a set of expansion operators {O_i} 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 {O_i |Ψ_ref⟩}.  As this basis may be non-orthogonal, 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

    H_ij = ⟨Ψ_ref| O_i^† H O_j |Ψ_ref⟩
    
        S_ij = ⟨Ψ_ref| O_i^† O_j |Ψ_ref⟩

with these matrices, one then uses canonical diagonalization to remove the approximately 0 eigenvalues of S and solve the generalized eigenvalue problem in the well conditioned subspace given by

    HC = SCE

where C is the matrix of eigenvectors in the basis {O_i |Ψ_ref⟩}, 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 <cit.>, it was suggested to use O_i approximating fermionic excitations of the form

    {a_p^† a_q | p, q ∈𝒜}

which when considered with the Hamiltonian, composed of only up to 2-particle operators, means that the matrix elements can be evaluated as sums over subsets of the 4-electron reduced density matrix (RDM) where the sum weights are determined by integrals in the Hamiltonian.  The 4-electron density matrix in the active space is 

    ^(4)D_𝒜  pqrs^tuvw = ⟨Ψ_ref|a_w^† a_v^† a_u^† a_t^† a_p a_q a_r a_s |Ψ_ref⟩.

Each of these elements can be evaluated on a quantum computer by repeated preparation of |Ψ_ref⟩ 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_A^8 where N_A is the number of active space orbitals.



§ 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 post-processing and measurements can be used to re-introduce 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 𝒜.  We now introduce a set of expansion operators

    {a_i^† a_p, a_μ^† a_q a_ν^† a_r | i ∈𝒜∪𝒱;  p, q, r ∈𝒜; μ, ν∈𝒱}

that act, in principle, on additional qubits that would define the virtual space.  However, by definition, the reference wavefunction 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 multi-reference configuration interaction (MRCI) methods here is that the reference wavefunction 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 quantum computing 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 have restricted 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_A^2(2+k) 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. (<ref>) are double excitations is

    M = ⟨Ψ_ref| a_ξ a_s^† a_η a_r^†  a_i^† a_j^† a_k a_l  a_μ^†  a_p a_ν^† a_q |Ψ_ref⟩ .

Using Wick's theorem where Greek indices denote virtual orbitals, we contract the operators on the virtual space,





    M     = Δ_ξη, μν ⟨ a_s^† a_r^† a_i^† a_j^† a_k a_l  a_p a_q ⟩ +∑_x= i, j ∑_y = l, k
          (-1)^δ_x j + δ_y k(δ_ξ xΔ_yη, μν- δ_η xΔ_yξ, μν)⟨ a_s^† a_r^† a_x̅^† a_y̅ a_p a_q ⟩
          +Δ_ξη, i jΔ_k l, μν ⟨ a_s^† a_r^† a_p a_q ⟩ ,

where Δ_ξη, μν = δ_ξνδ_ημ-δ_ξμδ_ην, x̅ (y̅) means x (y) changing to the other value in the sum (e.g. i̅ = j for x), and ⟨⋯⟩ denotes expectation values  with respect to the reference state |Ψ_ref⟩ in the active space.  The reduction in Eq. (<ref>) shows that the expectation value M can be derived by knowing the 4-election RDM of the reference state |Ψ_ref⟩ 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 post-processing and no additional qubits after the appropriate measurements have been performed.





§ 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 𝒜_v and a part of the active space which will be treated as correlated core orbitals 𝒜_c.  This reduces the scaling in the number of measurements to the cost of originally estimating the energy plus the size of the virtuals active space 𝒜_v, N_Av^8, 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 4 RDM using products of lower RDMs and perturbative corrections.  This may be denoted by

    ^(4)D =    ^(4)Δ + 4 ^(3)Δ∧^(1)Δ + 
        3 ^(2)Δ∧^(2)Δ + 6 ^(2)Δ∧^(1)Δ∧^(1)Δ + 
       ^(1)Δ∧^(1)Δ∧^(1)Δ∧^(1)Δ

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 greatly reduces the number of terms to measure back to N_A^4, but 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.






§ 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 were enabled by the OpenFermion software package for doing quantum chemistry on a quantum computer <cit.>.  The orbitals were obtained from a Hartree-Fock calculation, and the reference function |Ψ_ref⟩ was the exact solution within the active space consisting of the single occupied orbital and the lowest energy virtual orbital.  Only singles and doubles excitations from the active space to the additional virtual orbitals were included in these calculations. Correlation consistent basis set of double zeta quality (cc-pVDZ)  <cit.> are used as well as the 6-31G basis set <cit.> for reference.  We compare with the exact results with both basis sets.

The results of the numerical simulations are shown in Fig. <ref> and the energy differences with respect to Exact/cc-pVDZ are shown in Fig. <ref>. 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 4 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.





§ 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 wavefunction. Variants of this idea fall into the family of methods known as multiconfigurational self-consistent-field (MCSCF).  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  <cit.> and also motivated other embedding methods in the quantum computing context  <cit.>. Here we demonstrate how a single step of an MCSCF orbital relaxation can be used to further improve energies with a procedure that is designed for the limitations of NISQ simulations.  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. (<ref>) through one-body rotations U on the full space. Due to Thouless's theorem, the one-body rotations can be efficiently implemented as a rotation of the underlying basis <cit.>.  This can be formulated as the following non-linear optimization problem:

    min_U E = ∑_iju_i,i'u_j,j'^*h_ij⟨ a_i^†a_j⟩_𝒜 + 
    ∑_ijklu_i,i'u_j,j'u_k,k'^*u_l,l'^*h_ijkl⟨ a_i^†a_j^†a_ka_l⟩_𝒜
    s.t.  Ua_j^†U^† = ∑_j'u_j,j'a_j'^†
    
    U = e^X
    
    X = ∑_p,qt_p,qa_p^†a_q  ,   t_p,q = -t_q, p^* 
    
    U^†U = 𝕀

where ⟨ a_i^†a_j^†a_ka_l⟩_𝒜 and ⟨ a_i^†a_j⟩_𝒜 are the 2- and 1-RDM of the ground state active space wavefunction.

In the classical literature it is common to use the second order approximation of u

    u = e^t≈𝕀 + t + 1/2t^2

to minimize energy in a Newton-Raphson style <cit.>. 

A well known alternative to parameterizing the unitary as an exponentiated antihermitian matrix is the use of Givens rotations <cit.>.  This parameterization uses a set of angles {θ} associated with the set of non-redundant orbital rotation generators.  For 2-RDMs obtained from an exact diagonalization of the active space Hamiltonian the only non-redundant 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 (<ref>) can be expressed as a product of Givens matrices

    U = ∏_i ∈active∏_b∈core, virtualG_i,b(θ_i,b).

The optimal rotation of a single angle with respect to the input 2-RDM, one-, and two-electron integrals has been derived <cit.> along with a sweep procedure to find an energy minimizing U. This algorithm was used for developing 2-RDM MCSCF methods <cit.> and the first orbital-optimized coupled-cluster doubles methods <cit.>.

In the context of NISQ algorithms 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 single orbital relaxation procedure as post-processing. Because of the variational principle, relaxing the orbitals in a single step as post-processing is guaranteed to reduced the energy.  Furthermore, if the ground state of the active space is not achieved due to an approximate wavefunction 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  <cit.>.   Any circuit that contains one-body rotations at the end can be made shorter with classical post-processing by instead rotating the operator to be measured <cit.>.  

To demonstrate the utility of the orbital relaxation step for NISQ devices 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 post-processing.  In Figure <ref> we plot the energy of molecular hydrogen and observe a significant recovery of energy when performing a 4-qubit active space calculation with orbital relaxation on the full 8-qubit space.  The 6-qubit active space calculation is included to demonstrate a smooth interpolation of this method to exact diagonalization in the full 8-qubit space.  Figure <ref> shows the energy difference relative to an 8-qubit exact diagonalization, revealing we reach within 5× 10^-5 E_h of the exact solution with half the number of qubits.



§ DISCUSSION


As NISQ 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 approaches such as active space divisions of the orbitals will continue to play a dominant role.  While powerful, these techniques introduce a level of approximation that may prevent reaching the desired accuracy for problems of interest without additional qubits or gate depth.

Here we have introduced a method for going beyond the active space approximation with no additional qubits or gate complexity.  We showed how this method may be constructed and demonstrated 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 NISQ devices.

Additionally, in the long-term these techniques should useful for the competitiveness of fault-tolerant approaches to simulating quantum chemistry. Studies such as <cit.> have focused on performing error-corrected quantum computations of chemistry within an active space, and the methods developed here should help significantly to capture dynamic correlation in those simulations; in fact, an MCSCF procedure was even suggested in <cit.>. However, more research would be needed to make such schemes viable within error-correction. This is because both schemes presented 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.


apsrev4-1_with_title


