Statistical Approach to Quantum Phase Estimation

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.


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,9,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 controlledgate 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 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,18,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 ondemand 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 iterative PEA and introduces a statistical metric C for quantifying the proximity of any given input-state to its closest eigenstate. Section 3 describes the Statistical PEA 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 Appendix B and C. We conclude with a discussion on the performance of the SPEA and propose future directions and applications of the method in Section 5.

Phase Estimation Algorithms
Traditional PEA implementations, diagrammed in Figure 1, take any given unitaryÛ and any given eigenstate |ν ofÛ and return the corresponding eigenphase θ wherê The (approximate) eigenphaseθ ∈ [0, 1) (equivalently, ∈ [−.5, .5)) may be directly measured on the control qubits (or qudits, when the control is d c -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 iterative PEA (IPEA) was developed. The IPEA significantly reduces circuit depth requirements by approximating a particular θ one qubit (or d c -level qudit) at a time, starting from the least significant qubit (qudit). The IPEA requires a rotation gate -a linear phase across the control control register target register H . Figure 1: The traditional PEA using n control d c -dimensional qudits in the control register. The quantum gates are colored in orange and the measurement gates in red. The target register is highlighted in purple and control register in grey. Note that in the d c > 2 case, the H "Hadamard" gate acts as a discrete quantum Fourier transform (QFT) gate. Additionally, in the d c > 2 case the control gates acts as MVCGs, applying the gate to the q th power when the control qudit is in state |q . The circuit will estimate the phaseθ k of eigenstate |ν k to precision d Figure 2: The iterative PEA. See Fig. 1 for color-conventions and notes on the MVCG. To retrieve n dits of the eigenphase (i.e. phase precision ±1/d n c ), run the circuit with x = n − 1 and θ R = 0; the measured control state is the n th base-d c 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 Note that P 0 has infinite domain with period 1. Probability of the (d c -level) control register of an iterative PEA 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 P 0 (∆θ), 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 d c = 2 has no sidelobes (the probability is monotonic on either side of the central lobe). Also note P(∆θ) = 0 for ∆θ = d −1 c and the width of the central lobe is therefore ∆θ FWFM = 2d −1 c . register -to "subtract" off eigenphase information determined in previous iterations. (I.e. if the quantum circuit's state before R z (θ R ) is q α q |q |Φ q , then after the rotation gate the quantum circuit's state is q α q e −iq2πθ R |q |Φ q .) The iterative PEA, as the name suggests, requires a number of iterations equal to the number 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 highdimensional (d c −level). In this case, the Hadamard gates represent a d c -dimensional quantum Fourier transform gate and the control-Û gate is a multi-level control gate (MLCG) [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 d c -level control, the final state of the control qudit before measurement is 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 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 to |0 is WhereÛ is d t -by-d t -dimensional and the target register is d t -dimensional. Appreciate that C(|Φ , θ R ) = 1 if and only if |Φ is an eigenstate and θ R is its corresponding eigenphase.

Statistical Approach to PEAs
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 single-dit precision (n = 1) and the rotation gate standard to the IPEA (i.e. an iterative PEA 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 probabilityC ≈ 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 (makingC the negative cost function), the classical controller adjusts |Φ and θ R , until the quantum circuit near-deterministically returns |0 as the output state. Wheñ C(|Φ * , θ * R ) ≈ 1, the classical controller has found the (approximate) eigenstate |Φ * and the associated 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: 1. The classical controller chooses a |Φ at random 2. The classical controller constructs an orthogonal basis {|B m } including |Φ

Quantum System
State Preparation |Φ PEA-type Circuit: 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 orange and measurement gates in red. The (potentially high-dimensional) PEA-type circuit is simplified from the typical iterative PEA in that U need not be raised to high orders (d x c ) 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.

3.
• (Standard Method: viable when the quantum circuit can measure output bins |0 and not(|0 )) The quantum circuit evaluatesC(|Φ , θ R ) over a range of θ R and returns the maximum value C * • (Alternative Method: viable when the quantum circuit can measure all d c output bins: |0 through |d c − 1 .) The quantum circuit evaluatesC(|Φ , 0) and uses this result to approximate the eigenphase θ * . The quantum circuit then evaluatesC(|Φ , θ R = θ * ) and returns C * .
4. For all m = 0 to 2d t − 1, we set a=1 and the following occurs: • the classical controller generates the new state: • |Φ is fed to the quantum circuit, the maximum value returned is C * • if C * > C * , then |Φ = |Φ and C * = C * . Otherwise |Φ and C * are unchanged.
5. If |Φ was not updated during step 4, set a = a/2 and repeat step 4.
6. 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 2.
A few observations on the optimization process may be made. For each iteration, at least 2d t distinct input states are used. For each of these input states a set of {θ R } is applied (when using the 'standard approach' in step 3). 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 {|B m } (step 2) 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.

Statistical PEA 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 2-and 4-dimensional systems. The local computer simulations run a full spectral decomposition on a 16-dimensional system with various control levels d c .
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 3 whereas the IBM Q simulations follow the alternative method. The IBM Q simulation runs one measurement with R z (θ 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 R z (θ 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 R z 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 control level (|0 , |1 , ..., |d c − 1 ). The alternative approach (or some hybrid approach) is generally preferable if the hardware is available for d c detectors.

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 2-dimensional operator U 1 and the other two with 4-dimensional operators U 2 and U 3 , the matrix forms of which are shown in Appendix B. U 1 and U 2 are operators directly built with the default gates offered by the IBM Q and U 3 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 ST O − 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 U i 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 U 1 , U 2 , U 3 , 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 2. 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. Variations in mean phase error are likely a function of which particular eigenphase-eigenstate pair was converged to. During the simulation of U 3 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 (U 1 , U 2 ) and physically relevant (U 3 ) 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.

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 H H 2 O for the water molecule H 2 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, U H 2 O = e iH H 2 O , and the matrix's spectral decomposition is simulated with the statistical variational algorithm (SPEA) on a local computer.
The simulation is run for various control levels d c until 120 successful spectral decompositions are achieved. Each of the 16 eigenphases are retrieved in a random order. The optimization runs until C * ≥ C goal . If the optimizer is unable to reach C goal , the process for that eigenphase will conclude so long as C * ≥ C req . If C req is not met, the entire spectral decomposition is abandoned and the trial is classified as failed. Generally, C goal 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 )    Table 2: Statistics for 120 successful trials of complete spectral decomposition of U H2O matrix. C goal is theC value the optimizer attempts to reach, and generally does reach for at least the first 10 (of 16) eigen-estimates. C req is thẽ 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 d c = 2 was run on two different C goal , C req levels for comparison. were used to create the matrix, Letting M = U † H 2 O U retrieved , the fidelity is defined as fidelity = 1 n · (n + 1) T r(M M † ) + |T r(M )| 2 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, phase error = k |2πθ k − 2πθ k,true | n .
Note that the d c = 2 was run for two different sets of C goal and C req . 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 d c , the number of control levels, over increasing C * (analogous to C goal ), when possible. This is expected, as increasing d c 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.

Conclusion
In this work, we have proposed a novel statistical variational approach (SPEA) to the quantum phase estimation algorithm (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 Iterative PEA 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.
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.
A Quality of Eigenphase-Eigenstate retrieval Preliminary note: All phases θ ∈ [0, 1). The difference between two phases can be found via the function I.e. d(θ a , θ b ) ∈ [−.5, .5). In short, as the phase wraps around the phase difference also wraps around. E.g.
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 we are within the primary lobe of P θ . Let θ k * be the eigenvalue closest to θ R and define δ * ≡ d(θ k * , θ R ). When Equation 10 is true, then P 0 (δ * ) ≥ P 0 (θ k − θ R )∀k and As P 0 is symmetric and monotonic within the primary lobe, Therefore the estimated eigenphase θ * R is within ±P −1 0 (C * ) of the nearest eigenphase (whenever Equation 10 is met). Now to quantify the eigenstate estimate. Define a ∆-eigenstate |ν ∆ as any superposition of eigenstates where the corresponding eigenphases are within ±∆ of θ * (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 10 holds), Letting k∈{m} ∆ | ν k |Φ | 2 = f , The estimated eigenstate |Φ * matches some ∆-eigenstate (as defined by Equation 13) with fidelity f given by Equation 15 (whenever Equation 10 holds and |∆| ∈ [0, 1 dc ]).

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]: as well as the phase gate: and the Hadamard gate: The first operator U 1 is a single qubit rotation-Z gate with θ = π/2, 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 with the eigenstates The gate representations of the two operators can be found in Figure 6.
Then we take the exponential of the Hamiltonian to generate our unitary operators The exponential of the Hamiltonian is done with the scipy.linalg.expm in the scipy python package [26]. The eigenstates are 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 7. 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-U i 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 3 of Section 3), the algorithm runs twice: the first time the RZ 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 d c = 4 qudit system. The classical optimization process described in Section 3 is implemented using python. The basis set {|B m } 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 a factor is set to 1/2 in step 4 and is doubled Figure 7: Illustration of the 4-level control (d c = 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 4-dimensional control register. The bottom two rails represents the target register. The initialize block prepares the input state differently throughout the optimization process. U i represents the unitary the SPEA is running. For i = 2, 3 we implement the circuit as diagrammed. For i = 1, the unitary is 2-dimensional and we use only a single target rail. The three control gates together realize a multi(4)-level control gate as described in the main text. The two RZ-gates together realize a realize four-level R z (−θ 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). in step 5 (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.

C Full H 2 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 ST O − 3G minimal basis using PySCF [22] and chemistry package provided by the Qiskit [16]. The 16-by-16 Hamiltonian of the water molecule used for the local computer's spectral decomposition simulations is as follows [24]. The exponential of the Hamiltonian is done with the MATLAB expm funtion.