Iterative Qubit Coupled Cluster approach with efficient screening of generators
Ilya G. Ryabinkin Robert A. Lang Scott N. Genin Artur F. Izmaylov
=1




 

 

 
 
 
 
 
 
 
 
 
 
 
 
 



 
 

 








 



Iterative Qubit Coupled Cluster approach with efficient
  screening of generators 
Ilya G. Ryabinkin 
Robert A. Lang 
Scott N. Genin 
Artur F. Izmaylov 

Quantum computing, quantum chemistry 





 
 
 
 
 
 
⌈⌉ 
⌊⌋ 
 
 
 
 
 
 
 
 
 
 
 

    [number-unit-product =  ]cal
        
          range-phrase = –, 
        
        
        
        
        
        
        
        
        
        
        Y>X
        
        [longplural=degrees of freedom, firstplural=degrees of
          freedom (DOF), plural=DOF]DOFDOFdegree of freedom
        [longplural=equations of motion, firstplural=equations
          of motion (EOM), plural=EOM]EOMEOMequation of motion
        NISQNISQnoisy intermediate-scale quantum
        JWJWJordan–Wigner
        BKBKBravyi–Kitaev
        QPEQPEquantum phase estimation
        VQEVQEvariational quantum eigensolver
        QMFQMFqubit mean-field
        QCCQCCqubit coupled cluster
        iQCCiQCCiterative qubit coupled cluster
        IQCCIQCCiterative qubit coupled cluster
        PQAPQAparametrized quantum annealing
        DISDISdirect interaction set
        CICIconfiguration interaction
        FCIFCIfull configurational interaction
        CASCICASCIcomplete active space configurational
          interaction
        MCSCFMCSCFmulticonfigurational self-consistent
          field
        CASSCFCASSCFcomplete active space self-consistent
          field
        CCCCcoupled cluster
        UCCUCCunitary coupled cluster
        UCCSDUCCSDunitary coupled cluster singles and
          doubles
        CCSDCCSDcoupled-cluster singles and doubles
        CCSD-TCCSD(T)coupled-cluster singles and doubles and
          non-iterative triples
        RHFRHFrestricted Hartree–Fock
        UHFUHFunrestricted Hartree–Fock
        DMRGDMRGdensity-matrix renormalization group
        DFTDFTdensity-functional theory
        CASCAScomplete active space
        PESPESpotential energy surface
        PECPECpotential energy curve
        MOMOmolecular orbital
        SQPSQPsequential quadratic programming
        MMAMMAmethod of moving asymptotes
        
        
        
        ilya.ryabinkin@otilumionics.com
        OTI Lumionics Inc., 100 College St. #351, Toronto,
          Ontario  M5G 1L5, Canada
        Department of Physical and Environmental Sciences,
          University of Toronto Scarborough, Toronto, Ontario, M1C 1A4,
          Canada; and Chemical Physics Theory Group, Department of Chemistry,
          University of Toronto, Toronto, Ontario, M5S 3H6, Canada
        scott.genin@otilumionics.comOTI Lumionics Inc.,
          100 College St. #351, Toronto, Ontario  M5G 1L5, Canada
        artur.izmaylov@utoronto.caDepartment of Physical
          and Environmental Sciences, University of Toronto Scarborough,
          Toronto, Ontario, M1C 1A4, Canada; and Chemical Physics Theory
          Group, Department of Chemistry, University of Toronto, Toronto,
          Ontario, M5S 3H6, Canada
        
          An iterative version of the QCC
          method [I.G. Ryabinkin et al., J. Chem. Theory Comput.14, 6317 (2019)] is proposed. The new method seeks to find
          ground electronic energies of molecules on NISQ devices. Each
          iteration involves a canonical transformation of the Hamiltonian and
          employs constant-size quantum circuits at the expense of increasing
          the Hamiltonian size. We numerically studied the convergence of the
          method on ground-state calculations for LiH, H2O, and
          N2 molecules and found that the exact ground-state energies can
          be systematically approached only if the generators of the QCC
          ansatz are sampled from a specific set of operators. We report an
          algorithm for constructing this set that scales linearly with the
          size of a Hamiltonian.
        
            
            December 30, 2023
        =====================
        
        
        
        § INTRODUCTION
        
        
        The advent of commercial quantum computers greatly stimulated a desire
        to use them to solve practically relevant hard computational problems.
        One such problem is the electronic structure
        problem <cit.>. Current and near-future quantum
        computers are NISQ devices <cit.>, which
        are restricted in the number of available qubits, in qubit
        connectivity, and in the fidelity of single- and multi-qubit
        entangling gates. Algorithms for such hardware need to minimize the
        gate count and be able to withstand noise.
        
        VQE <cit.> is
        one such algorithmic framework. It engages both quantum and classical
        computers in an iterative optimization of the system wave function
        using the variational principle. The quantum computer constructs a
        wavefunction guess
        |Ψ(τ)⟩ = Û(τ)|0⟩ as
        a sequence of gates representing a parametrized unitary
        Û(τ) acting on some initial qubit wave
        function |0⟩; τ is a set of numerical
        parameters. To obtain the expectation value of energy,
        E(τ) = ⟨Ψ(τ) | Ĥ
          |Ψ(τ)|$⟩, the quantum computer performs a series of
        measurements, which involve|Ψ(τ)⟩and a qubit
        HamiltonianĤ.Ĥis derived from the second-quantized
        form of the electronic HamiltonianĤ_eof a problem using a
        fermion-to-qubit transformation. The classical computer accepts the
        energy estimateE(τ)and provides an updated set of
        parameters,τ', to start the next cycle of the
        algorithm.
        
        The essential component of VQE is a parametrized unitaryÛ(τ)which defines a form (ansatz) of a wave
        function and determines the accuracy of the method. The VQE does
        not specify it explicitly; the only practical constraint onÛ(τ)is its length when expressed in terms of
        universal qubit gates.
        
        One of the first ansätze explored for VQE was the UCCSD
        form <cit.>. The UCC parametrization has
        several advantages: 1) it is systematically improvable due to its
        clear fermionic excitation hierarchy, 2) it is
        size-consistent <cit.>, and 3) it is
        variational <cit.>. Besides, UCC is highly accurate
        already at the UCCSD level and rapidly convergent for molecules
        near equilibrium configurations <cit.>. However, due to general non-commutativity of
        involved operators, the UCC form cannot be directly translated
        into a sequence of quantum gates without an additional Trotter
        approximation <cit.>.
        Furthermore, fermionic excitation operators tend to produce redundant
        terms in the qubit representation <cit.>. This observation together with the
        NISQ hardware restrictions prompted a search for more efficient
        UCC forms <cit.>.
        
        On the other hand, experimental simulations of small molecules carried
        out on an existing NISQ device put forward a
        “hardware-efficient” ansatz <cit.>, which is a regular, periodic sequence
        of parametrized single-qubit and fixed-amplitude two-qubit
        gates. While closely matching the hardware requirements, this
        ansatz poses a problem of slow convergence with the number of gates.
        The latter creates an overhead for the classical computer because
        high-dimensional global minimization with parametrized circuits scales
        exponentially with the number of
        dimensions <cit.>.
        
        A hardware-oriented approach has stimulated an interest in
        methods <cit.> that operate directly in the space
        of multi-qubit operators:
        
            P̂_k = ∏_iσ̂_i^(k),  σ∈{1̂, x̂, ŷ, ẑ}.
        
        One such method, the QCC, was introduced in our early
        work <cit.> and used the following ansatz:
        
            Û(τ)  = ∏_kexp(-τ_k P̂_k/2).
        
        The total number ofP̂_ks grows exponentially with the number of
        qubits. Therefore, to select the most relevantP̂_ks, QCC
        uses a screening procedure. The screening is based on the energy
        derivative with respect to theP̂_kamplitude,τ_k, taken
        atτ = 0. Thus,P̂_ks with the largest energy
        derivative magnitudes are included in Eq. (<ref>) first.
        This ranking can be done efficiently on a classical computer. However,
        the procedure's simplicity is outweighed by the exponential number of
        operators that need to be tested, which limits the applicability of
        the QCC method.
        
        In Ref. , the exponential
        ranking problem of the QCC method was avoided by limitingP̂_ks to a polynomial number of those that are produced by all
        single and double fermionic excitations (an “operator pool”). Since
        this set cannot provide convergence to exact energy (otherwise the
        UCCSD method would be exact), an iterative scheme was employed:
        at each iteration,P̂_ks from the pool are ranked using partial
        derivatives of energy with the wave function generated at the previous
        iteration. Calculations of these derivatives are parallelizable, but
        require a quantum device. Unfortunately, the iterative refinement of
        the ansatz increases the size of the corresponding quantum circuit,
        and eventually exhausts the capacity of an NISQ device.
        Moreover, convergence towards the exact energy is warranted only if
        all parameters in the ansatz at each iteration are fully
        re-optimized. The computational cost of this optimization grows
        exponentially with the number of parameters <cit.>.
        
        In this work we address two problems. First, how to construct and
        characterize all operators that have significant energy derivatives
        without the exponential screening procedure. Second, how to formulate
        an iterative VQE-type procedure that avoids expansion of a
        quantum circuit by delegating additional work to the classical
        computer and performing more measurements. The use of fixed-size
        quantum circuits will enhance applicability of NISQ devices for
        solving the electronic structure problem.
        
        The rest of the paper is organized as follows. First, after
        introducing several prerequisites to the QCC approach, we
        formulate a new polynomially scaling generator screening procedure.
        Second, we formulate an iterative QCC scheme and estimate its
        resource requirements. Third, approaches for reducing the iterative
        QCC resource requirements are discussed. Fourth, numerical
        benchmarks for a few molecular systems (LiH, H2O, and
        N2) are presented.
        
        
        
        § THEORY
        
        
         §.§ A few definitions
        
        
        The QCC method starts at the second-quantized electronic
        Hamiltonian of a molecule:
        
            Ĥ_e = ∑_ij^N_so h_ijâ^†_i â_j + 1/2∑_ijkl^N_so
              g_ijklâ^†_i â^†_j â_l â_k,
        
        whereâ_i^†andâ_iare fermion creation and
        annihilation operators, whileh_ijandg_ijklare molecular
        one- and two-electrons integrals, written in a spin-orbital basis. The
        number of spin-orbitals,N_so, can be as large as2N_AO, twice the number of orbitals in the atomic basis set
        chosen for a molecule, but in the present work we consider
        active-space Hamiltonians withN_so < 2N_AO. To obtain
        an active-space Hamiltonian one has to prepare a basis of MO,
        typically by running the Hartree–Fock calculations and then
        transforming one- and two-electron integrals to that basis. After
        that, core (always occupied) and frozen virtual (always empty)
        orbitals must be specified; the contribution of the latter to the
        Hamiltonian (<ref>) may be simply dropped, while the
        contribution of the former must be re-calculated
        explicitly <cit.>. The resulting orbital count is:N_so = 2N_AO - 2N_core - 2N_frozen,
        where factors of 2 account for the doubling of the number of
        spin-orbitals as compared to the number of spatial orbitals.
        
        
        
        
        
        
        
        
        
        For VQE, the electronic Hamiltonian (<ref>) is converted to
        a qubit form by one of the fermion-to-qubit transformations,
        JW <cit.> or
        BK <cit.>
            Ĥ = ∑_k=1^M C_k P̂_k,
        
        whereC_kare coefficients, andP̂_kare Pauli strings
        [“words,” see Eq. (<ref>)]. The number of qubitsninĤis equal to the number of spin-orbitals in the active
        space,n = N_so. The total number of terms inĤ,MisO(n^4), because the fermion-to-qubit transformations map each
        term ofĤ_eto the constant number of Pauli
        words <cit.>.
        
        In QCC, an initial state|0⟩is parametrized as a direct product
        ofncoherent
        states <cit.>:
        
            |Ω⟩    =    ∏_j=1^n|Ω_j⟩, 
            |Ω_j⟩    =    cos(θ_j/2)|↑⟩_j + 
                                   ^ϕ_jsin(θ_j/2)|↓⟩_j,
        
        where|↑⟩_jand|↓⟩_jare eigenstates ofẑ_j, andΩ = {θ_1, ..., θ_n}∪{ϕ_1, ...,
        ϕ_n}are the corresponding Bloch angles taken as variational
        parameters. Such a parametrized product state is referred hereinafter
        as the QMF wave function <cit.>. The
        total QCC energy is
        
            E_QCC    = min_Ω, τ⟨Ω |
                Û(τ)^†ĤÛ(τ)|Ω|.⟩
        
        To estimate the contribution of eachP̂_ktoÛ(τ), one computes a modulus of the QCC
        energy derivative atτ = 0:
        
            | d E[P̂_k] /dτ_k|_τ=0| 
                =   |d/dτ_kmin_ΩΩ | ^τ_k P̂_k/2Ĥ^-τ_k P̂_k/2 | Ω|_τ_k =
                   0| 
                =   |QMF | -/2 [Ĥ,
                   P̂_k]|QMF|≥ 0,
        
        where|QMF⟩ = |Ω_min⟩is
        the QMF wave function at the QMF energy minimum.
        
        
        
         §.§ Efficient screening procedure
        
        
        Consider an arbitraryP̂_ithat satisfies the gradient
        condition (<ref>). In the absence of external
        magnetic fields, the electronic Hamiltonian in Eq. (<ref>)
        is real. This results in real coefficients ofĤand an even
        number ofŷterms inP̂_ks. Accounting for this, the
        energy gradient forP̂_ican be rewritten as:
        
            dE[P̂_i]/dτ = ∑_k C_k ImΩ_min| P̂_k P̂_i |
              Ω_min.
        
        For any|Ω_min⟩, non-vanishing
        contributions in Eq. (<ref>) can be produced
        only byP̂_k P̂_iwith purely imaginary matrix elements,
        requiringP_ito have odd powers ofŷterms.
        
        Frequently, the optimized QMF state|Ω_min⟩is an eigenstate of all{ẑ_i}_i=1^noperators and, hence, any product of them:
        
            ∏_iẑ_i |Ω_min⟩ = ±|Ω_min⟩.
        
        If the lowest-energy QMF solution does not satisfy this
        condition, one can define a “purified” mean-field state|Φ_0⟩that satisfies Eq. (<ref>) and
        has the maximum overlap with the lowest-energy QMF solution.
        Thus, for further analysis, we will use|Φ_0⟩.
        
        With a property of the independent-qubit reference
        Eq. (<ref>), the non-vanishing terms in
        Eq. (<ref>) are the productsP̂_k P̂_ithat contain onlyẑoperations, since
        
            ⟨Φ_0 | f̂_j | Φ_0|=⟩ 0,
        
        wheref̂_j ∈{x̂_j, ŷ_j }is a generalized
        flip operator acting on thej^ thqubit. We denote the
        flip indicesF(T̂)of a Pauli termP̂as
        
            F(P̂) = {j : f̂_j ∈P̂}.
        
        
        Using Eq. (<ref>) we find that the only non-zero
        contributions in[Ĥ, P̂_i]for a givenP̂_iare from
        thoseP̂_ks which have the same flip indices asP̂_i.
        This leads to partitioning of the original Hamiltonian as
        
            Ĥ    = ∑_k Ŝ_k,
        
        where
        
            Ŝ_k = ∑_j C_j P̂_j
        
        group the terms with the same flip indices,
        
            F(P̂_i)     = F(P̂_j),   ∀ (P_i, P̂_j) ∈Ŝ_k.
        
        Thus, Pauli words possessing the same flip indices introduce an
        equivalence relation on the set of Hamiltonian terms. As a result, all
        the generatorsP̂_iwith evenŷparity or withF(P̂_i) F(Ŝ_k),∀ kin
        Eq. (<ref>) have zero energy gradients.
        Furthermore, any two generatorsP̂_iandP̂_jwith the
        sameŷparity andF(P̂_i) = F(P̂_j)will have
        identical gradients up to a sign, and hence the same gradient
        magnitudes. Therefore, any generators obtained by replacements of1̂withẑoperators or permutations ofx̂andŷthat conserveŷparity have the same gradient
        magnitude. This leads to a set of2^n-1generators that are
        characterized by the same absolute energy gradients.
        
        The number of equivalence classes [terms in
        Eq. (<ref>)] is bound from above by the total
        number of terms inĤ,M = O(n^4). Thus, the set of all
        operators that satisfy the gradient
        condition (<ref>) has the sizeO(M 2^n-1)with a partitioning intoO(M)groups. We will refer to this set as
        the DIS.
        
        A representative operator from the DIS can be constructed as
        follows. First, partition the Hamiltonian as in
        Eq. (<ref>) by grouping its terms according to
        their flip indices. Second, for a givenŜ_ktake a product ofx̂operators for all but one index fromF(Ŝ_k),
        Eq. (<ref>), and multiply it by a singleŷoperator with the remaining index. The resulting Pauli wordP̂_khas the odd (1) number ofŷoperators and is characterized by
        the same flip set asŜ_k. Third, compute the energy gradient by
        Eq. (<ref>) by takingŜ_kasĤandP̂_kas a candidate; a modulus of the resulting value will characterize the
        gradient group corresponding toŜ_k. Finally, repeat these
        steps for eachŜ_kto find all representatives and their
        gradients, thus obtaining the full description of DIS. Since the
        number ofŜ_ks is bound by the polynomial number of terms inĤthe number of different gradient groups is also polynomial.
        Due to symmetries encoded in the Hamiltonian coefficientsC_j, some
        of the representatives can have gradients close to zero.
        
        
        
        
         §.§ The iQCC
        
        
        By rewriting the QCC energy expression (<ref>) as
        
            E_QCC = min_τ{min_Ω⟨Ω| Ĥ_d(τ) |Ω|⟩},
        
        one demonstrates that the QCC energy is the minimum of the
        QMF minima for a canonically transformed (“dressed”)
        Hamiltonian,
        
            Ĥ_d(τ) = U^†(τ) ĤÛ(τ),
        
        parametrized by the set of amplitudesτ.Ĥ_d(τ)can be evaluated recursively as
        
            Ĥ_d^(k)(τ_k, …, τ_1) =    ^τ_k P̂_k/2 Ĥ_d^(k-1)(τ_k-1, …, τ_1)  ^-τ_k P̂_k/2
             
              =    Ĥ_d^(k-1) + sinτ_k
                (-/2[ Ĥ_d^(k-1), P̂_k ]) 
                + 1/2(1 - cosτ_k) (P̂_k
                Ĥ_d^(k-1)P̂_k - Ĥ_d^(k-1)) ,
        
        wherek = 1, …, N_gandĤ_d^(0) = Ĥ. This
        procedure produces3^N_gdistinct operator terms and exposes the
        exponential complexity of the QCC form for a classical computer.
        However, as shown in Appendix <ref>, if amplitudesτare fixed, the complexity of the dressing
        Hamiltonian by Eq. (<ref>) is∼ M(3/2)^N_g.
        
        This observation suggests an iterative reformulation of the QCC
        procedure. Instead of a single-step optimization ofN_g > 1amplitudes, one can use multiple steps and optimizeN_gamplitudes
        sequentially. The number of operators introduced at each step
        is a constant that can be as low as 1, which means that a quantum
        circuit of a fixed size can be used at each iteration. However,
        this iterative formulation does not guarantee the convergence to the
        exact answer. We did not find rigorous conditions when such
        convergence was possible and resorted to numerical experiments (see
        Sec. <ref>).
        
        The iQCC algorithm is summarized below. The number of steps,N_steps, and the number of generators,N_g ≥ 1, which
        will be used at each step, are parameters of the scheme. The
        initialization step is the QMF energy minimization to determineE^(0)_QCC = E_QMFand the initial set of Bloch
        angles,Ω^(0). The iQCC loop is:
        
          *  Run a generator sampling algorithm using the current Hamiltonian
          Ĥ_d^(k-1) and Bloch angles from a previous iteration (to
          construct the mean-field reference state,
          |Ω^(k-1)⟩) to identify N_g generators with
          the largest absolute gradients in a given pool. If the highest
          gradient is lower than a threshold, terminate the loop.
        
        
          *  Minimize the QCC energy, Eq. (<ref>) with
          respect to N_g amplitudes and 2n Bloch angles starting from a
          random guess. If this search with a small (usually  10) number of
          guesses fails to locate the solution with lower energy than on a
          previous iteration, perform an additional minimization using
          τ = 0 and Bloch angles from the previous iteration
          as a guess. The last attempt is guaranteed to lower energy
          because the chosen generators have non-zero energy gradients by
          construction. The random-search stage is introduced to prevent
          sticking in local minima and saddle points. You may also terminate
          the procedure here if the energy difference between the current and
          previous iterations is below a threshold.
        
        
          *  Replace the current Hamiltonian by its dressed version
          calculated by Eq. (<ref>) using the N_g
          amplitudes optimized at the current iteration.
        
        
          *  Compress the resulting Hamiltonian using the techniques from
          Sec. <ref> (optional).
        
        
          *  If the number of steps exceeds N_steps, exit. Otherwise
          start a new iteration.
        
        The QCC energy at exit is the final result: the ground-state
        energy estimate for the HamiltonianĤ.
        
        
        
         §.§ Compression of intermediate Hamiltonians
        
        
        Since the size of the Hamiltonians increases in the course of the
        iterations, it is natural to seek a method to “compress” them in
        such a way as to guarantee that their ground-state energies differ
        less than a desired accuracyϵ,
        
            |E_0(Ĥ_d^(k)) - E_0(Ĥ_c^(k))| ≤ϵ.
        
        One can set, for example,ϵ = 1, which
        is better than the so-called “chemical accuracy,”1.
        
        The compression procedure that we propose is based on the Weyl's
        spectral perturbation theorem <cit.>:
        
            max_j |λ_j^↓(Ĥ) -
              λ_j^↓(Ĥ_c)| ≤Ĥ - Ĥ_c≤Ĥ - Ĥ_c_F,
        
        where{λ_j^↓}are eigenvalues of the
        corresponding operator arranged in decreasing order,Â = sup_Ψ = 1ÂΨis the operator
        norm, which is for a normal (diagonalizable) operator equal tomax_j{λ_j}, andÂ_F = √((Â^†Â))is
        the Frobenius norm ofÂ[To align our consideration
          with the Weyl's theorem, one should use the negate of the real
          Hamiltonian, -Ĥ.]. The Frobenius norm of the qubit
        Hamiltonian (<ref>) is easy to evaluate:
        
            Ĥ_F = 2^n/2√(∑_j |C_j|^2),
        
        sinceTr(P̂_i^†P̂_j) = 2^nδ_ij,
        whereδ_ijis the Kronecker symbol. Thus, if all|C_j|are
        sorted in descending order, we can define an approximate (compressed)
        Hamiltonian as
        
            Ĥ_c = ∑_j=1^J C_j P̂_j,
        
        whereJsatisfes
        
            √(∑_j = J+1 |C_j|^2)≤ϵ/2^n/2.
        
        At each iteration of the iQCC method one can replace the dressed
        HamiltonianĤ_d^(k)with its compressed version,Ĥ_c^(k), to use it as a starting operator for the next
        iteration. According to inequality (<ref>), this will
        change the spectrum by no more thanϵ.
        
        
        The suggested compression procedure is well-suited for use with the
        iQCC method. As we established in
        Appendix <ref>, the main reason for growing the
        dressed Hamiltonians is the commutator term in
        Eq. (<ref>). However, its average value on the
        QMF wave function is precisely the value of the gradient
        contribution of the corresponding generatorP̂, which is
        systematically reduced by the iQCC procedure. Thus, after
        the initial rapid growth in size of the intermediate Hamiltonians, one
        could expect progressively stronger compression when the commutator
        contributions start systematically falling below the compression
        thresholdϵ. We verify these expectations numerically in
        Sec. <ref>.
        
        
        § NUMERICAL EXAMPLES
        
        
         §.§ Electronic structure calculations and generation of qubit
          Hamiltonians
        
        
        Molecular configurations for all species used in this study, as well
        as the atomic basis sets and active spaces, are listed in
        Table <ref>. One- and two-electron integrals for
        Hamiltonians Eq. (<ref>) in the Hartree–Fock MO
        basis, were computed by a locally modified version [Feb. 14, 2018
        (R1)] of the gamess electronic structure
        package <cit.>. Qubit Hamiltonians were derived
        from their fermionic counterparts by applying the appropriate
        fermion-to-qubit transformation (see Table <ref>). In all
        cases, the qubit Hamiltonians had stationary
        qubits <cit.>: namely, qubit operators at
        positionsN_soandN_so/2, enter the Hamiltonian asẑor1̂. Thus,z-projections of these qubits are
        constant (stationary), and the corresponding operators can be replaced
        with their eigenvalues (± 1) to define qubit-reduced operators.
        This procedure has been applied to all Hamiltonians, and final qubit
        counts are listed in Table <ref>. Eigenvalues ofẑoperators for stationary qubits were chosen so as to ensure that the
        ground states of the full and reduced Hamiltonians are the same. The
        qubit reduction procedure has also been applied to operators other
        than a Hamiltonian; see Appendix <ref>.
        
        
        
         §.§ LiH and H2O near equilibrium geometry
        
        
        For small-qubit problems like LiH and H2O, increasing the
        size of intermediate Hamiltonians poses no difficulties, and no
        mitigation techniques are necessary. Moreover, rapid convergence of
        energy in this weak-correlation regime allowed us to gauge the ability
        of the iQCC method to approach the exact energy. We tested first
        the most challenging case,N_g =1, in which only one amplitude is
        optimized at each iteration. This regime highlights the importance of
        choosing appropriate generators for the QCC ansatz.
        
        We compare three operator pools: all two-qubit Pauli operators, the
        operators from spin-adapted single and double fermionic excitation
        operators transformed to a qubit space (like in
        Ref. ), and the DIS. In
        any case, a single top-gradient [according to Eq. (<ref>)]
        generator of entanglement was chosen for the QCC
        ansatz (<ref>) at each iteration.
        
        Convergence of the iQCC procedure with different operator pools
        is shown in Fig. <ref>.
        
        Energies for non-DIS pools saturate before reaching the exact energy,
        while DIS provides systematic and almost geometric convergence
        towards the exact energy. We conjecture that the success of the
        DIS is due to its additional variational flexibility as it grows
        together with the growth of intermediate Hamiltonians.
        
        
        
          §.§.§ Adding multiple generators
        
        
        A faster converging setup for the iQCC method is a scheme with
        the addition ofN_g > 1generators at each step to fully exploit the
        capabilities of an NISQ device. However, in this case, a problem
        of degeneracy appears. When generators are subjected to the
        condition (<ref>), some of them have identical
        gradient magnitudes and, hence, cannot be ordered. This is a
        consequence of the structure of the DIS, which is revealed in
        Sec.  <ref>, but to the extent in which the
        other operator pools overlap with the DIS, it is relevant to
        them as well. The simplest possible strategy in this situation is a
        stochastic sampling, in which we draw a single random representative
        fromN_ghighest-gradient groups. We believe that this procedure
        reduces the probability of being trapped in local minima.
        
        With multiple generators a faster convergence is observed; see
        Fig. <ref>. Only 5 iterations are needed
        forN_g = 6to bring the energy closer than1e-8to
        exact for the LiH molecule. A similar trend was observed for the
        H2O molecule (not shown).
        
        Sometimes grouping found less thanN_gdistinct gradient groups. In
        this caseN_gis dynamically adjusted to match the number of groups;
        these numbers are also shown in
        Fig. <ref>.
        
        We do not expect this problem to be important for larger-qubit
        problems, since the number of gradient groups will be much larger thanN_g.
        
        Larger QCC ansatz can also mitigate the issue of premature
        convergence when the limited operator pools are used. For example, we
        found (see Fig. <ref>) that withN_g = 6the iQCC
        procedure with sampling from the pool of two-qubit Pauli operators
        converges to a much lower (albeit not to the exact) energy than the
        same procedure withN_g = 1(cf. Fig. <ref>).
        Thus, with a more sophisticated QCC form, and a more capable
        NISQ device the use of limited operator pools may be justified.
        
        
         §.§ Stretched H2O
        
        
        To assess the iQCC method to treat strong correlation, we
        consider a stretched H2O molecule withR(O-H) = 2.35≈
        2.5 R_e[R_e = 0.95, is the O-H
          distance in the water molecule at equilibrium geometry]. To obtain
        the singlet QMF solution, we used the penalty
        method <cit.>, where a spin-penalized
        Hamiltonian
        
            Ĥ + μ/2Ŝ^2,  μ >0
        
        for sufficiently largeμ(≥0.25) is used.
        Noticeable deviations from the spin purity in the course of iterations
        are tolerable because the exact ground state is a singlet, and, as
        long as the convergence is established, spin deviations vanish.
        
        The convergence of the iQCC procedure with different operator
        pools is shown in Fig. <ref>.
        
        Compared to the weak-correlation case (Fig. <ref>) the
        convergence is slower. Three times more iterations are necessary to
        reach similar accuracy with the DIS. Moreover, the initial onset
        of geometric convergence lasts longer, and the convergence curve is
        non-linear in the logarithmic scale for the first∼ 70iterations. The other pools exhibit the same behaviour as in the
        weak-correlation case. They level off before reaching the exact
        energy. This behaviour is due to a sequential optimization of
        amplitudes in the ansatz. The amplitudes from earlier iterations
        cannot be adjusted later. Surprisingly, the number of terms in dressed
        Hamiltonians at intermediate iterations are comparable between
        different pools, suggesting that the DIS pool surpasses its
        competitors by providing more efficient dressing.
        
        
        
         §.§ Compression and extrapolation for large-qubit problems
        
        
        We consider a 14-qubit problem for N2 near the equilibrium
        geometry (i.e. in the weak-correlation regime) as a prototype
        for “large-scale” problems. In the tensor-product basis the
        Hamiltonian matrix has the size2^14× 2^14, and we expect
        a rapid growth of the size of intermediate dressed Hamiltonians. This
        is confirmed by the results in Fig. <ref>. After
        30 iterations, the uncompressed canonically transformed Hamiltonian
        has∼5e5terms.
        
        The onset stage, in which the growth follows the trend predicted in
        Sec. <ref>,M (3/2)^Iteration, continues for the
        first 10iterations. After that the proliferation of terms slows
        down.
        
        Compression allows one to have∼ 1.5times fewer terms with a
        negligible impact on energies (Fig. <ref>). The
        QCC estimates at iterations differ by less than
        3e-6 without apparent accumulation of errors. Thus, the
        spectral compression can be applied with even more aggressive
        thresholds due to a rather conservative estimate of ground-state
        disturbance.
        
        Having fewer terms in the Hamiltonian noticeably improves efficiency
        of the iQCC procedure. However, even after 40 iterations the
        deviation from the exact energy is more than 1
        (see Fig. <ref>), so that comparable accuracy in
        relative energies cannot be guaranteed. Noticing a regular
        (geometric) convergence pattern of the iQCC energies, we suggest
        using extrapolation to improve final estimates. This technique
        gained some popularity, for example, in the DMRG
        calculations <cit.>.
        
        As evident from Fig. <ref>, the convergence is close to
        geometric (linear in the logarithmic scale) after the first∼ 10iterations.
        
        If one assumes perfect geometric convergence, then
        
            log(E_QCC^(k) - E_exact) = ak + b.
        
        Of course, the exact energyE_exactis not known. However, by
        exponentiating, subtracting, and taking the logarithm again one can
        obtain
        
            log(E_QCC^(k-1) - E_QCC^(k)) = a'k + b',
        
        where
        
            a     = a', 
            
              b     = b' - log(10^-k' - 1)
        
        with an obvious requirementk' < 0.a'andb'can be determined
        via the least-square fitting of Eq. (<ref>), and the
        exact energy estimate at each iteration may be computed as
        
            E_exact = E_QCC^(k) - 10^(ak + b).
        
        We dropped the first 10 values from iQCC iterations and used the
        remaining ones to fit Eq. (<ref>) and to calculate
        the extrapolated values by Eq. (<ref>). For the last
        10 iterations, the extrapolated value was around
        -109.0340, which is 1 above the
        exact one -109.0350, a 5-fold increase in precision as
        compared to the raw values. Thus, extrapolation may be useful for
        large-scale iQCC iterations.
        
        
        
        § CONCLUSIONS
        
        
        
        We have presented and assessed a new technique for the ground-state
        electronic structure calculations dedicated specifically for use with
        NISQ devices – the iQCC method. It is an iterative
        modification of our previous technique, the QCC
        method <cit.>, with a few important
        amendments.
        
        First, a new selection procedure for generators of the QCC
        ansatz (<ref>) is proposed. It has computational complexity,
        which is linear in the number of terms (M) of the Hamiltonian, a
        significant improvement over the previous, exponentially difficult
        procedure. Apart from that, the new procedure unveiled a non-trivial
        structure of the set of operators satisfying the QCC selection
        condition, Eq. (<ref>). This DIS is a
        union ofO(M)groups, each of which contains precisely2^n-1operators with identical gradient magnitudes.
        
        Secondly, our numerical experiments indicate that the iQCC
        procedure converges to the exact energy even with a fixed-size
        QCC ansatz when the DIS is used as a pool of generators.
        Compared to competitive approaches, such as the ADAPT-VQE
        method <cit.>, the QCC method
        requires neither holding of the whole ansatz on a quantum device nor
        cumbersome full re-optimization of all its parameters, thus enhancing
        practical utility of NISQ devices.
        
        It must be admitted that these advantages are not free of cost. The
        canonical transformation [Eq. (<ref>)] increases the size of
        the Hamiltonian up to a factor of3/2at each step, which means an
        increased load on a classical computer and more measurements for a
        quantum one. Recently, multiple proposals for improving the
        measurement efficiency have appeared <cit.>. It
        would be interesting to apply them to the iQCC intermediate
        Hamiltonians. To partially compensate an increased load on a
        classical computer, we introduced two techniques, compression
        and extrapolation. The first is a rigorous scheme for removing
        “unimportant” terms from the canonically-transformed Hamiltonians,
        while the second is inspired by numerical behaviour of the iQCC
        energy estimates. Both techniques improve the efficiency of
        large-scale iQCC calculations.
        
        We believe, therefore, that the iQCC method is a viable
        alternative to other approaches for electronic structure calculations
        on NISQ devices with the unique feature of employing fixed-size
        quantum circuits, though its systematic convergence is shown only
        numerically.
        
        
        
        
        § COMPUTATION SCALING OF THE CANONICAL TRANSFORMATION,
          EQ. (<REF>)
        
        
        Consider dressing of the original HamiltonianĤ_d^(0) = Ĥ, which hasM = O(n^4)terms, withÛ(τ) = exp(-τP̂/2), for a fixed value ofτ. Fork = 1, Eq. (<ref>) contains three distinct
        operator terms,Ĥ,[Ĥ, P̂], andP̂ĤP̂multiplied by numerical coefficients. First, we
        notice that the operatorP̂ĤP̂has the same
        operator terms asĤ. This is a simple consequence of the
        anti-commutation relations for Pauli elementary operators:σ̂' σ̂+ σ̂σ̂' = 0(σ̂σ̂'). Multiplying this relation byσ̂from the left and using the involutory property,σ̂^2 = 1, we write:
        
            σ̂σ̂' σ̂=
              
                -σ̂',    σ̂' σ̂
            σ,    σ' = σ̂.
        
        Thus, every term ofP̂ĤP̂has either the same or the
        opposite sign compared to that inĤ. That is, the sum ofĤandP̂ĤP̂with an arbitrary numerical
        coefficient has no more terms thanĤitself. Therefore, new
        algebraically independent terms may only come from the commutator term[Ĥ, P̂].
        
        An arbitrary Pauli word commutes with a half and anti-commutes with
        the remaining half of the4^noperators of then-qubit Lie
        algebra <cit.>. Thus, “on average” the
        commutator[Ĥ, P̂]will contain only half the terms,≈ M/2. In total, we have∼3/2Mterms inĤ_d^(1)after the dressing with one generator, or∼ (3/2)^N_g Mafter a dressing withN_gof them.
        
        
        
        
        § QUBIT REDUCTION PROCEDURE FOR GENERAL OPERATORS
        
        
        For systems with a small number of qubits, the identification and
        removal of stationary qubits provide a noticeable reduction of
        computational complexity. Operators that commute with the Hamiltonian
        represent exact symmetries. Qubits in such operators can be reduced in
        the same manner as those in the Hamiltonian; see
        Sec. <ref>. We derived reduced qubit
        expressions for the total electron number operator,N̂, and the
        total spin operatorsŜ_zandŜ^2. The corresponding full
        qubit expressions were derived from the second-quantized counterparts
        by applying the chosen fermion-to-qubit transformations (see
        Table <ref>).
        
        However, not every operator can be reduced, as its symmetry may not be
        compatible with the symmetry of the Hamiltonian. In particular, the
        UCC ansatz, whose cluster operators are written in a
        spin-orbital basis, preserves the electron-number but not the spin
        symmetry. Some ofT̂operators do not commute withŜ^2.
        Such operators couple different spin sub-blocks of the Hamiltonian,
        and their qubit expressions have operators other than1̂orẑat the position of the stationary qubits. To find reducible
        combinations, we derived spin-adapted cluster amplitudes following a
        general scheme given by Eqs. 13.7.2 and 13.7.2 of
        Ref. . In particular, we solved equations
        
            [Ŝ_±, P̂_k]     = 0, 
            
              [Ŝ_z, P̂_k]     = 0,    k = 1 (singles), 
                                      2 (doubles), …
        
        for the qubit expressions ofP̂_k,Ŝ_z, andŜ_±(spin raising/lowering) operators. The solutions are
        spin-free cluster amplitudes that have stationary qubits at the same
        positions as the Hamiltonian, and thus can be reduced by the same
        procedure.
