Statistical approach to quantum phase estimation

You may also like PAPER • OPEN ACCESS Statistical approach to quantum phase estimation - Electron collisions with the CN radical: bound states and resonances Stephen Harrison and Jonathan Tennyson To cite this article: Alexandria J Moore et al 2021 New J. Phys. 23 113027 - Multichannel generalization of eigenphase-preserving supersymmetric transformations Andrey M Pupasov-Maksimov View the article online for updates and enhancements. - R-matrix calculations of low-energy electron collisions with methane Will J Brigg, Jonathan Tennyson and Martin Plummer This content was downloaded from IP address 128.210.107.25 on 29/11/2021 at 13:16 Statistical approach to quantum phase estimation Abstract We introduce a new statistical and variational approach to the phase estimation algorithm (PEA). Unlike the traditional and iterative PEAs which return only an eigenphase estimate, the proposed method can determine any unknown eigenstate–eigenphase pair from a given unitary matrix utilizing a simplified version of the hardware intended for the iterative PEA (IPEA). This is achieved by treating the probabilistic output of an IPEA-like circuit as an eigenstate–eigenphase proximity metric, using this metric to estimate the proximity of the input state and input phase to the nearest eigenstate–eigenphase pair and approaching this pair via a variational process on the input state and phase. This method may search over the entire computational space, or can efficiently search for eigenphases (eigenstates) within some specified range (directions), allowing those with some prior knowledge of their system to search for particular solutions. We show the simulation results of the method with the Qiskit package on the IBM Q platform and on a local computer. 1. Introduction Efficient spectral decomposition of large matrices is a key component to many optimization and machine learning algorithms, with applications ranging from factoring and searching algorithms to computational chemistry [1]. On classical computers, spectral decomposition scales super-linearly with the system dimension [2], making it intractable for large problems. Due to the utility of spectral decomposition and its classical limitations, quantum approaches to spectral decomposition and eigenvalue estimation have been pursued [3]. One significant approach is the quantum phase estimation algorithm (PEA) [4]—a means of determining unknown eigenphases of a unitary matrix—which is a key subroutine in a number of quantum algorithms including Shor’s factoring algorithm [5], quantum principal component analysis [6], the generalized Grover’s search algorithm [7], and quantum simulations [8–10]. Near-term quantum systems operate in the noisy intermediate-scale quantum (NISQ) regime [11], facing restrictions on both circuit depth and breadth due to decoherence and gate infidelity. Consequently, interest in the traditional PEA [4] and quantum principal component analysis [6] has been channeled toward developments in the iterative PEA (IPEA) [12]—a method which estimates an unknown phase over multiple circuit iterations—allowing for significant reduction in both qubit usage (circuit breadth) and controlled-gate operations (circuit depth). The IPEA has been demonstrated on photonic systems [13]. On the other hand, variational quantum algorithms (VQA)—which use a classical computer to control and optimize the parameters applied in a quantum circuit—have been developed for a variety of problems as they leverage the speedup of quantum algorithm with lower-depth circuits [14, 15]. Here, we introduce a quantum–classical hybrid algorithm combining the PEA with the VQA—which we call the statistical PEA (SPEA)—and show preliminary simulation results on the IBM Q platform with the Qiskit package [16] as well as simulations on a local computer. The method is able to determine any Figure 1. The traditional PEA using n control dc -dimensional qudits in the control register. The quantum gates are colored in brown and the measurement gates in blue. The target register is highlighted in purple with dashed outline and control register in gray with dotted outline. Note that in the dc > 2 case, the H ‘Hadamard’ gate acts as a discrete quantum Fourier transform (QFT) gate. Additionally, in the dc > 2 case the control gates acts as MVCGs, applying the gate to the qth power when the control qudit is in state |q. The circuit will estimate the phase θ̃k of eigenstate |ν k  to precision dc−(n) . The estimate of eigenphase θk is returned with probability |Φ|νk |2 . unknown eigenstate–eigenphase pair from a unitary matrix by utilizing hardware intended for the IPEA. Further, the SPEA can be applied repeatedly to obtain a full spectral decomposition. The SPEA may be compared to other variational quantum eigensolvers [17–19], the primary difference being other variational eigensolvers work directly on a (Hermitian) matrix encoded as a quantum state using specially designed quantum circuits. The SPEA assumes access to a gate representation of the unitary exponentiation of the state—or assumes simultaneous availability of several copies of the quantum state to approximate the quantum gate à la [6]. In return, the SPEA requires a polynomially-reduced number of (classical) optimization parameters—as it optimizes for a single eigenstate, rather than diagonalize the entire matrix simultaneously—and directly delivers eigenstate–eigenphase pairs (whereas other approaches may allow on-demand generation of eigenstates, but require tomography if knowledge of the state is needed). The SPEA is also able to search for eigenphases within specified ranges, allowing those with some prior knowledge of their system to search for particular solutions, whether ground state (near minimum eigenphase), principle (near-maximal eigenphase), or any other region of interest. This paper is organized as follows: section 2 reviews the traditional and IPEA and introduces a statistical metric C for quantifying the proximity of any given input-state to its closest eigenstate. Section 3 describes the SPEA and discusses the connections between the C factor and the quality (in terms of proximity) of the derived eigenstate–eigenphase pairs (with the derivation details in appendix A). Section 3 also outlines the optimization process for obtaining the eigenstate–eigenphase pairs. Simulation results on different platforms are reported and discussed in section 4; methodology details are provided in appendices B and C. Scaling for larger systems is considered in section 5. We conclude with a discussion on the performance of the SPEA and propose future directions and applications of the method in section 6. 2. Phase estimation algorithms Traditional PEA implementations, diagrammed in figure 1, take any given unitary Û and any given eigenstate |ν of Û and return the corresponding eigenphase θ where Û|ν = ei2πθ |ν. (1) The (approximate) eigenphase θ̃ ∈ [0, 1) (equivalently, ∈ [−.5, .5)) may be directly measured on the control qubits (or qudits, when the control is dc -dimensional) of the PEA. The target register is typically unmeasured during the process. For an arbitrary target register input |Φ, the probability of the circuit representing a particular eigenstate |ν k  and the associated eigenphase θk is |νk |Φ|2 . If |Φ is not itself an eigenstate, the eigenphase retrieved varies each time the PEA circuit is run. The prototypical PEA thus approximates a particular θ in a single trial. The traditional PEA requires large quantum circuits which are often unreliable in the NISQ regime. To overcome hardware constraints, the IPEA was developed. The IPEA significantly reduces circuit depth requirements by approximating a particular θ one qubit (or dc -level qudit) at a time, starting from the least significant qubit (qudit). The IPEA requires a rotation gate—a linear phase across the control register—to ‘subtract’ off eigenphase information determined in previous iterations (i.e. if the quantum circuit’s state Figure 2. The IPEA. See figure 1 for color-conventions and notes on the MVCG. To retrieve n dits of the eigenphase (i.e. phase precision ±1/dcn ), run the circuit with x = n − 1 and θR = 0; the measured control state is the nth base-dc dit of the unknown phase. Proceed to run the circuit for x = n − 2 and set θR according to the previous control results; the measured control state is the (n − 1)th dit. Continue the process iteratively until x = 0 and the entire phase is recovered. Note that the iterative method is diagrammed for a single qudit control, but may be realized with any number of control qudits, similar to the traditional PEA.  before Rz (θR ) is q αq |q|Φq , then after the rotation gate the quantum circuit’s state is  −iq2πθR |q|Φq ). The IPEA, as the name suggests, requires a number of iterations equal to the number q αq e of bits (dits) of precision desired from the eigenphase. Additionally, the input to the target register of an IPEA must either be an eigenstate (and identically prepared each iteration) or the previous iteration’s output must propagate forward and serve as the next iteration’s input. An IPEA circuit is diagrammed in figure 2. In the general case, the control qudit may be high-dimensional (dc -level). In this case, the Hadamard gates represent a dc -dimensional quantum Fourier transform gate and the control-Û gate is a multi-level control gate [20]: when the control state is |q, a Û q gate is applied to the target register. Consider the IPEA in its ‘last’ iteration’s settings (x = 0 in figure 2). When the target register is an eigenstate |Φ = |ν and the rotation gate is used to subtract off phase 2πθR = 2πθ, the control dits deterministically collapse to state |0. When either the target input is not an eigenstate and/or θR is not the corresponding eigenphase, the control dits will collapse to |0 with non-unity probability. Indeed, for eigenstate input |ν with a dc -level control, the final state of the control qudit before measurement is   dc −1 d c −1 q 1 2πin θ−θR − d c |q, |ΨC  = e (2) dc q=0 n=0 where θ is the eigenphase of |ν and −θR (θR ∈ [0, 1)) is the rotation applied by the rotation gate. The probability of measuring the system in output bin |0 is d −1 2 c  1   e2πin(θ−θR )  Pθ (−θR ) = |0|ΨC | = 2   dc  n=0 2 = P0 (θ − θR ). (3) Pθ (−θR ) goes to one as θ R approaches θ, as shown in figure 3. In the most general case, where the target register is an arbitrary (non-eigenstate) input state |Φ and the rotation gate subtracts off phase 2πθR , the probability that the control qudits will collapse to |0 is C(|Φ, θR ) = d t −1 |νk |Φ|2 Pθk (−θR ) k =0 d −1 2 dt −1 c   1 2 2πin(θk −θR )  = 2 |νk |Φ|  e .   dc n =0 (4) k =0 Where Û is dt -by-dt -dimensional and the target register is dt -dimensional. Appreciate that C(|Φ, θR ) = 1 if and only if |Φ is an eigenstate and θR is its corresponding eigenphase. 3. Statistical approach to PEA The non-deterministic nature of the IPEA (in the non-eigenstate case) disqualifies the circuit from use as an eigenphase estimator in the standard approach. The SPEA instead considers the probabilistic outputs of the IPEA (and PEA) as valuable information which—when coupled with a classical controller as in figure 4—allows quantum PEA-like hardware to be used in a variational approach to determine any unknown eigenphase–eigenstate pair. The quantum hardware required is that of a traditional PEA with Figure 3. Plot of P0 (Δθ) (equation (3)) from Δθ ∈ [−.5, .5) for various dc . Note that P0 has infinite domain with period 1. Probability of the (dc -level) control register of an IPEA collapsing to |0 as a function of difference between the eigenphase θ and the applied rotation θR , Δθ ≡ θ − θR for an eigenstate input. Note that when the applied rotation matches the eigenphase (Δθ = 0), the control collapses to |0 deterministically. Denote the region around Δθ = 0 (from dot to dot) as the central lobe of P0 (Δθ), and the small lobes with local maxima outside of it as the sidelobes. See that the higher the system’s dimensionality, the narrower the probability curve’s central lobe and the lower local maxima in the sidelobes. Note that dc = 2 has no sidelobes (the probability is monotonic on either side of the central lobe). Also note P(Δθ) = 0 for Δθ = dc−1 and the width of the central lobe is therefore ΔθFWFM = 2dc−1 . Figure 4. Diagram of variational classical–quantum system. Classical processes are indicated by double blue lines and quantum processes by single black lines. Quantum gates are shown in brown and measurement gates in blue. The (potentially high-dimensional) PEA-type circuit is simplified from the typical IPEA in that U need not be raised to high orders (dcx ) corresponding to desired eigenphase precision and the measurement gate need only distinguish between the |0-state and the not (|0)-state. The (negative) cost function (estimate) C˜ is returned after a predetermined number of trials of the quantum circuit, approximating a probability. The classical controller optimizes |Φ and θR based on the (negative) cost function returned by the quantum circuit. Process continues until C̃(|Φ, θR ) = 1 or improvement in C̃(|Φ, θR ) ceases. single-dit precision (n = 1) and the rotation gate standard to the IPEA (i.e. an IPEA with x set to 0). The classical controller determines |Φ and θR which are used in the PEA-type circuit. Multiple trials of the quantum circuit are run to approximate the probability C˜ ≈ C(|Φ, θR ) (of equation (4)). Note that the PEA-like circuit need only detect two measurement outcomes: |0 and not (|0), further reducing hardware ˜ requirements compared to typical high-dimensional PEAs. Treating the estimate −1 · C(|Φ, θR ) as a cost ˜ function in an optimization process (making C the negative cost function), the classical controller adjusts |Φ and θR , until the quantum circuit near-deterministically returns |0 as the output state. When ∗ ∗ ˜ , θR ) ≈ 1, the classical controller has found the (approximate) eigenstate |Φ∗ and the associated C(|Φ eigenphase θR∗ . The quality of the eigenstate |Φ∗ and eigenphase θR∗ retrieval can be quantified by C ∗ = C(|Φ∗ , θR∗ ). C ∗ can both (1) determine the maximum distance from the eigenphase θR to the nearest eigenphase θ k and (2) find the fidelity of |Φ∗ to actual eigenstate(s). Derivations of both are provided in appendix A. The (negative) cost function C acts as a metric for quality of eigenvalue–eigenstate retrieval as shown in appendix A; by finding |Φ and θR which maximize this metric, we arrive at good estimates for an eigenstate (|Φ) and eigenphase (θR ) pair. Following is the classical algorithm used to maximize C, which is similar to a gradient search algorithm: (a) The classical controller chooses a |Φ at random (b) The classical controller constructs an orthogonal basis {|Bm } including |Φ (c) 1 (Standard method: viable when the quantum circuit can measure output bins |0 and not (|0).) ˜ The quantum circuit evaluates C(|Φ, θR ) over a range of θR and returns the maximum value C ∗ 2 (Alternative method: viable when the quantum circuit can measure all dc output bins: |0 through |dc − 1.) ˜ The quantum circuit evaluates C(|Φ, 0) and uses this result to approximate the eigenphase θ∗ . ˜ The quantum circuit then evaluates C(|Φ, θR = θ∗ ) and returns C ∗ . (d) For all m = 0 to 2dt − 1, we set a = 1 and the following occurs: √ 1. If m  dt then z = −1. Otherwise z = 1. 2. The classical controller generates the new state: |A |Φ  =  , A|A (5) where |A = |Φ + z · a · (1 − C ∗ )|Bm mod dt . 3. |Φ  is fed to the quantum circuit, the maximum value returned is C ∗ 4. If C ∗ > C ∗ , then |Φ = |Φ  and C ∗ = C ∗ . Otherwise |Φ and C ∗ are unchanged. (e) If |Φ was not updated during step (d), set a = a/2 and repeat step (d). (f) If C ∗ is greater than the stopping condition or the maximum run-time has been exceeded, the classical controller concludes and returns |Φ, C ∗ , and θR∗ . Otherwise the process continues from step (b). A few observations on the optimization process may be made. For each iteration, at least 2dt distinct input states are used. For each of these input states a set of {θR } is applied (when using the ‘standard approach’ in step (c)). Initially, the {θR } range from 0 to 1 with coarse resolution; as the optimization proceeds, {θR } will become fine and include phases from a limited region. Notably, we can choose to run the optimization process limiting {θR } to a narrow range of space from the outset. In this way, we may choose to look only for ground state (small θ), principle (large θ), or any other particular solutions to equation (1). In addition, we may eliminate known eigenstates or directions not of interest by excluding them from {|Bm } (step (b)) each iteration. In this fashion, the SPEA may be used to determine a complete (or partial) spectral decomposition of Û. Finally, we note while the hardware conventional to a PEA is utilized, this system is superior to the original PEA, as it determines both the eigenstate and the eigenphase, given no prior knowledge. 4. SPEA simulation We test the proposed algorithm on the IBM Q platform and on a local computer. In both cases, a classical computer is used to simulate the C parameter (of equation (4)) delivered by a quantum circuit. These simulations of a quantum system are ideal: neither the IBM Q nor the local computer simulations include any noise terms i.e. all quantum gates are assumed to operate with perfect fidelity. The IBM Q trials study the convergence of the optimization algorithm to any single eigenstate on two- and four-dimensional systems. The local computer simulations run a full spectral decomposition on a 16-dimensional system with various control levels dc . The local computer simulations also consider the convergence of the optimization algorithm to any single eigenstate for various control levels dc on systems ranging from 8 to 512-dimensional, to provide some indication of the method’s scalability. Both simulations apply the variational algorithm as defined in section 3, with one primary difference: the local computer simulations follows the primary method of step (c) whereas the IBM Q simulations follow the alternative method. The IBM Q simulation runs one measurement with Rz (θR = 0) and applies the eigenphase estimation methodology introduced in the Discussion of [20]—under the (inaccurate first, but increasingly accurate) assumption that the input state is an eigenstate—to obtain a phase estimate θ∗ . Then, the measurement is run with Rz (θR = −θ∗ ) to obtain the metric C used for the optimization. By contrast, the local computer’s simulations follow the primary method, picking a representative sample of input phases to apply to the Rz gate and selecting the largest C that arises. The local computer’s simulations therefore require more runs of the quantum circuit per trial, but only require two control-qudit detectors: one for the |0 state and one for the not(|0), whereas the IBM Q methodology needs one detector for each Table 1. IBM Q Qiskit QASM simulator results. Three unitary operators are simulated on the IBM Q platform. The SPEA is run 20 times starting from each input state. The distance from the input state to the nearest eigenstate is quantified by the absolute inner product of the two vectors (inner product 1 being identical and smaller values indicting greater difference). The table records the mean and standard deviation (S.D.) of the number of iterations needed to reach the stopping condition (1 − C˜ = 10−4 ). The average absolute phase error is reported in radians. For each Ui the input states range from good approximations of one of the operator’s eigenstates to input states which are nearly equal-distance from every eigenstate. See that the iteration mean tends to increase with decreasing inner product but the ultimate phase error is generally agnostic to the input state difference. Iteration Unitary U1 U2 U3 Input state Abs. inner product Mean S.D. Mean phase error (0.1951, 0.9808) (0.3827,0.9239) (0.7071, 0.7071) (0, 0, 0.7432, 0.6690) (0, 0, 0.6690, 0.7432) (0, 0, 1, 0) (1, 0, 0, 0) (0.7071, 0, 0.7071, 0) (−0.1379, 0, 0, 0.9904) (0, 0.7807, 0.6247, 0) (0, 1, 0, 0) (0.7071, 0, 0, 0.7071) (0.5774, 0.5774, 0, 0.5774) 0.98 0.92 0.71 0.99 0.99 0.71 0.71 0.50 0.99 0.99 0.71 0.62 0.51 6.20 8.15 8.90 5.85 6.7 17.7 23.05 21.3 1.15 1.1 4.35 4.15 21.5 2.82 3.41 3.34 8.14 10.42 6.06 11.22 10.71 0.36 0.3 4.17 1.01 11.06 1.099 × 10−2 1.005 × 10−2 1.005 × 10−2 2.083 × 10−2 2.168 × 10−2 1.663 × 10−2 2.167 × 10−2 2.262 × 10−2 1.885 × 10−2 1.508 × 10−2 1.414 × 10−2 1.570 × 10−2 2.199 × 10−2 control level (|0, |1, . . . , |dc − 1). The alternative approach (or some hybrid approach) is generally preferable if the hardware is available for dc detectors. 4.1. Qiskit simulation On the IBM Q experience platform, we developed our quantum algorithms with Qiskit, the python-based programming package provided by IBM Q which offers all the facilities to design, simulate and execute quantum algorithms on IBM’s quantum computers [16]. In this section we present the simulation results of the SPEA on the Qiskit quantum simulator. Three sets of simulations are run on the IBM Q, one with two-dimensional operator U1 and the other two with four-dimensional operators U2 and U3 , the matrix forms of which are shown in appendix B. U1 and U2 are operators directly built with the default gates offered by the IBM Q and U3 is a unitary exponentiation based on the Hamiltonian of the hydrogen molecule generated with Bravyi–Kitaev transformation [21]. The second quantization Hamiltonian of a hydrogen molecule with a bond length 0.74 Å is calculated by the STO − 3G minimal basis using PySCF [22] and the transformation is done by OpenFermion [23]. We encode the matrix into the ‘operator’ class provided by Qiskit [16]. In the simulations of each unitary operation Ui where i = 1, 2, 3, we start with the input states that are good approximations of one of the operator’s eigenstates and then move to input states which are nearly equal-distance from every eigenstate. We quantify the distance of the input state |Φ to its nearest eigenstate |ν by taking the absolute inner product |Φ|ν| as reported in table 1. In each simulation we run the same input state 20 times and set the maximum iteration number to be 50 (to save the resources) and the stopping condition, which is the difference between the C factor and 1, to be 10−4 . The stopping condition is set so that when it is met we will have a reasonably good approximation of the eigenstate. We then calculate the average number of iterations and standard deviation of the number of iterations required to exceed the stopping condition. Most trials reach the stopping condition before exceeding the iteration limit and give a good approximation of one of the eigenstate–eigenphase pair, as indicated by the low mean phase error reported. The results are shown in table 1. For each operator U1 , U2 , U3 , input states which are initially close to an eigenstate (input states with a large absolute inner product) have lower required iteration number than those which are initially far from all eigenstates (low absolute inner product). Appreciate that the eigenstate converged to is non-deterministic, as the optimizer itself is non-deterministic due to randomness added by the random orthogonal basis in step (b). In other words, added randomness may converge the input state to an eigenstate other than the closest eigenstate. The mean phase error recorded in table 1 is calculated by taking the absolute value of the difference between the eigenphase of the converged input state and the true eigenphase of the eigenstate that the input state converged to. As the input state can converge to different eigenstates in the simulation, we report the absolute phase error rather than the error percentage. No correlation between the phase error and the absolute inner product is apparent, indicating the quality of the final eigenphase–eigenstate pair is agnostic to the proximity of the initial input state to any eigenstate. Table 2. Statistics for 120 successful trials of complete spectral decomposition of UH2 O matrix. Cgoal is the C̃ value the optimizer attempts to reach, and generally does reach for at least the first 10 (of 16) eigen-estimates. Creq is the C̃ value the optimizer is required to reach for all eigen-estimates, else the trial is abandoned. For successful trials, the mean Fidelity and standard deviation are reported, as well as the mean eigenphase error (in radians). Note that dc = 2 was run on two different Cgoal , Creq levels for comparison. Fidelity dc Cgoal Creq Trials Fails 2 2 3 4 5 6 7 8 0.999 0.995 0.995 0.995 0.995 0.995 0.995 0.995 0.95 0.9 0.9 0.9 0.9 0.9 0.9 0.9 120 120 120 120 120 120 120 120 77 13 17 32 30 26 60 132 Phase error Mean S.D. Mean 0.984 0.966 0.981 0.986 0.986 0.989 0.991 0.992 7.29 × 10−3 10.68 × 10−3 6.62 × 10−3 5.32 × 10−3 7.10 × 10−3 4.80 × 10−3 7.05 × 10−3 6.13 × 10−3 2.84 × 10−2 4.34 × 10−2 3.12 × 10−2 2.40 × 10−2 1.86 × 10−2 1.53 × 10−2 1.37 × 10−2 1.20 × 10−2 S.D. 10.7 × 10−2 17.3 × 10−2 9.14 × 10−2 7.36 × 10−2 4.98 × 10−2 3.80 × 10−2 3.65 × 10−2 3.54 × 10−2 Variations in mean phase error are likely a function of which particular eigenphase–eigenstate pair was converged to. During the simulation of U3 with an input state of equal weight combination of all the eigenstates—i.e. the hardest input state to converge to an eigenstate—there are few cases that the iteration limit is reached and the simulation did not reach the stopping condition. This can usually be fixed by increasing the iteration limit. In summary, these results indicate that the SPEA method is capable of delivering high-quality estimates of eigenphase–eigenstate pairs with no prior knowledge of the operator’s eigenstates, in the case of both arbitrary (U1 , U2 ) and physically relevant (U3 ) operators. The quality of the estimates is not influenced by prior system knowledge; however, the resources required to deliver an eigenstate–eigenphase pair may be reduced with prior knowledge. 4.2. Full spectral decomposition The statistical approach differs from some other variational approaches [17] in that it does not diagonalize the input state matrix, but solves for only one eigenphase–eigenstate pair. This allows for significant reduction in the number of parameters (and iterations) needed to perform the optimization. However, as shown below, a complete spectral decomposition is realizable. As a representative case, we consider the 16-by-16 Hamiltonian HH2 O for the water molecule H2 O with the H–O–H angle at 104.5◦ and the bond length at 1.0 a.u. given in appendix C [24]. The Hamiltonian is converted to a unitary exponentiation, UH2 O = eiHH2 O , and the matrix’s spectral decomposition is simulated with the SPEA on a local computer. The simulation is run for various control levels dc until 120 successful spectral decompositions are achieved. Each of the 16 eigenphases are retrieved in a random order. The optimization runs until C ∗  Cgoal . If the optimizer is unable to reach Cgoal , the process for that eigenphase will conclude so long as C ∗  Creq . If Creq is not met, the entire spectral decomposition is abandoned and the trial is classified as failed. Generally, Cgoal is achieved for the first 10 eigenvalues and the latter 2 to 6 eigenvalues must settle at a lower value (due to small cumulative errors). Results are recorded in table 2 and plotted in figure 5. To determine the fidelity of the spectral decomposition, the retrieved eigenphase–eigenstate pairs, (θk , |ν k ) were used to create the matrix,  ei2πθk |νk νk |. (6) Uretrieved = k Letting M = UH† 2 O Uretrieved , the fidelity is defined as fidelity = Tr(MM † ) + |Tr(M)|2 n · (n + 1) (7) following the average fidelity definition of [25] where n is matrix dimension (i.e. n = 16). The reported phase error is the average absolute phase error over all 16 phases,  |2πθk − 2πθk,true | . (8) phase error = k n Note that the dc = 2 was run for two different sets of Cgoal and Creq . Increasing these values increased the optimization failure rate, but also improved decomposition fidelity and reduced the average eigenphase error. The high failure rate suggests superior results will be achieved by increasing dc , the number of control Figure 5. (Left axis; black; circles) mean fidelity achieved for optimization simulation at control level dc , with error bars from the 25th to 75th percentile. (Right axis; blue; squares) mean phase error (per phase) for optimization simulation at control level dc , with error bars from the 25th to 75th percentile. levels, over increasing C ∗ (analogous to Cgoal ), when possible. This is expected, as increasing dc leads to a narrower cost function. Overall, these results indicate both the viability of the SPEA for full and partial eigenphase recovery and provides an example of a quantum algorithm which benefits from working with high-dimensional quantum states, i.e. qudits. 4.3. Scaling to higher dimensions The classical controller is a modified gradient descent algorithm and in this work has not been optimized for scaling to very high-dimensional systems. Here the controller serves to demonstrate that the function C —the output of the quantum circuit—is usable as a (negative) cost function for any classical optimizer. Nevertheless, we consider randomly generated systems from 8 to 512 dimensions under control levels dc = 2, 3, 4, 8 to briefly observe how the SPEA may scale. The classical controller was set to run with no value for the stopping condition. Instead, the value of C ∗ over the course of a single optimization was recorded for iterations 5, 20, 100, and 500 for randomly generated system U and a randomly selected starting state. One iteration of this process is considered one trial. Over all trials, the average values of 1 − C ∗ are reported in table 3. Each control level dc = 2, 3, 4, and 8 uses the same set of random systems and starting points. Different trials may drive the optimizer to different eigenstates, due to inherent randomness in the controller. After iteration 500, the absolute difference between the final estimated eigenphase and the true eigenphase is recorded. This average final absolute phase error is reported. Under an ideal optimization, C ∗ will approach 1 as rapidly as possible. It should be noted (see appendix A) that C ∗ should not be directly compared across different values of dc . The higher the value of dc the lower C ∗ will be for a given estimated eigenstate–eigenphase pair. See that, for a given du , the final phase error decreases as dc increases, whether or not C ∗ increases. Regardless, see that for du  64, higher values of dc may cause a lower initial value for C ∗ , but generally achieves a lower C ∗ near the end of the optimization process. The final C ∗ values for all dc for a given du are within an order of magnitude in all cases. All tested dimensions are able to achieve a final C ∗ better than 0.997. These results suggest the system is capable of scaling favorably to high-dimensional systems. The discussion section addresses improvements to the classical controller for realizing this scalability. 5. Discussion Section 4 demonstrates that C of equation (4) can be used effectively as a cost function to determine an eigenphase–eigenstate pair of an arbitrary unitary system despite the simplistic nature of the control algorithm (outlined in section 3). No prior system knowledge is required for successful convergence. The process was reported here for fundamental unitary gates, unitary systems relevant to computational chemistry, and entirely arbitrary choices of U. For large U, concern over the coherence time of the operations which need to take place can be reduced by systems which decompose U into a set of smaller Table 3. For du -by-du dimensional randomly generated unitary systems U, the quantity 1 − C ∗ is averaged over multiple trials at for SPEA systems with varying control levels dc . Quantity 1 − C ∗ is reported at iterations 5, 20, 100, 500 over the course of an optimization run with no stopping conditions. After iteration 500, the mean absolute phase error is reported in radians. Average 1 − C ∗ du 8 16 32 64 128 256 512 dc 2 3 4 8 2 3 4 8 2 3 4 8 2 3 4 8 2 3 4 8 2 3 4 8 2 3 4 8 Trials 100 100 100 100 100 100 26 5 20 100 500 Mean final phase error 3.84 × 10−2 4.36 × 10−2 4.71 × 10−2 6.37 × 10−2 8.39 × 10−3 6.63 × 10−3 3.97 × 10−3 2.84 × 10−3 0.69 × 10−4 0.38 × 10−4 0.23 × 10−4 0.17 × 10−4 3.43 × 10−5 1.63 × 10−5 1.01 × 10−5 0.48 × 10−5 0.04 × 10−2 0.02 × 10−2 0.01 × 10−2 0.005 × 10−2 4.16 × 10−2 6.38 × 10−2 8.93 × 10−2 14.95 × 10−2 8.04 × 10−3 7.45 × 10−3 6.51 × 10−3 6.93 × 10−3 1.29 × 10−4 0.67 × 10−4 0.48 × 10−4 0.17 × 10−4 11.76 × 10−5 3.24 × 10−5 1.73 × 10−5 0.82 × 10−5 0.33 × 10−2 0.04 × 10−2 0.02 × 10−2 0.01 × 10−2 5.03 × 10−2 11.12 × 10−2 18.44 × 10−2 31.71 × 10−2 9.11 × 10−3 10.46 × 10−3 12.17 × 10−3 22.15 × 10−3 1.48 × 10−4 1.4 × 10−4 1.16 × 10−4 0.95 × 10−4 12.79 × 10−5 18.31 × 10−5 10.86 × 10−5 2.34 × 10−5 0.54 × 10−2 0.45 × 10−2 0.17 × 10−2 0.01 × 10−2 8.74 × 10−2 23.27 × 10−2 33.65 × 10−2 49.44 × 10−2 16.45 × 10−2 35.69 × 10−2 44.98 × 10−2 58.69 × 10−2 10.48 × 10−3 25.76 × 10−3 96.27 × 10−3 154.09 × 10−3 24.17 × 10−3 175.24 × 10−3 135.54 × 10−3 137.72 × 10−3 1.79 × 10−4 2.27 × 10−4 2.32 × 10−4 1.54 × 10−4 2.26 × 10−4 4.59 · 10−4 4.68 × 10−4 4.78 × 10−4 27.74 × 10−5 20.49 × 10−5 38.2 × 10−5 5.97 × 10−5 34.45 × 10−5 65.51 × 10−5 54.59 × 10−5 59.57 × 10−5 1.12 × 10−2 0.50 × 10−2 0.63 × 10−2 0.04 × 10−2 1.12 × 10−2 0.98 × 10−2 0.78 × 10−2 0.41 × 10−2 25.29 × 10−2 43.44 × 10−2 51.71 × 10−2 65.86 × 10−2 32.12 × 10−2 49.53 × 10−2 58.17 × 10−2 73.14 × 10−2 156.88 × 10−3 124.12 × 10−3 193.64 × 10−3 400.17 × 10−3 146.61 × 10−3 288.44 × 10−3 398.08 × 10−3 411.68 × 10−3 5.37 × 10−4 9.71 × 10−4 11.82 × 10−4 20.37 × 10−4 19.61 × 10−4 22.36 × 10−4 51.21 × 10−4 26.01 × 10−4 60.39 × 10−5 110.67 × 10−5 136.93 × 10−5 119.06 × 10−5 117.08 × 10−5 227.46 × 10−5 239.06 × 10−5 173.41 × 10−5 0.68 × 10−2 0.64 × 10−2 0.59 × 10−2 0.48 × 10−2 0.44 × 10−2 0.44 × 10−2 0.37 × 10−2 0.28 × 10−2 gates whose operations may take place concurrently [15, 26, 27]. Given the controlled-nature of the PEA, a similar advantage may not be available to the SPEA, or its implementation may be more complex. However, the SPEA only requires a single control gate—differentiating it from the traditional PEA which has a number of operations scaling with 1/p where p is the desired eigenphase precision—allowing the SPEA to be more competitive with other circuit-depth sensitive variational approaches. Other advantages common to variational algorithms can be applied to the SPEA—notably the SPEA classical controller need not know the complete O(2d) parameters defining the input quantum state, only a (potentially) smaller number of parameters necessary to prepare the state, avoiding the N-representability problem. Like a more select number of variational solvers, the SPEA offers self-verification of its results (the proximity of C to one) [27]. The SPEA is not designed around systems which are representable with a polynomial number of terms [15, 27]; the SPEA may be more suitable for non-representable systems than other variational approaches. As C is the expectation value of a quantum circuit’s output, the algorithm as outlined in this work asks for an increasingly large number of runs of the quantum circuit in order to estimate C with sufficient accuracy as C approaches one. Such requirements can be mitigated with a more sophisticated classical controller, which may employ more efficient sampling techniques (such as a Metropolis–Hastings algorithm) or make better use of the information available from C (e.g. driving C(|ψ, θ) → 1 is closely related to forcing C(|ψ, θ + d1c ) → 0). The issue can also be ameliorated by substituting the control-U for control-Uk (for k > 1) as C approaches one, decreasing C for a given estimated eigenphase–eigenstate pair and thus allowing it to be estimated with fewer calls to the quantum circuit. The quality of the eigenstate–eigenphase retrieval may be more accurately defined by using Uk or (rather than following appendix A) by fitting the distribution of C(|ψ, θ + x) over x to the ideal distribution of Pθ (x). 6. Conclusion In this work, we have proposed a novel statistical variational approach (SPEA) to the quantum PEA. From the probabilistic output of a PEA circuit using non-eigen input states, we have defined a statistical metric C indicating the proximity of any given input state to the nearest eigenstate and develop an optimization process that can variationally retrieve all the eigenstate–eigenphase pairs of a given unitary operator. The SPEA takes advantage of the hardware intended for the IPEA and therefore requires no novel quantum hardware development. The main disadvantage of the SPEA is the non-deterministic nature of the measurements requires running the quantum circuit repeatedly for each measurement setting. However, in the near-term era, repeated runs of a quantum circuit per measurement is already the norm, due to noise and imperfect gate fidelity. One of the main advantages of the SPEA compared to the PEA and IPEA is the ability to systematically find both the eigenstates and associated eigenphases, rather than just the eigenphases. The simulations on the IBM Q platform with Qiskit proves the feasibility of applying the SPEA on standard quantum computation platforms. On the local computer, the full spectral decomposition of the operator generated from the water molecule Hamiltonian demonstrates the viability of the SPEA for applications in quantum chemistry. The ability to retrieve eigenstates and efficiency (in terms of low iterations requirement) of this method shows the improvement to the original PEA methods and offers the clear potential to work with larger physical and chemical systems. Further, SPEA is demonstrated to work with arbitrary U and does not require that the U be expanded into a sum of simpler operations. As compared to variational approaches that search only for the ground state, the SPEA is able to search for any eigenstate with no modifications to the quantum circuit. Future work includes improving the optimization process with a more sophisticated algorithm for the classical controller. In addition to improving efficiency and failure rate, this may also improve the accuracy of the eigenphase–eigenstate retrieval as well as the fidelity of the full spectral decomposition. The efficiency and the viability of our methods enable us to simulate more complex systems in quantum chemistry. Future work also includes implementing this method on real computational systems provided by the IBM Q and also on a photonic platform with high-dimensional control qudit capabilities. Data availability statement All data that support the findings of this study are included within the article (and any supplementary files). Appendix A. Quality of eigenphase–eigenstate retrieval Preliminary note: all phases θ ∈ [0, 1). The difference between two phases can be found via the function d(θa , θb ) ≡ arg(ei2π(θa −θb ) ) , 2π (A.1) where d(θa , θb ) ∈ [−.5, .5). In short, as the phase wraps around the phase difference also wraps around e.g. . And d( 34 , 18 ) = −3 , not 58 . Like usual differences, d(θa , θb ) = −d(θb , θa ). Appreciate d( 18 , 34 ) = 38 , not −5 8 8 that—for P0 of equation (3)—P0 (d(θa , θb )) = P0 (θa − θb )∀θa , θb . Here, we quantify the quality of the eigenstate |Φ∗ and eigenphase θR∗ using C ∗ = C(|Φ∗ , θR∗ ). When ∗ C is greater than the largest (non-global) local maximum of Pθ (of equation (3)), then θR∗ must be within the primary lobe of Pθ (examples shown in figure 3). That is, when C∗ > max P0 (ζ)  ζ∈ d1 ,0.5 (A.2) c we are within the primary lobe of Pθ . Let θk∗ be the eigenvalue closest to θR and define δ ∗ ≡ d(θk∗ , θR ). When equation (A.2) is true, then P0 (δ ∗ )  P0 (θk − θR )∀k and C∗ = d t −1 |θk |Φ|2 P0 (θk − θR )  k =0 = P0 (δ ∗ ) d t −1 k =0 d t −1 k =0 |θk |Φ|2 = P0 (δ ∗ ) |θk |Φ|2 P0 (δ ∗ ) C ∗  P0 (δ ∗ ). (A.3) As P0 is symmetric and monotonic within the primary lobe, |δ ∗ |  P0−1 (C ∗ ). (A.4) Therefore the estimated eigenphase θR∗ is within ±P0−1 (C ∗ ) of the nearest eigenphase (whenever equation (A.2) is met). Now to quantify the eigenstate estimate. Define a Δ-eigenstate |ν Δ  as any superposition of eigenstates where the corresponding eigenphases are within ±Δ of θR∗ |νΔ  =  αm |νm , (A.5) m where |d(θm , θR∗ )|  Δ∀m  (and where m |αm |2 = 1). That is, |ν Δ  is a superposition of eigenstates (indexed {m}Δ ) that are nearly degenerate: the corresponding eigenphases are all within 2Δ of one another. Proceeding from equation (4) (whenever equation (A.2) holds), C ∗ = C(|Φ∗ , θR∗ ) =   d t −1 |νk |Φ|2 P0 (θk − θR∗ ) k =0 |νk |Φ|2 (1) +   k∈{m}Δ |νk |Φ| 2  |νk |Φ|2 + |νk |Φ|2 P0 (Δ). (A.6) k∈{m} / Δ k∈{m}Δ Letting |νk |Φ|2 P0 (θk − θR∗ ) k∈{m} / Δ k∈{m}Δ   = f, C ∗  f + (1 − f )P0 (Δ) ∴ f  C ∗ − P0 (Δ) 1 − P0 (Δ) (A.7) The estimated eigenstate |Φ∗ matches some Δ-eigenstate (as defined by equation (A.5)) with fidelity f given by equation (A.7) (whenever equation (A.2) holds and |Δ| ∈ [0, d1c ]). Appendix B. Details for the IBM Q SPEA calculations On IBM Q we realize an SPEA with a four-dimensional control register by using two qubits (the top two rails) as controls. The target is either two- or four-dimensional, using the bottom one or two rails, respectively. We list out the three operators in matrix form used in our simulations on the IBM Q accompanied by the eigenstates of each matrix as well as showing how we achieved these matrices with the Qiskit. In the following we use the rotation-Z gate as defined by Qiskit [16]:  −iθ/2 e RZ(θ) = 0 as well as the phase gate: P(θ) = and the Hadamard gate: 1 H=√ 2 0  (B.1) eiθ/2   1 0 0 eiθ  1 1  1 . −1 (B.2) (B.3) The first operator U1 is a single qubit rotation-Z gate with θ = π/2, U1 = RZ(π/2) =  −iπ/4 e 0 0 eiπ/4  (B.4) Figure B1. Gate representation of the unitary operators applied in the simulation. (A) Represents operator U1 , a two-dimensional operator applied to a single qubit. (B) Represents operator U2 , a four-dimensional operator applied to two qubits. RZ is the rotation-Z gate, P is the phase gate and H is the Hadamard gate. When a control qubit is present the RZ gate and the P gate become controlled gates. The Hadamard gates will serve as their own inverse and therefore do not need to be implemented as controlled gates. with eigenstates v1 =   0 , 1 v2 =   1 . 0 The second operator is a two qubit operation achieved by a phase gate P(θ = π/4) acting on the first qubit and a rotation-Z gate RZ(θ = π/2) sandwiched between two Hadamard gates H, acting on the second qubit. The matrix form is  U2 = P(π/4) ⊗ H · RZ(π/2) · H ⎛ ⎞ 1 −i 0 0 0 0 ⎟ 1 ⎜−i 1 ⎟ (B.5) =√ ⎜ i π4 −i π4 ⎠ ⎝ 0 e e 2 0 π π 0 0 e−i 4 ei 4 with eigenstates ⎛ ⎞ 0 ⎜0⎟ ⎟ v3 = ⎜ ⎝1⎠ , 1 ⎛ ⎞ 0 ⎜0⎟ ⎟ v4 = ⎜ ⎝ 1 ⎠, −1 ⎛ ⎞ 1 ⎜1⎟ ⎟ v5 = ⎜ ⎝0⎠ , 0 ⎛ ⎞ 1 ⎜−1⎟ ⎟ v6 = ⎜ ⎝ 0 ⎠. 0 (B.6) The gate representations of the two operators can be found in figure B1. For the third operators we start with the lowest energy Hamiltonian of the hydrogen molecule generated with Bravyi–Kitaev transformation in STO − 3G basis, ⎛ ⎞ 0.487 049 0 0 0.180 653 ⎜ ⎟ 0 −0.337610 0.180653 0 ⎟. H=⎜ (B.7) ⎝ ⎠ 0. 0.180 653 −0.337 610 0 0.180 653 0 0 −1.117 194 Then we take the exponential of the Hamiltonian to generate our unitary operator ⎞ ⎛ 0.8686 + 0.4687i 0 0 0.0499 + 0.1531i ⎜ ⎟ 0 0.9282 − 0.3259i 0.0595 + 0.1695i 0 ⎟. U3 = eiH = ⎜ ⎝ ⎠ 0 0.0595 + 0.1695i 0.9282 − 0.3259i 0 0.0499 + 0.1531i 0 0 0.4256 − 0.8905i (B.8) The exponential of the Hamiltonian is done with the scipy.linalg.expm in the scipy python package [28]. The eigenstates are ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ −0.1105 0 0 −0.9939 ⎜ 0 ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎟ , v8 = ⎜ 0.7071 ⎟ , v9 = ⎜0.7071⎟ , v10 = ⎜ 0 ⎟ . (B.9) v7 = ⎜ ⎝ 0 ⎠ ⎝−0.7071⎠ ⎝0.7071⎠ ⎝ 0 ⎠ 0.9939 0 0 0.1105 We design the phase estimation algorithms that work with up to two qubits in the control register and four qubits in total as shown in figure B2. Due to the restrictions in the qubit numbers, the simulations on the IBM Q are focused on the lower dimensional systems, i.e. one or two qubits in the target register. The rotation RZ gates are applied to each qubit in the control register after the controlled-Ui operations but before inverse Fourier transform. Every time the quantum algorithm is called to generate a new C factor (following the ‘alternative method’ in step (c) of section 3), the algorithm runs twice: the first time the RZ Figure B2. Illustration of the four-level control (dc = 4) SPEA circuit implemented using Qiskit. H is the Hadamard gate, RZ is the rotation-Z gate, P is the phase gate, and M represents a measurement. The top two rails are the control qubits and realize a four-dimensional control register. The bottom two rails represents the target register. The initialize block prepares the input state differently throughout the optimization process. Ui represents the unitary the SPEA is running. For i = 2, 3 we implement the circuit as diagrammed. For i = 1, the unitary is two-dimensional and we use only a single target rail. The three control gates together realize a multi(four)-level control gate as described in the main text. The two RZ-gates together realize a realize four-level Rz (−θR ) as described in the main text. The double line at the bottom represents the classical information retrieved from the measurement gates (a total of 2 bits). gates are set to zero and the phase factor θR is calculated statistically; the second time the upper RZ gate applies phase −θR and the lower gate applies phase −2θR , together acting as a −θR -rotation would on a dc = 4 qudit system. The classical optimization process described in section 3 is implemented using python. The basis set {|Bm } is generated using the Gram–Schmidt methods with the input vector plus a set of linearly independent vectors obtained from a randomly generated unitary matrix. As a deviation from how the algorithm is described in the main text, the search step (d) a factor is set to 1/2 in step (d) and is doubled in step (e) (rather than halved) if C-factor is not updated, up to seven times. The optimization concludes when the C factor meets the stopping condition which means the input state is converged to an eigenstate, or when the maximum iteration time is exceeded. Appendix C. Full H2 O matrix The Hamiltonian of the water molecule with the H–O–H angle at 104.5◦ and the bond length at 1.0 a.u. is calculated by STO − 3G minimal basis using PySCF [22] and chemistry package provided by the Qiskit [16]. The 16-by-16 Hamiltonian of the water molecule [24] used for the local computer’s spectral decomposition simulations is sparse with nonzero entries, HH2 O (1, 1) = 1.027 · 10−15 , HH2 O (2, 2) = −2.594, HH2 O (3, 3) = −2.654, HH2 O (4, 4) = −4.583, HH2 O (5, 5) = −2.594, HH2 O (6, 6) = −4.427, HH2 O (7, 7) = −4.529, HH2 O (8, 8) = −5.696, HH2 O (9, 9) = −2.654, HH2 O (10, 10) = −4.529, HH2 O (11, 11) = −4.428, HH2 O (12, 12) = −5.637, HH2 O (13, 13) = −4.583, HH2 O (14, 14) = −5.696, HH2 O (15, 15) = −5.637, HH2 O (16, 16) = −6.085, HH2 O (6, 11) = HH2 O (11, 6) = 0.054, HH2 O (7, 10) = HH2 O (10, 7) = 0.054. (C.1) The exponential of the Hamiltonian is done with the MATLAB expm function. References [1] Kais S, Whaley K B, Dinner A R and Rice S A 2014 Quantum information and computation for chemistry Advances in Chemical Physics (New York: Wiley) [2] Liberty E, Woolfe F, Martinsson P-G, Rokhlin V and Tygert M 2007 Randomized algorithms for the low-rank approximation of matrices Proc. Natl Acad. Sci. 104 20167–72 [3] Daskin A and Kais S 2011 Decomposition of unitary matrices for finding quantum circuits: application to molecular Hamiltonians J. Chem. Phys. 134 144112 [4] Nielsen M A and Chuang I 2002 Quantum Computation and Quantum Information (Cambridge: Cambridge University Press) [5] Shor P W 1994 Algorithms for quantum computation: discrete logarithms and factoring Proc. 35th Annual Symp. Foundations of Computer Science (Piscataway, NJ: IEEE) pp 124–34 [6] Lloyd S, Mohseni M and Rebentrost P 2014 Quantum principal component analysis Nat. Phys. 10 631–3 [7] Byrnes T, Forster G and Tessler L 2018 Generalized Grover’s algorithm for multiple phase inversion states Phys. Rev. Lett. 120 060501 [8] Wang H, Kais S, Aspuru-Guzik A and Hoffmann M R 2008 Quantum algorithm for obtaining the energy spectrum of molecular systems Phys. Chem. Chem. Phys. 10 5388–93 [9] Aspuru-Guzik A, Dutoi A D, Love P J and Head-Gordon M 2005 Simulated quantum computation of molecular energies Science 309 1704–7 [10] Daskin A and Kais S 2018 Direct application of the phase estimation algorithm to find the eigenvalues of the Hamiltonians Chem. Phys. 514 87–94 [11] Preskill J 2018 Quantum computing in the NISQ era and beyond Quantum 2 79 [12] Dobšíček M, Johansson G, Shumeiko V and Wendin G 2007 Arbitrary accuracy iterative quantum phase estimation algorithm using a single ancillary qubit: a two-qubit benchmark Phys. Rev. A 76 030306(R) [13] Paesani S, Gentile A A, Santagati R, Wang J, Wiebe N, Tew D P, O’Brien J L and Thompson M G 2017 Experimental Bayesian quantum phase estimation on a silicon photonic chip Phys. Rev. Lett. 118 100503 [14] Wang D, Higgott O and Brierley S 2019 Accelerated variational quantum eigensolver Phys. Rev. Lett. 122 140504 [15] Peruzzo A, McClean J, Shadbolt P, Yung M-H, Zhou X-Q, Love P J, Aspuru-Guzik A and O’Brien J L 2014 A variational eigenvalue solver on a photonic quantum processor Nat. Commun. 5 4213 [16] Aleksandrowicz G et al 2019 Qiskit: an open-source framework for quantum computing Zenodo https://doi.org/10.5281/ zenodo.2562111 [17] LaRose R, Tikku A, O’Neel-Judy É, Cincio L and Coles P J 2019 Variational quantum state diagonalization npj Quantum Inf. 5 57 [18] Cerezo M, Sharma K, Arrasmith A and Coles P J 2020 Variational quantum state eigensolver (arXiv:2004.01372) [19] O’Malley P J J et al 2016 Scalable quantum simulation of molecular energies Phys. Rev. X 6 031007 [20] Lu H H, Hu Z, Alshaykh M S, Moore A J, Wang Y, Imany P, Weiner A M and Kais S 2020 Quantum phase estimation with time-frequency qudits in a single photon Adv. Quantum Technol. 3 1900074 [21] Seeley J T, Richard M J and Love P J 2012 The Bravyi–Kitaev transformation for quantum computation of electronic structure J. Chem. Phys. 137 224109 [22] Sun Q et al 2018 Pyscf: the Python-based simulations of chemistry framework WIREs Comput. Mol. Sci. 8 e1340 [23] McClean J R et al 2020 OpenFermion: the electronic structure package for quantum computers Quantum Sci. Technol. 5 034014 [24] Bian T, Murphy D, Xia R, Daskin A and Kais S 2019 Quantum computing methods for electronic states of the water molecule Mol. Phys. 117 2069–82 [25] Pedersen L H, Møller N M and Mølmer K 2007 Fidelity of quantum operations Phys. Lett. A 367 47–51 [26] McClean J R, Romero J, Babbush R and Aspuru-Guzik A 2016 The theory of variational hybrid quantum–classical algorithms New J. Phys. 18 023023 [27] Kokail C et al 2019 Self-verifying variational quantum simulation of lattice models Nature 569 355–60 [28] Virtanen P et al ((SciPy 1.0 Contributors) 2020 SciPy 1.0: fundamental algorithms for scientific computing in Python Nat. Methods 17 261–72