Quantum-assisted Monte Carlo algorithms for fermions

Quantum computing is a promising way to systematically solve the longstanding computational problem, the ground state of a many-body fermion system. Many efforts have been made to realise certain forms of quantum advantage in this problem, for instance, the development of variational quantum algorithms. A recent work by Huggins et al. reports a novel candidate, i.e. a quantum-classical hybrid Monte Carlo algorithm with a reduced bias in comparison to its fully-classical counterpart. In this paper, we propose a family of scalable quantum-assisted Monte Carlo algorithms where the quantum computer is used at its minimal cost and still can reduce the bias. By incorporating a Bayesian inference approach, we can achieve this quantum-facilitated bias reduction with a much smaller quantum-computing cost than taking empirical mean in amplitude estimation. Besides, we show that the hybrid Monte Carlo framework is a general way to suppress errors in the ground state obtained from classical algorithms. Our work provides a Monte Carlo toolkit for achieving quantum-enhanced calculation of fermion systems on near-term quantum devices.


Introduction
The many-body Schrödinger equation is a longstanding computational challenge in quantum chemistry [2,3], condensed matter physics [4,5], nuclear physics [6,7], quantum chromodynamics [8,9], etc. Over the years, many numerical methods are developed to address this problem, for instance, mean-field theory [10,11], Monte Carlo method [12] and tensor-network formalism [13,14]. These methods have played essential Ying Li: yli@gscaep.ac.cn roles in basic science and have also been applied to various fields outside physics, such as economic modelling [15,16] and machine learning [17,18]. However, because of the infamous difficulty that the Hilbert space dimension grows exponentially with the particle number, all known methods are limited to small-size systems, approximate solutions or specific models.
Among conventional numerical methods, quantum Monte Carlo (QMC) is a group of classical algorithms that can bypass the difficulty of exponentially growing space dimensions [19]. Specifically, QMC applies probabilistic sampling in a subset (e.g. Slater determinants) of the entire Hilbert space, resulting in a polynomial scaling of the memory size. However, for a generic Hamiltonian, QMC suffers from the sign problem where the sign of the sampling integrands oscillates alternatively from positive to negative [20]. The sign problem can lead to an impermissible statistical error. Commonly used methods to constrain the sign problem, like the phaseless approximation [21] and fixed-node approximation [22], usually bias the result. The bias is largely determined by the so-called trial state.
In recent years, quantum computing has been experiencing fast development and has shown the potential to solve classically intractable problems. Quantum computing offers a scalable way to solve the many-body Schrödinger equation given a good guiding wavefunction [23][24][25]. An established quantum-computing (QC) approach combines Hamiltonian simulation and quantum phase estimation [26,27]. However, this approach requires fault-tolerant quantum technologies, which are still years away. At present, much attention has been put into quantum algorithms suitable for noisy intermediate-scale quantum (NISQ) devices [28]. The variational quantum eigensolver (VQE) and its variants [29], as the currently most popular near-term algorithms, have proved successful in applications of a wide area. Despite this, VQE has challenges such as barren plateaus [30,31], which motivates the recent development of optimiser and variational quantum circuit design [32,33]. Meanwhile, feasible and non-variational approaches with potential quantum advantages are urged.
A promising hybrid approach is quantum computing incorporating QMC. Recently, Huggins et al. proposed a new way to calculate the ground state of fermion systems, a quantum-classical hybrid auxiliary-field Monte Carlo (AFMC) algorithm, with an experimental demonstration that it can reduce the bias [1]. The algorithm utilises a trial state prepared on a quantum computer replacing the conventional trial state. An alternative way to incorporate QMC is sampling states not from a subset in favour of classical computing but those generated by quantum circuits [34]. In this way, one can simulate the time evolution with a trade-off between the sign problem and the quantum circuit complexity. Either way, one of the central issues is whether a quantum advantage, i.e. a smaller bias or sample size compared with classical algorithms, is attainable on NISQ devices.
This paper presents a framework of assisting QMC algorithms with states prepared on a quantum computer. We generalise the original work of Huggins et al. [1] by utilizing quantum computing at flexible levels, and the QC state is used in various ways other than guiding the QMC sampling. Although shadow tomography has been proposed as the interface that passes the state from a quantum computer to a classical computer, we focus on the Hadamard test as an alternative way of measuring amplitudes of the QC state. We notice that the idea of flexible-level quantum-assisted computing works both ways. One of the purposes of this work is to address the measurement cost issue, which has raised debates about the potential quantum advantage [35,36]. By carrying out QMC calculation in a bi-trial-state manner, i.e. taking use of a classical-computing (CC) trial state in addition to a QC trial state, we can reduce the measurement cost required for demonstrating quantum advantage. Here we list the main results as follows: 1. We give a framework of quantum-assisted variational Monte Carlo (VMC), Green's function Monte Carlo (GFMC) and AFMC algorithms; 2. Depending on the extent a quantum com-puter is involved in computation, we name two strategies, the consistent quantumassisted (CQA, or consistent quantumclassical) algorithms and the quantumassisted energy evaluation (QAEE) algorithms, and the latter use minimal QC resources; 3. We propose a general Bayesian inference method to reduce quantum measurements and achieve a stable quantum advantage; 4. We show that quantum-assisted QMC is error-resilient due to an inherent symmetry projection; 5. Lastly, we propose the quantum-classical Monte Carlo subspace diagonalisation (QCMCSD) algorithm and show that hybrid QMC is a general way to reduce errors in state-of-the-art classical algorithms.
In QMC algorithms, evaluating amplitudes of the trial state is the elementary operation. Each amplitude evaluation on a quantum computer requires a number of quantum measurements (which are called circuit shots in this paper). When the circuit shot number is finite, the amplitude evaluation has a statistical error. In what follows, we use QC variance to denote this error source. We can reduce the impact of this error in various ways. In CQA algorithms, the output energy is variational in the sense that it is always higher than the true ground-state energy, even if the amplitude evaluation has a significant statistical error. Because of the variational property, we can enhance the quantum advantage by postselection, i.e. we run the algorithm a few times and select the lowest energy. Using the Bayesian inference method, we find that a bias smaller than the classical algorithm occurs with a high probability even with a small circuit shot number. In QAEE algorithms, we only utilise the quantum computer to evaluate the local energy at the last stage, such that the cost of circuit shots is minimised.
This paper is structured as follows. In section 2 we give an introduction to the Monte Carlo methods related to this work. We present these methods in a way appropriate for the generalisation to quantum computing. The main algorithms and numerical results are presented in section 3.
In section 4 we introduce three scalable methods of evaluating amplitudes on a quantum computer. Section 5 introduces a method based on the Bayesian inference that can reduce the variance due to the finite QC resources. Section 6 discusses the symmetries in fermion systems and how we can make use of symmetries to reduce errors. In section 7, we show an application of the hybrid QMC approach, i.e. the QCMCSD algorithm, to obtain a solution closer to the true ground state given optimal solutions from classical and quantum algorithms. Finally, this paper is concluded in section 8.

Quantum Monte Carlo algorithms
In VMC [37][38][39], GFMC [22,40] and AFMC [21,41] algorithms formalised with the mixed estimator, the ground-state energy E g is computed according to Here, |Ψ⟩ is an approximation to the ground state |Ψ g ⟩, and |Ψ T ⟩ is the trial state. If |Ψ⟩ = |Ψ g ⟩, E = E g regardless of the trial state. The trial state is important in these algorithms as we will show later.
QMC algorithms use Monte Carlo integration/summation to evaluate the properties of a physical system. Specifically, to compute the ground-state energy according to Eq. (1), we generate a linear expression of |Ψ⟩ in the form where l is the label of the random walker, w l and θ l are the weight and phase of the corresponding state ϕ l , respectively, and ϕ l is generated in a stochastic process. On a classical computer, the memory size for storing a generic quantum state usually increases exponentially with the number of particles. In QMC algorithms, ϕ can be stored without the exponential memory cost, for example in the form of Slater determinants for fermion systems. In the following, we call ϕ l the walker state. Given the linear expression, the energy becomes where is called the local energy. Without using a quantum computer, we need to take a trial state that quantities ⟨Ψ T |ϕ⟩ and ⟨Ψ T | H |ϕ⟩ can be efficiently computed on the classical computer.
In VMC, the walker states |ϕ⟩ are taken from an orthonormal basis of the Hilbert space R = {|R⟩}. For fermion systems, usually we take Slater determinants (SDs) as the basis. A usual way of constructing the trial state is to apply a Jastrow factor to a mean-field wavefunction [37], and VMC works for all trial states that the amplitudes ⟨Ψ T |R⟩ can be computed. Coupled-cluster states [38] and tensor-network states [39] have been investigated for this purpose in the literature. Given the basis, we can decompose the trial state as The summation in the decomposition is evaluated using the Monte Carlo method: N random basis states {|R l ⟩ |l = 1, 2, . . . , N } are generated according to the probability distribution |⟨Ψ T |R⟩| 2 . Usually the Metropolis-Hastings algorithm or some variant is used in this procedure [42]. Taking ϕ l = R l , w l = N −1 and θ l = 0 in Eq. (2), we have This state |Ψ⟩ converges to |Ψ T ⟩ in the limit N → ∞. Substituting |Ψ⟩ in Eq. (6) for |Ψ T ⟩ (not for ⟨Ψ T |) in Eq. (1), we have the final expression of the energy which coincides with Eq. (3). Essential steps of the energy evaluation in VMC can be found in Algorithm 1. Given the energy evaluation, VMC minimizes the expected energy of a parameterised trial state |Ψ T (λ)⟩ variationally, where λ denotes parameters. The expected energy is GFMC is often referred as a projector QMC algorithm. For any initial state |Ψ I ⟩ with a nonzero overlap with the ground state, the final state e −n∆βH |Ψ I ⟩ ≈ (1 − ∆βH) n |Ψ I ⟩ converges to the ground state in the limit n → ∞.
Here, e −n∆βH is the imaginary-time propagator, and ∆β is a real parameter that has to be properly chosen. Similar to VMC, we generate random basis states to represent the final state: First, we generate an initial random basis state; then, it evolves in a stochastic process driven by the operator 1 − ∆βH for n steps. In general, phases θ l are not zero in GFMC, and they can cause a large statistical error, which is the notorious sign problem. We can deal with this problem by introducing approximations in fixednode GFMC (fn-GFMC). In the fixed-node approximation, the Hamiltonian is replaced with an oscillatory-sign-free Hamiltonian H f n (called the fixed-node Hamiltonian) that depends on the trial state. In the approach introduced by van Bemmel et al. [22,40], E in the limit N, n → ∞ converges to E f n g , the ground-state energy of H f n , and E f n g is variational, i.e. it always holds that E f n g ≥ E g . We can optimise the trial state to minimise the difference E f n g − E g according to the variational principle. When |Ψ T ⟩ = |Ψ g ⟩, the fixed-node ground-state energy is exact, i.e. E f n g = E g . See Appendix A for a detailed formalism of GFMC.
AFMC is another projector QMC algorithm. Similar to GFMC, the ground-state energy is computed by simulating the imaginary-time evolution e −n∆βH |Ψ I ⟩. In AFMC, ϕ states are general SDs instead of basis states, and the trial and initial states are also SDs, e.g. worked out with the Hartree-Fock (HF) or density functional theory methods. Using the Trotter-Suzuki formula and Hubbard-Stratonovich transformation, we can approximate the imaginary-time propagator with an integral over a product of oneparticle propagators. This procedure is crucial because the evolution of general SDs under oneparticle propagators is tractable on a classical computer. These one-particle propagators depend on an auxiliary field. To generate ϕ l and the corresponding weights and phases, the auxiliary field is randomly sampled for each time step, and the initial SD evolves driven by corresponding one-particle propagators. For a generic fermion Hamiltonian, AFMC also has the sign problem (more specifically, the phase problem). Using the phaseless approximation, which is a different approach from the fixed-node approximation, we can implement AFMC without an oscillatory sign/phase, and this algorithm is called phaseless AFMC (ph-AFMC) [41]. With the phaseless approximation, the energy does not converge to the exact ground-state energy in the limit N, n → ∞. The bias to the ground-state energy depends on the trial state and vanishes when the trial state is exactly the ground state. See Appendix B for a detailed formalism of AFMC.
In QMC algorithms, there are two sources of errors, bias and variance. In VMC, E in Eq. (7) converges to E(λ) in the limit N → ∞. The bias is the difference between E(λ) and the true ground-state energy E g . When N is finite, the Monte Carlo method has a finite accuracy due to the statistical error, i.e. an error to the true mean E(λ). Both of them contribute to the final error in E. In GFMC and AFMC, we can systematically reduce the bias to zero; however, the variance may increase exponentially with time steps because of the sign problem. We can prevent the sign problem by modifying the exact GFMC and AFMC formalisms and introducing approximations. In this paper, we focus on oscillatory-signfree GFMC and AFMC algorithms. Specifically, we consider fn-GFMC and ph-AFMC. These algorithms have a residual bias because of approximations.
To summarise, some features of these algorithms can be found in Table 1.

Consistent and energy evaluation algorithms
The bias of the ground-state energy in VMC, fn-GFMC and ph-AFMC is determined by the trial state. Conventionally, options of the trial state are limited in each of the algorithms, which must permit the efficient computing of amplitudes ⟨Ψ T |ϕ⟩ on a classical computer. See Table 1 for examples. In quantum-assisted algorithms, one can choose more general trial states, for example the variant coupled-cluster states used for ph-AFMC in Ref. [1]. In this paper, we generalise this quantum-assisted approach to VMC and fn-GFMC and propose the CQA algorithms. In these algorithms, the QC trial state prepared on a quantum computer is used in two ways: generating walkers (ϕ l , w l ) (θ l = 0 for these oscillatory- Table 1: Variational Monte Carlo (VMC), Green's function Monte Carlo (VMC) and auxiliary-field Monte Carlo (AFMC) algorithms with a classical-computing trial state or a quantum-computing trial state. Unitary-coupledcluster and Hardware-efficient trial states can be prepared and optimised using the variational quantum eigensolver algorithm [43].
sign-free algorithms) and evaluating the local energy in Eq. (3). As the walker generation is costly, to minimise the QC cost, we further propose QAEE algorithms, in which the quantum computer is only used for evaluating the local energy.
Before giving the algorithms, we note that all state-related quantities in VMC, GFMC and AFMC can be derived from amplitudes. As introduced previously, the oscillatory-sign-free energy obtained by QMC methods can generally be expressed as In VMC and GFMC, states ϕ = R are basis states, and ⟨Ψ T |R⟩ are amplitudes in the basis. The other quantity we need to compute is the numerator of the local energy, i.e. ⟨Ψ T | H |R⟩. This quantity can be rewritten as where H R ′ ,R = ⟨R ′ | H |R⟩ are matrix elements of the Hamiltonian. It is required that H is sparse in the basis. Assuming the number of nonzero H R ′ ,R for the given R is L, we can compute ⟨Ψ T | H |R⟩ with L amplitudes. In AFMC, ϕ are general SDs instead of basis SDs. In addition to amplitudes, quantities used in AFMC are ⟨Ψ T | H |ϕ⟩ and ⟨Ψ T | A |ϕ⟩, where A are one-particle operators, i.e. a linear combination of fermion operators a † p a q up to a constant term, see Appendix B. The Hamiltonian in the second quantisation form is usually formed of one-particle and two-particle terms like a † p a q and a † p a † p ′ a q ′ a q . Because a † p a q |ϕ⟩ and a † p a † p ′ a q ′ a q |ϕ⟩ are also SDs, quantities ⟨Ψ T | H |ϕ⟩ and ⟨Ψ T | A |ϕ⟩ can be expressed as linear combinations of amplitudes in the form ⟨Ψ T |ϕ⟩. This paper reports three sets of algorithms. CQA and QAEE algorithms and their numerical results are presented in this section. In section 7, we present the QCMCSD algorithm, which is a generic hybrid way to reduce errors in classical algorithms in the framework of quantum-assisted QMC. Sections 4 and 5 are subroutines of the algorithms. In Section 4, we present three methods to evaluate the amplitude ⟨Ψ T |ϕ⟩ with a quantum computer. Section 5 introduces a Bayesian inference method to reduce the variance due to quantum computing. In section 6 we discuss the symmetries in fermion systems and propose a method that can reduce the errors based on symmetries.

Consistent quantum-assisted VMC and fn-GFMC
In CQA algorithms, we take a state prepared on the quantum computer as the trial state. For clarity, we use |Ψ Q ⟩ to denote such a QC trial state, and we take |Ψ T ⟩ = |Ψ Q ⟩ in VMC and fn-GFMC (see Algorithms 1 and 2). However, in quantum computing, we cannot exactly evaluate amplitudes ⟨Ψ Q |R⟩ because of the finite computational resources. Specifically, with a finite number of circuit shots to measure the amplitude, the result has a finite variance. We show that in CQA algorithms there is an effective trial state Ψ Q that determines the final result of computation, and the variational property of VMC and fn-GFMC are preserved. We can take advantage of the variational property to reduce the impact of the QC variance.
Letŷ R be the estimate of ⟨Ψ Q |R⟩ obtained in quantum computing. We measure amplitudes in a consistent way as follows: When ⟨Ψ Q |R⟩ is required for the first time, we evaluate it on the quantum computer and recordŷ R in a classical register; If the same amplitude is required later, instead of reevaluating it, we read the estimatê y R from the register. Compared with reevaluating amplitudes whenever required, the consistent approach reduces the overall QC cost by increasing the CC complexity.
To introduce the effective trial state, we think of the limit of infinite walkers, i.e. N → ∞. In this limit, all relevant amplitudes (irrelevant amplitudes are those forbidden due to certain symmetries) have been evaluated and recorded. Then, the effective trial state reads Ψ Q = Rŷ R |R⟩. The bias in VMC and fn-GFMC is determined by this trial state. We remark that in practice the number of walkers is always finite, and only amplitudes queried in the finite walkers are recorded. Notice that if amplitudes are reevaluated every time when it is required, estimates are different every time, and in this case there is not a consistent effective trial state.
The consistent effective trial state Ψ Q can also be constructed using shadow tomography [1]. Our method is proposed in the spirit of importance sampling, i.e. only those amplitudes queried in VMC or fn-GFMC are evaluated. In this work, we will not benchmark which method is more efficient. We would like to focus on the advantage of using a consistent effective trial state, which is reducing circuit shots according to the variational principle.
In CQA algorithms, the energy is variational, i.e. higher than the ground-state energy (up to a statistical error due to finite N ). In VMC, the energy in the N → ∞ limit is the expected energy of the effective trial state, i.e. ⟨ΨQ|H|ΨQ⟩ ⟨ΨQ|ΨQ⟩ .
In fn-GFMC, the fixed-node Hamiltonian is constructed according to the effective trial state, i.e. replacing Ψ T withΨ Q in Appendix A, and the corresponding fixed-node energy is always higher than the ground-state energy (using the approach in Refs. [22,40]). Therefore, according to the variational principle, we can reduce the bias by repeating the VMC or fn-GFMC calculations a few times and taking the one with the minimum energy as the final result.
With the variational principle, we can achieve the quantum improvement/advantage at a smaller QC cost. The bias in CQA algorithms is determined by the effective trial state, which is stochastic due to the variance of amplitudes. Effective trial states leading to a bias reduction occurs probabilistically. By post-selecting the lowest energy, a larger QC variance is tolerable in observing the bias reduction. Without the post-selection, we have to reduce the variance to a level that the bias reduction occurs with certainty. With the post-selection, we only need to reduce the variance to the level that bias reduction occurs with a finite probability, e.g. 50%, and the overall probability of bias reduction after the post-selection in a few trials can be high. In Sec. 5, we give a protocol based on Bayesian inference to attain a robust success probability.
We remark that the effective trial state must be consistent in the entire algorithm. In VMC, we need to sample R according to the distribution |⟨Ψ Q |R⟩| 2 . This step is straightforward in quantum computing: We can prepare and measure the state |Ψ Q ⟩ in the basis, then the probability of the outcome R is naturally |⟨Ψ Q |R⟩| 2 . However, to properly implement the consistent approach such that the variational property holds, all quantities must be derived from the same trial state Ψ Q . Therefore, we need to generate samples with a probability proportional to |ŷ R | 2 instead of |⟨Ψ Q |R⟩| 2 .
We would like to remark that VMC is discussed here because of its simplicity in concept, such that one can grasp the general idea without getting into details of more complicated algorithms, i.e. GFMC and AFMC. One may notice that the CQA VMC algorithm is likely not to outperform the VQE algorithm, however, we keep it in this work to demonstrate how the CQA approach works.

Quantum-assisted energy evaluation algorithms
The trial state is used in two phases in QMC algorithms. First, it is used to generate random walkers and yield an expression of the (approximate) ground state shown in Eq. (2). Second, given the expression and trial state, the energy in Eq. (1) is evaluated. In Ref. [1] and our CQA algorithms, the quantum trial state is used in both phases. As we see, the walker generation phase is costly. In fn-GFMC and ph-AFMC, the state |ϕ l ⟩ of each walker evolves in a stochastic process for many time steps. In each step, we need to evalu-ate a set of amplitudes ⟨Ψ T |ϕ⟩. With a QC trial state, i.e. |Ψ T ⟩ = |Ψ Q ⟩, amplitudes are evaluated on a quantum computer, and a certain number of circuit shots are needed for each amplitude. Note that tens of thousands of walkers and time steps may be required in ph-AFMC [21], which leads to a large total number of amplitude queries.
In this section, we propose an alternative way of using the QC trial state, i.e. QAEE algorithms, to minimize the QC cost. For clarity, we use |Ψ Q ⟩ to denote the QC trial state, which is prepared and measured on the quantum computer, and |Ψ C ⟩ to denote a CC trial state, which can be efficiently evaluated on a classical computer. Some examples of QC and CC trial states of different QMC algorithms can be found in Table 1.
In QAEE algorithms, we only use the QC trial state in the energy evaluation phase. First, we take the trial state |Ψ T ⟩ = |Ψ C ⟩ to generate the expression of |Ψ⟩ in Eq. (2) entirely on the classical computer. Second, we evaluate the energy in Eq. (1) with the quantum computer, in which we take |Ψ T ⟩ = |Ψ Q ⟩. Therefore, the final energy becomes Here we have taken θ l = 0 in oscillatory-signfree algorithms. The factor ⟨ΨQ|ϕl⟩ ⟨Ψ C |ϕ l ⟩ may have a nonzero phase, which can potentially cause the phase problem. However, because both |Ψ Q ⟩ and |Ψ C ⟩ are approximations to the ground state, i.e. the difference between them is small, we find that the phase problem is negligible in numerical simulations of a H 4 linear chain: in VMC and fn-GFMC the average phases are 0.999995 and 0.999996, respectively 1 . In the case that the problem is severe, we can use the phaseless approximation similar to ph-AFMC to eliminate the phase problem.
The optimal way to evaluate energy is using the importance sampling method. The denominator of Eq. (10) is a linear combination of amplitudes ⟨Ψ Q |ϕ l ⟩. Instead of estimating each amplitude with the same QC resources, we assign more circuit shots to an amplitude if the absolute value 1 The average phase is of its coefficient is larger. We give details of the importance sampling protocol for computing the real part of the denominator using the Hadamardtest circuit (see Sec. 4), and it is similar for the imaginary part, numerator and other quantum circuits: 1. Generate N pairs of {w l , ϕ l } on a classical computer with |Ψ T ⟩ = |Ψ C ⟩ according to VMC, fn-GFMC or ph-AFMC; 2. Randomly draw a walker l with the prob- Compose the Hadamard-test circuit to measure the real part of the amplitude . Implement the circuit for one shot to obtain an outcome η = ±1 of the X measurement. Start again from step 2 and repeat M tot times.
Note that in practice, one does not follow steps 2 to 3 for M tot times, but can repeat step 2 for M tot times and record all the walkers w l , before going to step 3 to evaluate the amplitudes on a quantum computer. In the above protocol, we have followed the Monte Carlo summation method and rewritten the denominator real . For clarity, let η (l) be a random variable that represents the outcome of the Hadamard-test circuit for measuring Re(⟨Ψ Q |ϕ ′ l ⟩), in contrast with η, which represents the outcome for a randomly chosen l. As a property of Hadamardtest circuits, η (l) is an unbiased estimator of , the distribution ofB is binomial (up to a transformation); and the probability of η = 1 is (B/C B + 1)/2. Therefore, the variance ofB is If |Ψ Q ⟩ is closer to the ground state than |Ψ C ⟩, we expect that the error in E with respect to the ground-state energy is smaller after replacing Besides, in some cases when |⟨Ψ Q |Φ 0 ⟩| 2 and |⟨Ψ C |Φ 0 ⟩| 2 are close to each other, where |Φ 0 ⟩ is the ground state, the final energy from |Ψ Q ⟩ can still be better. This is probably originated from the inherent symmetry projection on the trial state. In Sec. 6, we will discuss the symmetries in fermion systems and give such an example. We will show that in quantum-assisted algorithms, the knowledge of symmetry can be used to reduce errors.
Next, we will demonstrate that the QAEE algorithm can reduce the bias in numerical simulations. Besides the numerical evidence, there are limited theoretical arguments supporting the bias reduction other than the intuitive conjecture: Replacing the trial state in the energy evaluation with a state closer to the true ground state may reduce the bias. This problem will be solved in the QCMCSD algorithm, in which the bias reduction is a rigorous theoretical result under the assumption without statistical error, i.e. the bias is always smaller than in the classical algorithm. We remark that QCMCSD can also reduce the bias in the quantum algorithm that prepares the QC trial state, such as VQE. In general, QAEE and QCMCSD reduce the requirement on quantum computing for observing quantum advantage, i.e. we only need to prepare a QC trial state better than the CC trial state to reduce the bias.

Simulation results
We demonstrate quantum-assisted algorithms numerically with the H 4 linear chain. The molecule has four hydrogen atoms with uniform spacing, and the interatomic distance is 0.74Å. In the STO-3G basis, the number of spin orbitals is eight, which can be encoded into eight qubits using the Jordan-Wigner transformation. We use the J-KSzGHF state [37] as the CC trial state in VMC and fn-GFMC and the UCCSD state (with one Trotter step [44]  UCCSD state, the expected energy can be measured on a quantum computer. In our case, we directly compute the energy by exact simulation. Because the molecule is small, both trial states can reach an accurate ground-state energy after optimization. The purpose of this numerical simulation is to study the impact of QC variance on the potential quantum advantage. Therefore, we don't intend to maximize the performance of classical algorithms. In each case, we only need a pair of CC and QC trial states, such that without QC variance the bias is smaller if we use the QC trial state. With such two trial states, we can illustrate how the QC variance harms the quantum advantage and to what extent our algorithms can reduce its impact. Therefore, we deliberately stop the optimisation of trial states before the actual minimum is reached. We take the UCCSD trial state with an error of 2.096 mE h (milli hartree) and the J-KSzGHF trial state with an error of 6.084 mE h in VMC and 2.068 mE h in fn-GFMC.
We illustrate the numerical result of the classical and quantum-assisted fn-GFMC with an increasing β in Fig. 1. A similar result of ph-AFMC can be found in Appendix B.2. We can find that the energy obtained by the classical algorithm has a larger bias and variance than those obtained by quantum-assisted algorithms, even though J-KSzGHF and UCCSD trial states (re-  Table 2: Errors in the ground-state energy (mE h ) of a H 4 linear chain computed with the classical and consistent quantum-assisted (CQA) algorithms. When circuit shots per amplitude query M is infinite, the energy is calculated with the quantum-computing trial state |Ψ T ⟩ = |Ψ Q ⟩. When M is finite, a Hadamard-test circuit is used when generating the effective trial states, and we calculate the energy with |Ψ T ⟩ = Ψ Q . In each case, the average energy calculated from 100 effective trial states is given in the table, associated with the corresponding standard deviations shown in the bracket. spectively used in classical and quantum-assisted algorithms) have approximately the same errors in energy. Particularly as N, β → ∞, as shown by the dashed (orange) and dash-dotted (green) curves, the QAEE algorithm leads to an energy closer to the exact ground-state energy than the CQA algorithm.
The impact of QC variance in CQA algorithms is summarised in Table 2. Here, we generate the effective trial states with a finite number of circuit shots M , then we compute the corresponding expected energy in VMC and fixed-node energy in fn-GFMC, which correspond to the limit N, β → ∞. Two approaches of amplitude estimation are considered. In the first approach, we use the Hadamard-test circuit to directly estimate the amplitude using the empirical mean estimation, see Sec. 4.1. In the second approach, the Bayesian inference is used, see Sec. 5. For both approaches, we can find that CQA algorithms outperform classical algorithms given a sufficiently large M , and the quantum advantage is observed in the Bayesian inference approach when the empirical mean estimation fails.
In the numerical simulation of QAEE algorithms, we take into account the statistical fluctuation due to both finite number of walkers N and number of circuit shots. The results are shown in Table 3, where we take N = 10 6 and β = 10  Table 3: Errors in the ground-state energy (mE h ) of a H 4 linear chain computed with the classical QMC and quantum-assisted energy evaluation (QAEE) algorithms. We take N = 10 6 walkers in Monte Carlo algorithms. When M tot is finite, we compute the energy with a Hadamard-test circuit and the importance sampling protocol; in this case the energy is random because of the variance in quantum computing. In fn-GFMC, we take ∆β = 10 −3 and the total number of steps n = 10 4 (where the energy stops descending and can be regarded as β → ∞), i.e. we get an energy at β = 10. For each case shown in the table, we compute the average of 100 energy samples and the standard deviation is given in the bracket.
(β can be regarded as ∞ as the energy stops descending). Besides, an infinite and a finite total number of circuit shots M tot are also considered. As the results indicate, with a finite N , the classical algorithm generates an energy close to the N → ∞ case, see Table 2. Besides, for QAEE algorithms, a finite M tot leads to roughly the same energies as in the case with an infinite M tot .
So far, we have given two ways of using a quantum computer in QMC algorithms, i.e. CQA algorithms and QAEE algorithms. In QAEE algorithms, the cost of quantum computing is minimized but still can reduce the bias. In general, there are various ways to tune the QC cost. For example, in fn-GFMC, we can take |Ψ T ⟩ = |Ψ C ⟩ and simulate the imaginary-time evolution on a classical computer, then we switch the trial state to |Ψ T ⟩ = |Ψ Q ⟩ and make use of a quantum computer to simulate the evolution for a relatively short time. In this picture, the QAEE algorithm is the extreme case that the quantum computer is only used in the last step.

Quantum computing of amplitudes
Computing amplitudes ⟨Ψ T |ϕ⟩ is essential in VMC, GFMC and AFMC. As discussed before, all quantities related to the quantum state used in these QMC algorithms can be derived from amplitudes. There are two approaches to evaluate the amplitude with a quantum computer. In the · · · · · · · · · · · · np nq-np first approach, one prepares the trial state |Ψ T ⟩ on the quantum computer, implements shadow tomography [45] and computes amplitudes on the classical computer using tomography data. This approach has been demonstrated in the experiment reported in Ref. [1], which has potential scaling issues as pointed out by the authors. The second approach, also mentioned in the same work, is scalable with the system size: Instead of representing the trial state with tomography data, one estimates the specific amplitude ⟨Ψ T |ϕ⟩ when it is queried. We focus on this nontomography approach, and techniques developed in this paper can be generalised to the tomography approach.
In quantum computing, we can accelerate the amplitude estimation using amplitude amplification and quantum phase estimation, in which the error depends on the circuit depth [46]. Because the accelerated estimation usually uses deep circuits, we focus on the direct amplitude estimation, which is more practical for near-term quantum technologies. In this section, we consider three types of quantum circuits for the amplitude estimation, including the conventional Hadamard test [47], the vacuum-reference protocol [48,49] and the Hartree-Fock-reference protocol. In the following, we assume that states |Ψ T ⟩ and |ϕ⟩ are all normalised. We will show that circuits for VMC and GFMC could be simpler than AFMC.

Hadamard test
Hadamard test is the conventional method for amplitude estimation [47]. The quantum circuit is shown in Fig. 2(a), with n q data qubits and an ancillary qubit. We suppose that unitary operators U T and U ϕ transform the initial state into the trial state and |ϕ⟩, respectively, i.e. |Ψ T ⟩ = U T |0⟩ ⊗nq and |ϕ⟩ = U ϕ |0⟩ ⊗nq . The S µ gate on the ancillary qubit controls the measurement in the X or Y basis by taking µ = 0, 3, respectively. The amplitude is ⟨Ψ T |ϕ⟩ = ⟨X⟩ + i⟨Y ⟩.
To measure the expected value of a Pauli operator P = X, Y , we run the corresponding circuit for M shots. Each circuit run has a measurement outcome η i = ±1. For the empirical mean, the estimator of ⟨P ⟩ is 1

Vacuum reference
We can measure the amplitude without the ancillary qubit in models with the particle number conservation, as proposed in Ref. [48,49]. We consider the Jordan-Wigner transformation for encoding fermions into qubits, and the number of ones (or zeros depending on the scheme) in the qubit state is the particle number [20]. We suppose that the HF state encoded into qubits is |HF ⟩ = |1⟩ ⊗np ⊗ |0⟩ ⊗(nq−np) , where n p is the particle number. The trial state |Ψ T ⟩ and the walker state |ϕ⟩ have the same particle number, then there exist particle-number-preserving unitary operators V T and V ϕ (i.e. [V T , N p ] = [V ϕ , N p ] = 0, where N p is the particle number operator) that can transform the HF state into the trial state and ϕ state, respectively (i.e. |Ψ T ⟩ = V T |HF ⟩ and |ϕ⟩ = V ϕ |HF ⟩). V T depends on the trial state, and V ϕ is formed of Givens rotation gates [50] when ϕ is a Slater determinant. The circuit is shown in Fig. 2(b). By measuring the first qubit in X or Y basis and all the other qubits in |0⟩, we effectively measure X and Y , where P = P ⊗ |0⟩ ⟨0| ⊗nq−1 . Then ⟨Ψ T |ϕ⟩ = ⟨ X⟩ + i⟨ Y ⟩. The detailed derivation is given in Appendix C.1. This protocol can be used for all three QMC algorithms.

Hartree-Fock reference
Now, we propose a protocol that only works for VMC and GFMC but uses shallower quantum cir-cuits. In VMC and GFMC, |ϕ⟩ = |R⟩, where |R⟩ is a basis state. Suppose that the HF state is one of the basis states in R, and each basis state is encoded as a qubit state b 1 , b 2 , . . . , b nq , where b i are binary numbers. This condition is satisfied in the Jordan-Wigner transformation [20]. We take the HF state as the reference state, as its overlap with the trial state is usually large. The trial state is prepared by a unitary |Ψ T ⟩ = U T |0⟩ ⊗nq . Without loss of generality, we assume ⟨Ψ T |HF ⟩ is positive: We can always add a phase factor to |Ψ T ⟩, i.e. take |Ψ T ⟩ ← e iθ |Ψ T ⟩, such that ⟨Ψ T |HF ⟩ is positive, noticing that the quantum state is the same after changing the phase factor.
Besides, we construct a unitary operator that realises the transformation |HF ⟩ = U R |0⟩ ⊗nq and |R⟩ = U R |1⟩ ⊗ |0⟩ ⊗(nq−1) . We suppose that |HF ⟩ = b HF 1 , b HF 2 , . . . , b HF nq and |R⟩ = b 1 , b 2 , . . . , b nq . Without loss of generality, we suppose that b HF i 0 ̸ = b i 0 , and S 1,i 0 is a swap gate on the first and i 0 th qubits. Then, the operator reads where Λ i 0 ,i is a controlled-NOT gate with the control qubit i 0 and target qubit i. We remark that the swap gate is unnecessary in practice, i.e. we only need to swap operations on the first qubit and qubit i 0 when measuring qubits.
The circuit in Fig. 2(c) is used to measure the real and imaginary parts of the amplitude ⟨Ψ T |R⟩ by measuring the first qubit (or qubit i 0 if the swap gate is removed) in X and Y bases and the rest in state |0⟩, i.e. effectively measure ⟨ X⟩ and ⟨ Y ⟩, respectively: where P r = | ⟨Ψ T |HF ⟩ | 2 . For a detailed derivation, refer to Appendix C.2. Compared with the vacuum reference protocol, the circuit used in the Hartree-Fock reference protocol is usually simpler in two aspects. First, it does not need the GHZ state preparation. Second, the transformation V ϕ used in the vacuum-reference protocol for preparing a general SD can be realised with O(n 2 q ) Givens rotation gates. In the Hartree-Fock-reference protocol, the additional transformation U R is formed of O(n q ) controlled-NOT gates. Therefore, the Hartree-Fock reference protocol has potential advantages in the number of multi-qubit gates.

Bayesian inference amplitude estimation
As we have shown, in quantum computing the estimator of ⟨Ψ T |ϕ⟩ has a finite variance, which could have a significant impact on quantumassisted QMC algorithms. First, the variance contributes to the bias. For example, in fn-GFMC the bias is the difference between groundstate energies of H f n and H [22,40]. The sign of the amplitude ⟨Ψ T |R⟩ is crucial for constructing H f n . It is difficult to determine the sign if the variance is large in comparison with the absolute value of the amplitude, and CQA fn-GFMC may lose its advantage because of the randomness in signs. Similarly, the phaseless approximation in ph-AFMC depends on the phase of ⟨Ψ T |ϕ⟩. Second, the amplitude estimation is intensively queried in QMC algorithms, thus reducing the variance with a large number of circuit shots will notably increase the QC cost.
Here we propose a method for evaluating the amplitude ⟨Ψ T |ϕ⟩ based on the Bayesian inference. With this method, we can significantly reduce the QC variance (and the measurement cost) and have a stable advantage over the classical algorithm. We consider the trial states |Ψ Q ⟩ and |Ψ C ⟩ as approximations to the ground state. Therefore, |Ψ C ⟩ contains information about |Ψ Q ⟩. Given ⟨Ψ C |ϕ⟩, we have some prior knowledge about ⟨Ψ Q |ϕ⟩ before evaluating it on a quantum computer, and we can make use of such information to reduce the QC variance.
As different quantum circuits are used to measure the real and imaginary parts of the amplitude, the real and imaginary parts are independent and thus can be estimated separately with the Bayesian inference method. In the following, we focus on the real part, and it is the same for the imaginary part. To simplify the expressions, we introduce notations y 0 ≡ Re(⟨Ψ C |ϕ⟩), y ≡ Re(⟨Ψ Q |ϕ⟩) and x, where x is the measurement outcome of Re(⟨Ψ Q |ϕ⟩) obtained directly on a quantum computer (e.g. with Hadamardtest circuits for M shots). In the spirit of the Bayesian method, both x and y are treated as values of random variables, denoted by X and Y, respectively. To carry out the Bayesian method, we need a prior distribution of Y, which is determined by y 0 .
Given a prior distribution p Y (y), the posterior distribution is For many methods of measuring Re(⟨Ψ Q |ϕ⟩) on a quantum computer, the distribution X|Y only depends on y. Considering the standard Hadamard-test circuit, p X|Y (x|y) is essentially a binomial distribution, i.e. p X|Y (x|y) = M k q k (1− q) M −k , where q = (1−y)/2, k = M (1−x)/2, and M is the number of circuit shots.
The eventual performance of the Bayesian inference depends on the prior. To demonstrate our method, we assume that the prior is a normal distribution centered at y 0 with the standard deviation σ 0 : Y ∼ N (y 0 , σ 0 ). We approximate the binomial distribution p X|Y (x|y) with the normal distribution (which can be done when M is M . For simplification, we can further approximate σ Q with 1−y 2 0 M or its upper bound 1 M . Under these approximations (specifically, we take the upper bound in the following), the standard deviation of X|Y is independent of y, and the posterior distribution Y|X is normal, i.e. Y|X ∼ N (ŷ,σ), whereŷ Note that x → y when the number of samples M → ∞. We takeŷ as the final estimate of y. The benefit of the Bayesian inference method comes from two sides. On the one hand, quantum computing provides a correction to the classical amplitude y 0 . The correction is more evident when the quantum computing result is more certain, i.e. M is larger. Eventually, in the limit M → ∞, the final estimateŷ converges to its true value. On the other hand, when x deviates from y, e.g. in the case that M is small, the classical amplitude y 0 can instead serve as a correction to the quantum estimate to stabilise the final estimateŷ. Consequently, even if M is small, we In the BI estimation, we take σ 2 0 = 10 −5 . The horizontal lines denote the fixed-node energies of the convectional J-KSzGHF state (solid), quantum-computing UCCSD state (dashed), and the true ground-state energy (dotted). (b) The percentage of samples below the fixed-node energy of the J-KSzGHF state in these 100 samples.
expect that quantum-assisted algorithms can still outperform classical algorithms up to some controllable fluctuation.
Numerical simulations are conducted to compare the Bayesian inference estimation [i.e., computeŷ according to Eq. (15)] with the empirical mean estimation (i.e., directly takeŷ = x) when evaluating the amplitude. We implement the CQA fn-GFMC algorithm here as the example, with parameter settings the same as in Table 2.
The results are shown in Fig. 3. We find that the energy from the Bayesian inference method is much more stable than that from the empiricalmean method. When M is small (i.e. the QC vari-ance is large), the quantum-assisted algorithm using the empirical-mean method has a larger bias than the classical algorithm; using the Bayesian inference method, the bias is smaller. When M is large, the energy converges faster to the limit set by the exact QC trial state (i.e. M → ∞) using the Bayesian inference method. Results of VMC are similar as summarised in Table 2.
Using the Bayesian inference method, we can achieve the quantum advantage, i.e. a reduced bias, with a relatively small M . When M = 100 (the minimum number taken in the simulation), the bias in the CQA fn-GFMC is smaller than the classical algorithm with a probability of about 50%. When M = 3200, the bias is reduced by 10% with a probability of 50%. Notice that the energy is variational, and we can repeat the computation to generate a set of energies and choose the lowest energy, as discussed in Sec. 3.1. Then, the 50% success probability means that we can achieve the quantum advantage by repeating the computation two times on average.
Using ⟨Ψ C |ϕ⟩ as the prior guess of ⟨Ψ Q |ϕ⟩, we have assumed that |Ψ C ⟩ is approximately normalised, as the state |Ψ Q ⟩ prepared on the quantum computer is normalised. In AFMC, we can take |Ψ C ⟩ as a normalised SD (e.g. HF state). In VMC and GFMC, we can construct |Ψ C ⟩ by applying a Jastrow factor to a normalised SD. However, the normalisation assumption generally does not always hold, e.g. when |Ψ C ⟩ is a J-KSzGHF state. In this case, we can process the state as follows. We measure ⟨Ψ Q |ϕ⟩ for a selected state ϕ, for instance the HF state, and compute ⟨Ψ C |ϕ⟩. Then, we multiply |Ψ C ⟩ by a factor and take |Ψ ′ C ⟩ = ⟨ΨQ|ϕ⟩ ⟨Ψ C |ϕ⟩ |Ψ C ⟩. Then, we use |Ψ ′ C ⟩ as the prior instead of |Ψ C ⟩ in the Bayesian inference method.

Symmetry projected quantumcomputing trial state
The ground state is usually in a subspace due to the symmetries of the Hamiltonian. Applying a projection operator of the subspace on a state removes its components orthogonal to the subspace, which is irrelevant to the ground state. In general the projection improves the approximation of a trial state to the ground state. In classical computing, we can introduce correlations to mean-field states by applying the symmetry projection, which is used when generating the trial states in VMC and fn-GFMC, e.g. KSzGHF states [37]. In this paper, we propose that in (fully or partially) quantum-assisted QMC algorithms, we can project the QC trial state according to certain symmetries in order to reduce errors.
Verifying the symmetry in quantum computing can be used to reduce error, and it is a commonlyused method in error mitigation [51,52]; however, such a technique usually increases the circuit complexity. Specifically, in a quantum computer, we can directly measure qubits in the computation basis, i.e. the basis of Pauli Z operators. To verify the symmetry in quantum computing, we need to transform the symmetry (i.e. the corresponding observable) to an observable that can be directly measured. Such a transformation is realised using quantum gates, which increases the circuit complexity [32,33]. Depending on the measurement outcome, the state is probabilistically projected onto the subspace or a state orthogonal to the subspace. On the other hand, in quantum-assisted QMC algorithms, we can apply the symmetry projection without additional gates. The proof is straightforward. Because |ϕ⟩ ∈ H S for all ϕ, P |Ψ⟩ = |Ψ⟩. Then According to Theorem 1, the QC trial state is effectively projected if we restrict walker states in the subspace. In the following, we show this restriction on walker states is for free in VMC, GFMC and AFMC. For fermion systems, the symmetries usually include the particle number N p , spin component S z , etc. In VMC and GFMC, it is natural that we choose the basis according to these symmetries. We take the H 2 molecule as an example. In the STO-3G basis, the model has four spin orbitals filled with two electrons. We use A ↑, A ↓, B ↑ and B ↓ to label the four spin orbitals, in which A, B denote two orbitals, and ↑, ↓ denote spins, respectively. The ground state has the S z symmetry, and the component of the total spin in the z direction is zero. Therefore, one electron is in A ↑ and B ↑, and the other electron is in A ↓ and B ↓. The dimension of the entire Hilbert space of four spin orbitals is 16, and the subspace dimension with N p and S z symmetries is four. The four states are |1 A↑ , 0 B↑ , 1 A↓ , 0 B↓ ⟩, |1 A↑ , 0 B↑ , 0 A↓ , 1 B↓ ⟩, |0 A↑ , 1 B↑ , 1 A↓ , 0 B↓ ⟩ and |0 A↑ , 1 B↑ , 0 A↓ , 1 B↓ ⟩, which respect the two symmetries and are sufficient for expressing the ground state. Therefore, in VMC and GFMC, we can take |R⟩ only in these four states. In AFMC, using the Hubbard-Stratonovich transformation, we express the imaginary-time propagator as an integral of one-particle propagators, and we can let the one-particle propagators preserve the S z symmetry of the Hamiltonian, i.e. A operators (see Appendix B) have the symmetry between two spins. With |1 A↑ , 0 B↑ , 1 A↓ , 0 B↓ ⟩ as the initial state (we suppose that it is the HF state), the initial state is in the four-dimensional subspace, then the walker state is always in the subspace because of the symmetry of A operators, i.e. walker states ϕ l are SDs with the N p and S z symmetries.
We can generalise this approach to symmetries other than N p and S z . For example, for the H 2 molecule, we can take |R⟩ only from the three states |1 A↑ , 0 B↑ , 1 A↓ , 0 B↓ ⟩, |0 A↑ , 1 B↑ , 0 A↓ , 1 B↓ ⟩ and 1 √ 2 (|1 A↑ , 0 B↑ , 0 A↓ , 1 B↓ ⟩ + |0 A↑ , 1 B↑ , 1 A↓ , 0 B↓ ⟩), such that the |Ψ Q ⟩ is effectively projected onto the subspace with the N p , S z and S 2 (total spin) symmetries. In many models, e.g. molecules without an external magnetic field, the Hamiltonian is real, and the ground state is real, which is called the complex conjugation symmetry. If |Ψ⟩ is real, such as in VMC and fn-GFMC, we can impose the complex conjugation symmetry by modifying Eq. (1) to The inherent symmetry projection effectively improves the quality of the QC trial state. It turns out that even in the case CC and QC trial states have roughly the same error in energy (or have approximately the same state fidelity with respect to the ground state), the latter can lead to a smaller error in the final ground state energy, see Fig. 1 (and Appendix B.2) for the numerical example. We remark that symmetry may not be the only reason for this phenomenon. To illustrate the impact of symmetry projection, we decompose CC and QC trial states with approximately the same state fidelity (which are used in the simulation presented in Appendix B.2) into eigenstates |Φ i ⟩ of the Hamiltonian of a H 2 molecule and plot amplitudes in Fig. 4. For simplification and better illustration, we plot the absolute value of the amplitude | ⟨Ψ T |Φ i ⟩ |. From the left side, we see that the CC trial state has non-zero components only on two eigenstates, while the QC trial state has a more uniform distribution except for the ground state (the top bar). The right side shows amplitudes | ⟨ϕ|Φ i ⟩ | of five walker states which we picked randomly in the simulations. We find walker states have nonzero components only on three eigenstates due to symmetries. Because of the symmetries of walker states, the QC trial state is effectively projected onto the same three eigenstates, and components on other eigenstates are effectively removed. After the projection, the fidelity of the QC trial state is increased. This may explain why with a similar state fidelity the QC trial state can lead to a smaller error in the ground-state energy.
7 Quantum-classical Monte Carlo subspace diagonalisation In CQA algorithms, the energy as the result is variational. Because of this reason, the quantum advantage, i.e. a reduced bias, is verifiable: If the energy is lower than those from classical VMC and fn-GFMC algorithms, the bias of quantum-assisted algorithms is smaller. However, ph-AFMC and QAEE algorithms do not have this property, and it could be the case that the energy is even lower than the true groundstate energy; in this case a lower energy may not be better. Therefore, we would like to ask the question, is there a way to reduce the bias that is systematic, verifiable and generalisable to other advanced classical algorithms? Diagonalisation in a subspace is one of the most important methods for computing extreme eigenvalues. In classical computing, it is standard to choose the Krylov subspace generated by applying the Hamiltonian power on a reference state [53]. This approach has been generalised to quantum computing, in which a set of states prepared on the quantum computer spans the subspace [54][55][56]. A crucial property of subspace diagonalisation is that, according to the Rayleigh quotient theorem, the lowest eigenenergy is the lowest expected energy in the subspace. In other words, given the best solution to the ground state from classical computing, if we construct a subspace containing this classical solution and an optimal state from quantum computing, the diagonalisation in this subspace always yields an energy closer to the true ground-state energy. This property is exactly what we need.
To carry out the subspace diagonalisation involving optimal classical and quantum solutions, we need to prepare the classical-solution state on the quantum computer. For applications on nearterm quantum devices, it is essential to reduce the gate number. In the following, we show that using the hybrid QMC method the state preparation requires O(n q ) gates, where n q is the qubit number.
In VMC, GFMC and AFMC algorithms, we express the approximation to the ground state |Ψ⟩ as a linear combination of walker states |ϕ⟩, see Eq.
(2). As we have discussed in this paper, we can evaluate quantities ⟨Ψ Q |Ψ⟩ and ⟨Ψ Q | H |Ψ⟩ without physically preparing the state |Ψ⟩, i.e. we only need to estimate amplitudes ⟨Ψ Q |ϕ⟩. This feature is particularly useful for problems where |Ψ⟩ is harder to prepare than |ϕ⟩. Such states include solutions of state-of-the-art methods, for instance coupled cluster [38] and tensor network [13,39], which can be prepared on a quantum computer: Preparing a coupled-cluster state requires O(n 4 ) gates [57], and it is similar for a tensor-network state [58]. In VMC and fn-GFMC, to express |Ψ⟩ as a linear combination of basis states to estimate ⟨Ψ Q |Ψ⟩ [the denominator of Eq. (10)], preparing ϕ only requires at most n q single-qubit gates when ϕ is a basis state. Therefore, in the framework of QMC, the circuit complexity can be reduced. It is similar for In the hybrid Monte Carlo subspace diagonalisation, we construct the subspace as follows. Let |Ψ C ⟩ and |Ψ Q ⟩ be optimal solutions from classical and quantum computing, respectively. We can use VQE algorithms to find a quantum solution. The classical solution is from, for instance, QMC, coupled-cluster or tensor-network algorithms. For example, thinking of VMC, we have found the optimal parameters λ ⋆ , and the corresponding final solution is |Ψ C ⟩ = |Ψ T (λ ⋆ )⟩; Together with the exact expression, VMC also yields a walker-state representation of |Ψ T (λ ⋆ )⟩, see Eq. (6). In GFMC, the state closest to the true ground state is e −n∆βH |Ψ I ⟩. Instead of an exact expression of this state, fn-GFMC outputs a walker-state representation |Ψ⟩ as an approximation to e −n∆βH |Ψ I ⟩: ϕ l and w l in Eq. (2) are generated according to Algorithm 2, where θ l = 0. It is similar for ph-AFMC. For coupledcluster and tensor-network states, we can also work out the corresponding walker-state representation following Refs. [38,39]. In general, any classical-solution state with amplitudes in a basis can be efficiently calculated, and we can work out a walker-state representation according to VMC. Overall, we suppose that there is a walker-state representation |Ψ⟩ of |Ψ C ⟩ in the form of Eq. (2). The subspace is the span of two states |Ψ⟩ and |Ψ Q ⟩. This approach can be directly generalised to multiple classical and quantum solutions.
The subspace diagonalisation works as follows. Let |Φ 1 ⟩ = |Ψ⟩ and |Φ 2 ⟩ = |Ψ Q ⟩. We need to work out the matrix elements H i,j = ⟨Φ i | H |Φ j ⟩ and S i,j = ⟨Φ i |Φ j ⟩, where i, j = 1, 2. Elements H 1,1 and S 1,1 can be calculated on the classical VMC fn-GFMC Classical algorithms (N = 10 6 ) 6.042 1.345 Subspace diagonalisation 1.761 0.861 Table 4: Errors in the ground-state energy (mE h ) of a H 4 linear chain computed with the classical Monte Carlo and quantum-classical Monte Carlo subspace diagonalisation algorithms. In the classical algorithms, we take N = 10 6 and the same trial states as in Table 3.
computer. Because |Ψ Q ⟩ is physically prepared on the quantum computer, it is normalized, and S 2,2 = 1; H 2,2 can be evaluated on the quantum computer. For off-diagonal elements, they can be evaluated as the same as ⟨Ψ Q | H |Ψ⟩ and ⟨Ψ Q |Ψ⟩ in QAEE algorithms. Given two matrices, we solve the generalised eigenvalue problem Hx = ESx, where x is a column vector [59,60]. The minimum eigenvalue E is the lowest expected energy of all states (including |Ψ C ⟩ and |Ψ Q ⟩) in the subspace.
To demonstrate the QCMCSD algorithm, we consider VMC and fn-GFMC for classical solutions and VQE with UCCSD for the quantum solution, see Sec. 3.3. We remark that the optimisation in VMC and VQE can be further pushed to produce better solutions than what we use in the numerical simulation, and these pseudo-optimal solutions are sufficient for the preliminary demonstration. The result is summarised in Table 4. Expressed as a linear combination of 10 6 walkers according to VMC, the energy of |Ψ⟩ has an error of 6.042 mE h . |Ψ Q ⟩ is a UCCSD state with an error of 2.096 mE h . The diagonalisation in the subspace Span({|Ψ⟩ , |Ψ Q ⟩}) results in a minimum eigenvalue with the error of 1.761 mE h . The result is similar for fn-GFMC, where we take |Ψ⟩ as a walker-state representation of the fixed-node ground state. This approach can also be applied to the output |Ψ⟩ of ph-AFMC.
In the numerical result, we find that the QCMCSD algorithm is less accurate than other quantum-assisted algorithms. However, it is of particular interest. As far as we know, there is no rigorous theory showing that utilising a QC trial state can reduce the bias with certainty, although the bias reduction is observed in experiments [1] and numerics. According to the theory of subspace diagonalisation, the QCMCSD algo-rithm always results in a smaller bias than the classical algorithm, under the assumption of no statistical error.

Conclusion
In summary, this work introduces a set of quantum-assisted Monte Carlo algorithms. Combining quantum computing with QMC, the advanced classical computational techniques, one can reduce the bias when constraining the sign problem in computing the ground-state energy. This form of quantum advantage is achieved by utilising a trial state prepared on the quantum computer. The potential challenge for hybrid Monte Carlo algorithms of this kind is the cost for measuring trial state amplitudes. We propose these algorithms with the purpose to ease this issue. Our methodology is to control the involvement of the quantum trial state, specifically in two ways. First, we use both the quantum trial state and a trial state of classical computing in different stages of the computation. We use the quantum trial state in the entire computation in CQA algorithms and only at the last stage in QAEE algorithms, and intermediate cases are possible. Second, we take the quantum trial state as an adaptable correction (depending on the cost budget) to the classical trial state in the Bayesian inference amplitude estimation. Additionally, applying post-selections in CQA algorithms lowers the requirement for the precision of amplitude measurement. In this work, we numerically demonstrate the performance of our algorithms with hydrogen molecules. Overall, theoretical and numerical results suggest that we can avoid the high-precision amplitude measurement to attain a bias reduction.
The hybrid Monte Carlo framework provides new tools for quantum computing. We can reduce errors in a quantum trial state with the inherent symmetry projection, which is almost free in Monte Carlo algorithms. We conjecture that this symmetry projection can also mitigate machine noise in quantum computing. The hybrid Monte Carlo is a gate-efficient interface between quantum computing and classical computing to leverage both computation paradigms. In this paper, we propose the QCMCSD algorithm as an example. QMC is a numerical method developed for decades and encompasses a family of variants.
In this view, hybrid Monte Carlo offers a powerful suite of techniques to explore to demonstrate practical quantum advantage.

Algorithm 1 Energy evaluation in the variational Monte Carlo algorithm.
1: Input H, Ψ T , N . 2: for l = 1 to N do 3: Generate R according to the probability distribution |⟨Ψ T |R⟩| 2 .

4:
Compute the local energy E l = E loc (R) according to Eq. (4).
Generate R according to the probability distribution A −1 ⟨Ψ T |R⟩ ⟨R|Ψ I ⟩.

9:
Generate R according to the probability distribution A −1 S f n (R, R ′ ).

A Green's function Monte Carlo
Given the basis, the imaginary-time Green's function reads Considering the first-order expansion, we define In the limit n → ∞, e −n∆βH |Ψ I ⟩ converges to the ground state if |Ψ I ⟩ has a finite overlap with the ground state. It is similar for the operator F = 1 − ∆βH. Eigenstates of the Hamiltonian are eigenvectors of F , and the eigenvalue of the ground state is 1 − ∆βE g . If 1 − ∆βE g is the largest absolute eigenvalue, F n |Ψ I ⟩ also converges to the ground state in the limit n → ∞. If ∆β∥H∥ 2 ≤ 1, F is positive semi-definite, and the condition always holds. Note that sometimes it is helpful to add a proper constant to the Hamiltonian, i.e. replace H with H − E 0 , where E 0 is the constant. According to F n |Ψ I ⟩, we compute the ground-state energy by evaluating In the last line, a similarity transformation is applied to F (R, R ′ ), and Algorithm 3 Phaseless auxiliary-filed Monte Carlo algorithm. w l,0 ← 1 ▷ Initialise the walker.
B Auxiliary-filed Monte Carlo

B.1 Theory
For many fermion models, such as molecules, the Hamiltonian can be expressed in the form where A j = A j,1 1 + k,q A j,p,q a † p a q are one particle operators, and a q is the fermion annihilation operator of the qth spin orbital. The evolution time is divided into n small time steps, i.e. β = n∆β. Using the Trotter-Suzuki formula and Hubbard-Stratonovich transformation, the time evolution operator of each time step is rewritten as where x = (x 1 , x 2 , . . . , x L ) denotes the auxiliary-filed, is the normal distribution, are one-particle propagators, and E 0 is a constant taken a value close to the ground-state energy. Then we compute the ground-state energy by evaluating = dx 1 · · · dx n p(x 1 ) · · · p(x n )E loc (Ψ n )I(x n ,x n , Ψ n−1 ) · · · I(x 1 ,x 1 , Ψ 0 ) dx 1 · · · dx n p(x 1 ) · · · p(x n )I(x n ,x n , Ψ n−1 ) · · · I(x 1 ,x 1 , Ψ 0 ) + O(n∆β 2 ). (30) where the importance function is and |Ψ 0 ⟩ ≡ |Ψ I ⟩. We take eachx according tō and |ϕ⟩ = |Ψ k−1 ⟩ forx k . In ph-AFMC, we approximate the importance function with a non-negative number. Specifically, in the phaseless approximation, we take I(x,x, ϕ) ≈ e −∆β[ReE loc (ϕ)−E 0] × max 0, cos arg ⟨ψ T | B(x −x) |ϕ⟩ ⟨ψ T |ϕ⟩ .
Essential steps of ph-AFMC are given in Algorithm 3.  Classical QAEE-better QAEE-worse Exact Figure 5: The energy (mE h ) of a H 2 molecule computed using the classical and the quantum-assisted energy evaluation ph-AFMC algorithms. We take N = 2 × 10 5 and ∆β = 0.01. The yellow and green curves represent the cases with lower and higher state fidelities (with respect to the true ground state) of quantum computing, respectively. The exact ground-state energy is plotted as the red dotted curve. A zoomed-in image of the range of β from around 2 to 4 is shown as the inserted graph for better illustration.

B.2 Numerical simulation results
state of H 2 and H 4 , corresponding to an error of 20.525mE h and 40.99mE h in the ground state energy, respectively. For the quantum-assisted algorithm, we use two UCCSD trial states for H 2 , with a state fidelity of 98.72% and 99.41%, respectively, which correspond to an error of 12.001mE h and 4.733mE h ; for H 4 , we use a UCCSD trial state with a fidelity of 98.38% and an error of 19.885mE h . The phaseless approximation is used in the simulation to eliminate the sign problem. The result of the energy of a H 2 molecule is shown in Fig. 5. As β increases, all the three curves stabilise to approach the exact ground state, while the blue curve (representing the result from the classical algorithm) is observed to have a larger bias and variance than the two other curves (representing results from the quantum-assisted algorithm with different trial states), indicating that the quantum-assisted algorithm generates a better result. Similar results are observed for H 4 , as summarised in Table 5, which shows the error in the ground state energy. Compared with the results for VMC and fn-GFMC, we see the error from both the classical and quantum-assisted ph-AFMC algorithms is smaller, as is usually the case for classical algorithms; that's because in AFMC, the walker state is a general Slater determinant close to the ground state and the local energy fluctuates around the actual energy. On the other hand, such a feature also makes AFMC more challenging to implement for both classical and quantum-assisted algorithms.

C.1 The vacuum reference method
In this section we show the circuit in Fig. 2(b) can compute the real and imaginary parts of the amplitude ⟨Ψ T |ϕ⟩.