Implement numerical simulation of the Ising Hamiltonian evolution

[SimoneGasperini]

It would be interesting to have the possibility to numerically simulate the full Ising Hamiltonian (as an Hermitian complex matrix) time evolution during the annealing process. This would allow to study for instance the behaviour over time of the minimum energy gap (i.e. difference between the lowest eigenvalue $\lambda_0$, corresponding to the Hamiltonian ground state, and the second-lowest eigenvalue $\lambda_1$, corresponding to the Hamiltonian first excited state).

The superconducting QPU at the heart of the D‑Wave quantum annealing system is a controllable, physical realization of the quantum Ising spin system in a transverse field. Analog control lines enables to implement the following time-dependent Hamiltonian:

In the equation above, $\hat{\sigma}_{x}$ and $\hat{\sigma}_{z}$ represent the Pauli-X and Pauli-Z matrix operators, $h_{i}$ and $J_{i,j}$ are the linear and quadratic coefficients defining the final Hamiltonian (encoding the optimization problem), and the functions $A(s)$ and $B(s)$ control the annealing schedule ($s$ is time rescaled in $[0,1]$). See D-Wave Systems docs for an extensive description of the quantum annealing implementation and controls.

[SimoneGasperini]

Here is a sketch of a naive implementation just to give an idea:

import functools
import numpy as np

identity = np.array([[1, 0],
                     [0, 1]])

sigma_x = np.array([[0, 1],
                    [1, 0]])

sigma_z = np.array([[1, 0],
                    [0, -1]])

def get_H_init(num_qubits):
    H_init = np.zeros(shape=(2**num_qubits, 2**num_qubits))
    for i in range(num_qubits):
        terms = [identity] * num_qubits
        terms[i] = sigma_x
        H_init -= functools.reduce(np.kron, terms)
    return H_init

def get_H_final(num_qubits, h, J):
    H_final = np.zeros(shape=(2**num_qubits, 2**num_qubits))
    for i in range(num_qubits):
        terms = [identity] * num_qubits
        terms[i] = sigma_z
        H_final += h[i] * functools.reduce(np.kron, terms)
    for i in range(num_qubits):
        for j in range(i + 1, num_qubits):
            terms = [identity] * num_qubits
            terms[i] = sigma_z
            terms[j] = sigma_z
            H_final += J[i, j] * functools.reduce(np.kron, terms)
    return H_final

def anneal_schedule(s):
    A = 1 - s
    B = s
    return A, B

def get_ising_H(H_init, H_final, s):
    A, B = anneal_schedule(s=s)
    ising_H = A/2 * H_init + B/2 * H_final
    return ising_H

To make it more efficient, a sparse representation of the Hamiltonian matrix could be used (see scipy.sparse module). Note that the linalg sub-module provides for instance the function eigsh(A, k=2, which='SA') to compute the first $k=2$ smallest eigenvalues of an Hermitian sparse matrix $A$.