Adiabatic evolution

[scarrazza]

We believe that the inclusion of an annealing simulator in the first release of Qibo is a good idea. One possibility is to implement the QUBO graph and then proceed with the following algorithms:

• Monte-Carlo based simulated quantum annealer, for approximate solutions with large number of qubits.
• Brute-force search, for smaller number of qubits.

The DWave Ocean SDK provides a C++ library for simulated annealing called neal, which uses samples with QUBO and Ising. There are other libraries such as sqaod which may be relevant for benchmarking.

[scarrazza]

Let me try to summarize the discussion we had this morning.

We are planning to implement in Qibo an Adiabatic Quantum Computation algorithm based on the Adiabatic Evolution presented by JI.

This is one of the possible approaches for quantum annealing, it is very close to what quantum hardware is designed to do. This approach is not public available in other libraries. Public libraries implement classical annealing algorithms (see post above).

Implementation

For this first implementation in Qibo our target goal is:

1. simulate the time evolution of a given Hamiltonian, using Schrödinger equation, and numerical ODE solvers (with multi-threading or GPU support)
2. attach a scheduling optimization algorithm which given a target ground state, tunes the adiabatic steps by re-evaluating point 1.
3. the scheduling could be used to achieve the trotterization of the Hamiltonian evolution and thus obtain local terms that can be translated into gates.

While point 1 is simulation specific, points 2 and 3 can be reused on real quantum hardware, i.e. for large number of qubits we use the quantum hardware to perform time evolution while the classical computer optimizes the scheduling. The point 3 is also interesting because allows us to perform an adiabatic evolution using quantum circuit hardware (instead of annealing), so this seems more like shifting the "simulation" to a real device emulation (i.e. circuits simulating time evolution).

Practical implementation suggestion

1. we should create an AdiabaticEvolution class in models.py that takes an input hamiltonian, fixed scheduling, and performs a numerical ODE solution, e.g. using rk4.
2. re-factor the VQE minimization and loss methods to determine the best steps. We can start with numpy basic optimizers and eventually move to GAs and NN.
3. after trotterization convert the local operators to known gates (using the gates decomposition idea).

Please correct me if you see some misinterpretation or wrong assumption.

[stavros11]

Thank you for the great summary. I agree with both the general points and the practical implementation for Qibo. As a start, I did a simple integration of Schrodinger's equation using RK4 in numpy. I hide the code in a spoiler because it is a bit long but it should be possible to execute it if you copy/paste:

Adiabatic evolution using Runge-Kutta in numpy

```Python """Implements adiabatic evolution using Runge-Kutta integrator on numpy.""" import argparse import numpy as np from scipy.linalg import expm import matplotlib.pyplot as plt from qibo import matrices def multikron(matrices): h = 1 for m in matrices: h = np.kron(h, m) return h class TimeHamiltonian: """Hamiltonian used for adiabatic evolution. Args: h0 (np.ndarray): Easy Hamiltonian matrix of shape (2^nqubits, 2^nqubits). This Hamiltonian should be easy to prepare (typically sum(X)). h1 (np.ndarray): Hard Hamiltonian matrix of shape (2^nqubits, 2^nqubits). This Hamiltonian encodes the problem. s (Callable): Function s(t) that defines the scheduling. """ def __init__(self, h0, h1, s): self.h0 = h0 self.h1 = h1 self.s = s def __call__(self, t): """Returns the Hamiltonian matrix at a given time. Calculated as H(t) = (1 - s(t)) * H0 + s(t) * H1 Args: t (float): Time value. """ return (1 - self.s(t)) * self.h0 + self.s(t) * self.h1 def initial_state(self): """Returns the ground state of ``h0`` to use as initial state.""" _, psi = np.linalg.eigh(self.h0) return psi[:, 0] def energy(self, psi): """Calculates the

expectation value on a given state.""" return psi.conj().dot(self.h1.dot(psi)) def eigenvalues(self): """Calculates the eigenvalues of the hard Hamiltonian ``h1``. Used for checking the final solution. """ return np.linalg.eigvalsh(self.h1) class StateEvolution: """Keeps track of energies and norms during the state evolution. Every time the state vector is updated for the next time step its energy and norm is logged. """ def __init__(self, hamiltonian): self.hamiltonian = hamiltonian self._vector = hamiltonian.initial_state() self._energies = [] self._norms = [] @property def vector(self): return self._vector @property def energies(self): return np.array(self._energies) @property def norms(self): return np.array(self._norms) @vector.setter def vector(self, vector): self._vector = vector self._energies.append(self.hamiltonian.energy(vector)) self._norms.append(np.sqrt((np.abs(vector) ** 2).sum())) def sx_sum(nqubits): """Creates sum(X) Hamiltonian matrix.""" return - sum(multikron((matrices.X if i == j % nqubits else matrices.I for j in range(nqubits))) for i in range(nqubits)) def ising(nqubits): """Creates Ising sum(ZZ) Hamiltonian matrix.""" return -sum(multikron((matrices.Z if i in {j % nqubits, (j+1) % nqubits} else matrices.I for j in range(nqubits))) for i in range(nqubits)) def main(nqubits, T, dt): # Create Hamiltonian using sum(X) as the easy Hamiltonian and Ising (sum(ZZ)) as the problem Hamiltonian # Use linear scheduling s(t) = t / T ham = TimeHamiltonian(sx_sum(nqubits), ising(nqubits), lambda t: t / T) nsteps = int(T / dt) # Runge-Kutta steps overlaps = [] state = StateEvolution(ham) # initialize state to evolve with RK state_exp = StateEvolution(ham) # initialize state to evolve with exp(-1j dt H) for n in range(nsteps): # Evolution using RK t = n * dt k1 = ham(t).dot(state.vector) k2 = ham(t + dt / 2.0).dot(state.vector + dt * k1 / 2.0) k3 = ham(t + dt / 2.0).dot(state.vector + dt * k2 / 2.0) k4 = ham(t + dt).dot(state.vector + dt * k3) state.vector += -1j * dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6.0 # Evolution using exponential state_exp.vector = expm(-1j * dt * ham(t)).dot(state_exp.vector) # Overlap between the two evolutions overlaps.append(state_exp.vector.conj().dot(state.vector)) target_energy = ham.eigenvalues()[0] tt = np.linspace(0, T, nsteps) plt.figure() plt.subplot(131) # normalize overlap overlaps = np.array(overlaps) / (state.norms * state_exp.norms) plt.plot(tt, 1 - overlaps) plt.xlabel("$t$") plt.ylabel("1 - Overlap") plt.subplot(132) plt.plot(tt, state.norms - 1, color="blue") plt.plot(tt, state_exp.norms - 1, color="red") plt.xlabel("$t$") plt.ylabel("Norm - 1") plt.subplot(133) plt.plot(tt, state.energies, color="blue") plt.plot(tt, state_exp.energies, color="red", linestyle="--") plt.axhline(y=target_energy, color="green") plt.xlabel("$t$") plt.ylabel(r"$\left \langle H_0\right\rangle $") plt.show() if __name__ == "__main__": parser = argparse.ArgumentParser() parser.add_argument("--nqubits", default=3, help="Number of qubits", type=int) parser.add_argument("--T", default=1, help="Total evolution time", type=float) parser.add_argument("--dt", default=1e-4, help="Runge-Kutta time step", type=float) args = parser.parse_args() main(**vars(args)) ```

This follows JI's slides and defines H(s) = (1 - s)H0 + sH1 where the "easy" Hamiltonian H0 is the common sum of X and the "problem" Hamiltonian H1 is the classical Ising (sum of ZZ). I use a linear scheduling s(t) = t / T but this can be easily changed.

I integrate the time-dependent Schrodinger equation using two different ways (RK4 and multiplying with the propagator exp(-1j dt H)) and keep track of the norm and H1 energy for both evolutions as well as the overlap between the two evolutions. The main observations so far are:

• The two integration methods agree relatively well (overlap is close to 1), at least for small nqubits. RK's accuracy may drop as nqubits is increased though, requiring smaller dt.
• RK doesn't conserve the norm exactly (should be conserved since evolution is unitary). Integration is expected to have some numerical errors but I am not sure how big of an issue this will be for more qubits.
• H1 energy converges to the Ising ground state energy with some tuning of T (the total time of adiabatic evolution). For example 3 qubits need T > 3. I guess a better choice of scheduling can reduce the required T.

If you agree, I could start by implementing a Qibo model based on this simple RK implementation written in Tensorflow instead of numpy. Then we can continue with some benchmark and check other integrators that are more accurate or easier to parallelize for GPU.

One last note is that in the script I calculate the Hamiltonian as a (2^N, 2^N) array and take dot products to the state vector. For many qubits it will probably be more efficient (and will use less memory) to do this as a sum of local terms to avoid creating huge matrices. For example the H1.dot(psi) for the Ising Hamiltonian could be calculated as sum(ZZ.dot(psi)) where our custom operators can be used for Z.dot(psi). This can be generalized for any local Hamiltonian and may also be useful for the VQE.

[scarrazza]

Thanks for the snipped, the structure is pretty much clear and follows what we have discussed, so I agree in porting this interface to QIbo.

Concerning performance, your script requires just 0.5s on my CPU to evaluate the rk4 and 2s for the exponential solution, so we may provide as option both solutions.

The main difficulty of this algorithm is to find the proper ODE solver, which maintains quality and is parallelizable.

[scarrazza]

Here some references collected by Artur:

• https://arxiv.org/pdf/1901.05824.pdf

Fundamentals on QA:

• https://arxiv.org/abs/quant-ph/0001106
• https://arxiv.org/abs/cond-mat/9804280

QSA vs CA:

• https://arxiv.org/abs/1601.03030

Implementations:

• E. M. Inack and S. Pilati, Phys. Rev. E 92, 053304 (2015)
• https://projecteuclid.org/euclid.ss/1475001234