Efficient and Noise Resilient Measurements for Quantum Chemistry on Near-Term Quantum Computers

TL;DR

This paper introduces a measurement strategy for quantum chemistry on near-term quantum computers that significantly reduces measurement time and improves noise resilience through low rank tensor factorization and postselection techniques.

Contribution

A novel low rank tensor factorization method that reduces measurement complexity and enhances noise resilience in variational quantum algorithms for molecular systems.

Findings

Measurement times reduced by three orders of magnitude.

Elimination of sampling challenges with non-local operators.

Effective error mitigation via postselection.

Abstract

Variational algorithms are a promising paradigm for utilizing near-term quantum devices for modeling electronic states of molecular systems. However, previous bounds on the measurement time required have suggested that the application of these techniques to larger molecules might be infeasible. We present a measurement strategy based on a low rank factorization of the two-electron integral tensor. Our approach provides a cubic reduction in term groupings over prior state-of-the-art and enables measurement times three orders of magnitude smaller than those suggested by commonly referenced bounds for the largest systems we consider. Although our technique requires execution of a linear-depth circuit prior to measurement, this is compensated for by eliminating challenges associated with sampling non-local Jordan-Wigner transformed operators in the presence of measurement error, while enabling a powerful form of error mitigation based on efficient postselection. We numerically characterize these benefits with noisy quantum circuit simulations for ground state energies of strongly correlated electronic systems.

Tables

Table 1. Table 1: A history of ideas reducing the measurements required for estimating the energy of arbitrary basis chemistry Hamiltonians with the variational quantum eigensolver. Here N 𝑁 N represents the number of spin-orbitals in the basis. Gate counts and depths are given in terms of arbitrary 1- or 2-qubit gates restricted to the geometry of 2-qubit gates specified in the connectivity column. What we mean by “compatible” Pauli groupings is that the terms can be measured at the same time with only single qubit rotations prior to measurement. We report whether terms are measured in a diagonal representation as this is important for enabling strategies of error-mitigation by postselection. The number of partitions refers to the number of unique term groupings which can each be measured with a single circuit - thus, this reflects the number of unique circuits required to generate at least one sample of each term in the Hamiltonian. However, we caution that one cannot infer the total number of measurements required from the number of partitions, and often this metric is highly misleading. The overall number of measurements required is also critically determined by the variance of the estimator of the energy. As explained in the first entry of this table, when terms are measured simultaneously one must also consider the covariance of those terms. In some cases, a grouping strategy can decrease the number of partitions but increase the total number of measurements required by grouping terms with positive covariances. Alternatively, strategies such as the third entry in this table actually increase the number of partitions while reducing the number of measurements required overall by lowering the variance of the estimator.

Ref.	Partitioning Method	Circuit Description	# of Partitions	Gate Count	Depth	Connectivity	Diagonal
McClean et al. (2016)	commuting Pauli heuristic	-	$O\left(N^4\right)$	-	-	-	-
Kandala et al. (2017)	compatible Pauli heuristic	single rotations	$O\left(N^4\right)$	$N$	1	any	no
Rubin et al. (2018)	$n$-representability constraints	single rotations	$O\left(N^4\right)$	$N$	1	any	no
Izmaylov et al. (2019)	mean-field partitioning	fast feed-forward	$O\left(N^4\right)$	$O\left(N\right)$	$O\left(N\right)$	full	no
Verteletskyi et al. (2020)	compatible Pauli clique cov.	single rotations	$O\left(N^4\right)$	$N$	$1$	any	no
O’Gorman et al. (2019)	counting argument	swap networks	$O\left(N^3\right)$	$O\left(N^2\right)$	$O\left(N\right)$	linear	no
Jena et al. (2019)	commuting Pauli graph color.	stabilizer formalism	$O\left(N^3\right)$	-	-	full	no
Izmaylov et al. (2020)	anticommuting Pauli clique cov.	Pauli evolutions	$O\left(N^3\right)$	$O\left(N^2\log \left(N\right)\right)$	-	full	no
Yen et al. (2019)	commuting Pauli clique cover	symplectic subspaces	$O\left(N^3\right)$	$O\left(N^2/\log N\right)$	-	full	no
Gokhale et al. (2019)	commuting Pauli clique cover	stabilizer formalism	$O\left(N^3\right)$	$O\left(N^2\right)$	-	full	no
here	integral tensor factorization	Givens rotations	$O\left(N\right)$	$N^2/4$	$N/2$	linear	yes

Table 2. Table 2: List of the molecular systems considered in this work, displayed in order of increasing number of qubits, for each type of system. The hydrogen systems consist of a chain of atoms arranged in a line, with equal interatomic spacing. The interatomic spacing for the water molecules refers to the length of the symmetrically stretched bonds O-H bonds, which are separated by a fixed angle of 104.5 ​ deg 104.5 degree 104.5\deg . The active space used for each system has one spatial orbital for every two qubits. A non-zero number of frozen orbitals indicates the number of molecular orbitals fixed in a totally occupied state.

System	Interatomic	Basis	Frozen	Number of
	Spacings (Å)	Set	Orbitals	Qubits
H2	.6, .7, …1.3	STO-3G	None	4
H2	.6, .7, …1.3	6-31G	None	8
H4	.6, .7, …1.3	STO-3G	None	8
H6	.6, .7, …1.3	STO-3G	None	12
H4	.6, .7, …1.3	6-31G	None	16
H8	.6, .7, …1.3	STO-3G	None	16
H2	.6, .7, …1.3	cc-pVDZ	None	20
H10	.6, .7, …1.3	STO-3G	None	20
H6	.6, .7, …1.3	6-31G	None	24
H2O	.8, .9, …1.5	STO-3G	1	12
H2O	.8, .9, …1.5	STO-3G	None	14
H2O	.8, .9, …1.5	6-31G	1	24
N2	.9, 1.0, …1.6	STO-3G	2	16
N2	.9, 1.0, …1.6	STO-3G	None	20

Table 3. Table 3: Bounds and uncertainties result from Bayesian inference using a Monte-Carlo approximation with 10 6 superscript 10 6 10^{6} particles for all Hydrogen FCI data  Granade  et al. ( 2012 ) . We assume log ⁡ ( N meas ) = log ⁡ ( a ) + x ^ + b ​ log ⁡ ( N ) subscript 𝑁 meas 𝑎 ^ 𝑥 𝑏 𝑁 \log(N_{\rm meas})=\log(a)+\hat{x}+b\log(N) where x ^ ∼ 𝒩 ​ ( 0 , 0.1 ) similar-to ^ 𝑥 𝒩 0 0.1 \hat{x}\sim\mathcal{N}(0,0.1)   end . The prior distributions are uniform for log ⁡ ( a ) 𝑎 \log(a) and b 𝑏 b over [ − 15 , 1 ] 15 1 [-15,1] and [ 1 , 20 ] 1 20 [1,20] respectively. Here σ ​ ( b ) 𝜎 𝑏 \sigma(b) is the posterior standard deviation for b 𝑏 b and Δ ​ ( a ) Δ 𝑎 \Delta(a) is the posterior standard deviation of log ⁡ ( a ) + x ^ 𝑎 ^ 𝑥 \log(a)+\hat{x} . ”RDM Constraints” refers to the Pauli Word Grouping approach with the RDM constraints applied, as in the text.

Measurement Strategy	$⟨\log \left(a\right)⟩$	$\Delta \left(a\right)$	$⟨\text{b}⟩$	$\sigma \left(b\right)$
Bound from Qubit Hamiltonian	-6.0	0.3	4.90	0.02
Separate Measurements	-9.3	0.4	5.70	0.06
Pauli Word Grouping	-8.9	0.4	4.88	0.06
RDM Constraints	-10.8	0.4	5.63	0.06
Basis Rotation Grouping	-6.0	0.3	2.75	0.01

Table 4. Table 4: Consider the normal ordered second quantized quantum chemistry Hamiltonian of Eq.  ( 11 ) , calculated for a chain of eight hydrogen atoms equally spaced 1.0 1.0 1.0 Å apart in an STO-3G basis. We group the terms in this Hamiltonian into five partitions. Partitions I and II contain the one particle terms from the first summation. Partition I consists of those terms where p = q 𝑝 𝑞 p=q , while II consists of those where p ≠ q 𝑝 𝑞 p\neq q . Partitions III, IV, and V contain the two particle terms from the second summation. Partition III consists of those where there are two unique values among p, q, r, and s, while IV consists of those with three unique values and V consists of those with four eigenvalues. For each partition, we report the sum of the absolute values of the coefficients of these terms in the fermionic representation of the Hamiltonian (counting the coefficient of a term and its Hermitian conjugate only once). We also report the same quantity calculated in the qubit representation after applying the Jordan-Wigner transformation. We drop any constant terms which appear as a result of the Jordan-Wigner transformation, since these do not contribute to the variance. Additionally, we report the sum of the absolute value of the coefficients for the entire Hamiltonian calculated in both ways in the final column.

Hamiltonian Partition	I	II	III	IV	V	Whole
						Hamiltonian
Example Term	$a_p^{†}{a}_p$	$a_p^{†}{a}_q+h.c.$	$a_p^{†}{a}_q^{†}{a}_p{a}_q$	$a_p^{†}{a}_q^{†}{a}_q{a}_r+h.c.$	$a_p^{†}{a}_q^{†}{a}_r{a}_s+h.c.$	-
Fermionic ${\sum }_{\ell }\left|w_{\ell }\right|$	32.288	2.852	34.214	7.436	35.579	112.368
Qubit ${\sum }_{\ell }\left|w_{\ell }\right|$	16.144	2.852	25.660	6.794	17.790	33.500

Table 5. Table 5: Variances for a symmetrically stretched chain of 8 8 8 hydrogen atoms in an STO-3G basis. Rows 1-5 show values of the variance bounds and approximations to these that are described in the text. The variances are presented in units of 100 E 2 h superscript subscript absent ℎ 2 {}_{h}^{2} . The bound in row 1 is calculated using Eq.  ( 9 ) , while row 2 uses Eq.  ( 16 ) , and row 3 uses Eq.  ( 12 ) . The approximations in rows 4-5 are calculated using the methodology of Ref.  3 , which amounts to using Eq.  ( 9 ) or Eq.  ( 12 ) but neglecting some of the terms in the Hamiltonian as described in the text below. The last two rows, 6-7, present the actual variance of an estimator that measures each term in the Jordan-Wigner transformed Hamiltonian separately, for the Hartree-Fock state and for the ground state, respectively.

Interatomic Spacing (Angstrom)	.6	.7	.8	.9	1.0	1.1	1.2	1.3
Qubit Variance Bound	22.499	17.657	14.680	12.675	11.222	10.134	9.297	8.629
Fermionic Variance Bound	99.393	86.271	76.352	68.683	62.596	57.721	53.768	50.462
Naive Fermionic Variance Bound	205.251	177.997	156.704	139.847	126.267	115.234	106.168	98.545
Fermionic Variance Approximation	31.501	27.404	24.580	22.557	21.038	19.893	19.023	18.321
Qubit Variance Approximation	7.212	6.225	5.567	5.118	4.786	4.540	4.356	4.209
Hartree-Fock Variance	7.211	6.224	5.565	5.118	4.785	4.539	4.355	4.208
Ground State Variance	9.206	8.194	7.568	7.181	6.929	6.779	6.695	6.641

Equations

$$M\le {\left(\frac{{\sum }_{\ell }\left|{\omega }_{\ell }\right|}{ϵ}\right)}^2,\text{where}H=\sum _{\ell }{\omega }_{\ell }{P}_{\ell }$$

$$H=U_0\left(\sum _p{g}_p{n}_p\right)U_0^{†}+\sum _{\ell =1}^L{U}_{\ell }\left(\sum _{pq}{g}_{pq}^{(\ell )}{n}_p{n}_q\right)U_{\ell }^{†},$$

$$U=\text{exp}\left(\sum _{pq}{\kappa }_{pq}{a}_p^{†}{a}_q\right),U{a}_p^{†}{U}^{†}$$

$$⟨H⟩=\sum _p{g}_p{⟨n_p⟩}_0+\sum _{\ell =1}^L\sum _{pq}{g}_{pq}^{(\ell )}{⟨n_p{n}_q⟩}_{\ell },$$

$$f_{\ell }=\frac{\left|{\omega }_{\ell }\right|\sqrt{1-{⟨H_{\ell }⟩}^2}}{{\sum }_j\left|{\omega }_j\right|\sqrt{1-{⟨H_j⟩}^2}}.$$

$${⟨P_{\ell }⟩}_{\mathrm{bitflip}}=(1-2p{)}^K{⟨P_{\ell }⟩}_{\mathrm{true}},$$

$${⟨H⟩}_{\mathrm{proj}}=\frac{\mathrm{Tr}\left(P\rho H\right)}{\mathrm{Tr}\left(P\rho \right)},$$

$${⟨H⟩}_{\mathrm{proj}}=\frac{\mathrm{Tr}\left(P\rho H\right)}{\mathrm{Tr}\left(P\rho \right)}=\frac{\mathrm{Tr}\left(\rho H\right)+\mathrm{Tr}\left(\rho \Pi H\right)}{1+\mathrm{Tr}\left(\rho \Pi \right)}.$$

$$M_q\le {\left(\frac{{\sum }_{\ell }\left|{\omega }_{\ell }\right|}{ϵ}\right)}^2,\text{where}H=\sum _{\ell }{\omega }_{\ell }{P}_{\ell }.$$

$$M={\left(\frac{{\sum }_{\ell }\left|{\omega }_{\ell }\right|{\sigma }_{\ell }}{ϵ}\right)}^2.$$

$$H=\sum _{pq}{h}_{pq}{a}_p^{†}{a}_q+\sum _{pqrs}{h}_{pqrs}{a}_p^{†}{a}_q^{†}{a}_r{a}_s,$$

$$M_f\le {\left(\frac{{\sum }_{p,q\le p}|h_{pq}|+{\sum }_{pqr,s<q}|h_{pqrs}|+{\sum }_{pq,r\le p,s=q}|h_{pqrs}|}{ϵ}\right)}^2.$$

$$h_{pqrs}=h_{qpsr}=h_{rspq}=h_{srqp}=h_{rqps}=h_{qrsp}=h_{psrq}=h_{spqr}=h.$$

$$2h{a}_4^{†}{a}_2^{†}{a}_3{a}_1+h.c.+2h{a}_3^{†}{a}_2^{†}{a}_4{a}_1+h.c.,$$

$$\scriptsize\begin{split}\frac{-h}{4}\big{(}X_{1}X_{2}X_{3}X_{4}+X_{1}X_{2}Y_{3}Y_{4}-X_{1}Y_{2}X_{3}Y_{4}+X_{1}Y_{2}Y_{3}X_{4}+Y_{1}X_{2}X_{3}Y_{4}-Y_{1}X_{2}Y_{3}X_{4}+Y_{1}Y_{2}X_{3}X_{4}+Y_{1}Y_{2}Y_{3}Y_{4}\big{)}+\\ \frac{-h}{4}\big{(}X_{1}X_{2}X_{3}X_{4}+X_{1}X_{2}Y_{3}Y_{4}+X_{1}Y_{2}X_{3}Y_{4}-X_{1}Y_{2}Y_{3}X_{4}-Y_{1}X_{2}X_{3}Y_{4}+Y_{1}X_{2}Y_{3}X_{4}+Y_{1}Y_{2}X_{3}X_{4}+Y_{1}Y_{2}Y_{3}Y_{4}\big{)}=\\ \frac{-h}{4}\big{(}X_{1}X_{2}X_{3}X_{4}+X_{1}X_{2}Y_{3}Y_{4}+Y_{1}Y_{2}X_{3}X_{4}+Y_{1}Y_{2}Y_{3}Y_{4}+X_{1}X_{2}X_{3}X_{4}+X_{1}X_{2}Y_{3}Y_{4}+Y_{1}Y_{2}X_{3}X_{4}+Y_{1}Y_{2}Y_{3}Y_{4}\big{)}.\end{split}$$

$$M_f\le {\left(\frac{\frac{1}{2}{\sum }_p|h_p|+{\sum }_{p,q<p}|h_{pq}|+{\sum }_{pqr,s<q}|h_{pqrs}|+{\sum }_{pq,r<p,q}|h_{pqrq}|+\frac{1}{2}{\sum }_{pq}|h_{pqpq}|}{ϵ}\right)}^2,$$

$$H=U_0\left(\sum _p{g}_p{n}_p\right)U_0^{†}+\sum _{\ell =1}^L{U}_{\ell }\left(\sum _{pq}{g}_{pq}^{(\ell )}{n}_p{n}_q\right)U_{\ell }^{†},$$

$$H=\sum _{\sigma \in \{\uparrow ,\downarrow \}}\sum _{pq}{T}_{pq}{a}_{p,\sigma }^{†}{a}_{q,\sigma }+\frac{1}{2}\sum _{\alpha ,\beta \in \{\uparrow ,\downarrow \}}\sum _{pqrs}{V}_{pqrs}{a}_{p,\alpha }^{†}{a}_{q,\alpha }{a}_{r,\beta }^{†}{a}_{s,\beta }.$$

$$V_{pqrs}=V_{srqp}=V_{psqr}=V_{qprs}=V_{qpsr}=V_{rsqp}=V_{rspq}=V_{srpq},$$

$$V_{pqrs}=\sum _{\ell }{w}_{\ell }{v}_{pq}^{(\ell )}{v}_{rs}^{(\ell )}$$

$$\frac{1}{2}\sum_{\alpha,\beta\in\{\uparrow,\downarrow\}}\sum_{pqrs}V_{pqrs}a^{\dagger}_{p,\alpha}a_{q,\alpha}a^{\dagger}_{r,\beta}a_{s,\beta}=\frac{1}{2}\sum_{\ell}w_{\ell}\Big{(}\sum_{\sigma\in\{\uparrow,\downarrow\}}\sum_{pq}v^{(\ell)}_{pq}a^{\dagger}_{p,\sigma}a_{q,\sigma}\Big{)}^{2},$$

Full text

I Introduction

Given the recent progress in quantum computing hardware, it is natural to ask where the first demonstration of a quantum advantage for a practical problem will occur. Since the first experimental demonstration by Peruzzo et al. Peruzzo et al. (2014), the variational quantum eigensolver (VQE) framework has offered a promising path towards utilizing small and noisy quantum devices for simulating quantum chemistry. The essence of the VQE approach is the use of the quantum device as a co-processor which prepares a parameterized quantum wavefunction and measures the expectation value of observables. In conjunction with a classical optimization algorithm, it is possible to then minimize the expectation value of the Hamiltonian as a function of the parameters, arriving at approximations for the wavefunction, energy, and other properties of the ground state Peruzzo et al. (2014); McClean et al. (2016); Wecker et al. (2015); O’Malley et al. (2016); Kandala et al. (2017); Lee et al. (2018); Parrish et al. (2019); O’Brien et al. (2019). A growing body of work attempting to understand and ameliorate the challenges associated with using VQE to target non-trivial systems has emerged in recent years Bonet-Monroig et al. (2018); McClean et al. (2020); McArdle et al. (2019); Temme et al. (2017); Sagastizabal et al. (2019); McClean et al. (2017a); Otten and Gray (2019); O’Gorman et al. (2019); Verteletskyi et al. (2020); Jena et al. (2019); Yen et al. (2019); Izmaylov et al. (2020, 2019); Gokhale et al. (2019). In this article we address the challenge posed by the large number of circuit repetitions needed to perform accurate measurements and propose a new scheme that dramatically reduces this cost. Additionally, we explain how our approach to measurement has reduced sensitivity to readout errors and also enables a powerful form of error mitigation.

Within VQE, expectation values are typically estimated by Hamiltonian averaging. Under this approach, the Hamiltonian is decomposed into a sum of operators that are tensor products of single-qubit Pauli operators, commonly referred to as Pauli strings. The expectation values of the Pauli strings are determined independently by repeated measurement. When measurements are distributed optimally between the Pauli strings $P_{\ell}$, the total number of measurements $M$ is upper bounded by

[TABLE]

is the Hamiltonian whose expectation value we estimate as $\sum_{\ell}\omega_{\ell}\mathinner{\langle{P_{\ell}}\rangle}$, the $\omega_{\ell}$ are scalars, and $\epsilon$ is the target precision Wecker et al. (2015); Rubin et al. (2018). Prior work assessing the viability of VQE has used bounds of this form and concluded that chemistry applications require “a number of measurements which is astronomically large” (quoting from Ref. 3).

Several recent proposals attempt to address this obstacle by developing more sophisticated strategies for partitioning the Hamiltonian into sets of simultaneously measurable operators Izmaylov et al. (2019); O’Gorman et al. (2019); Verteletskyi et al. (2020); Jena et al. (2019); Yen et al. (2019); Izmaylov et al. (2020). We summarize their key findings in Table 1. This work has a similar aim, but we take an approach rooted in a decomposition of the two-electron integral tensor rather than focusing on properties of Pauli strings. We quantify the performance of our proposal by numerically simulating the variances of our term groupings to more accurately determine the number of circuit repetitions required for accurate measurement of the ground state energy. This contrasts with the analysis in other recent papers that have instead focused on using the number of separate terms which must be measured as a proxy for this quantity. By that metric, our approach requires a number of term groupings that is linear in the number of qubits - a quartic improvement over the naive strategy and a cubic improvement relative to these recent papers. However, we argue that the number of distinct term groupings alone is not generally predictive of the total number of circuit repetitions required, because it does not consider how the covariances of the different terms in these groupings can collude to either reduce or increase the overall variance. We will show below that our approach benefits from having these covariances conspire in our favor; for the systems considered here our approach gives up to three orders of magnitude reduction in the total number of measurements, while also providing an empirically observed asymptotic improvement.

Although there are a variety of approaches to simulating indistinguishable fermions with distinguishable qubits Setia and Whitfield (2018); Jiang et al. (2018); Bravyi and Kitaev (2002), the Jordan-Wigner transformation is the most widely used. This is due to its simplicity and to the fact that it allows for the explicit construction of a number of useful circuit primitives not available under more sophisticated encodings. These include the Givens rotation network that exactly implements a change of single-particle basis Babbush et al. (2018); O’Gorman et al. (2019); Kivlichan et al. (2018); Motta et al. (2018). A disadvantage of using the Jordan-Wigner transformation is the fact that it maps operators acting on a constant number of fermionic modes to qubit operators with support on up to all $N$ qubits. In the context of measurement, the impact of this non-locality can be seen by considering a simple model of readout error such as a symmetric bitflip channel. Under this model, a Pauli string with support on $N$ qubits has $N$ opportunities for an error that reverses the sign of the measured value, leading to estimates of expectation values that are exponentially suppressed in $N$ (see Section II.3). It has recently been shown that techniques based on fermionic swap networks can avoid the overheads and disadvantages imposed by the non-locality of the Jordan-Wigner encoding in a variety of contexts, including during measurement Kivlichan et al. (2018); Motta et al. (2018); O’Gorman et al. (2019). Our work will likewise avoid this challenge without leaving the Jordan-Wigner framework, allowing estimation of 1- and 2-particle fermionic operator expectation values by the measurement of only 1-local and 2-local qubit operators, respectively.

In addition to this reduction in the support of the operators that we measure, our work offers another opportunity for mitigating errors. It has been observed that when one is interested in states with a definite eigenvalue of a symmetry operator, such as the total particle number, $\eta$, or the $z$-component of spin, $S_{z}$, it can be useful to have a method which removes the components of some experimentally prepared state with support on the wrong symmetry manifold McArdle et al. (2019); Bonet-Monroig et al. (2018); O’Brien et al. (2019); McClean et al. (2020). Two basic strategies to accomplish this have been proposed. The first of these strategies is to directly and “non-destructively” measure the symmetry operator and discard those outcomes where the undesired eigenvalue is observed, projecting into the proper symmetry sector by postselection. In order to construct efficient measurement schemes, prior work in this direction has focused on measuring the parities of $\eta$ and $S_{z}$, rather than the full symmetry operators, McArdle et al. (2019); Bonet-Monroig et al. (2018). These proposals involve non-local operations that usually require $O(N)$ depth, which may induce further errors during their implementation. The second class of strategies builds upon the foundation of Ref. 14 and uses additional measurements together with classical post-processing to calculate expectation values of the projected state without requiring additional circuit depth Bonet-Monroig et al. (2018); O’Brien et al. (2019); McClean et al. (2020), a procedure which can be efficiently applied to the parity of the number operator in each spin sector. In this work, we show how our proposal for measurement naturally leads to the ability to postselect directly on the proper eigenvalues of the operators $\eta$ and $S_{z}$, rather than on their parities.

II Results

II.1 Using Hamiltonian Factorization for Measurements

The crux of our strategy for improving the efficiency and error resilience of Hamiltonian averaging is the application of tensor factorization techniques to the measurement problem. Using a representation discussed in the context of quantum computing in Refs. 30; 31; 29, we begin with the factorized form of the electronic structure Hamiltonian in second quantization:

[TABLE]

where the values $g_{p}$ and $g_{pq}^{(\ell)}$ are scalars, $n_{p}=a^{\dagger}_{p}a_{p}$, and the $U_{\ell}$ are unitary operators which implement a single particle change of orbital basis. Specifically,

[TABLE]

where $[e^{\kappa}]_{pq}$ is the $p,q$ entry of the matrix exponential of the anti-Hermitian matrix $\kappa$ that characterizes $U$.

Numerous approaches that accomplish this goal exist, including the density fitting approximation Whitten (1973); Aquilante et al. (2010), and a double factorization which begins with a Cholesky decomposition or eigendecomposition of the two-electron integral tensor Pedersen et al. (2009); Motta et al. (2018); Beebe and Linderberg (1977); Koch et al. (2003); Aquilante et al. (2010); Purwanto et al. (2011); Mardirossian et al. (2018); Peng and Kowalski (2017). In this work, we use such an eigendecomposition and refer readers to Appendix C and to Refs. 29 and 31 for further details. The eigendecomposition step permits discarding small eigenvalues to yield a controllable approximation to the original Hamiltonian. While such low rank truncations are not central to our approach and would not significantly reduce the number of measurements, doing so would asymptotically reduce $L$ (and thus ultimately, the number of distinct measurement term groupings). Such decompositions have been explored extensively in the context of electronic structure on classical computers on a far wider range of systems than those considered here Røeggen and Wisløff-Nilssen (1986); Koch et al. (2003); Røeggen and Johansen (2008); Boman et al. (2008); Pedersen et al. (2009); Peng and Kowalski (2017). It has been found that $L=O(N)$ is sufficient for the case of arbitrary basis quantum chemistry, both in the large system and large basis set limits Pedersen et al. (2009). Furthermore, specific basis sets exist where $L=1$, such as the plane wave basis or dual basis of Ref. 27.

Our measurement strategy, which we shall refer to as “Basis Rotation Grouping,” is to apply the $U_{\ell}$ circuit directly to the quantum state prior to measurement. This allows us to simultaneously sample all of the $\mathinner{\langle{n_{p}}\rangle}$ and $\mathinner{\langle{n_{p}n_{q}}\rangle}$ expectation values in the rotated basis. We can then estimate the energy as

[TABLE]

where the subscript $\ell$ on the expectation values denotes that they are sampled after applying the basis transformation $U_{\ell}$. The reason that the $\mathinner{\langle{n_{p}}\rangle}_{\ell}$ and $\mathinner{\langle{n_{p}n_{q}}\rangle}_{\ell}$ expectation values can be sampled simultaneously is because under the Jordan-Wigner transformation, $n_{p}=(\openone+Z_{p})/2$, which is a diagonal qubit operator. In practice, we assume a standard measurement in the computational basis, giving us access to measurement outcomes for all diagonal qubit operators simultaneously. Thus, our approach is able to sample all terms in the Hamiltonian with only $L+1=O(N)$ distinct term groups.

Fortunately, the $U_{\ell}$ are exceptionally efficient to implement, even on hardware with minimal connectivity. Following the strategy described in Ref. 28, and assuming that the system is an eigenstate of the total spin operator, any change of single-particle basis can be performed using $N^{2}/4-N/2$ two qubit gates and gate depth of exactly $N$, even with the connectivity of only a linear array of qubits Kivlichan et al. (2018). This gate depth can actually be improved to $N/2$ by further parallelizing the approach of Kivlichan et al. (2018), making using ideas that are explained in the context of multiport interferometry in Clements et al. (2016). In fact, a further optimization is possible by performing the second matrix factorization discussed in Ref. 29. This would result in only $O(\log^{2}N)$ distict values of the $g_{pq}^{(\ell)}$ and a gate complexity for implementing the $U_{\ell}$ which is reduced to $O(N\log N)$; however, we note that this scaling is only realized in fairly large systems when $N$ is growing towards the thermodynamic (large system) rather than continuum (large basis) limit.

The primary objective of our measurement strategy is to reduce the time required to measure the energy to within a fixed accuracy. Because different hardware platforms have different repetition rates, we focus on quantifying the time required in terms of the number of circuit repetitions. We shall present data for electronic ground states that demonstrates the effectiveness of our Basis Rotation Grouping approach in comparison to three other measurement strategies and the upper bound of Eq. (9). All calculations were performed using the open source software packages OpenFermion and Psi4 McClean et al. (2017b); Parrish et al. (2017). Specifically, we used exact calculations of the variance of expectation values with respect to the full configuration interaction (FCI) ground state to determine the number of circuit repetitions required. The calculations presented here are performed for symmetrically stretched Hydrogen chains with various bond lengths and numbers of atoms, for a symmetrically stretched water molecule, and for a stretched Nitrogen dimer, all in multiple basis sets. We justify our focus on the electronic ground states here by noting that most variational algorithms for chemistry attempt to optimize ansatze that are already initialized near the ground state. For reference, we provide analogous data calculated with respect to the Hartree-Fock state in the supplementary materials (see Appendix D).

In order to calculate the variance of the estimator of the expectation value of the energy, it is necessary to determine the distribution of measurements between the different term groupings. Refs. 3 and 23 provide a prescription for the optimal choice. They demand that (in the notation of Eq. (9)) each term $H_{\ell}$ is measured a fraction of the time $f_{\ell}$ equal to

[TABLE]

In practice, the expectation values in the above expression are not known ahead of time and so the optimal measurement fractions $f_{\ell}$ cannot be efficiently and exactly determined a priori. For the purposes of this paper, we approximated the ideal distribution of measurements by first performing a classically tractable configuration interaction singles and doubles (CISD) calculation of the quantities in Eq. (5). We shall show that this approximation introduces a negligible overhead in measurement time for all systems considered in this work. One could also envisage using an adaptive measurement scheme that makes additional measurements based on the observed sample variance, in order to approximate the ideal partitioning of measurement time, such as the one described in Ref 46.

II.2 Circuit Repetitions Required for Energy Measurement

In Figure 1 we plot the number of circuit repetitions for our proposed “Basis Rotation Grouping” measurement approach (black circles), together with three other measurement strategies and the upper bound based on Eq. (9) for the systems listed in Table 2. The first and most basic alternative strategy is simply to apply no term groupings and measure each Pauli string independently, a strategy we refer to as “Separate Measurements” (lime green circles). A more sophisticated approach, similar to the one described in Ref. 17, is to partition the Pauli strings into groups of terms that can be measured simultaneously. In the context of a near-term device we consider two Pauli strings $P_{j}$ and $P_{k}$ simultaneously measurable if and only if they act with the same Pauli operator on all qubits on which they both act non-trivially. Pauli Strings that satisfy this condition can be simultaneously measured using only single qubit rotations and measurement. In order to efficiently partition the Pauli strings into groups we choose to take all of the terms which only contain $Z$ operators as one partition and then account for the remaining Pauli words heuristically by adding them at random to a group until no more valid choices remain before beginning a new group. We refer to this approach as “Pauli Word Grouping” (teal circles). The final strategy that we compare with preprocesses the Hamiltonian by applying the techniques based on the fermionic marginal (RDM) constraints described in Ref. 23, before applying the Jordan-Wigner transformation and using the same heuristic grouping strategy to group simultaneously measurable Pauli strings together 111In the process of preparing this manuscript we have become aware of several recent works that employ more sophisticated strategies for grouping Pauli words together or employing a different family of unitary transformations than those we consider to enhance the measurement process Verteletskyi et al. (2020); Jena et al. (2019); Yen et al. (2019); Izmaylov et al. (2020). It would be an interesting subject of future work to calculate and compare the number of circuit repetitions required by these approaches.. We call this latter strategy “Pauli Word Grouping, RDM Constraints” (dark blue circles).

We refer to the bound of Eq. (9) as being based on the Hamiltonian coefficients and calculate it from the Jordan-Wigner transformed Hamiltonian, (meaning that the $\omega_{\ell}$ in Eq. (9) are the coefficients of Pauli strings). This bound is indicated by salmon-colored circles in Figure 1. We note that attempting to calculate a similar bound directly from the fermionic Hamiltonian (meaning that the $\omega_{\ell}$ in Eq. (9) would be the coefficients of the terms $a^{\dagger}_{p}a_{q}$ or $a^{\dagger}_{p}a_{q}^{\dagger}a_{r}a_{s}$) leads to different bounds. These are derived in Appendix A, where they are shown to be substantially looser for the systems we consider in this work. While one would not measure the fermion operators directly, it is surprising that these bounds would be significantly different. We refer the interested reader to Appendix A for an analysis and discussion of this phenomenon.

Considering first the Hydrogen chain systems in Figure 1 (left panel, A), we note that our Basis Rotation Grouping approach consistently outperforms the other strategies for simulations with more than four fermionic modes, requiring significantly fewer measurements. Interestingly, while the bounds from the qubit Hamiltonian and other three methods appear to have relative performances that are stable across a variety of system sizes, the Basis Rotation Grouping method appears to have a different asymptotic scaling, at least for Hydrogen chains of increasing length and basis set size. This is likely due to large scale effects that only manifest when approaching a system’s thermodynamic limit (which one approaches particularly quickly for Hydrogen chains) on the Many-Electron Problem et al. (2017). In Table II we quantify this asymptotic scaling by assuming that the dependence of the variance on the number of qubits $N$ in the Hydrogen chain’s Hamiltonian can be modeled by the functional form $aN^{b}$ for some constants $a$ and $b$ which we fit using a Bayesian analysis described in the table caption Granade et al. (2012). By contrast, the data from the minimal basis water molecule (panel B in Figure 1) shows no benefit in measurement time from our method compared to the heuristic grouping strategies. However, the advantage of our approach becomes significant for that system in larger basis sets, a trend which is also apparent to a lesser extent for the Nitrogen dimer (panel C in Figure 1).

We find that applying the RDM constraints of Ref. 23 to our Pauli Word Grouping strategy (the combination is plotted with dark blue circles in Figure 1) does not significantly reduce the observed variance, despite the fact that the use of the RDM constraints have been previously shown to dramatically reduce the bounds on the number of circuit repetitions required Rubin et al. (2018). In Appendix B, we explore the possibility that this is due to the fact that these constraints were applied to minimize a bound of the same form as Eq. (9) that is however formulated using the fermionic representation of the Hamiltonian. We present evidence in Appendix A that, in the context of such bounds, the use of the Jordan-Wigner transformed operators leads to surprisingly different results. However, as we show in Appendix B, we find that the actual variance with respect to the ground state is not substantially changed by applying the same constraints and performing the minimization using the qubit representation of the Hamiltonian.

Earlier we explained that the data presented in Figure 1 was calculated by distributing the measurements between different term groupings according Eq. (5) using the variance of each term calculated with a classically efficient CISD approximation to the ground state. Any deviation from the ideal allocation of measurement cycles (obtained by evaluating Eq. (5) with respect to the true ground state) must increase the time required for measurement. In Figure 2 we present the ratio between the time required with the approximate distribution and the time required under the optimal one for each of the systems treated in the work. We find that impact from this approximation is negligible, with the largest observed increase in measurement time being below $3\%$.

Overall, Figure 1 speaks for itself in showing that in most cases there is a very significant reduction in the number of measurements required when using our strategy - sometimes by up to three orders of magnitude for even modestly sized systems. Furthermore, these improvements become more significant as system size grows.

II.3 Error Mitigation

Beyond the reduction in measurement time, our approach also provides two distinct forms of error mitigation. First, it reduces the susceptibility to readout errors by replacing the measurement of $O(N)$ qubit operators with 1 and 2 qubit operators. Second, it allows us to perform postselection based on the eigenvalues of the particle number operators in each spin sector. Both properties stem from measuring the Hamiltonian only in terms of density operators in different basis sets.

The first benefit, the reduction in readout errors, is a consequence of only needing to measure expectation values of operators that have support on one or two qubits. Direct measurement of the Jordan-Wigner transformed Hamiltonian using only single-qubit rotations and measurement involves measuring operators with support on $O(N)$ qubits. To demonstrate how reducing the support of the operators helps to mitigate errors, we consider a simple model of measurement error: the independent, single-qubit symmetric bitflip channel. When estimating the expectation value of a Pauli string $P_{\ell}$ acting on $K$ qubits with a single-qubit bitflip error rate $p$, a simple Kraus operator analysis shows that $P_{\ell}$ is modified to

[TABLE]

which means that the noise channel will bias the estimator of the expectation value towards zero by a factor exponential in $K$. Thus, the determination of expectation values is highly sensitive to the extent of locality of the $P_{\ell}$, a behavior that we expect to persist under more realistic models of readout errors.

One could also accomplish the reduction in the support of the operators that our method achieves by other means. For example, one could measure each of the $O(N^{4})$ terms separately, localizing each one to a single qubit operator by applying $O(N)$ two-qubit gates. Other schemes have been proposed which simultaneously allow generic two-electron terms to be measured using $O(1)$ qubits each while simultaneously accomplishing the parallel measurement of $O(N)$ terms at a time, at the cost of using $O(N^{2})$ or $O(N^{2}log(N))$ two-qubit gates Izmaylov et al. (2020); Gokhale et al. (2019); O’Gorman et al. (2019). One unique advantage of our approach is that we achieve this reduction in operator support at the same time as the large reduction in the number of measurement repetitions presented in Sec. II.2 above.

Our approach also enables a second form of error mitigation. Each measurement we prescribe is also simultaneously a measurement of the total particle number operator, $\eta$, and of the $z$ component of spin, $S_{z}$. We can therefore reduce the impact of circuit and measurement errors by performing postselection conditioned on a desired combination of quantum numbers for each of these operators. Let $P$ denote the projector onto the corresponding subspace and let $\rho$ denote the density matrix of our state. We obtain access to the projected expectation value,

[TABLE]

directly from the experimental measurement record by discarding those data points which fall outside of the desired subspace. The remaining data points are used to evaluate the expectation values of the desired Pauli strings.

This postselection is efficient in the sense that it require no additional machinery beyond what we have already proposed. The only cost is a factor of $\approx 1/{\rm Tr}(P\rho)$ additional measurements. This factor is approximate because discarding measurements with the wrong particle number is likely to lead to a lower observed variance. Specifically, by removing measurements in the wrong particle number sector, we avoid having to average over large fluctuations caused by the energetic effects of adding or removing particles. This therefore presents an additional route by which our Basis Rotation Grouping scheme will reduce the number of measurements in practice.

Several recent works have proposed error mitigation strategies which allow for the targeting of specific symmetry sectors. We make a brief comparative review of these here in order to place our work in context. One class of strategies focuses on non-destructively measuring one or more symmetry operators McArdle et al. (2019); Bonet-Monroig et al. (2018). After performing the measurements and conditioning on the desired eigenvalues, the post-measurement state becomes $P\rho P/{\rm Tr}(P\rho)$ and the usual Hamiltonian averaging can be performed. These approaches share some features with our strategy in that they also require an additional number of measurements that scale as $1/{\rm Tr}(P\rho)$ and an increased circuit depth. However, they also have some drawbacks that we avoid. Because they separate the measurement of the symmetry operator from the measurement of the Hamiltonian they require the implementation of relatively complicated non-destructive measurements. As a consequence, existing proposals focus on measuring only the parity of the $\eta$ and $S_{z}$ operators, leading to a strictly less powerful form of error mitigation than the approach we propose. Additionally, most errors that occur during or after the symmetry operator measurement are undetectable, including errors incurred during readout.

A different class of approaches avoids the need for additional circuit depth at the expense of requiring more measurements Bonet-Monroig et al. (2018); O’Brien et al. (2019); McClean et al. (2020). To understand this, let $\Pi$ denote the fermionic parity operator and $P=(1+\Pi)/2$ the projector onto the $+1$ parity subspace. Then,

[TABLE]

To construct the projected energy it then suffices to measure the expectation values of the Hamiltonian, the parity operator, and the product of the Hamiltonian and parity operators. A stochastic sampling scheme and a careful analysis of the cost of such an approach reveals that it is possible to use post-processing to estimate the projection onto the subspace with the correct particle number parity in each spin sector at a cost of roughly $1/{\rm Tr}(P_{\uparrow}P_{\downarrow}\rho)^{2}$ (where $P_{\uparrow}$ and $P_{\downarrow}$ are the parity projectors for the two spin sectors) McClean et al. (2020). Unlike our approach, this class of error mitigation techniques does not easily allow for the projection onto the correct eigenvalues of $\eta$ and $S_{z}$, owing to the large number of terms required to construct these projection operators. Furthermore, the scaling in the number of additional measurements we described above, already more costly than our approach, is also too generous. This is because the product of the parity operators and the Hamiltonian will contain a larger number non-simultaneously measurable terms than the same Hamiltonian on its own. Maximum efficiency may require grouping schemes that consider this larger number of term groupings.

The most significant drawback of our method in the context of error mitigation is that the additional time and gates required for the basis transformation circuit lead to additional opportunities for errors. We believe the reduction in circuit repetitions we have shown makes our method the most attractive choice when it is feasible to use an additional $O(N^{2})$ two-qubit gates during the measurement process. We therefore focus on comparing the performance of our strategy with a strategy that requires no additional gates and uses a quantum subspace error mitigation approach that effectively projects onto the correct parity of the number operator on each spin sector Bonet-Monroig et al. (2018); McClean et al. (2020). In order to do so, we use the open source software package Cirq Developers (2019) to simulate the performance of both strategies for measuring the ground state energy of a chain of six Hydrogen atoms symmetrically stretched to 1.3Å in an STO-3G basis. We take an error model consisting of i) applying a single qubit depolarizing channel with some probability to both qubits following each two qubit gate, and ii) applying a bitflip channel during the measurement process with some other probability. We report results for a wide range of gate and readout noise levels inspired by the capabilities of state of the art superconducting and ion trap quantum computers Kjaergaard et al. (2020); Bruzewicz et al. (2019); Heinsoo et al. (2018); Preskill (2018). Specfically, we consider single qubit depolarizing noise with probabilities ranging from $2.5\times 10^{4}$ to $8\times 10^{3}$ and single qubit bitflip error probabilities between $6.25\times 10^{4}$ and $1\times 10^{2}$. Here we do not consider the effect of a finite number of measurements and instead report the expectation values from the final density matrix.

Figure 3 shows the error in the measurement of the ground state energy for the error-mitigated Basis Rotation Grouping (far right panel) and Pauli Word Grouping (second panel from right) approaches together with the expectation values for both measurement strategies without error mitigation (two left panels). In these calculations we assumed that the ground state wavefunction under the Jordan-Wigner transformation is prepared without error. Circuit level noise is considered only during the execution of the Givens rotation required for our Basis Rotation Grouping approach. In order to include the impact of our proposed error mitigation strategy on state preparation as well as measurement, we have also carried out calculations including circuit noise during state preparation. The results of these calculations are presented in Figure 4. Here we have approximated a realistic state preparation circuit by applying three random basis rotations which compose to the identity to the ground state wavefunction. These state preparation circuits are simulated with the same gate noise as the measurement circuits. This choice is motivated by the assumption that low-depth circuits will be required for the successful application of VQE and the expectation that $90$ two-qubit gates represents a reasonable lower bound to the size of circuit for a strongly-correlated problem on $12$ qubits.

Figures 3 and 4 show that the Pauli Word Grouping and Basis Rotation Grouping approaches to measurement benefit significantly from their respective error mitigation strategies. Despite the fact that our proposed Basis Rotation Grouping technique requires $30$ additional two-qubit gates compared to the Pauli Word Grouping approach, we see that the errors remaining after mitigation are comparable in some regimes and are lower for our strategy when noise during measurement is the dominant error channel (compare the bottom right corners of the two rightmost panels in both figures). Focusing first on Figure 3, we can see that this is true even when the errors during state preparation are not taken into account. Examining the left two panels of both figures, we can see that even without applying postselection,the locality of our Jordan-Wigner transformed operators leads to a considerable benefit in suppressing the impact of readout errors.

We note that the absolute errors we find when including noise during state preparation (Figure 4), even at the lowest noise levels considered here, are larger than the usual target of “chemical accuracy” ($\sim$ 1 mHa). In practice, an experimental implementation of VQE on non-trivial systems will require the combination of multiple forms of error mitigation. Prior work has shown that error mitigation by symmetry projection combines favorably with proposals to extrapolate expectation values to the zero noise limit McArdle et al. (2019). We expect that such an extrapolation procedure could significantly improve the numbers we present here. Other avenues for potential improvements are also available. For example, one could rely on the error mitigation and efficiency provided by our measurement strategy during the outer loop optimization procedure, before utilizing a richer quantum subspace expansion in an attempt to reduce errors in the ground state energy after determining the optimal ansatz parameters.

III Discussion

We have presented an improved strategy for measuring the expectation value of the quantum chemical Hamiltonian on near-term quantum computers. Our approach makes use of well studied factorizations of the two-electron integral tensor, in order to rewrite the Hamiltonian in a form which is especially convenient for measuring under the Jordan-Wigner transformation. By doing so, we obtain $O(N)$ distinct sets of terms which must be measured separately, instead of the $O(N^{4})$ required by a naive counting of terms approach. Application to specific molecular systems show that in practice, we require a much smaller number of repetitions to measure the ground state energy to within a fixed accuracy target. For example, assuming an experimental repetition rate of 10 kHz (consistent with the capabilities of commercial superconducting qubit platforms), a commonly referenced bound based on the Hamiltonian coefficients suggests that approximately $55$ days are required to estimate the ground state energy to of a symmetrically stretched chain of 6 Hydrogen atoms encoded as a wavefunction on 24 qubits to within chemical accuracy, while our approach requires only $44$ minutes. Our proposed measurement approach also removes the susceptibility to readout error caused by long Jordan-Wigner strings and allows for postselection by simultaneously measuring the total particle number and $S_{z}$ operators with each measurement shot.

The tensor factorization that we used to realize our measurement strategy is only one of a family of such factorizations. Future work might explore the use of different factorizations, or even tailor the choice of single particle bases for measurement to a particular system, by choosing them with some knowledge of the variances and covariances between terms in the Hamiltonian. As a more concrete direction for future work, the data we show in Appendix A regarding the difference between the bounds when calculated directly from the fermionic operators and the same approach applied to the Jordan-Wigner transformed operators, suggests that the cost estimates for error-corrected quantum algorithms should be recalculated using the qubit Hamiltonian.

For the largest systems we consider in this work, the 24 qubit Hydrogen chain and water simulations, and the 20 qubit Nitrogen calculations, our numerical results indicate that using our approach results in a speedup of more than an order of magnitude when compared to recent state-of-the-art measurement strategies. Furthermore, we observe a speedup of more than three orders of magnitude compared to the bounds commonly used to perform estimates in the literature. We also present strong evidence for an asymptotic improvement in our data on Hydrogen chains of various sizes. We performed detailed circuit simulations that show that reduction in readout errors combined with the error mitigation enabled by our work largely balances out the requirement for deeper circuits, even when compared against a moderately expensive error mitigation strategy based on the quantum subspace expansion Bonet-Monroig et al. (2018). We expect that the balance of reduced measurement time and efficient error mitigation provided by our approach will be useful in the application of variational quantum algorithms to more complex molecular systems.

Finally, we note that these techniques will generally be useful for quantum simulating any fermionic system, even those for which the tensor factorization cannot be truncated, such as the Sachdev-Ye-Kitaev model of many-body chaotic dynamics Kitaev (2015); Babbush et al. (2019). In that case, $L$ will attain its maximal value of $N^{2}$, and our scheme will require $N^{2}+1$ partitions. Likewise, if the goal is to use the basis rotation grouping technique to estimate the fermionic 2-particle reduced density matrix rather than just the energy, one would need to measure in all $O(N^{2})$ bases.

Acknowledgements

The authors thank Dominic Berry for his insight that the techniques of Clements et al. (2016) can be used to halve the gate depth of the basis rotation circuits presented in Kivlichan et al. (2018). KBW was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Quantum algorithm Teams Program, under contract number DE-AC02-05CH11231.

Ethics Declarations

The authors declare no competing interests.

Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code Availability

Much of the code that support the findings of this study is already available in the OpenFermion library McClean et al. (2017b). The remainder is available from the corresponding author upon reasonable request.

Appendix A Variance Bounds

In the text of the main article, we reviewed the standard approach to upper bounding the number of measurements $(M)$ required to measure the energy of a Hamiltonian $(H)$ to within a desired target precision $\epsilon$ Wecker et al. (2015); Rubin et al. (2018). This bound is calculated in a straightforward way from the qubit representation that is obtained by making the Jordan-Wigner transformation on a second quantized representation of the fermionic states, by expressing $H$ as the sum of Pauli strings (products of single-qubit Pauli operators) $P_{\ell}$ acting on the qubit representation of the state. Then we have

[TABLE]

More generally, if the Hamiltonian is expressed as a linear combination $H=\sum_{\ell}\omega_{\ell}O_{\ell}$, one can work out the optimal way of distributing independent measurements between these terms and the overall number of measurements required for the resulting estimator to attain a target precision. We refer the reader to Ref. 23 for more details and simply recall the expression here,

[TABLE]

The notation is the same as above, except that $\sigma_{\ell}$ is the positive square root of the variance of the operator $O_{\ell}$. The upper bound of Eq. (9) is derived by noting that the variance of a Pauli operator measurement is at most one and by performing the appropriate substitutions. Our primary concern here is to show how the calculation of such a bound for fermionic Hamiltonians depends in a subtle manner on the representation of the Hamiltonian. We have denoted the bounds above $M_{q}$, to refer to their evaluation in the qubit representation of the Hamiltonian.

When calculating an upper bound of this type for a quantum chemical Hamiltonian,

[TABLE]

it might seem natural to work directly with the coefficients $h_{pq}$ and $h_{pqrs}$ and the fermionic representation rather than performing the Jordan-Wigner transformation to the qubit representation. Provided that one is careful to count the coefficient for each term and its Hermitian conjugate only once, it is possible to obtain an upper bound to the number of measurements required directly from the coefficients $h_{pq}$ and $h_{pqrs}$. Provided that the Hamiltonian is normal-ordered first, this bound can be expressed as

[TABLE]

We denote this bound by $M_{f}$ to indicate it was derived in the fermionic representation of the Hamiltonian. However, this bound is looser than necessary in multiple ways. In Table 4 we show a breakdown of the calculation of the sum of the absolute values of the coefficients for the Hamiltonian of a chain of eight equally spaced hydrogen atoms. We consider five different types of terms from Eq. (11) and calculate the sum of the absolute values of the coefficients for all terms of each type in both the qubit and fermionic representations. By comparing the two approaches in this way we show below that we can shed light on the difference in the resulting bounds.

In column I of Table 4 we begin with this analysis for the ‘number operator’ terms, $a^{\dagger}_{p}a_{p}$ of partition I. We see that the value reported in the qubit representation is exactly half of that reported in the fermionic one. This is because the Jordan-Wigner transformation applied to $a^{\dagger}_{p}a_{p}$ yields $\frac{1}{2}\mathcal{I}+\frac{1}{2}Z_{p}$. However, we may neglect the first contribution since this is a constant and so does not affect the variance. Another way of understanding this difference between the qubit and fermionic contributions is to realize that the bound of Eq. (9) is derived with the assumption that the maximum variance of each term is $1$. The number operator, however, has eigenvalues [math] and $1$ rather than $-1$ and $1$ like a Pauli operator, and so its maximum variance is lower. We shall present an alternative to Eq. (12) below in Eq. (16) that accounts for this lower variance. The analysis of the other one-body term, $a_{p}^{\dagger}a_{q}+h.c.$, from partition II, is simpler. Column II shows that the part of the total magnitude of the coefficients of these terms is the same, regardless of which representation is used for the calculation.

The two-body terms in the Hamiltonian display more varied behavior. Taking a term from partition III and applying the Jordan-Wigner transformation results in exactly $\frac{1}{4}$ of the weight being assigned to a constant term, explaining the difference between the values for this partition. Because terms such as $a^{\dagger}_{p}a^{\dagger}_{q}a_{p}a_{q}$ can be rewritten as the product of two number operators, they must have eigenvalues [math] and $1$ and thus, a maximum variance smaller than one. This improved bound would actually be lower than the one suggested by the analysis of the Jordan-Wigner transformed terms, because it would account for the covariances between the $Z_{p}$, $Z_{q}$, and $Z_{p}Z_{q}$ terms that emerge. We shall incorporate this tighter bound on the variance of the individual terms in this class into the alternative to Eq. (12) presented below as Eq. (16).

The disparity for partitions IV and V has a different source. If one performs the Jordan-Wigner transformation of a term from either class individually, there is no difference between the total magnitudes of the coefficients for the fermionic operators and their qubit counterparts. However, when one sums all such terms together there is some cancellation between the qubit terms that reduces the overall total magnitude. Specifically, the terms in class V benefit from some cancellation due to the eight-fold symmetries of the two-electron integral tensor for real orbitals Szabo2012-ag,

[TABLE]

We claim that the cancellations between symmetric terms result in a value for the sum of the magnitudes of the coefficients that is exactly half as large in the qubit representation as it is in the fermionic one. As an example, consider the case with $p=4,q=2,r=3,s=1$. After normal-ordering, the eight terms become four and we have

[TABLE]

where $h$ is given by Eq. (13) and denotes the value of the coefficients before normal-ordering. The Jordan-Wigner transformation leads to terms from $a^{\dagger}_{4}a^{\dagger}_{2}a_{3}a_{1}+h.c.$ that cancel with terms from $a^{\dagger}_{3}a^{\dagger}_{2}a_{4}a_{1}$ as shown below:

[TABLE]

It is straightforward, although tedious, to prove that the same cancellation occurs generically for terms in class V as a consequence of this eight-fold symmetry. As a result, for this class of terms the sum of the magnitudes of the coefficients is exactly half as large in the qubit representation as it is in the fermionic one. Analogous cancellations in the sum of the class IV terms do not show an obvious symmetry but they are also the source of the difference between the contributions from the two representations in partition IV.

Further cancellation is also apparent when one combines all five classes and calculates the sum of the absolute values of the coefficients (dropping the constant terms) for both representations of the Hamiltonian. The sum of magnitudes of the individual classes in the fermionic representation is the same as the magnitude of the sum. However, in the qubit representation, this value calculated for the entire Hamiltonian is roughly half the size of the sum of the individual partitions. One substantial reason for this difference is the fact that the terms in partitions I and III naturally give rise to terms proportional to products of single qubit $Z$ operators having opposite signs. This is behavior that we should expect for any molecular Hamiltonian, where the single number operator terms arise from the Coulomb attraction between nuclei and individual electrons (negative sign), while the terms containing two number operators arise from the Coulomb repulsion between pairs of electrons (positive sign). Furthermore, unlike the tighter bounds achievable by accounting for the smaller variance of the terms in partitions I and III, this cancellation derives from the underlying form of the Hamiltonian and can not be accounted for in a straightforward way by using a better upper bound for the $\sigma_{\ell}$ values in Eq. (10) when deriving a fermionic bound like Eq. (12).

In the first three rows of Table 5, we now tabulate the bounds on the variance of the energy estimator, in units of 100 $E_{h}^{2}$, that arise from the sums of the absolute values of the coefficients. We perform these calculations for a chain of eight hydrogen atoms at various symmetric stretched interatomic spacings, including the $1.0$Å distance explored in Table 4. In this table, ‘Qubit Variance Bound’ refers to the bound of Eq. (9) calculated using the qubit form of the Hamiltonian. ‘Naive Fermionic Variance Bound’ is calculated in a similar way, except using the sums of the fermionic coefficients, as in Eq. (12). As noted above, the terms in the Fermionic Hamiltonian which consist of number operators (class I in Table 4) or products of number operators (class II in Table 4) actually have a variance which is upper-bounded by $\frac{1}{4}$ rather than one. One can substitute this tighter bound in Eq. (10) to yield the expression

[TABLE]

where we have assumed that the Hamiltonian is normal-ordered to simplify the expression. We present calculation based on this improved bound in the row titled ‘Fermionic Variance Bound.’ However, it is clear that the bounds obtained directly from the qubit representation of the Hamiltonian are considerably tighter than either of the bounds obtained using the fermionic representation. This difference is explained by the cancellation effects that we have described above.

In addition to these bounds on the variance of the estimator for the Hamiltonian, we also consider two approximations to this variance that are not guaranteed to be upper bounds (rows 4 and 5 of Table 5). These approximations, which we refer to for brevity as FVA and QVA, are calculated using the methodology of Ref. 3 using the fermionic and qubit Hamiltonians respectively. In that work Wecker et al. reasoned that, in a typical quantum chemical calculation, the orbitals would have occupation numbers near [math] or $1$. Therefore, the number and number-number terms in the Hamiltonian (partition I and partition III in Table 4) would have a variance that is close to zero. This assumption is satisfied exactly for the Hartree-Fock state when the appropriate single-particle basis is used, and should be approximately true when Hartree-Fock is qualitatively correct. Based on this assumption, Ref. 3 neglected these terms and then approximated the variance of the remaining terms in the Hamiltonian using the type of bounds we have already discussed.

Rows 4 and 5 of Table 5 show that that there is still a substantial difference between the variances calculated under this approximation using the two different representations of the Hamiltonian, i.e., between FVA and QVA. This is primarily due to the reduction caused by the cancellations among the double-excitation terms (class V in Table 4). Interestingly, the numbers presented for the ‘Qubit Variance Approximation’ (QVA) are nearly identical to those for the actual variance expected when measuring the Hartree-Fock state (row 6). In fact, any differences between these values are found to arise purely from numerical precision issues in the data. One can examine the Pauli operators that arise from performing the Jordan-Wigner transformation on the Hamiltonian after deleting the diagonal terms and see that each of them has an expectation value of exactly zero on the Hartree-Fock state. Measurements of these terms therefore achieve the maximum possible variance of $1$, while measurements of the deleted diagonal terms would have a variance that is exactly [math]. Thus, the calculation of the actual variance (when measuring each term in the Hamiltonian separately) using Eq. (10) for the Hartree-Fock state then yields the same value as the calculation of the bound of Eq. (9) under the approximation proposed by Ref. 3.

Appendix B Applying the Fermionic RDM Constraints to the Qubit Hamiltonian

In the previous section we saw a substantial difference between the bounds calculated from the fermionic operators and those calculated from the qubit operators after applying the Jordan-Wigner transformation. In light of this, it is natural to ask how the reduced density matrix (RDM) approach of Ref. 23 might perform when formulated using the qubit representation of the Hamiltonian. In Ref, Rubin et al. (2018), Rubin et al. proposed that one could use known $n$-representability constraints on the expectation values of few-fermion operators, in order to construct estimators for the expectation value of the Hamiltonian that will have lower variance. They showed how one could take a collection of algebraic equalities from these fermionic n-representability constraints and use them to construct a new Hamiltonian $\tilde{H}$ from the original $H$. According to their approach, $\tilde{H}$ is constructed to have the same expectation value as $H$, but a lower maximum variance according to the bounds discussed above. They performed this minimization of the upper bound using standard linear programming techniques.

We are primarily focused here on the impact of these techniques for a real-world experiment. Therefore, we shall compare the impact of performing this minimization on the fermionic and qubit representations of the Hamiltonian, using the actual observed variance with respect to the ground state as the figure of merit, rather than employing the bounds or approximations discussed above. We take the same Pauli Word Grouping strategy described in the main text and apply it to the Hamiltonians $\tilde{H}_{fermionic}$ and $\tilde{H}_{qubit}$. We define $\tilde{H}_{fermionic}$ and $\tilde{H}_{qubit}$ as the Hamiltonians that result from performing the upper bound minimization procedure of Ref. 23 in the fermionic and qubit representations respectively.

In Figure 5 we plot the ratio between the variances of $\tilde{H}_{qubit}$ and $\tilde{H}_{fermionic}$ for the ground state of each of the systems considered in the main text of this work. We list these systems in Table 2 below, duplicating Table II of the main text for convenience. Despite the substantial differences in the variance bounds formulated in the two representations, the impact of applying the RDM constraints to the qubit Hamiltonian rather than the fermionic one is found to be marginal, at best. For the majority of the systems it appears that the qubit-based bounds perform slightly better, but there are also a number of cases where this pattern is reversed.

Appendix C Low Rank Decomposition

In the main text, we explained that our strategy for measurement is based on rewriting the standard quantum chemical Hamiltonian in the following form:

[TABLE]

where the values $g_{p}$ and $g_{pq}^{(\ell)}$ are scalars, $n_{p}=a^{\dagger}_{p}a_{p}$, and the $U_{\ell}$ are unitary operators which implement a single particle change of orbital basis. Here we shall explain how one obtains that factorization starting from a standard representation. We follow the presentation of Berry et al. with minor deviations and refer the reader to Ref. Berry2019-ux for more details.

First, it is necessary to obtain the Hamiltonian in the chemist’s standard form,

[TABLE]

This differs from the physicist’s convention of Eq. (11), where the operators in the two-electron component of the Hamiltonian have both creation operators to the left of both annihilation operators. We assume the use of purely real spatial orbitals, and therefore the tensor $V_{pqrs}$ inherits the eight-fold symmetry,

[TABLE]

from the definition of the two-electron integrals Szabo2012-ag.

Now we can perform the decomposition. We treat the tensor $V$ as a matrix indexed by the collective indices $pq$ and $rs$. We can eigendecompose this matrix to yield

[TABLE]

In the above equation, $w_{\ell}$ are the eigenvalues of $V$, $v^{(\ell)}$ are the eigenvalues. We proceed by using this equality to rewrite the two-electron component of the Hamiltonian,

[TABLE]

with $v^{(\ell)}_{pq}$ inheriting the symmetry between the $p$ and $q$ indices from $V$.

The final remaining step is to transform Eq. (21), as well as the one-electron component of the Hamiltonian by diagonalizing each of the one-body operators. It is straightforward to express each of the one-body operators as diagonal operators in a rotated single-particle basis. The appropriate change of basis matrices can be obtained from the eigenvalues of the coefficient tensors, $T$ and the $g^{(\ell)}$s in our case. We can therefore express the Hamiltonian in the form of Eq. (17), dropping the spin indices for simplicity. The coefficients $g_{p}$ come from rotating to a basis where $T_{pq}$ is diagonal. The coefficients $g^{(\ell)}_{pq}$ likewise come from rotating to a series of bases where the tensors $v^{\ell}$ are diagonal between their $p$ and $q$ indices. The operators $U_{\ell}$ are the inverse of the operators that diagonalize the one-body operators $\sum_{\sigma\in\{\uparrow,\downarrow\}}\sum_{pq}T_{pq}a^{\dagger}_{p,\sigma}a_{q,\sigma}$ and $\sum_{\sigma\in\{\uparrow,\downarrow\}}\sum_{pq}v^{(\ell)}_{pq}a^{\dagger}_{p,\sigma}a_{q,\sigma}$. Note that the $p$ and $q$ indices of Eq. (17) represent new dummy indices and that the $w_{\ell}$ terms have been absorbed into $g^{(\ell)}_{pq}$, together with the contributions from the squares of the diagonalized $v^{\ell}_{pq}$ terms.

Appendix D Description of Data

In addition to this appendix, we also include the raw data we have generated through our numerical calculations which does not already appear in tables throughout the manuscript. The data is provided as a table in the csv file format with each row corresponding to one of the particular systems listed below in Table 2 (identical to Table II in the main text) and each column corresponding to the variance of a different estimator (or an ancillary piece of data, such as the energy of the system). Energies are provided in units of $E_{h}$ and variances in units of $E_{h}^{2}$.