Julia SoC '19 : Quantum Algorithms for Differential Equations

"$VERSION"

0.9s

Julia

"1.1.0"

Hi. I am Divyanshu Gupta. I am one of the participants in Julia's Season of Contribution, 2019.

An Introduction

Over the next couple of months, my focus will be on implementing quantum algorithms for solving differential equations. These include algorithms for first-order linear and non-linear (of quadratic nature) differential equations, and partial differential equations such as the wave equation and heat equation.

I am using Yao.jl, a quantum computer simulator that allows you to design quantum algorithms in Julia.

QuDiffEq.jl will serve as the the repository for the aforementioned algorithms.

QuantumBFS/QuDiffEq.jl

Public

A Quantum Algorithm for Linear Differential Equation

Background

A Quantum algorithm can heavily cut down on the time is takes to solve a specific problem, in comparison to its classical counterpart. Though it may seem appealing in theory, a physically feasible algorithm is difficult to come by. There are many caveats to keep in mind - there is no easy way to input data or read the output, and there are questions about how effectively evolution of states can be brought about (Hamiltonian evolution).

An earlier attempt at solving linear differential equations has been presented in Berry's paper, where HHLhas been employed. In this, a multi-step forward difference method is modelled as a system of linear equations which in turn is solved using an HHL circuit. HHL is a linear equations solver that uses the phase estimation algorithm. The HHL circuit has very strict requirements on the kind of Hamiltonian it can simulate. In the months before JSoC, I was involved with implementing these methods.

The algorithm, I shall be implementing, is derived from Berry's paper on 'Hamiltonian simulation by truncated Taylor series'. Berry's approach shows many clear advantages over the former - the complexity has better dependence on the precision and it allows for a greater set of Hamiltonians to be simulated. The article describes it theoretically as well as gives an instance where it has been successfully executed on a quantum computer.

Note: The output from these algorithms is a superposition of computational basis with resultant vector's component as the probability amplitude.So we can't just read out whatever values the registers have for the final answer. In addition, inputing data requires what is called a QRAM. I shall be looking past all of this during my development.

Mathematical details

The algorithm solves problems of the form

\frac{dx}{dt} = Mx + b ,

is an arbitrary matrix,andare dimensional vectors. The exact solution for is give by -

x(t) = e^{Mt}x(0) + (e^{Mt} - I)M^{-1}b\;.

This can be Taylor expanded up to theorder as -

x(t) \approx \sum^{k}_{m=0}\frac{(Mt)^{m}}{m!}x(0) + \sum^{k-1}_{n=1}\frac{(Mt)^{n-1}t}{n!}b\;.

We need to prescribe a way to express the vectors and as states, which then will be acted upon by quantum gates to leave the desired output. As stated earlier, the vectors are encoded in a way that their components are probability amplitudes for computational basis states: and , with as the computational basis states. We can also write as:. is now:

|x(t)\rangle \approx \sum^{k}_{m=0}\frac{||x(0)||(||M||t)^{m}}{m!}\mathcal{M}^{m}|x(0)\rangle + \sum^{k-1}_{n=1}\frac{||b||(||M||t)^{n-1}t}{n!}\mathcal{M}^{n-1}|b\rangle\;.

The task at hand is to simulate an operation on the encoded vector states to collectively obtain the above. We have to consider two cases:

1.is unitary. In addition to the vector state register, there are two ancillary registers of 1 qubit and qubits.

The input to the circuit is an all zero ancilla and work registerthat is control-rotated (by and ) to the required vectors. The circuit for the algorithm can be seen below. The block decomposes of the first register into a superposition of and that couples withand respectively.

Essentially, the first (size 1) register controls the multiplication of and through control-rotation gates (and ). The second register controls the multiplication of (power ofis controlled byin the second register). After decoding the ancilla, the measurement is done in the subspace where the ancilla are in state.

2.is non-unitary. may be expressed as a linear combination of unitaryi.e.. For an arbitrary , this can be done definitively by breaking down into a hermitian and non-hermitian matrices and constructing 2 unitary matrices, from each of the components, which then add up to give . Here.

The quantum circuit consists of three ancillary register:qubit, qubits and-level qudits. The first register has the same role as it does in the previous case.

The second register and the third register together control the multiplication ofon the vector states. In the second register, states with '1's where are raised to probability amplitudes equal to the term with the power in the summation above, while rest of the states are given zero probability. The mappings used is , corresponds to the state in the second register. This register governs the number of times thehas to be operated. blocks generate a superposition with state having amplitude . The third register helps in generating the linear combination - it ensures that are multiplied with the correct coefficient. The ancillary registers are decoded and the measurement is again done in the subspace where the ancilla are in the state.

Usage and Example

Constructing an arbitrary linear differential equation problem:

using QuDiffEq
using OrdinaryDiffEq, Test
N = 1 #size of the work register
k = 3 #order of approximation
siz = 1<<N
b = (rand(ComplexF64, siz)) 
u0 = (rand(ComplexF64, siz)) #initial condition
M = rand(ComplexF64,siz,siz)
tspan = (0.0, 0.4)
qprob = QuLDEProblem(M, b, u0, tspan) #problem is defined as a QuLDEProblem

327.6s

Julia

QuLDEProblem with uType Array{Complex{Float64},1} and tType Float64. In-place: true timespan: (0.0, 0.4) u0: Complex{Float64}[0.0841947+0.900611im, 0.640969+0.869549im]

The algorithm is named QuLDE(). It takes the argument k to build a circuit of theorder approximation. The algorithm has the following two usage:

For linear differential equations:

solve(::QuLDEProblem,::QuLDE)

For non-linear differential equations through linearisation (covered below):

solve(::ODEProblem,::QuLDE)

out = solve(qprob, QuLDE(k))

19.1s

Julia

2-element Array{Complex{Float64},1}: -0.323781+1.8603im 0.331014+1.30299im

To compare the result, the problem is solved using Tsit5() from OrdinaryDiffEq.jl

f(u,p,t) = M*u + b;
prob = ODEProblem(f, u0, tspan)
sol = solve(prob, Tsit5())
s = sol.u[end]

@test isapprox.(s,out, atol = 0.05) |> all

7.8s

Julia

Test Passed

We shall use this algorithm to solve a system of non-linear equation.

\begin{array}{rcl} \frac{dx}{dt} & = & f(x,y) \;, \\ \frac{dy}{dt} &=& g(x,y) \;.\end{array}

We proceed by linearising the equations about the point and plugging the new matrix and vector into the circuit. This is repeated multiple times, with being updated every iteration.

J = \begin{pmatrix} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \\ \frac{\partial g}{\partial x} &\frac{\partial g}{\partial y} \end{pmatrix}\:,

b = \begin{pmatrix} f(x^{*},y^{*})\\ g(x^{*},y^{*}) \end{pmatrix}\;.

The differential equation to be solved will be of the form

\frac{d\Delta u}{dt} = J * \Delta u + b\:,

whereis naturally zero. We then have, at end the of every iteration.

using QuDiffEq
using OrdinaryDiffEq, Test
"""
Linearizing non linear ODE and solving it using QuLDE
du/dt = f
u0 -> inital condition for the ode
Δu -> difference for fixed point
h -> time step
k -> order of Taylor Expansion in QuLDE circuit
Equation input to QuLDE circuit : d(Δu)/dt = J * Δu + b
"""
# Arbitary non-linear equation
function fo(du,u,p,t)
  du[1] = -2*u[2]^2*u[1]
  du[2] = 3*u[1]^(1.5) - 0.1*u[2]
end
u0 = [0.2,0.1]
h = 0.1 
k = 3
tspan = (0.0,0.8)

prob = ODEProblem(fo,u0,tspan)
qsol = solve(prob,QuLDE(k),dt = h)
sol = solve(prob,Tsit5(),dt = h,adaptive = false)
v = transpose(hcat(sol.u...))

@test isapprox.(v,qsol, atol = 1e-3) |> all

9.2s

ExampleJulia

Test Passed

References

High-order quantum algorithm for solving linear differential equations: https://arxiv.org/abs/1010.2745v2
Quantum algorithm for solving linear systems of equations: https://arxiv.org/abs/0811.3171
A Quantum Algorithm for Solving Linear Differential Equations: Theory and Experiment: arxiv.org/abs/1807.04553
Simulating Hamiltonian dynamics with a truncated Taylor series: https://arxiv.org/abs/1412.4687
Quantum Machine Learning Algorithms: Read the Fine Print