Quantum simulation and circuit design for solving multidimensional Poisson equations
QQuantum simulation and circuit design for solvingmultidimensional Poisson equations
Michael Holzmann and Harald K¨ostler Friedrich-Alexander-Universit¨at Erlangen-N¨urnberg, Germany
Abstract —Many methods solve Poisson equations by using gridtechniques which discretize the problem in each dimension. Mostof these algorithms are subject to the curse of dimensionality,so that they need exponential runtime. In the paper ”Quantumalgorithm and circuit design solving the Poisson equation” aquantum algorithm is shown running in polylog time to producea quantum state representing the solution of the Poisson equation.In this paper a quantum simulation of an extended circuit designbased on this algorithm is made on a classical computer. Ourpurpose is to test an efficient circuit design which can breakthe curse of dimensionality on a quantum computer. Due to theexponential rise of the Hilbert space this design is optimizedon a small number of qubits. We use Microsoft’s QuantumDevelopment Kit and its simulator of an ideal quantum computerto validate the correctness of this algorithm.
I. I
NTRODUCTION
The goal of this work is to implement a quantum algorithmsolving the d-dimensional Poisson equation with Dirichletboundary conditions: − ∆ u ( x ) = f ( x ) , x in Ω u ( x ) = g ( x ) , x on δ Ω Ω = (0 , d One way to solve this problem is to discretize Ω in M+1 gridpoints in each dimension. M is an exponent of base 2 in thiswork. The solution u(x) is a vector of ( M − d entries. Tocalculate a problem of one dimension a linear equation systemhas to be solved: h A d =1 u ... u m − = h − − ... ... ... ... − u ... u m − = f + h u ... f m + h u m = (cid:126)b For multidimensional problems the linear equation uses matrix A d , which can described as [1]: A d = d (cid:88) i =1 i − (cid:79) I ⊗ A d =1 ⊗ d − i (cid:79) I (1)Matrix A d has the dimension of ( M − d × ( M − d .The best numerical algorithms for solving this problem runpolynomially to matrix size [3], so the runtime increasesexponentially with the dimension of the problem. In this papera quantum algorithm is used to produce a quantum staterepresenting the normalized solution of the problem. Since thistechnique runs in polylog time the curse of dimensionality canbe broken. II. Q UANTUM ALGORITHM AND CIRCUIT DESIGN
The quantum algorithm for solving linear equations is basedon the HHL09 algorithm [4]. If successfull this algorithmreplaces the input (cid:126)b with the normalized solution (cid:126)u of u ( x ) .The algorithm used in this work follows several steps: • Produce a quantum state (cid:80) j β j | j (cid:105) in RegC. β j ∈ (cid:126)b • Use Phase Estimation Algorithm (PEA) on RegC andRegB. Register B consists of n qubits. In the processof PEA several Hamiltonian simulations of U = e h A d t with t = πi n k k = 0 , ..., n − are applied toRegC. PEA entangles the eigenvalues λ j of A d in RegBwith the eigenstates (cid:126)u j in RegC: (cid:80) j b j |≈ λ j (cid:105) | (cid:126)u j (cid:105) with b j = (cid:104) (cid:126)b | (cid:126)u j (cid:105) • Calculate the reciprocal of the eigenvalues in RegA. Thesystem has now the state: (cid:80) j b j | k j (cid:105) | k j (cid:105) | (cid:126)u j (cid:105) with k j ≈ λ j • Apply a controlled rotation on an ancilla qubit to producefollowing system: (cid:80) j b j | k j (cid:105) | k j (cid:105) | (cid:126)u j (cid:105) ( (cid:113) − α k j ) | (cid:105) + αk j | (cid:105) ) with k j ≈ λ j and amplitude factor α • Uncompute RegA and RegB • Measure the ancilla qubit. If the measurement of the qubitresults in state | (cid:105) , the algorithm successfully transformsRegC into the normalized solution (cid:126)u = C (cid:80) j b j αk j | (cid:126)u j (cid:105) = CA − (cid:126)b with C as a normalization constant. Otherwise thealgorithm has to be restarted. A. Encoding of register C
For the first step the right side of the linear equationhas to be encoded in register C. If β j can be calculatedefficiently, a quantum oracle can be applied to a zero stateregister to produce (cid:80) j β j | j (cid:105) in polylog time [5]. Since log ( M ) qubits can encode M states, the zero state is notused and has amplitude zero in the one-dimensional case. Formultidimensional problems d log ( M ) qubits are needed. Tosolve these problems Hamiltonian simulations of A d =1 eachacting on log ( M ) qubits are used. Because every Hamiltoniansimulation does not allow a zero state, every log ( M ) qubitblock of register C must not be in zero state. This yieldsM d − ( M − d invalid states. i.e for d = 1 and M = 4 β is the amplitude of state | (cid:105) , while for the two-dimensionalcase it is encoded in the amplitude of state | (cid:105) . Since registerC is used for the input and output of the algorithm, the sameencoding applies for the normalized solution. a r X i v : . [ c s . ET ] J un . Phase Estimation Algorithm of A d (PEA) One subroutine of the algorithm is the one-dimensionalHamiltonian simulation U = e h A t . Using the spectraltheorem [1]: e h A t = Se t ∆ S , ∆ = λ λ . . .
00 0 0 λ M − (2)S desribes the discrete sine transform acting on log ( M ) + 1 qubits, which can be implemented with the quantum fouriertransform (QFT), two transformation circuits and one ancillaqubit [6]. For the Hamiltonian simulation of a 1-sparse matrix e t ∆ two oracle calls are needed. The first one calculates theeigenvalues λ j of the one-dimensional Poisson matrix, whichare used to produce the right rotation in the secondary oraclecall [7]. Since A is a toeplitz matrix, λ j are [8]: λ j = 4 M sin ( jπ M ) j = 1 , ..., M − (3)The sine function can be approximated with a quantum oracleusing repeated squaring [1]. Even though it can be wellcalculated with respect to M, quantum arithmetics use manyancilla qubits to run. Since the referenced quantum algorithmof the repeated squaring process results in multiple ancillaregisters, the Hilbert space rises exponentially. In the paper“Quantum Fast Poisson Solver: the algorithm and modularcircuit design” [2] a module is shown calculating theseeigenvalues with O ( n M ) qubits where n describes the sizeof register B. Since n depends on the matrix size M theirsimulations need to run on a supercomputer. We instead usea quantum circuit design which runs linear to M and can besimulated easily on a classical computer. This gives us thepossiblity to simulate multidimensional problems. To entanglethe eigenvalues with their eigenstates PEA needs to run severalHamiltonian simulations with the power up to n − . U k = Se t ∆2 k S for k = 0 , · · · , n − (4)Thus the application of multiple U can be directly encoded inthe phase shifts of e t ∆ , so that the runtime of U k does notincrease. The Hamiltonian simulation of a diagonal matrix canbe split into: e t ∆2 k = (cid:81) M − j =1 diag (1 , ..., j − , e λ j t k , j +1 , ..., M − ) (5)This can be implemented with M − phase shift gates R j = (cid:18) e iλ j t k (cid:19) . Because we do not use any oracle producingthese phases, they have to be calculated classicly. This isthe reason why this typ of circuit design does not give anyexponential speedup in the one-dimensional case. It can bedemonstrated how to simulate multidimensional Poisson ma-trices by using U. The splitting formula (1) and its Hamiltoniansimulation show that the simulation of a d-dimensional prob-lem can be implemented with d one-dimensional Hamiltoniansimulations [1]. In fact this breaks the curse of dimensionality Fig. 1: Hamiltonian simulation of a 1-sparse matrix with seveneigenvalues using M − phase shift gatesand brings the exponential quantum speed up with respect tod. Figure (2) shows the quantum circuit for PEA of a two-dimensional problem. The ancilla qubit, which is needed forFig. 2: Circuit for PEA with M = 4 and d = 2 applying the discrete sine transform, can be reused for everyHamiltonian simulation, because it is kept in state | (cid:105) . Thiscircuit design is the reason for the specific encoding of registerC, since every log ( M ) qubit block uses the subroutine U sep-arately. The phase estimation algorithm entangles eigenstateswith eigenvalues of A d and produces the following system: | ψ (cid:105) = (cid:88) j b j |≈ λ j (cid:105) | (cid:126)u j (cid:105) with b j = (cid:104) (cid:126)b | (cid:126)u j (cid:105) C. Reciprocal calculation of λ j (INV) For the next step the eigenvalues have to be entangled withits reciprocals. We use one iteration of the Newton-Raphson-Division method to calculate their values: λ − ≈ x − λx (6)For the start value x we use an approximation of thereciprocal: x = 12 p , p ∈ R with | p − λ j | min (7)Let k = k , ..., k n − be the binary representation of λ j as aninteger. The following circuit describes a method to store x in a new n-sized register A for λ j ≥ . The circuit iteratesthrough all qubits in register B starting with the qubit with thehighest value. If it reads a | (cid:105) block, x will be rounded off,otherwise x is greater than the reciprocal of λ j . The iterationstops after the first one is read. The ancilla qubit indicatesat the end whether the inversion happened or not. To get amore accurate reciprocal of λ j the first iteration of the newtonmethod gets calculated: x = 2 − p +1 − λ − p (8)ig. 3: Quantum circuit to transform | (cid:105) | k j (cid:105) into | x j (cid:105) | k j (cid:105) This could be implemented with one quantum subtraction, onemultiplication and one squaring circuit. Most of the knownquantum multiplication circuits store the result in an additionaltwo times sized ancilla register [9]. Even if truncation of theresult is used, the simulation cost of the system would increaserapidly with respect to n. In our work the calculation of x isperformed in the rotational part of the algorithm. Register Aand register B together contain the whole information about x . Instead calculating x as a fixed precision number, ourinversion modul produces a floating point number, where RegBcan be unterstood as the mantissa and RegA as the binaryexponent. The quantum state can now be described as: IN V | ψ (cid:105) = (cid:88) j b j | (cid:126)u j (cid:105) | k j (cid:105) | k − j (cid:105) with k − j = x j D. Rotation of the reciprocal (ROT)
The goal of this module is to rotate the α -fold reciprocal x into the amplitude of an ancilla qubit. Let φ , ..., φ n − φ i ∈ (0 , be the binary representation of a number with the valueorder − , ..., − n . The rotaton of a n-qubit large number canbe implemented with n controlled rotations about the y-axesof the Bloch sphere [1]. R y ( αφ ) = φ R y ( α − ) · φ R y ( α − ) · ... · φ n − R y ( α − n ) with R y ( αφ ) = (cid:18) cos αφ − sin αφ sin αφ cos αφ (cid:19) (9)Every scalar factor of the angle φ can be implemented directlyin the rotation gates. To rotate the ancilla qubit about theangle αx we need two rotations. The first part applies theangle α − p +1 . Since x is not known before runtime, wehave n different cases. This can be solved with n differentcontrolled rotations. The second part rotates the floating pointnumber − αλ − p . The mantissa αλ also can be implementedwith n different controlled rotations described in formula(9). The multiplication with the binary exponent could bedone by shifting the mantissa. This yields bigger registersand larger simulation costs. We instead rotate the mantissan times shifted by the exponential factor each one belongs.Let a , ..., a n − a i ∈ (0 , be the binary representation of x . The rotation of the floating point number αλ − p can besplit into: R y ( αλ − p ) = n − (cid:89) i =0 a i R y ( αλ − i +1) ) (10)The subtraction itself is implemented by rotating in the otherdirection. This method keeps the number of qubits smalland allows us to simulate more difficult Poisson problems,however it needs much more runtime ( O ( n ) rotation gates).The quantum state of the system after using this module lookslike: (cid:88) j b j | (cid:126)u j (cid:105) | k j (cid:105) | k − j (cid:105) ( (cid:114) (1 − sin ( αk j )) | (cid:105) + sin( αk j ) | (cid:105) ) (11)Since the amplitude sin( αk j ) is not a multiple of the recip-rocal, we have to adapt the rotation angle to arcsin( αk j ) .The referenced paper shows an algorithm for approximatingthis by using bisection and their sine function module [1].Because this method uses repeated squaring as a subroutinein the sine function, the simulation cost of the system wouldincrease exponentially for the same reason as the eigenvaluecalculation. Thus we skip this part of the algorithm and use thesmall angle approximation of the sine function: sin( αk j ) ≈ αk j for αk j ∈ [0; 0 . . Skipping the arcsin module gives us theopportunity to shift the calculation of x in the rotational partof the algorithm. Thus we do not need any quantum arithmeticcircuits. E. Amplitude factor α After uncomputing registers A and B a measurement of therotated ancilla qubit leads with success probability Ω succ tostate | (cid:105) . Ω succ = (cid:88) j b j sin ( α ≈ λ j ) (12)If the measurement leads to one, register C collapses to thenormalized solution (cid:126)u of the problem: (cid:126)u = C (cid:88) j b j α ≈ λ j | (cid:126)u j (cid:105) with C = ( (cid:115)(cid:88) j ( b j α ≈ λ j ) ) − Due to the normalization of C small amplitude factors do notinfluence the solution, but they do boost the success probabilityof the algorithm. Doubling the amplitude of the qubit meansfour times better success rates. Because we do not want toproduce big errors, α has a upper boundary of λ . Thus theamplitude for state | (cid:105) is at least κ , so that the runtime ofthe algorithm is O (4 κ ) . By using an arcsin module α can bedoubled, which results in the same runtime described in theHHL09 algorithm. Other methods like variable time amplitudeamplification may be used to achieve closely linear runtimesto κ [10].ig. 4: Rotation circuit to transform state | λ j (cid:105) | x j (cid:105) | (cid:105) into | λ j (cid:105) | x j (cid:105) ( (cid:113) − sin ( αx j ) | (cid:105) + sin( αx j ) | (cid:105) ) III. R
ESULTS
To validate our circuit design for the algorithm we usefollowing example for a problem with M = 4 and d = 1:Ex. 1: A = 16 · − − − − (cid:126)b = with eigenvalues λ = 9 . , λ = 32 , λ = 54 . ,eigenvectors: (cid:126)u = . √ . , (cid:126)u = √ − √ , (cid:126)u = . − √ . and b = , b = √ , b = . First we check if the phaseestimation algorithm of the Hamiltonian simulation producesthe right quantum state | ψ (cid:105) . For the register size of RegAand RegB we choose that λ j can be quantified as fixedprecision integer. Based on formula (3) we can estimate aupper boundary for the eigenvalues: λ max = 4 d M sin ( ( M − π M ) ≤ d M (13)Thus we use n = 2 + log d + 2 log ( M ) qubits for registersA and B. One way to varify the achieved quantum state | ψ (cid:105) is to measure register B. Consequently the wave functioncollapses with probability b j to | λ j (cid:105) | (cid:126)u j (cid:105) . MeasurementsFig. 5: Eigenvalue destribution for 2000 measurements withregister size n-1, n and n+1show the expected result of PEA. The state for the highesteigenvalue is exactly represented with its probability. Sinceeigenvalue λ and λ can not be stored as an integer, PEAbuilds a distribution around its eigenvalue. The probability toget the nearest eigenvalue approximation is at least π [11]. We can use more qubits to increase the resolution whichleads to a more precise eigenvalue distribution. It is possibleto reduce the register size to n − k, where k is the bestsmallest eigenvalue approximation. This leads to a resolutionerror of k − . Thus k less qubits are needed and the registersize is O ( κ ) . We can use this method to reduce the simulationcost of the system. By utilizing the register dump methodof Microsoft’s quantum simulator we can check the state ofregister C. By subtracting the global phase we can see, thatthe amplitudes and phases correspond to the right eigenstatewhich belongs to its eigenvalue.In the following we show the influence of the amplitude factor α to the solution and the success probability of the algorithm.The algorithm transforms register C with probability ofsuccess Ω succ into a state of the form (cid:80) x u ( x ) | x (cid:105) . We useMicrosoft’s quantum simulator to read out the amplitudesof register C. The output encoding for register C has tobe considered. The success probability depends on b j , α and the condition number of A d . By increasing α we canrotate a multiple of the reciprocal. For a valid small angleapproximation doubling α results in a four times higher Ω succ . By choosing α above the upper boundary λ on theFig. 6: Measurements of register C for example one withdifferent α one hand we violate the small angle approximation, on theother the amplitude may be over rotated. This is the reasonwhy high amplitude factors produce high errors in the output.If b is close to zero, we can set the upper boundary for α tothe next higher eigenvalue λ , since the error of sin( αλ ) ismultiplied with b . This concept may be continued to higher a) solution (b) α = 30 , Ω ≈ (c) α = 300 , Ω ≈ Fig. 7: Measured output for different α with M = 8 and d = 2 upper boundaries. Because b j is not known, one methodcould be to increase α stepwise until the algorithm showssuccess with a certain precision.It is posssible to solve more complex problems with higherdimensions. The described algorithm solves multidimensionalproblems by simulating multiple Hamiltonian simulations inthe PEA step. When we increase the dimension of the problemby one, at least log ( M ) additional qubits are needed to run,which leads to a linear qubit count to d. All together we use d + (4 + d ) log ( M ) qubits. i.e. for M = 8 and d = 2
27 qubits are needed. Figure (7) shows the measuredoutput for a two dimensional problem. By changing α wecan choose between accuracy and success probability. While(b) violates the small angle approximation, (c) overrotatesits amplitude in addition. This results in higher errors, butbetter success rates. Thus this circuit design can be used toachieve an approximation of the normalized solution withhigh probability. To get more accurate results the algorithmhas to repeated O (4 κ ) with an amplitude factor less than λ . Amplitude amplification techniques may be used to getbetter runtimes to κ .Since the values of (cid:126)u are encoded in the amplitudes ofregister C, we can not read them out directly. To get the fullsolution of the equation the algorithm has to be successfullyrepeated at least O (( M − d ) times. This does not give us anyquantum advantage over classical methods, though it does ifone is interested in an expectation value of an operator actingon (cid:126)u . For example a measurement of register C may collapsewith high probability into a state | x (cid:105) which represents theposition of a maximum in the solution of the problem. Thusit is possible to get some correlated information about thesolution by using several quantum operators.We showed a modular circuit design which can solvemultidimensional Poisson equations and can be simulatedon a simple classical computer. Since we only use a smallnumber of qubits the memory consumption is kept low. Itmay also be possible to test this design on a real quantumcomputer. By adding some error correcting qubits the today’smost advanced quantum computers could run simple problemswith a similar circuit design. Rwhich represents theposition of a maximum in the solution of the problem. Thusit is possible to get some correlated information about thesolution by using several quantum operators.We showed a modular circuit design which can solvemultidimensional Poisson equations and can be simulatedon a simple classical computer. Since we only use a smallnumber of qubits the memory consumption is kept low. Itmay also be possible to test this design on a real quantumcomputer. By adding some error correcting qubits the today’smost advanced quantum computers could run simple problemswith a similar circuit design. R