Quantum Machine Learning for Electronic Structure Calculations
Rongxin Xia Sabre Kais


1]Rongxin Xia
1,2,3]Sabre Kais kais@purdue.edu
[1]Department of Physics and Astronomy, Purdue University, West Lafayette, IN, 47907 USA
[2]Department of Chemistry and Birck Nanotechnology Center, Purdue University,
West Lafayette, IN 47907 USA
[3]Santa Fe Institute, 1399 Hyde Park Rd, Santa Fe, NM 87501




Quantum Machine Learning for Electronic Structure Calculations
    [
    
==============================================================






Considering recent advancements and successes in the development of efficient quantum algorithms for electronic structure calculations — alongside impressive results using machine learning techniques for computation — hybridizing quantum computing with machine learning for the intent of performing electronic structure calculations is a natural progression. Here we report a hybrid quantum algorithm employing a restricted Boltzmann machine to obtain accurate molecular potential energy surfaces. By exploiting a quantum algorithm to help optimize the underlying objective function, we obtained an efficient procedure for the calculation of the electronic ground state energy for a small molecule system. Our approach achieves high accuracy for the ground state energy for H_2, LiH, H_2O at a specific location on its potential energy surface with a finite basis set.  With the future availability of larger-scale quantum computers, quantum machine learning techniques are set to become powerful tools to obtain accurate values for electronic structures. 




§ INTRODUCTION

Machine learning techniques are demonstrably powerful tools displaying remarkable success in compressing high dimensional data <cit.>. These methods have been applied to a variety of fields in both science and engineering, from computing excitonic dynamics  <cit.>, energy transfer in light-harvesting systems  <cit.>, molecular electronic properties <cit.>, surface reaction network  <cit.>, learning density functional models <cit.> to classify phases of matter, and the simulation of classical and complex quantum systems <cit.>. Modern machine learning techniques have been used in the state space of complex condensed-matter systems for their abilities to analyze and interpret exponentially large data sets <cit.> and to speed-up searches for novel energy generation/storage materials <cit.>. 

Quantum machine learning <cit.> - hybridization of classical machine learning techniques with quantum computation – is emerging as a powerful approach allowing quantum speed-ups and improving classical machine learning algorithms 
 <cit.>. Recently, Wiebe et. al. <cit.> have shown that quantum computing is capable of reducing the time required to train a restricted Boltzmann machine (RBM), while also providing a richer framework for deep learning than its classical analogue. The standard RBM models the probability of a given configuration of visible and hidden units by the Gibbs distribution with interactions restricted between different layers. Here, we focus on an RBM where the visible and hidden units assume {+1,-1} forms <cit.>. 
 
 Accurate electronic structure calculations for large systems continue to be a challenging problem in the field of chemistry and material science. Toward this goal — in addition to the impressive progress in developing classical algorithms based on ab initio and density functional methods — quantum computing based simulation have been explored <cit.>. Recently, Kivlichan et. al. <cit.> show that using a particular arrangement of gates (a fermionic swap network) it is possible to simulate electronic structure Hamiltonian with linear depth and connectivity. These results present significant improvement on the cost of quantum simulation for both variational and phase estimation based quantum chemistry simulation methods. 

Recently, Troyer and coworkers proposed using a restricted Boltzmann machine to solve quantum many-body problems, for both stationary states and time evolution of the quantum Ising and Heisenberg models <cit.>. However, this simple approach has to be modified for cases where the wave function's phase is required for accurate 
calculations <cit.>.

Herein, we propose a three-layered RBM structure that includes the visible and hidden layers, plus 
a new layer correction for the signs of coefficients for basis functions of the wave function.
We will show that this model has the potential to solve complex quantum many-body problems and to obtain very accurate results for simple molecules as compared with the results calculated by a finite minimal basis set, STO-3G. We also employed a quantum algorithm to help the optimization of training procedure.




§ RESULTS



Three-layers restricted Boltzmann machine. We will begin by briefly outlining the original RBM structure as described by  <cit.>. For a given Hamiltonian, H, and a trial state, |ϕ⟩ = ∑_xϕ(x)|x⟩, the expectation value can be written as<cit.>:

    ⟨ H⟩=⟨ϕ|H|ϕ⟩/⟨ϕ|ϕ⟩=∑_x,x'⟨ϕ|x⟩⟨ x|H|x'⟩⟨ x'|ϕ⟩/∑_x⟨ϕ|x⟩⟨ x|ϕ⟩=∑_x,x'ϕ(x)⟨ x|H|x'⟩ϕ(x')/∑_x|ϕ(x)|^2

where ϕ(x) = ⟨ x|ϕ⟩ will be used throughout this letter to express the overlap of the complete wave function with the basis function |x⟩, ϕ(x) is the complex conjugate of ϕ(x).

We can map the above to a RBM model with visible layer units 1, 2... n and hidden layer units h_1, h_2... h_m with i, h_j∈{-1,1}. We use visible units i to represent the spin state of a qubit i – up or down. The total spin state of n qubits is represented by the basis |x⟩=|12...n⟩. ϕ(x) = √(P(x)) where P(x) is the probability for x from the distribution determined by the RBM. The probability of a specific set x={1,2...n} is:


    P(x) = ∑_{h}e^(∑_ia_ii+∑_jb_jh_j+∑_i,jw_ijih_j)/∑_x'∑_{h}e^(∑_ia_ii+∑_jb_jh_j+∑_i,jw_ijih_j)


Within the above a_i and b_j are trainable weights for units i and h_j. w_ij are trainable weights describing the connections between i and h_j (see Figure 1.)

By setting ⟨ H⟩ as the objective function of this RBM, we can use the standard gradient decent method to update parameters, effectively minimizing ⟨ H⟩ to obtain the ground state energy. 




However, previous prescriptions considering the use of RBMs for electronic structure problems have found difficulty as ϕ(x_i) can only be non-negative values. 
We have thus appended an additional layer to the neural network architecture to compensate for the lack of sign features specific to electronic structure problems.

We propose an RBM with three layers. The first layer, σ^z, describes the parameters building the wave function. The h's within the second layer are parameters for the coefficients for the wave functions and the third layer s, represents the signs associated |x⟩:


    s(x) = s(1, 2...n) = tanh(∑_id_ii+c)





The s uses a non-linear function tanh to classify whether the sign should be positive or negative.  Because we have added another function for the coefficients, the distribution is not solely decided by RBM. We also need to add our sign function into the distribution.  Within this scheme, c is a regulation and d_i are weights for i. (see Figure 1). Our final objective function, now with |ϕ⟩ = ∑_xϕ(x)s(x)|x⟩, becomes:


    ⟨ H⟩ =∑_x,x'ϕ(x)s(x)⟨ x|H|x'⟩ϕ(x')s(x')/∑_x|ϕ(x) s(x)|^2


After setting the objective function, the learning procedure is performed by sampling to get the distribution of ϕ(x) and calculating to get s(x). We then proceed to calculate the joint distribution determined by ϕ(x) and s(x). The gradients are determined by the joint distribution and we use gradient decent method to optimize ⟨ H ⟩ (see Supplementary
Note 1). Calculating the the joint distribution is efficient because s(x) is only related to x.


Electronic Structure Hamiltonian Preparation. The electronic structure is represented by N single-particle orbitals which can be empty or occupied by a spinless electron<cit.>:


    Ĥ= ∑_i,jh_ija_i^† a_j+1/2∑_i,j,k,l h_ijkla_i^† a_j^† a_ka_l


where h_ij and h_ijkl are one and two-electron integrals. In this study we use the minimal basis (STO-3G) to calculate them. a_j^† and a_j are creation and annihilation operators for the orbital j.

Equation (5) is then transformed to Pauli matrices representation, which is achieved by the Jordan-Wigner transformation<cit.>. The final electronic structure Hamiltonian takes the general form with σ_α^i ∈{σ_x,σ_y,σ_z,I} where σ_x,σ_y,σ_z are Pauli matrices and I is the identity matrix<cit.>:


    H=∑_i,αh_α^iσ_α^i+∑_i,j,α,βh_αβ^ijσ_α^iσ_β^j+∑_i,j,k,α,β,γh_αβγ^ijkσ_α^iσ_β^jσ_γ^k+...




Quantum algorithm to sample Gibbs distribution. We propose a quantum algorithm to sample the distribution determined by RBM. The probability for each combination y ={σ^z, h} can be written as:


    P(y) = e^∑_ia_ii+∑_jb_jh_j+∑_i,jw_ijih_j/∑_y'e^∑_ia_ii+∑_jb_jh_j'+∑_i,jw_ijih_j'


Instead of P(y), we try to sample the distribution Q(y) as:


    Q(y) = e^1/k(∑_ia_ii+∑_jb_jh_j+∑_i,jw_ijih_j)/∑_y'e^1/k(∑_ia_ii+∑_jb_jh_j'+∑_i,jw_ijih_j')


where k is an adjustable constant with different values for each iteration and is chosen to increase the probability of successful sampling. In our simulation, it is chosen as O(∑_i,j|w_ij|).

We employed a quantum algorithm to sample the Gibbs distribution from the quantum computer. This algorithm is based on sequential applications of controlled-rotation operations, which tries to calculate a distribution Q'(y) ≥ Q(y) with an ancilla qubit showing whether the sampling for Q(y) is successful<cit.>.

This two-step algorithm uses one system register (with n+m qubits in use) and one scratchpad register (with one qubit in use) as shown in Figure 2. 

All qubits are initialized as |0⟩ at the beginning. The first step is to use R_y gates to get a superposition of all combinations of {σ^z,h} with θ_i = 2arcsin(√(e^a_i/k/e^a_i/k+e^-a_i/k)) and γ_j=2arcsin(√(e^b_j/k/e^b_j/k+e^-b_j/k)):



    ⊗_iR_y(θ_i)|0_i⟩⊗_jR_y(γ_j)|0_j⟩|0⟩ =
    ∑_y√(O(y))|y⟩|0⟩


where O(y) = e^∑_ia_ii/k+∑_jb_jh_j/k/∑_y'e^∑_ia_ii/k+∑_jb_jh_j'/k and |ϕ_y⟩ corresponds to the combination |y⟩ =|1..nh_1...h_m⟩.

The second step is to calculate e^w_ijih_j. We use controlled-rotation gates to achieve this. The idea of sequential controlled-rotation gates is to check whether the target qubit is in state |0⟩ or state |1⟩ and then rotate the corresponding angle (Figure 2). If qubits i and h_j are in |00⟩ or |11⟩, the ancilla qubit is rotated by R_y(θ_ij,1) and otherwise by R_y(θ_ij,2), with θ_ij,1 = 2arcsin(√(e^w_ij/k/e^|w_ij|/k)) and θ_ij,2 = 2arcsin(√(e^-w_ij/k/e^|w_ij|/k)). Each time after one e^w_ijih_j is calculated, we do a measurement on the ancilla qubit. If it is in |1⟩ we continue with a new ancilla qubit initialized in |0⟩ , otherwise we start over from the beginning (details in  Supplementary Note 2).



After we finish all measurements the final states of the first m+n qubits follow the distribution Q(y).  We just measure the first n+m qubits of the system register to obtain the probability distribution. After we get the distribution, we calculate all probabilities to the power of k and normalize to get the Gibbs distribution.



The complexity of gates comes to O(mn) for one sampling and the qubits requirement comes to O(mn). If considering the reuse of ancilla qubits, the qubits requirements reduce to O(m+n) (see Supplementary Note 4). The probability of one successful sampling has a lower bound e^-1/k∑_i,j2|w_ij| and if k is set to O(∑_i,j|w_ij|) it has constant lower bound (see Supplementary Note 3). If N_s is the number of successful sampling to get the distribution, the complexity for one iteration should be O(N_smn) due to the constant lower bound of successful sampling as well as processing distribution taking O(N_s). In the meantime, the exact calculation for the distribution has complexity as O(2^m+n).  The only error comes from the error of sampling if not considering noise in the quantum computer.



Summary of numerical results. We now present the results derived from our RBM for H_2, LiH and H_2O molecules. It can clearly be seen from Figure 4 that our three layer RBM yields very accurate results comparing to the disorganization of transformed Hamiltonian which is calculated by a finite minimal basis set, STO-3G. Points deviating from the ideal curve are likely due to local minima trapping during the optimization procedure. This can be avoided in the future by implementing optimization methods which include momentum or excitation, increasing the escape probability from any local features of the potential energy surface. 

Further discussion about our results should mention instances of transfer learning. Transfer learning is a unique facet of neural network machine learning algorithms describing an instance (engineered or otherwise) where the solution to a problem can inform or assist in the solution to another similar subsequent problem. Given a diatomic Hamiltonian at a specific intermolecular separation, the solution yielding the variational parameters — which are the weighting coefficients of the basis functions — are adequate first approximations to those parameters at a subsequent calculation where the intermolecular separation is a small perturbation to the previous value.

Except for the last point in the Figure 4 d, we use  1/40 of the iterations  for the last point in calculations initiated with transferred parameters from previous iterations of each points and still achieve a good result. We also see that the local minimum is avoided if the starting point achieve global minimum.




Discussion. In conclusion, we present a combined quantum machine learning approach to perform electronic structure calculations. Here, we have a proof of concept and show results for small molecular systems. Screening molecules to accelerate the discovery of new materials for specific application is demanding since the chemical space is very large! For example, it was reported that the total number of possible
small organic molecules that populate the ‘chemical space’ 
exceed 10^60<cit.>. Such an enormous size makes a thorough exploration of chemical
space using the traditional electronic structure methods impossible. Moreover, in a recent perspective<cit.>in Nature Reviews Materials the potential of machine learning algorithms to accelerate the discovery of materials was pointed out.  Machine learning algorithms have been used for material screening.
For example, out of the GDB-17 data base, consisting of about 166 billion molecular graphs, one can make organic and drug-like molecules with up to 17 atoms and 134 thousand smallest molecules with up to 9 heavy atoms were calculated using hybrid density functional (B3LYP/6-31G(2df,p). Machine learning algorithms trained on these data, were found to predict molecular properties of subsets of these molecules <cit.>.  

In the current simulation, H_2 requires 13 qubits with the number of visible units n=4, the number of hidden units m=8 and additional 1 reusing ancilla qubits . LiH requires 13 qubits with the number of visible units n=4, the number of hidden units m=8 and additional 1 reusing ancilla qubits. H_2O requires 13 qubits with the number of visible units n=6, the number of hidden units m=6 and additional 1 reusing ancilla qubits. The order of scaling of qubits for the system should be O(m+n) with reusing ancilla qubits. The number of visible units n is equal to the number of spin orbitals. The choice of the number of hidden units m is normally integer times of n which gives us a scaling of O(n) with reusing ancilla qubits . Thus, the scaling of the qubits increases polynomially with the number of spin orbitals. Also, the complexity of gates O(n^2) scales polynomially with the number of spin orbitals while the scaling of classical Machine Learning approaches calculating exact Gibbs distribution is exponential. With the rapid development of larger-scale quantum computers and the possible training of some machine units with the simple dimensional scaling results for electronic structure, quantum machine learning techniques are set to become powerful tools to perform electronic structure calculations and assist in designing new materials for specific applications.



§ METHODS


Preparation of the Hamiltonian of H_2, LiH and H_2O. We treat H_2 molecule with 2-electrons in a minimal basis STO-3G and use the Jordan-Wigner transformation<cit.>. The final Hamiltonian is of 4 qubits.  We treat LiH molecule with 4-electrons in a minimal basis STO-3G and use the Jordan-Wigner transformation<cit.>. We assumed the first two lowest orbitals are occupied by electrons and the the final Hamiltonian is of 4 qubits.  We treat H_2O molecule with 10-electrons in a minimal basis STO-3G, we use Jordan-Wigner transformation<cit.>. We assume the first four lowest energy orbitals are occupied by electrons and first two highest energy orbitals are not occupied all time. We also use the spin symmetry in <cit.> to reduce another two qubits. With the reduction of  the number of qubits, finally we have 6 qubits Hamiltonian <cit.>. All calculations of integrals in second quantization and transformations of electronic structure are done by OpenFermion<cit.> and Psi4<cit.>.

Gradient estimation.
The two functions ϕ(x) and s(x) are both real function. Thus, the gradient for parameter p_k can be estimated as  2(⟨ E_locD_p_k⟩ -⟨ E_loc⟩⟨ D_p_k⟩) where E_loc(x)=⟨ x|H|ϕ⟩/ϕ(x)s(x) is so called local energy, D_p_k(x) = ∂_p_k(ϕ(x)s(x))/ϕ(x)s(x). ⟨...⟩ represents the expectation value of joint distribution determined by ϕ(x) and s(x) (details in  Supplementary Note 1).

Implementation Details.
In our simulation we choose small constant learning rate 0.01 to avoid trapping in local minimum. All parameter are initialized as a random number between (-0.02, 0.02). The range of initial random parameter is to avoid gradient vanishing of tanh. For each calculation we just need 1 reusing ancilla qubit all the time.  Thus, in the simulation, the number of required qubits is m+n+1.  All calculations do not consider the noise and system error (details in  Supplementary Note 5).

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

unsrt





§ ACKNOWLEDGEMENT

We would thank to Dr. Ross Hoehn, Dr. Zixuan Hu and Teng Bian for critical reading and useful discussions. S.K and R.X are grateful for the support from Integrated Data Science Initiative Grants, Purdue University.



§ AUTHOR CONTRIBUTION

S.K designed the research. R.X performed the calculations. Both discussed the results and wrote the paper. 



§ ADDITIONAL INFORMATION

Supplementary Information. Supplementary Materials are available.

Competing Interests. The authors declare no competing interest.






§ SUPPLEMENTARY NOTE 1



 §.§ Derivation of the gradient

For an electronic structure Hamiltonian prepared by second quantization and Jordan-Wigner transformation<cit.>, H, and a trial wave function, |ϕ⟩= ∑_xϕ(x)s(x)|x⟩, the expectation value can be written as<cit.>:


    ⟨ H⟩=⟨ϕ|H|ϕ⟩/⟨ϕ|ϕ⟩=∑_x,x'ϕ(x)s(x)⟨ x|H|x'⟩ϕ(x')s(x')/∑_x|ϕ(x)s(x)|^2



x is a combination of  {1,2...n} and |x⟩=|12...n⟩.

If we set Φ(x) = ϕ(x)s(x), because ϕ(x) and s(x) are all real value functions, then the gradient can be calculated as<cit.>:


    ∂_p_k⟨ H⟩   =∑_x(∂_p_kΦ(x))⟨ x|H|ϕ⟩+∑_x⟨ϕ|H|x⟩(∂_p_kΦ(x))/∑_x|Φ(x)|^2
       -∑_xΦ(x)⟨ x|H|ϕ⟩/∑_x|Φ(x)|^2∑_x((∂_p_kΦ(x))Φ(x)+Φ(x)∂_p_kΦ(x))/∑_x|Φ(x)|^2


If we set E_loc(x)=⟨ x|H|ϕ⟩/Φ(x) and D_p_k(x)=∂_p_kΦ(x)/Φ(x), the gradient can be written as<cit.>:


    ∂_p_k⟨ H⟩   =∑_xD_p_k(x)E_loc(x)|Φ(x)|^2+∑_xE_loc(x)D_p_k(x)|Φ(x)|^2/∑_x|Φ(x)|^2
       -∑_x|Φ(x)|^2E_loc(x)/∑_x|Φ(x)|^2∑_x(D_p_k(x)+D_p_k(x))|Φ(x)|^2/∑_x|Φ(x)|^2
       =2⟨ E_locD_p_k⟩ - 2⟨ E_loc⟩⟨ D_p_k⟩


where ⟨...⟩ represent the expectation value of distribution determined by Φ(x). ⟨ x|H|ϕ⟩ = ⟨ϕ| H|x⟩ for that  H is a real symmetric matrix due to Jordan-Wigner transformation.

p_k is the parameters a_i, b_j, w_ij, d_i, c for k_th iterations. Thus we have<cit.>:

    D_a_i(x)=1/2i -1/2⟨i⟩_RBM,
       D_b_j(x)=1/2tanh(θ_j)-1/2⟨ h_j⟩_RBM,
       D_w_ij(x)=1/2tanh(θ_j)i-1/2⟨ih_j⟩_RBM,
       D_c(x)=1/s(x)-s(x),
       D_d_i(x)=i(1/s(x)-s(x)),


where θ_j=∑_iw_ijσ^z_i+b_j. ⟨...⟩_RBM represents the distribution determined solely by RBM. We do not need to calculate the second term of D_a_i, D_b_i and D_w_ij for that they will be cancelled when calculating the gradient ∂_p_k⟨ H⟩. We use the gradient decent method to optimize our RBM, yielding the global minimum corresponding to the ground energy.


    p_k+1=p_k-α_k∂_p_k⟨ H⟩


Where α_k is the learning rate for k_th iteration, controlling the convergence rate. We can continue iterating until we reach the maximum number of iterations. The gradient is estimated by the distribution calculated by sampling.




§ SUPPLEMENTARY NOTE 2




 §.§ Sequential applications of controlled-rotation algorithm


The probability for each combination y ={σ^z, h} can be written as:


    P(y) = e^∑_ia_ii+∑_jb_jh_j+∑_i,jw_ijih_j/∑_y'e^∑_ia_ii+∑_jb_jh_j'+∑_i,jw_ijih_j'



However, we do not directly calculate the P(y)  but we do some modification on P(y) to increase the successful probability of our algorithm. We calculate Q(y) =e^∑_ia_ii/k+∑_jb_jh_j/k+∑_i,jw_ijih_j/k/∑_y'e^∑_ia_ii/k+∑_jb_jh_j'/k+∑_i,jw_ijih_j'/k where k is a large number to increase the successful probability of our measurements.

First we use R_y gate to achieve a superposition of all possible σ^z and h. The system qubits and the ancilla qubit are initialized at state |0⟩.


    ⊗_iR_y(2arcsin(√(e^a_i/k/e^a_i/k+e^-a_i/k)))|0_i⟩⊗_jR_y(2arcsin(√(e^b_j/k/e^b_j/k+e^-b_j/k)))|0_j⟩|0⟩=∑_y√(O(y))|y⟩|0⟩


where O(y) = e^∑_ia_ii/k+∑_jb_jh_j/k/∑_y'e^∑_ia_ii/k+∑_jb_jh_j'/k and |y⟩ =|1..nh_1...h_m⟩.

The next step is to calculate each term of e^∑_i,jw_ijih_j, which is achieved by controlled rotations gates. The idea is, for each time controlled rotation, we calculate two angles θ_ij,1 = 2arcsin(√(e^w_ij/ke^-|w_ij|/k)) and θ_ij,2 = 2arcsin(√(e^-w_ij/ke^-|w_ij|/k)). We use controlled-rotation CR_y(θ_ij,1) and CR_y(θ_ij,2) which are controlled by combination of i, h_j as working qubits  to do rotation on the ancilla qubit. The controlled rotation is to check the working qubits and then do the corresponding rotation θ_ij,1 or θ_ij,2.

All controlled rotation gates can be expressed as below:


    CR_w_ij,1=C_i,h_j⊗ R_y(2arcsin(√(e^w_ij/ke^-|w_ij|/k)))+(D_i,h_j+E_i,h_j+F_i,h_j)⊗ I 
       CR_w_ij,2=D_i,h_j⊗ R_y(2arcsin(√(e^-w_ij/ke^-|w_ij|/k)))+(C_i,h_j+E_i,h_j+F_i,h_j)⊗ I
       CR_w_ij,3=E_i,h_j⊗ R_y(2arcsin(√(e^-w_ij/ke^-|w_ij|/k)))+(C_i,h_j+D_i,h_j+F_i,h_j)⊗ I
       CR_w_ij,4=F_i,h_j⊗ R_y(2arcsin(√(e^w_ij/ke^-|w_ij|/k)))+(C_i,h_j+D_i,h_j+E_i,h_j)⊗ I


where C_i,h_j = B_i⊗ B_h_j, D_i,h_j = A_i⊗ B_h_j, E_i,h_j = B_i⊗ A_h_j, F_i,h_j = A_i⊗ A_h_j and  


    A=
    [ 1 0; 0 0;   ] B=
    [ 0 0; 0 1;   ]
 

Between the calculation of two w_ij, we need to do a measurement on the ancilla qubit to make sure the state of system qubits collapse to the wanted state. Measuring ancilla qubit in |1⟩ means the state of system qubits collapse to the wanted state as we initialize the ancilla qubit in |0⟩. 

We then do the measurement, if and only if the ancilla qubit is in |1⟩ we continue with a new ancilla qubit initialized in |0⟩, otherwise we start from beginning. The probability of success is very large since we choose k as a large number.

After we finish all measurements, the distribution should be Q(y). We just measure the first n+m qubits of the system register to obtain the probability distribution. After we get the distribution, we calculate all probabilities to the power of k and normalize to get the Gibbs distribution. 




§ SUPPLEMENTARY NOTE 3



 §.§ Lower bound of successful sampling

The successful probability P can be written as:


    P = ∑_σ^z,he^1/k(∑_ia_ii+∑_jb_jh_j+∑_i,jw_ijih_j)/∑_σ^z,he^1/k(∑_ia_ii+∑_jb_jh_j)e^1/k(∑_i,j|w_ij|)≥e^-1/k(∑_i,j|w_ij|)/e^1/k(∑_i,j|w_ij|)=1/e^1/k(∑_i,j2|w_ij|)


If we choose k = O(∑_i,j|w_ij|), we have P≥1/e^O(1) which means the lower bound of probability of successful sampling is a constant. In the simulation, we choose k=max(1/2∑_i,j|w_ij|,1) because larger k  introduces larger sampling errors. This particular choice of k gives us lower bound of success as e^-4. But in numerical simulation, the probability is much larger than e^-4, see Supplementary Note 5.




§ SUPPLEMENTARY NOTE 4



 §.§ Complexity and error of the algorithm

For a C^2(U) conditioned by |11⟩, it can be decomposed as the below<cit.>:



where V^2 = U. In our algorithm, U=R_y(θ), thus we can choose V=R_y(θ/2) to achieve the decomposition. C^2(U) conditioned  by |00⟩, |10⟩ or |01⟩ can be achieved by adding X gates on controlling qubits. For each w_ij, we have 4 C^2(U) which means the gates complexity scales to O(mn) and the number of qubits for our algorithms scales to O(mn), which can be reduced to O(m+n) if considering qubit reuse. Because the lower bound of probability of successful sampling is constant, if the number of successful sampling is N_s, the complexity for each iteration is O(N_smn). The only error comes from the error of sampling if not considering noise in the quantum computer.



§ SUPPLEMENTARY NOTE 5



 §.§ Implementation details for H_2, LiH and H_2O


Here we present the probabilities of successful sampling when calculating H_2, LiH and H_2O for bond length equals to 1.75 Angstrom.




Here we present the the changes of energy during optimization when calculating H_2, LiH and H_2O for bond length equals to 1.75 Angstrom..





The distribution we want to sampling for the quantum algorithm is:


    Q(y) =e^∑_ia_ii/k+∑_jb_jh_j/k+∑_i,jw_ijih_j/k/∑_y'e^∑_ia_ii/k+∑_jb_jh_j'/k+∑_i,jw_ijih_j'/k


In our controlled-rotation algorithm, we use a k as regulation to increase the probability of success as the proof in the Supplementary Note 3, the lower bound of probability of success would become :


    1/e^1/k(∑_i,j2|w_ij|)


Thus, if no regulation (k=1), the probability of success would become 1/e^∑_i,j2|w_ij| which means we need exponential number of measurements to get enough successful sampling, making no speedup in quantum algorithm.


If we add a regulation of k, in simulation we use k=1/2∑_i,j|w_ij|, the probability becomes 1/e^4 which needs constant number of measurements to get enough successful sampling. 


After we get the distribution, we need to calculate all distribution to the power of k and normalize to get the wanted distribution. k is a large number, which is around 5 at final for H_2, LiH and H_2O in our simulation. To decrease the errors in calculating power of k, we have to increase the number of sampling for our quantum algorithm when k is large, which requires large number of sampling and may not be efficient. In the Supplementary Figure 3, we can see that at the final procedure of optimization, the fluctuation is very large due to large k, which can be decreased by increasing the number of sampling. Because we investigated small molecule system H_2, LiH and H_2O, k is not very large and the quantum algorithm is efficient. For large k, our quantum algorithm may require large sampling.


