Heisenberg spin model

import numpy as np
import scipy.sparse.linalg as spla
from qiskit.transpiler import CouplingMap
import fulqrum as fq

Construct operator

Here we construct a 16-qubit Heisenberg spin model over a square-grid using a Qiskit CouplingMap to define the grid

# Build 16-qubit coupling map
cmap = CouplingMap.from_grid(4, 4)
num_qubits = cmap.size()

# Generate Hamiltonian
H = fq.QubitOperator(num_qubits, [])
touched_edges = set({})
coeffs = [-1 / 2, -1 / 2, -1]
for edge in cmap.get_edges():
    if edge[::-1] not in touched_edges:  # Only add edge once to Hamiltonian
        touched_edges.add(edge)
        H += fq.QubitOperator(
            num_qubits,
            [("XX", edge, coeffs[0]), ("YY", edge, coeffs[1]), ("ZZ", edge, coeffs[2])],
        )

Construct subspace

In lieu of actually running on hardware or a simulator, here we generate a fake set of bit-strings:

counts = []
for kk in range(2**num_qubits):
    counts.append(bin(kk)[2:].zfill(num_qubits))

S = fq.Subspace([counts])

Create subspace Hamiltonian and initial vector

The Hamiltonian and subspace are combined into a SubspaceHamiltonian object for eigensolving. Here we also set the initial vector for reproducibility.

Hsub = fq.SubspaceHamiltonian(H, S)
v0 = np.ones(S.size(), dtype=Hsub.dtype)

Perform diagonalization

Here we perform the diagonalization twice; first matrix-free and then using a CSR matrix

Matrix-free

evals, evecs = spla.eigsh(Hsub, k=1, which="SA", v0=v0)

CSR-matrix solution

M = Hsub.to_csr_linearoperator_fast()
evals2, evecs2 = spla.eigsh(M, k=1, which="SA", v0=v0)

Compare solutions

First compare relative error in found eigenvalues

abs((evals[0] - evals2[0]) / evals[0])

np.float64(2.960594732333747e-16)

Show that the two eigenstates agree:

Hsub.interpret_vector(evecs)

{'0000000000000000': 0.7071067811865476,
 '1111111111111111': 0.7071067811865476}

Hsub.interpret_vector(evecs2)

{'0000000000000000': 0.7071067811865476,
 '1111111111111111': 0.7071067811865476}