A Hybrid Classical/Quantum Approach for Large-Scale Studies of Quantum Systems with Density Matrix Embedding Theory

A Hybrid Classical/Quantum Approach for Large-Scale Studies of Quantum Systems with Density Matrix Embedding Theory

Nicholas C. Rubin Rigetti Quantum Computing, 775 Heinz Ave., Berkeley, CA, 94710 nick@rigetti.com

Determining ground state energies of quantum systems by hybrid classical/quantum methods has emerged as a promising candidate application for near-term quantum computational resources. Short of large-scale fault-tolerant quantum computers, small-scale devices can be leveraged with current computational techniques to identify important subspaces of relatively large Hamiltonians. Inspired by the work that described the merging of dynamical mean-field theory (DMFT) with a small-scale quantum computational resource as an impurity solver [Bauer et al., arXiv:1510.03859v2], we describe an alternative embedding scheme, density matrix embedding theory (DMET), that naturally fits with the output from the variational quantum eigensolver and other hybrid approaches. This approach is validated using a quantum abstract machine simulator [Smith et al., arXiv:1608.03355] that reproduces the ground state energy of the Hubbard model converged to the infinite limit.

[inlinelist,1]label=(), \pdfstringdefDisableCommands

I Introduction

Many current computational studies of strongly correlated chemical or material systems involve a high accuracy simulation of their electronic structure. A ubiquitous theme in classical computational methods is to identify a relevant subspace of the full Hamiltonian such that a high fidelity simulation in this subspace results in a wavefunction with a high degree of overlap with the true full system wavefunction. Examples of this methodology are active space methods, popular in quantum chemistry, and embedding or impurity problems, prevalent in studies of strongly correlated materials. The embedded problem commonly requires a full-configuration interaction-like solution (FCI, an exact diagonalization of the embedded Hamiltonian) to significantly improve the description of the electronic structure. Unfortunately, FCI scales exponentially in the number of orbitals, restricting its application to small problems.

The integration of active space or impurity methods with quantum computation offers a route to improving wavefunction approximations by increasing the size of the active space to include relevant low-energy states. Bauer et al. Bauer et al. (2015) recently proposed the integration of dynamical mean-field theory (DMFT) and the quantum phase-estimation algorithm which was demonstrated in a proof-of-concept study by Kreula et al. Kreula et al. (2016). In a similar vein, Bravyi and Gosset studied the complexity of impurity methods integrated with quantum computation Bravyi and Gosset (2016). An integration of complete active space self-consistent field (CASSCF) and the phase-estimation algorithm is an example of one such integration that has already been suggested Reiher et al. (2016).

The variational quantum eigensolver (VQE) is a hybrid classical/quantum algorithm to approximate the ground state eigenvalues and density matrices of a Hamiltonian Peruzzo et al. (2014); McClean et al. (2016a); McClean (2015); McClean et al. (2016b); Wang et al. (2015); Wecker et al. (2015). The VQE is a classical optimization loop invoking a quantum abstract machine (QAM) Smith et al. (2016) used for state preparation and operator measurement. Recent work investigating this approach yielded quantum circuits for ansatz state preparation that are short and demonstrated that the algorithm is potentially robust to standard noise models McClean et al. (2016b); Sawaya et al. (2016).

The VQE using a wavefunction ansatz that is exponentially expensive on a classical computer, such as unitary coupled-cluster McClean et al. (2016a) or approximate Hamiltonian evolution Wecker et al. (2015), offers a potentially advantageous alternative to other high fidelity solvers (FCI, density matrix renormalization group, etc.) for active space calculations Reiher et al. (2016) and embedding schemes Bauer et al. (2015); Dallaire-Demers and Wilhelm (2016). In general, the VQE can be integrated into any computational chemistry methodology that requires a high-fidelity solver with an interface requiring the one-particle and two-particle reduced density matrices (- and -RDM, respectively) Helgaker et al. (2014); Werner and Knowles (1985); Gidofalvi and Mazziotti (2007).

Density matrix embedding theory (DMET) is a powerful alternative to DMFT that replaces the problem of finding the impurity’s local Green’s function with determining the embedded problem’s -RDM Knizia and Chan (2012, 2013); Wouters et al. (2016); Bulik et al. (2014a, b); Bulik (2015). DMET maps the problem of finding the ground state of a large -particle system onto many interacting impurity problems. This is achieved by self-consistently matching the -RDM from a low-level wavefunction for the entire system with the -RDMs of the impurities calculated with a high-level technique. Though it is not necessarily clear that such an iterative scheme would result in accurate local RDMs, DMET has been demonstrated to be highly accurate for one- and two-dimensional lattice models Knizia and Chan (2012); Bulik et al. (2014a), extended chemical systems Zheng and Chan (2016); Bulik et al. (2014b); Knizia and Chan (2013), and for ground state and transition state molecules Wouters et al. (2016). Furthermore, spectral functions (on the real-frequency axis) are accessible from DMET without any bath discretization error allowing for the construction of arbitrary dynamic correlation functions Booth and Chan (2015). These studies demonstrate the amazing quality of DMET with modestly sized impurities.

A key component of DMET is the embedded problem solver. Ideally, this solver produces a high-accuracy -RDM for the embedded problem such that the fragment piece of the mean-field -RDM can be matched to the correlated -RDM for the fragment. In this work we propose the use of VQE for the embedded problem solver as an alternative to FCI or DMRG routines. This change potentially allows for {enumerate*}

the study of a larger embedding Hamiltonian, allowing access to longer range correlation functions, and

generally accelerates the DMET algorithm by quickly finding the solution to the embedded problem on a quantum computer. As a proof of principle, the ground state energy of the repulsive- Hubbard model is determined for small rings and the thermodynamic limit.

Ii Density Matrix Embedding

In this section we review the salient features of density matrix embedding and discuss how the method naturally incorporates the use of the variational quantum eigensolver as an embedded Hamiltonian solver. References Wouters et al. (2016); Bulik (2015) provide a more detailed review on density matrix embedding while Bulik et al. (2014a, b) present alternatives to the self-consistent procedure that appears in DMET.

ii.1 Preliminaries

For a large system , its wavefunction can be arbitrarily bi-partitioned into a fragment (or impurity) and bath (or environment). Some examples of a fragment piece include a single or multiple sites of a lattice model or localized atomic orbitals corresponding to a piece of a larger molecular basis set. The total wavefunction is naturally expressed in the tensor product of basis states of the fragment and the bath, , which has a linear dimension of where (resp. ) is the dimension of the Hilbert space of the fragment (resp.  bath). Utilizing the Schmidt decomposition, an eigenstate of can be written as a sum over the tensor products of the Schmidt basis states111This Schmidt basis is constructed by the singular value decomposition of the coefficient tensor. Wouters et al. (2016); Peschel (2012); Nielsen and Chuang (2010).


This reformulation demonstrates that there is a local fragment basis and environment basis of the same size that produces an equivalent representation of . Without loss of generality, we can assume the Hilbert space of the fragment is smaller than the environment. In the Schmidt basis, no matter how large the bath, a fragment can only entangle with bath states. The Schmidt states can be used to project the Hamiltonian into a combined impurity/bath basis that necessarily has the same ground state as the original Hamiltonian but is significantly smaller in size:


For a general large-scale quantum system the wavefunction can only be computed approximately. This constraint naturally leads to the fundamental approximation in DMET: the embedded Hamiltonian is approximated by constructing bath states from the Schmidt decomposition of the ground state of an approximate quadratic Hamiltonian for the total system. As a consequence, the embedded Hamiltonian now contains an interacting fragment embedded in a non-interacting bath. The exact embedding Hamiltonian is approximated by matching the fragment’s -RDM with the low-level mean-field density matrix of the system by varying the embedding potential that appears in the quadratic Hamiltonian Knizia and Chan (2012).

ii.2 The Low-Level Quadratic Hamiltonian

For an arbitrary chemical system, the Hamiltonian contains terms that are at most four creation/annihilation operators corresponding to one- and two-particle interactions:


In order to vary the form of the bath states used when constructing the embedded Hamiltonian, a potential is added to the full chemical Hamiltonian. The augmented total system Hamiltonian is then approximated using the Hartree–Fock method to find an optimal single-particle basis:


ii.3 Bath Oribtals from a Slater Determinant

An alternative way to construct the fragment and bath embedding states that is similar to the Schmidt decomposition in (II.1) is a contraction over the bath basis:


Orthogonalizing into a set of states gives:


Note that while each , there are no more than of them. With this form for the embedding basis, one simply needs to find the form of the bath states constructed from the mean-field wavefunction. This is achieved by projecting the occupied states with overlap on the fragment onto the environment orbitals and normalizing.

There are no formal restrictions on the low-level wavefunction representing the total system , though it is desirable that one can quickly construct an approximate eigenfunction of . Alternative approximate wavefunctions that have been used are the antisymmeterized geminal power wavefunctions Tsuchimochi et al. (2015) that introduce some correlation into the bath, products of Bogoliubov quasi-particle states Zheng and Chan (2016), and a Slater determinant Knizia and Chan (2012). The following is a derivation of the embedding basis from a Slater determinant wavefunction. This low-level wavefunction for the full system can be expressed as the product of fermionic operators in the occupied subspace acting on the true vacuum:


The single particle orbitals can be constructed from linear combinations of site orbitals, , in the lattice:


The matrix corresponds with the Hartree–Fock transform from local lattice spin-orbitals to the basis that minimizes the Hartree–Fock Hamiltonian of the system. A rotation to the fragment and bath basis can be constructed in a similar manner to the canonical orthogonalization of atomic orbitals that appears in the Hartree–Fock procedure applied to molecules Szabo and Ostlund (1989); Helgaker et al. (2014).

The overlaps of the occupied orbitals projected onto the fragment space are calculated using the fragment projection operator acting on the occupied set of orbitals indexed by and :


The eigenvectors of this overlap matrix correspond to an orbital rotation that forms the fragment and bath basis. Diagonalizing with a unitary gives


Naturally, the eigenvalues of the overlap matrix are between and . Zero eigenvalues correspond to occupied states that have zero overlap on the fragment. The bath states can be constructed by projecting the occupied orbitals with overlap on the fragment into the bath:


Finally, the rotation to the Schmidt basis is defined by the direct sum of the bath transformation and an identity matrix of dimension equal to the number of orbitals on the fragment, :


The core states with zero overlap on the fragment are eliminated from the embedded Hamiltonian by including their interaction at the mean-field level similar to generating an active space interacting with a frozen set of core states Hosteny et al. (1975).

ii.4 Embedded Hamiltonian

The embedded Hamiltonian is constructed by projecting into the fragment and bath states


which, for a Slater determinant wavefunction, corresponds to a single- and double-particle integral transform Knizia and Chan (2013):


In the above transformation is the interaction with the non-entangled non-overlapping core states, is a matrix from (14) whose columns are the transformation vectors to the embedding basis, and and are the one- and two-particle integral tensors defined in the lattice basis.

ii.5 Determining the Optimal Embedding Potential

The key step in the DMET procedure requires improving the bath states so they approximate the bath states of the true embedded Hamiltonian. Schmidt bath state tuning is achieved by varying the potential that was added to the system Hamiltonian in Sec. II.2. Variation in this potential is the connection between the embedded Hamiltonian and the quadratic Hamiltonian representing the entire system. The total DMET procedure can be summarized as follows:

  1. Find the ground state Slater determinant wavefunction for the entire system with the Hartree-Fock Hamiltonian augmented with the embedding potential.

  2. For all fragments, use to construct the embedding basis and project the system Hamiltonian into each embedded Hamiltonian. Solve for the -RDM and -RDM with a high-level method.

  3. Adjust the embedding potential in the mean-field Hamiltonian such that the -RDM of the impurity calculated with the high-level wavefunction matches the -RDM calculated from the quadratic Hamiltonian.

  4. Repeat until particle conservation and the embedding potential do not change between iterations.

Step 3 is carried out by minimizing the squared norm between the fragment piece of the mean-field -RDM and the fragment -RDM computed from the ground state wavefunction for the embedded Hamiltonian:


In Wouters et al. (2016) this cost function is known as the fragment-only density matrix fitting. Alternatives such as full matching of the -RDMs in the embedding basis or matching the electron densities on the fragment have also been studied Bulik et al. (2014a, b).

The exponentially scaling FCI has prompted studying alternative methods for finding the correlated ground state -RDM Wouters et al. (2016); Zheng et al. (2016). Though this may not lead to an optimal approximation to the embedded Hamiltonian fragment and bath states it does allow for the study of significantly larger fragments–important for assessing longer range correlation functions and more accurate local expectation values Zheng and Chan (2016). The ability to study larger fragments at higher fidelity has many implications for accurate simulations of materials and chemistry. For example, if one were studying a reaction center in a metal-organic framework a fragment would likely need to contain the metal site and ligand orbitals for an accurate description of the energy and properties. The VQE enters in step when solving for the -RDM and -RDM of the embedded Hamiltonian. DMET’s need for the -RDM of the embedded Hamiltonian makes VQE an ideal quantum algorithm for this use case.

Iii The Unitary Cluster Ansatz and Trotterization in VQE

The VQE is a functional minimization scheme that leverages fast construction of the wavefunction to find expectation values of operators O’Malley et al. (2016); McClean et al. (2016a). The energy of the system is minimized by varying over parameters for the wavefunction ansatz (21):


The expectation value is determined by summing the expectation of each term in the Hamiltonian:


For a general quantum chemical Hamiltonian, scales quadratically with respect to basis set size. In order to evaluate the Hamiltonian expectation value, the wavefunction is prepared many times with a subsequent measurement of the operator. In this work the unitary coupled cluster (UCC) state is used as the function for the energy functional though others could be considered such as the variational adiabtic ansatz.

UCC is a many-body expansion wavefunction ansatz parameterized by the cluster coefficients associated with each generator. We use the UCC-singles-doubles ansatz corresponding to anti-Hermitian generators performing single- and double-particle excitations and de-excitations:


In quantum computing, the antisymmetric fermionic creation/annihilation operators are represented by distinguishable qubits. We utilize the Jordan–Wigner transformation (JW) for mapping fermionic creation/annihilation operators to Pauli spin-operators that preserve the anti-commutation relations and parity of the second quantized operators (27Whitfield et al. (2011). It should be noted that the representation of second quantized operators in Pauli matrices by the JW transform scales linearly in the number of Pauli terms222There are alternatives such as the Bravyi–Kitaev transformation that have logarithmic scaling when representing second quantized operators Tranter et al. (2015); Seeley et al. (2012). as the basis index is increased. The maps are defined as




For each shot of the VQE algorithm the wavefunction of the system—according to a UCC ansatz—is constructed by first preparing the Hartree–Fock reference, in the number of gates required to prepare the state where is the number of orbitals, then evolving the initial wavefunction according to (23), using in the number of gates for the exponentiation of the unitary generators.

The - and -RDM are measured after an optimal set of parameters are determined by measuring expectation values corresponding to the -RDM:


The -RDM is obtained by contraction from the -RDM:


The fragment -RDM of the total system is matched to the fragment piece of the -RDM obtained from the VQE algorithm by varying the embedding potential in (3).

iii.1 Commutation Relations for Generators in Unitary Coupled-Cluster

Unlike the classical cluster operators that commute at all orders Crawford and Schaefer (2000) the unitary coupled-cluster operators (25)–(26), henceforth denoted as , do not commute with each other at all orders. Therefore approximating the exponential by a particular order of the Suzuki–Trotter decomposition Suzuki (1994); Trotter (1959) may be important for constructing accurate distributions. The anti-Hermitian generators do not commute because the de-excitation operators do not commute with the excitation operators. This can be seen for the first order operators. Here we use to denote the particle space and to denote the hole space.


The first and last commutators vanish while the middle two terms survive when there is an overlapping index in and :


As a result, Trotter error must be considered when generating circuits for UCC state preparation. In Sec. IV we discuss the performance of UCCSD with Trotter orders and with Trotter steps beyond .333The term Trotter order refers to the order of the series approximation to the exponential of two non-commuting operators. The term Trotter steps refers to the number of slices that each order is broken into in order to minimize the the correction term.

Iv Results and Discussion

Table 1: Energy-per-site of the half-filled Hubbard model for a four-site ring with anti-periodic boundary conditions evaluated with DMET using one- and two-site fragments. The exact solution is computed with an exact diagonalization of the Hamiltonian. Both Trotter order and Trotter steps were equal to one for all UCCSD calculations

We considered the one-dimensional Hubbard model (35) with repulsive- interactions and anti-periodic boundary conditions as a representative test system for DMET.


The DMET self-consistency loops were run using the QC-DMET Wouters (2015) code with an interface to PySCF Su (2016) for exact diagonalization of the embedded Hamiltonian. The VQE algorithm for solving the embedded Hamiltonian used Quil and Rigetti Computing’s pyQuil package to construct the UCCSD ansatz and evolved the wavefunction on a noiseless quantum virtual machine Smith et al. (2016). We describe the programming environment in more detail in Appendix VII.1. All UCCSD-VQE runs started with a second-order Møller–Plesset perturbation theory guess for the cluster amplitudes and required about 20 iterations of BFGS Nocedal and Wright (2006) with the gradient numerically approximated. For single-site DMET calculations the low-level and high-level fragment density matrices were equivalent and thus the DMET loop only involved setting the chemical potential such that the number of electrons in each fragment summed to the total system’s electron count. For two-site models the fragment -RDM was generally not equal to the fragment Slater determinant -RDM and thus required a numerical search for the optimal embedding potential.

The non-commuting nature of the cluster operators do not drastically effect the accuracy of the UCCSD ansatz. Figure 1 plots the convergence of the energy-per-site as the Trotter order and Trotter steps are increased for . Trotter order and Trotter steps indicate the structure of the wave function ansatz for UCCSD. For example, the first order Trotter with -Trotter steps produces an ansatz of form:


where and correspond to the single and double anit-Hermitian generators with their coefficients which are described in Eq. (25) and Eq. (26). The persistence of the energy gap as Trotter order and Trotter steps are increased demonstrates that the single reference nature and the truncated set of generators are the main sources of error.

Figure 1: Energy-per-site of the four-site Hubbard model with with anti-periodic boundary conditions calculated with Trotter orders 1 and 2 with Trotter steps 1, 2, 3, and 4. The error between Trotter steps of Suzuki–Trotter orders is generally hundredths of a milli-Hartree indicating that the significant gap between FCI and UCCSD is based on the wavefunction ansatz.

The DMET methodology introduces three sources of error: {enumerate*}

the fragment size and thus the number of non-interacting bath states that can entangle with the fragment,

the embedded Hamiltonian solver, and

the pieces of the embedded Hamiltonian -RDM that are chosen to match for updating the embedding potential. The effects of (i) and (ii) are determined by studying small lattices and comparing the ground state energy against the FCI solution and the full UCCSD solution. Table 1 compares the ground state energy-per-site of a four-site ring at various interaction strengths between the FCI solution, the UCCSD solution, and various fragment sizes for DMET. The DMET solutions are labeled by DMET(), where is the number of sites considered in the fragment, followed by an acronym for the type of embedded solver–i.e. ED corresponds to FCI. The accuracy of the UCCSD solution to the four-site-ring lattice generally decreases as is increased demonstrating the expected result that singles and doubles generators cannot parametrize the full unitary group. A four-site lattice requires eight-spin-orbitals and thus eight qubits. The UCCSD ansatz contains 14 cluster amplitudes when restricting the wavefunction to the singlet subspace Bulik et al. (2015); Scuseria et al. (1988). Under a greedy parallelization of commuting instructions for breaking the unitary evolution into time steps, a single state evolution at first order Trotter involves time steps where multiple qubit operations are performed at each step. Each time step has an average of 0.8692 one-qubit gates and 0.8101 two-qubit gates. The coupled-cluster reference is usually chosen in a basis where the one-particle piece of the Hamiltonian is diagonal. For the Hubbard model this involves a Fourier transform on the one-particle and two-particle integral tensors resulting in one-body terms and two-body terms where is the number of lattice sites. The two-body terms are proportional to .

A single fragment site results in a trivial two-orbital embedded Hamiltonian. For this system, UCCSD + VQE is equivalent to FCI, CISD, and CCSD solvers. Considering a two-site DMET at and produced energies-per-site of and , respectively. Though not exactly equal to the DMET(2)-ED result, the error is on the order of the difference between the UCCSD solution for the full lattice and the FCI solution for the full lattice.

Figure 2: top: Energy-per-site in the 100-site Hubbard model with a one-, two-, and four-impurity sites in the DMET scheme along with the Bethe solution described in Shiba (1972). FCI is the exact diagonalization solution to the single impurity site problem. bottom: The energy difference between the Bethe ansatz solution and the DMET solutions. As expected, the error decreases significantly as the number of sites increases. The number of qubits required for each embedded Hamiltonian eigenvalue problem is four times the number of sites included in the fragment.

In Figure 2, we plot the energy-per-site for a -site lattice calculated with DMET and the Bethe ansatz Lieb and Wu (2003); Shiba (1972). The top subfigure demonstrates that DMET(1)-UCCSD method converges to the correct solution at high-. The bottom subfigure describes the error for DMET(1)-ED, DMET(1)-UCCSD, DMET(2)-ED, and DMET(4)-ED relative to the Bethe ansatz. As expected the quality of the solution improves as fragment size is increased. An site fragment requires qubits and thus we are restricted to small fragment sizes when performing a classical simulation.

V Conclusion

The hybrid classical/quantum computation model naturally fits with current state-of-the-art methods for studying the electronic structure of molecules and materials. Many of the current methods involve selecting an active space or embedded piece of the Hamiltonian to treat at higher fidelity. These techniques are especially important when describing the electronic structure of correlated materials and molecules such as metal-oxides, heterogeneous catalysis, multi-electron redox reactions, and catalysts. In contrast to exact diagonalization solvers commonly used in impurity or active space problems, quantum computation offers a technique to treat these Hamiltonians with high accuracy in polynomial time.

The first steps towards integrating a quantum computational resource with embedding schemes was outlined in Bauer et al. (2015). That work integrated quantum phase-estimation with DMFT as the impurity solver. Their results indicate that a modestly sized quantum computational resource of about 100 qubits can be used to simulate a relevant physical problem. Though their work describes DMFT as the impurity methodology they do discuss the use of alternative embedding schemes along with alternative algorithms to use on the quantum computation side.

In this work we have elaborated on one alternative in particular: DMET integrated with the VQE. DMET works by mapping the problem of finding the -particle ground state of a large system to many smaller interacting impurity problems. Therefore, DMET is applicable to very large scale problems be it either molecules or materials. Unlike DMFT that requires the frequency dependent two-particle Green’s function, DMET requires the -RDM of the embedded problem with high accuracy. This requirement makes DMET and the VQE a natural pair. The -RDM and -RDM are computed at each iteration of the VQE optimization.

The VQE algorithm is a general functional minimization routine that allows for the use of any quantum state preparation method. In this work we utilized the UCC ansatz and demonstrated that despite the non-commuting nature of the generators, the error for the ansatz comes from the truncation of the series and not Trotter error. Therefore, at least in the case of simulating local interaction models, low-order Trotter approximation and Trotter number are sufficient for accurate state preparation.

Current classical techniques for describing quantum chemistry or correlated materials are naturally accelerated with a hardware quantum processing unit (QPU). Many of the classical methods that require a -RDM or -RDM of an active space can be integrated with the VQE algorithm or other hybrid classical/quantum techniques. Further integration can lead to methodologies to study large quantum systems with higher resolution than what is currently possible with state of the art classical methods.

Vi Acknowledgments

The author acknowledges useful discussions of VQE and DMET with Jarrod McClean and members of the Rigetti Computing software team: Will Zeng and Robert Smith. A special thanks to Robert Smith for writing the Appendix on quantum programming and detailed editing of the draft.

Vii Appendix

vii.1 Details of the Quantum Programming Environment

In this work, execution of quantum programs was done using the quantum abstract machine and Quil, its quantum instruction language Smith et al. (2016). In this section, we outline elements of this quantum programming environment, including examples of essential techniques used to compute the - and -RDMs.

The VQE is, at its core, a classical optimizer around a quantum execution unit. In particular, a classical optimization loop produces the next parameters of (21) for minimization. Mathematically, parameterizes a wavefunction ansatz, but computationally, they are regarded as parameters to generate a new quantum program. These programs are run on a quantum execution unit modeled as a restricted quantum abstract machine444The machine is one with quantum state , classical state , static gates , parametric gates , program , and program counter . with


where the gates are understood to be the collection of operators acting on each qubit combination of the quantum state . The classical memory is as large as —the number of qubits—to hold the -basis measurements for each loop of the VQE, which are collected to approximate . Alternatively, by using a quantum virtual machine, the amplitudes of can be observed in the computational basis and can be computed directly.

The object of interest to simulate is usually specified in terms of sums of fermionic operations, such as the Hamiltonian in (16). After a fermionic transform such as (27), said object is represented as a sum of products of Pauli spin operators. We represent these Pauli sums algebraically as first-class objects in a canonical form. We wish to compute a program which acts identically to the exponentiation of this Pauli representation—as in (23)—on the quantum abstract machine’s . Exponentiating a Pauli sum requires two steps: factorization via Trotterization, and exponentiation of Pauli terms. Consider the Pauli terms


These are represented as objects in Python using the pyQuil library:

> A = 2.0*PauliTerm(’X’, 0)*PauliTerm(’X’, 1)
> B = -0.5*PauliTerm(’X’, 0)*PauliTerm(’Z’, 2)
> print "A =", A, "\nB =", B
A = 2.0*X0*X1
B = -0.5*X0*Z2

Exponentiation of a Pauli term has three parts: a change to the -basis, an entanglement of the qubits on which the operator non-trivially acts, and a rotation about . For example, we can see this by looking at a Quil program which allows one to compute the action of .

> print exponentiate(A)
H 0
H 1
CNOT 0 1
RZ(4.0) 1
CNOT 0 1
H 0
H 1

This Quil program itself is also a first-class object. Indeed, the result of a Suzuki–Trotter approximation of, say, must combine programs which compute and . We can see this by looking at the Quil program which computes in a first-order decomposition:

> print trotterize(A, B)
H 0
H 1
CNOT 0 1
RZ(4.0) 1
CNOT 0 1
H 0
H 1
H 0
CNOT 0 2
RZ(-1.0) 2
CNOT 0 2
H 0

Quil programs are written as a serial list of instructions, but in fact the semantics are not changed if the order between successive commuting instructions is changed. As such, we can think of such instructions occurring in a single time slice555If the th commuting instruction executes in time, then the time slice itself executes in time and is the worst case for arbitrary with maximal parallelization. If all are known for each time slice, then it is possible to do aggressive parallelization by overlapping commuting instructions from subsequent time slices, despite the time slices as a whole (i.e., its action on the total Hilbert space) being non-commuting.. A straight-line Quil program itself is parallelized if it is transformed according to a maximal semantics-preserving parallelization. Consider the following 16-instruction sample from a state evolution according to the UCCSD ansatz with first-order Trotter and one time slice for the four-site Hubbard model:

X 2
X 3
H 3
RX(1.5707963267948966) 5
CNOT 3 4
CNOT 4 5
RZ(0.00001) 5
CNOT 4 5
CNOT 3 4
H 3
RX(-1.5707963267948966) 5
RX(1.5707963267948966) 1
H 5
CNOT 1 2
CNOT 2 3
CNOT 3 4

Instruction sequences like X 2 and X 3 commute and may execute in the same time slice. However, X 3 and H 3 do not commute and must execute in different time slices. The parallelized Quil program is as follows:

Time Slice #1:  X 2
                X 3
                RX(1.5707963267948966) 5
                RX(1.5707963267948966) 1
Time Slice #2:  H 3
                CNOT 1 2
Time Slice #3:  CNOT 3 4
Time Slice #4:  CNOT 4 5
Time Slice #5:  RZ(0.00001) 5
Time Slice #6:  CNOT 4 5
Time Slice #7:  CNOT 3 4
                RX(-1.5707963267948966) 5
Time Slice #8:  H 3
                H 5
Time Slice #9:  CNOT 2 3
Time Slice #10: CNOT 3 4

The 16-instruction straight-line program has been reduced to a 10-time-slice program averaging one-qubit gates and two-qubit gates per time slice. All instructions being equal, this gives a 38% improvement in timing.

When a final Quil program is prepared through the various means described, it is dispatched to a representation of a QAM: either a hardware QPU or a software quantum virtual machine. In this work, a remotely deployed quantum virtual machine was used for all Quil execution.


  • Bauer et al. (2015) B. Bauer, D. Wecker, A. J. Millis, M. B. Hastings,  and M. Troyer, arXiv preprint arXiv:1510.03859  (2015).
  • Kreula et al. (2016) J. M. Kreula, L. García-Álvarez, L. Lamata, S. R. Clark, E. Solano,  and D. Jaksch, EPJ Quantum Technology 3, 11 (2016).
  • Bravyi and Gosset (2016) S. Bravyi and D. Gosset, arXiv preprint arXiv:1609.00735  (2016).
  • Reiher et al. (2016) M. Reiher, N. Wiebe, K. M. Svore, D. Wecker,  and M. Troyer, arXiv preprint arXiv:1605.03590  (2016).
  • Peruzzo et al. (2014) A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik,  and J. L. O’Brien, Nature communications 5 (2014).
  • McClean et al. (2016a) J. R. McClean, J. Romero, R. Babbush,  and A. Aspuru-Guzik, New Journal of Physics 18, 023023 (2016a).
  • McClean (2015) J. R. McClean, Algorithms Bridging Quantum Computation and Chemistry, Ph.D. thesis, Harvard University (2015).
  • McClean et al. (2016b) J. R. McClean, M. E. Schwartz, J. Carter,  and W. A. de Jong, arXiv preprint arXiv:1603.05681  (2016b).
  • Wang et al. (2015) Y. Wang, F. Dolde, J. Biamonte, R. Babbush, V. Bergholm, S. Yang, I. Jakobi, P. Neumann, A. Aspuru-Guzik, J. D. Whitfield, et al., ACS nano 9, 7769 (2015).
  • Wecker et al. (2015) D. Wecker, M. B. Hastings,  and M. Troyer, Phys. Rev. A 92, 042303 (2015).
  • Smith et al. (2016) R. S. Smith, M. J. Curtis,  and W. J. Zeng, arXiv preprint arXiv:1608.03355  (2016).
  • Sawaya et al. (2016) N. P. D. Sawaya, M. Smelyanskiy, J. R. McClean,  and A. Aspuru-Guzik, Journal of Chemical Theory and Computation 12, 3097 (2016), pMID: 27254482, http://dx.doi.org/10.1021/acs.jctc.6b00220 .
  • Dallaire-Demers and Wilhelm (2016) P.-L. Dallaire-Demers and F. K. Wilhelm, Phys. Rev. A 93, 032303 (2016).
  • Helgaker et al. (2014) T. Helgaker, P. Jorgensen,  and J. Olsen, Molecular electronic-structure theory (John Wiley & Sons, 2014).
  • Werner and Knowles (1985) H. Werner and P. J. Knowles, The Journal of Chemical Physics 82, 5053 (1985).
  • Gidofalvi and Mazziotti (2007) G. Gidofalvi and D. A. Mazziotti, The Journal of Chemical Physics 127, 244105 (2007), http://dx.doi.org/10.1063/1.2817602.
  • Knizia and Chan (2012) G. Knizia and G. K.-L. Chan, Phys. Rev. Lett. 109, 186404 (2012).
  • Knizia and Chan (2013) G. Knizia and G. K.-L. Chan, Journal of Chemical Theory and Computation 9, 1428 (2013), pMID: 26587604, http://dx.doi.org/10.1021/ct301044e .
  • Wouters et al. (2016) S. Wouters, C. A. Jiménez-Hoyos, Q. Sun,  and G. K.-L. Chan, Journal of Chemical Theory and Computation 12, 2706 (2016), pMID: 27159268, http://dx.doi.org/10.1021/acs.jctc.6b00316 .
  • Bulik et al. (2014a) I. W. Bulik, G. E. Scuseria,  and J. Dukelsky, Phys. Rev. B 89, 035140 (2014a).
  • Bulik et al. (2014b) I. W. Bulik, W. Chen,  and G. E. Scuseria, The Journal of Chemical Physics 141, 054113 (2014b), http://dx.doi.org/10.1063/1.4891861.
  • Bulik (2015) I. W. Bulik, Electron correlation in extended systems via quantum embedding, Ph.D. thesis, Rice University (2015).
  • Zheng and Chan (2016) B.-X. Zheng and G. K.-L. Chan, Phys. Rev. B 93, 035126 (2016).
  • Booth and Chan (2015) G. H. Booth and G. K.-L. Chan, Physical Review B 91, 155107 (2015).
  • Peschel (2012) I. Peschel, Brazilian Journal of Physics 42, 267 (2012).
  • Nielsen and Chuang (2010) M. A. Nielsen and I. L. Chuang, Quantum computation and quantum information (Cambridge university press, 2010).
  • Tsuchimochi et al. (2015) T. Tsuchimochi, M. Welborn,  and T. Van Voorhis, The Journal of Chemical Physics 143, 024107 (2015), http://dx.doi.org/10.1063/1.4926650.
  • Szabo and Ostlund (1989) A. Szabo and N. S. Ostlund, Modern quantum chemistry: introduction to advanced electronic structure theory (Courier Corporation, 1989).
  • Hosteny et al. (1975) R. P. Hosteny, T. H. Dunning, R. R. Gilman, A. Pipano,  and I. Shavitt, The Journal of Chemical Physics 62, 4764 (1975).
  • Zheng et al. (2016) B.-X. Zheng, J. S. Kretchmer, H. Shi, S. Zhang,  and G. K. Chan, arXiv preprint arXiv:1608.03316  (2016).
  • O’Malley et al. (2016) P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik,  and J. M. Martinis, Phys. Rev. X 6, 031007 (2016).
  • Whitfield et al. (2011) J. D. Whitfield, J. Biamonte,  and A. Aspuru-Guzik, Molecular Physics 109, 735 (2011).
  • Tranter et al. (2015) A. Tranter, S. Sofia, J. Seeley, M. Kaicher, J. McClean, R. Babbush, P. V. Coveney, F. Mintert, F. Wilhelm,  and P. J. Love, International Journal of Quantum Chemistry 115, 1431 (2015).
  • Seeley et al. (2012) J. T. Seeley, M. J. Richard,  and P. J. Love, The Journal of chemical physics 137, 224109 (2012).
  • Crawford and Schaefer (2000) T. D. Crawford and H. F. Schaefer, Reviews in computational chemistry 14, 33 (2000).
  • Suzuki (1994) M. Suzuki, Communications in mathematical physics 163, 491 (1994).
  • Trotter (1959) H. F. Trotter, Proceedings of the American Mathematical Society 10, 545 (1959).
  • Wouters (2015) S. Wouters, “QC-DMET,” https://github.com/SebWouters/QC-DMET (2015).
  • Su (2016) Q. Su, “PySCF,” https://github.com/sunqm/pyscf (2016).
  • Nocedal and Wright (2006) J. Nocedal and S. Wright, Numerical optimization (Springer Science & Business Media, 2006).
  • Bulik et al. (2015) I. W. Bulik, T. M. Henderson,  and G. E. Scuseria, Journal of Chemical Theory and Computation 11, 3171 (2015), pMID: 26575754.
  • Scuseria et al. (1988) G. E. Scuseria, C. L. Janssen,  and H. F. Schaefer, The Journal of Chemical Physics 89, 7382 (1988).
  • Shiba (1972) H. Shiba, Phys. Rev. B 6, 930 (1972).
  • Lieb and Wu (2003) E. H. Lieb and F. Wu, Physica A: Statistical Mechanics and its Applications 321, 1 (2003).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description