Experimental realization of quantum algorithm for solving linear systems of equations
Abstract
Quantum computers have the potential of solving certain problems exponentially faster than classical computers. Recently, Harrow, Hassidim and Lloyd proposed a quantum algorithm for solving linear systems of equations: given an matrix and a vector , find the vector that satisfies . It has been shown that using the algorithm one could obtain the solution encoded in a quantum state using quantum operations, while classical algorithms require at least steps. If one is not interested in the solution itself but certain statistical feature of the solution ( is some quantum mechanical operator), the quantum algorithm will be able to achieve exponential speedup over the best classical algorithm as grows. Here we report a proofofconcept experimental demonstration of the quantum algorithm using a 4qubit nuclear magnetic resonance (NMR) quantum information processor. For all the three sets of experiments with different choices of , we obtain the solutions with over fidelity. This experiment is a first implementation of the algorithm. Because solving linear systems is a common problem in nearly all fields of science and engineering, we will also discuss the implication of our results on the potential of using quantum computers for solving practical linear systems.
pacs:
03.65.Ta, 03.65.Xp, 76.60.kThe proposition of quantum computer dates back to 1980s feynman86 (), but it was not until the late 80’s and early 90’s that quantum computers are shown to be more powerful than classical computers on various specialized problems DJ92 (); Grover96 (); Lloyd96 (); Shor94 (). For example, the DeutschJozsa algorithm DJ92 (), Shor’s quantum algorithm for factoring integers Shor94 (), Grover’s quantum search algorithm Grover96 () and algorithms for Hamiltonian simulation of quantum systems Lloyd96 () have been found to require significantly less computational steps than their classical counterparts and thus render many classically intractable problems realistically solvable with a quantum computer. There has been experimental demonstrations of important quantum algorithms such as Shor’s algorithm shor (), Grover search algorithm grover (), optimization problems chuang_opt (), quantum simulation of molecular systems white_quan_chem (); walther_quan_chem (); baugh_quan_sim (); cory_quan_sim (), all using smallscale quantum information processors. In this work we implement an algorithm for solving linear systems originalpaper () using NMR.
Linear systems of equations play an important role in nearly all fields of science and engineering. In quantum reactive scattering, the Kohn variational calculation involves the inversion of the augmented stiffness matrix levine_scattering (); wyatt_scattering (), which is equivalent to solving a linear system in certain occasions. In chemistry, linear equations arise commonly in problems such as electrostatic calculation in density functional theory, where the discretized Poisson equation takes a linear form parr_DFT (). Recently the Finite Element method starts to be adopted for solving electronic structure problems in quantum chemistry martinez_FEA (), where Schrödinger’s equation is recast in form of a linear equation. Also solving linear systems of equations often play a role as intermediate step in many algorithms such as quantum algorithm for data fitting datafitting ().
Here we consider the particular type of linear system where we are given an sparse Hermitian matrix with condition number and unit vector and we are interested not in the solution itself but certain feature of that can be written in form of ( is some linear quantum mechanical operator). As the size of matrix grows, the size of the data sets which define the equations increase rapidly over time. For classical algorithms such as the Conjugate Gradient Method shewchuk94 (), it takes about total runtime of to get the solution when is positive definite or when is not. For all classical algorithms, it is shown that the linear lower bound in runtime scaling cannot be broken even we are only interested in rather than itself originalpaper ().
The quantum algorithm proposed by Harrow, Hassidim and Lloyd originalpaper () for solving linear systems of equations is shown to improve the runtime scaling to the regime. The algorithm starts with a quantum state and a few ancilla qubits in state . Here the state vector of in the computational basis represents the unit vector . Applying the wellknown phase estimation subroutine nielsen00 () on with as the unitary operation, we obtain the final state of the tworegister system which is approximately up to a normalization constant (When the eigenvalues of can be exactly encoded using the ancilla bits, which is the case of our experiment as we will see later in the paper, the final state of the phase estimation is exactly proportional to ). Here is the eigenbasis of and . When the phase estimation is completed, an ancilla bit is added to the system and a controlled rotation is performed on it with the register as the control register. The state of the system then becomes where is a normalization constant. By inverting all the previous quantum operations except for the controlled rotation on the ancilla bit, we uncompute the register back to the state . Conditioning on measuring in the ancilla qubit, The algorithm probabilistically outputs a state in the register that is initially in the state . The state vector of in the computational basis is proportional to the solution of the linear system . The total number of quantum operations needed for the algorithm is . This implies an exponential speedup over the Conjugate Gradient Method, which is the best classical algorithm for solving the problem.
Let us consider a system to find such that where
(1) 
.
where is normalized. The spectral decomposition has , and , . The quantum algorithm to solve this problem can be summarized as the following six steps originalpaper () (see Fig. 8).

Input the vector as a quantum state stored in a quantum register (termed register ) and prepare the register in the state . Here and represents the computational bases of registers and that encode the value of with the qubit and qubits respectively.

Apply the conditional Hamiltonian evolution on both registers and up to error . Here ‘’ in the notation marks the control register and .

Apply quantum Fourier transform to the register . Denote the Fourier basis as . In general, at this stage in the superposition state of both registers, the amplitudes of the basis states are concentrated on values that approximately satisfy , where is the th eigenvalue of the matrix . However, is defined as in Eq. (1), the eigenvalues , are exactly powers of 2. Therefore if we further let , , the phase estimation subroutine will produce exactly .

Add an ancilla qubit and apply conditional rotation on it, controlled by the register . This rotation transforms the qubit to , where is a nonzero normalization constant . This is a key step of the algorithm and it involves finding the reciprocal of the eigenvalue quantum mechanically and rotate the qubit with an angle such that . To achieve both, we first use a SWAP gate (Fig. 8) to accomplish the transformation
(2)
in register . Then we use the register with states as the control register to apply the controlled rotation on the ancilla qubit. As shown in Fig. 8, the controlled gates applies rotation with where is a parameter. In this experiment we let and use the approximation . With this approximation, after the controlled rotation gates, the state of the system is close to
(3) and is 0.736 according to the calculation based on Ref.jmp ().

Uncompute register and by inverting the operation in the steps 13.

Measure the ancilla bit. If it returns 1, the register of the system is in the state up to a normalization factor, which is equal to the solution of the linear system .
The experiment was carried out on a Bruker AV400 spectrometer () at 303.0K. We chose iodotrifiuoroethylene dissolved in chloroform, where a nucleus and three nuclei consist of a fourqubit quantum system. We label as the first qubit, , and as the second, third and forth qubit. Fig. 2 shows the measured properties of this fourqubit quantum system Luozhihuang (). We control the fluctuation of the temperature to be within 0.1K such that the effect of the chemical shift variations is suppressed. This system is first prepared into a pseudopure state (PPS) with representing the unity operator and the polarization, using the lineselectivetransition method PPS (); Luozhihuang (). It needs two gradient ascent pulse engineering (GRAPE) pulses grape () and two gradient field pulses. We apply quantum state tomography tomography () to reconstruct its experimental density matrix of the prepared PPS and the experimental fidelity is around 98.7. The state fidelity is calculated by , where and represent experimentally measured density matrices and the theoretical expectation, respectively. Then, we perform a rotation operation (i.e.,a rotation along the axis with an angle to the ) quan_adia (); quan_sim (); cloning (), to obtain the initial state , where the normalized state with the state vector .
Now, we implement the quantum circuit of the algorithm shown in Fig. 8 on the prepared input state . The quantum circuit is realized as a whole block using a shaped radiofrequency (r.f.) pulse that is optimized by the gradient ascent pulse engineering (GRAPE) algorithm grape (); nmr_qc (); factor143 (). The GRAPE pulse is characterized by 1500 segments, the pulse duration of 22.5 ms, and is robust to r.f. inhomogeneities, with a theoretical fidelity 0.995. Finally, we measure the final state to obtain the solution of the linear equations. The desired final state is with , the solution is encoded into the state of qubit 3 in the subspace labeled by . Here we consider three different by preparing three input states with different . We perform a partial state tomography to get the information of and (see Method). The experimental C spectra are shown in the left side of Fig. 4. are the ratio of the probabilities of projecting the final state to the states and , and the relative phase between and can be obtained by the coherence term or of qubit 3. Since A and both are real, is real and the relative phase is either 0 or by the sign of or . There may exist a global phase which can not be determined in the solution , but this global phase is often not important when one estimates the expectation value of some operator associated with . We list all results of the three experiments in Fig. 4. The experimental errors are about which is measured by the formula . Furthermore, we perform the complete quantum state tomography tomography () for the final state in the Hilbert space spanned by these four qubits which is described in supplementary materials. The real parts of the reconstructed density matrices in the subspace labelled by are shown in the right side of Fig. 3, with the experimental fidelities are all above .
The errors of the experiment mainly come from the inhomogeneity of magnetic fields, the imperfection of the GRAPE pulses and the chemical shift variations. The imperfection of the GRAPE pulses produces about error to the final state which can cause about error to the intensities of some peaks. The experimental time is about 50 ms, about of the coherence time ( in Fig. 2). Hence, the decay of the signals due to the limitation of coherence time is an important source of errors. The chemical shift of fluorine changes by about Hz when the temperature varies by 1K material (), so we control the fluctuation of the temperature below 0.1K to reduce the error to about . We obtain the intensity of the peaks by fitting the spectra using the Lorentzian curves to reduce the error by the noise. The total derivation between the experimental final state and the theoretical algorithm outcome is about 3.
The theoretical errors of algorithm concentrate on the phase estimation in the second step and the controlrotation operation in the forth step. In our cases, for the reason that the eigenvalues of which we chose satisfy , the phase estimation is perfect and produces no error in theory originalpaper (). Most theoretical errors come from the controlrotation operation which depends on the positive integer parameter . This error essentially depends on the quality of the approximation , where is the rotation angle which is integral multiple of acted on the ancillary bit. The theoretical error in this step decreases as increases when the intensities of the counterpart peaks of the final states are inversely proportional to jmp (). Therefore, we take a balanced choice in the experiments, which causes that the theoretical error is contributing to and the intensities of the counterpart peaks are about of the intensity of PPS. As presented in Fig. 4, we know that the theoretical error caused by this approximation is one of the major errors.
In conclusion, we experimentally demonstrate for the first time the quantum algorithm for solving linear systems of equations in a 4qubit NMR system. We acquire the solutions with errors about for all the three group experiments, which indicates fine experimental accuracy and the validity of the algorithm. This quantum experiment can be a good evidence that the quantum computer can help solving common and important problems. For solving linear equation plays important role in many classical algorithm and classical computer, this experiment is the first step to create and realize more quantum algorithm. It gives us the hope to broaden the application of quantum computer.
The application of the quantum linear equation solver could be extended to a range of applications in many fields. For example, the ability of using the quantum algorithm to solve Poisson equation poisson () would allow quantum chemists to speed up electrostatic calculation in density function theory. Furthermore, since the quantum algorithm could be used for efficiently solving linear systems of differential equations berry (), quantum computer might prove useful for solving the differential equations systems that arise in technical application.
Methods
Hamiltonian simulation. In this experiment, since the unitary evolution is a singlequbit gate, when realizing the operations we are not directly implementing any Hamiltonian simulation algorithms in the literature (such as BACS07 ()) which are responsible for the exponential speedup of the quantum algorithm originalpaper (). Instead, we decompose the gates into elementary quantum gates using Group Leader optimization algorithm groupleader ().
Finding the reciprocal of the eigenvalue . The inversion technique which we use in Eq. (2) to find the reciprocals of is applicable only to some particular type of matrices such as the one defined in Eq. (1), where a SWAP gate combined with a proper choice of can give us a sufficiently good approximation of the amplitude in the final state Eq. (3) of the ancilla qubit. For the general case where the qubit register holds a superposition of states that are approximately , one can resort to techniques such as Newton iteration for finding the reciprocals of . The quantum algorithm implementing Newton iteration has also been proven efficient poisson (). Furthermore, in step 4 we also used the approximation . In general this approximation does not hold but one could use bisection method poisson () to efficiently find .
Partial Tomography. The natural abundance of the sample in which just one carbon is is about . To distinguish those molecules against the large background, we read out all three qubits via the channel, by applying SWAP gates and reading out the qubit. The solution can be obtained from the subspace labeled by , i.e., the information of and . This can be achieved by 5 readout pulses (,, , , ). Here represents the unity operator. and represent, respectively, a rotation operation along and axis. denotes a SWAP operation between th and th qubits nielsen00 (). The first four readout pulses are to get the probabilities of and , while the last readout pulse is to get the relative phase of .
Acknowledgements.
The authors would like to thank the support of National Key Basic Research Program of China (Grant No. 2013CB921800) and the Strategic Priority Research Program (B) of the Chinese Academy of Sciences (Grant No. XDB01030400)¡£ This work is also supported by the National Natural Science Foundation of China (Grant no. 91021005, 21073171, 11275183, 11104262, 10975124) and the Ministry of Education of China (Grant no. 20113402110044 and the Scientific Research Foundation for the Returned Overseas Chinese Scholars). Prof. Sabre Kais thanks the NSF Center for Quantum Information and Computation for Chemistry, Award number CHE1037992.author contributions
J. Pan, X. Yao, Z. Li, C. Ju, X. Peng and J. Du designed and implemented the experiment for realizing the quantum circuit. Y. Cao and S. Kais devised the quantum circuit based on the original paper by Harrow et al. All authors discussed the results and cowrote the manuscript.
Competing financial interests
The authors declare no competing financial interests.
References
 (1) Feynman, R. P. Quantum mechanical computers. Foundations of Physics, 16, (1986).
 (2) Deutsch, D. & Jozsa, R. Rapid solution of problems by quantum computation. Proc. R. Soc. Lond. A, 439, 553 (1992).
 (3) Shor, P. W. Algorithms for Quantum Computation: Discrete Logarithm and Factoring. (IEEE Computer Society Press, New York, 1994).
 (4) Grover, L. A Fast Quantum Mechanical Algorithm for Database Search. (ACM Press, New York, 1996).
 (5) Lloyd, S. Universal quantum simulators. Science, 273, 10731078 (1996).
 (6) Lieven M. K. et al. Experimental realization of Shor’s quantum factoring algorithm using nuclear magnetic resonance. Nature 414, 883887 ( 2001).
 (7) Dasa, R., Mahesha, T. S. & Kumar, A. Experimental implementation of Grover search algorithm using efficient quantum state tomography. Chemical Physics Letters 369, 815 (2003).
 (8) Steffen, M., Dam, W. V., Hogg, T., Breyta, G. & Chuang, I. Experimental implementation of an adiabatic quantum optimization algorithm. Physical Review Letters 90, 067903 (2003).
 (9) Lanyon, B. P. et al. Towards quantum chemistry on a quantum computer. Nature Chemistry, 2, 106111 (2010).
 (10) AspuruGuzik, A. & Walther, P. Photonic quantum simulators. Nature Physics, 8, 285 (2012).
 (11) Zhang, J., Yung, M., Laflamme, R., AspuruGuzik, A. & Baugh. J. Digital quantum simulation of the statistical mechanics of a frustrated magnet. Nature Communication 3, 880 (2012).
 (12) Somaroo, S., Tseng, C. H., Havel, T. F., Laflamme, R. & Cory, D. G. Quantum simulations on a quantum computer. Phys. Rev. Lett. 82, 53815384 (1999).
 (13) Harrow, A. W., Hassidim, A. & Lloyd, S. Quantum algorithm for linear systems of equations. Phys. Rev. Lett. 15 150502 (2009).
 (14) Levine, R. D. Quantum Mechanics of Molecular Rate Processes. (Dover Publications, 2012).
 (15) Manolopoulos, D. E., D’Mello, M., & Wyatt, R. E. Quantum reactive scattering via the log derivative version of the Kohn variational principle: General theory for bimolecular chemical reactions. J. Chem. Phys. 91, 60966102 (1989).
 (16) Parr, R. G., & Yang, W. DensityFunctional Theory of Atoms and Molecules. (Oxford Science Publications, 1989).
 (17) Alizagedan, R., Hsia, K. J., & Martinez, T. J. A divide and conquer real space finiteelement HartreeFock method. The Journal of Chemical Physics 132, 034101 (2010).
 (18) Alizagedan, R., Hsia, K. J., & Martinez, T. J. A divide and conquer real space finiteelement HartreeFock method. Phys. Rev. Lett 109, 050505 (2012).
 (19) Wilson. D. J. Solution of systems of linear equations in analytical chemistry. Anal. Chem. 30, 15781579, (1958).
 (20) Shewchuk, J. R. An Introduction to the Conjugate Gradient Method Without the Agonizing Pain. (Carnegie Mellon University Pittsburgh, PA, USA, 1994).
 (21) Nielsen, M. A. & Chuang, I. L. Quantum Computation and Quantum Information. (Cambridge University Press, 2000).
 (22) Cao, Y., Daskin, A., Frankel, S. & Kais, S. Quantum circuit design for solving linear systems of equations. Molecular Physics. Special Issue: Scaling Mount Impossible: A Festschrift for Dudley Herschbach 110, 16751680 (2012).
 (23) Peng, X. et al. Quantum simulation of topologicalordertransition in the WenPlaquette Model. In preparation.
 (24) Peng, X. et al. Preparation of pseudopure states by lineselective pulses in nuclear magnetic resonance. Chem. Phys. Lett. 340, 509516 (2001)
 (25) Peng, X. et al. Quantum adiabatic algorithm for factorization and its experimental implementation. Phys. Rev. Lett. 101, 220405 (2008).
 (26) Peng, X., Zhang, J., Du, J. & Suter, D. Quantum simulation of a system with competing two and threebody interactions. Physical Review Letters 103, 140501 (2009).
 (27) Chen, H. et al. NMR experimental demonstration of probabilistic quantum cloning. Phys. Rev. Lett. 106, 180404 (2011).
 (28) Lu, D. et al. Simulation of chemical isomerization reaction dynamics on a NMR quantum simulator. Physical Review Letters 107, 020501 (2011).
 (29) Xu, N., Zhu, J., Lu, D., Zhou, X., Peng, X. & Du, J. Quantum factorization of 143 on a dipolarcoupling nuclear magnetic resonance system. Phys. Rev. Lett. 108, 130501 (2012).
 (30) Khaneja, N. et al. Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms. J. Magn. Reson. 172, 296 (2005).
 (31) Lee, J. S. The quantum state tomography on an NMR system. Phys. Lett. A 305, 349 (2002).
 (32) Supplementary materials.
 (33) Cao, Y., Papageorgiou, A., Petras, I., Traub, J. & Kais, S. Quantum algorithm and circuit design solving the Poisson equation. arXiv:1207.2485v2 (2012).
 (34) Berry, D. W. Quantum algorithms for solving linear differential equations. arXiv:1010.2745 (2010).
 (35) Berry, D. W., Ahokas, G., Cleve, R. & Sanders, B. C. Efficient quantum algorithm for simulating sparse Hamiltonians. arXiv: quantph/0508139 (2006).
 (36) Daskin, A. & Kais, S. Decomposition of unitary matrices for finding quantum circuits. J. Chem. Phys. 134, 144112 (2011).
Supplementary Materials
Appendix A. Pulse sequence of the experiment
Figure 5 is the pulse sequence of the experiment. The first two GRAPE pulses and two gradient field pulses were used to prepare a pseudopure state (PPS) . The third GRAPE pulse is to perform the quantum circuit. The time of the readout pulses differs from 0.4ms to 25ms. The implementation time of circuit and the readout pulse is about 0.05ms in total.
Appendix B. Change of chemical shift caused by temperature fluctuation
The chemical shifts are given with respect to reference frequencies of 376.47 MHz (fluorines) and 100.64 MHz (carbons) at 303.0K. The chemical shifts of fluorines change by about 3 Hz when the temperature varies by 1K while the chemical shift of carbon is almost the same. The expression for the relationship between the chemical shifts of fluorines and the temperature is obtained by fitting the experimental data:
where , , represent the the chemical shifts of , , respectively and represents the temperature. The units are Hz and K. We controlled the temperature in in the experiment.
Appendix C. State tomography.
The state tomography for the whole 4bit quantum state needs 44 readout pulses: , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , . ( represents the unity operator. represents a rotation operation with the angle along axis while along axis. means a swap operation between the th and th qubits.). As discussed in the methods, we can pick up the accurate and complete information of using only 5 readout pulses: , , , , to get the solution . Combining the first four readout pulses , we can obtain the diagonal elements of the density matrix. Using the fifth readout pulse, we can get and . The state tomography results using 44 readout pulses for PPS and the experimental final states are showed in Figure 6 and Figure 7.
Appendix D. Relationship between peaks and elements of density matrix
As shown in Figure 8, there are eight peaks for carbon. Four of them are almost zero, and the intensity of the four significant peaks quantify the probabilities of the states , , , from left to right respectively. In fact, defining the probabilities of the states , , … , are , , …, , the four peaks which can be seen obviously are proportional to , , , from left to right. Since for all in experiment, the intensity of the four peaks are approximate proportional to the probabilities of the states , , , respectively.