Using QubitOperators

One of the two core data structures used by Fulqrum is the QubitOperator that represents Hamiltonians defined over qubit subsystems. Here we will go through examples of the core QubitOperator functionality that is pertinent to most use cases. To begin, we will import fulqrum:

import fulqrum as fq

Single-term operators

Create an empty operator over 5 qubits:

fq.QubitOperator(5)

<QubitOperator[], width=5>

Create a N qubit single-term operator from a dense (including identity operators) string, where N is equal to the length of the string:

fq.QubitOperator.from_label("IIXYI")

<QubitOperator[('Y:1 X:2', (1+0j))], width=5, extended=0, group=-1>

Create a 5 qubit, single-term operator using sparse notion specifying an X operator on qubit 0 and a 1 operator on qubit 3, and with a coefficient of 2.2:

H = fq.QubitOperator(5, [("X1", [0, 3], 2.2)])
H

<QubitOperator[('X:0 1:3', (2.2+0j))], width=5, extended=1, group=-1>

Check the properties of the QubitOperator. Number of qubits:

H.width

5

Number of terms:

H.num_terms

1

or the number of terms can be found as:

len(H)

1

Multi-term operators

Create a multi-term QubitOperator using the sparse notion:

H = fq.QubitOperator(5, [("X1", [0, 3], 2.2), ("ZZY", [0, 1, 4], -1)])
H

<QubitOperator[('X:0 1:3', (2.2+0j)), ('Z:0 Z:1 Y:4', (-1+0j))], width=5>

H.num_terms

2

Create a multi-term operator by adding one operator to another, returning a new operator:

H1 = fq.QubitOperator.from_label("ZIXYI", 2)
H2 = fq.QubitOperator(5, [("ZZIZZ", range(5), -1.0)])
H1 + H2

<QubitOperator[('Y:1 X:2 Z:4', (2+0j)), ('Z:0 Z:1 Z:3 Z:4', (-1+0j))], width=5>

Create a multi-term operator by adding an operator in-place to an existing operator:

H1 += H2
H1

<QubitOperator[('Y:1 X:2 Z:4', (2+0j)), ('Z:0 Z:1 Z:3 Z:4', (-1+0j))], width=5>

Operator splitting

Often times it is of interest to decompose an operator into its diagonal and off-diagonal components. For example, consider th following lattice Hamiltonian:

from qiskit.transpiler import CouplingMap

# Build 16-qubit coupling map
N = 4
cmap = CouplingMap.from_grid(N, N)
num_qubits = cmap.size()
num_edges = 2 * N * (N - 1)

# 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:
        touched_edges.add(edge)
        H += fq.QubitOperator(
            num_qubits,
            [("XX", edge, coeffs[0]), ("YY", edge, coeffs[1]), ("ZZ", edge, coeffs[2])],
        )

H.num_terms

72

The operator consists of a collection of diagonal terms (ZZ) and off-diagonal terms (XX or YY). We can split these comtributions into their respective pieces:

diag_H, offdiag_H = H.split_diagonal()

The number of diagonal terms should be half that of the off-diagonal terms:

diag_H.num_terms

24

offdiag_H.num_terms

48

We can of course inspect the operators to see that the splitting works, e.g. for the diagonal component:

diag_H

<QubitOperator[('Z:0 Z:4', (-1+0j)), ('Z:0 Z:1', (-1+0j)), ('Z:1 Z:5', (-1+0j)), ('Z:1 Z:2', (-1+0j)), ('Z:2 Z:6', (-1+0j)), ('Z:2 Z:3', (-1+0j)), ('Z:3 Z:7', (-1+0j)), ('Z:4 Z:8', (-1+0j)), ('Z:4 Z:5', (-1+0j)), ('Z:5 Z:9', (-1+0j)), ('Z:5 Z:6', (-1+0j)), ('Z:6 Z:10', (-1+0j)), ('Z:6 Z:7', (-1+0j)), ('Z:7 Z:11', (-1+0j)), ('Z:8 Z:12', (-1+0j)), ('Z:8 Z:9', (-1+0j)), ('Z:9 Z:13', (-1+0j)), ('Z:9 Z:10', (-1+0j)), ('Z:10 Z:14', (-1+0j)), ('Z:10 Z:11', (-1+0j)), ('Z:11 Z:15', (-1+0j)), ('Z:12 Z:13', (-1+0j)), ('Z:13 Z:14', (-1+0j)), ('Z:14 Z:15', (-1+0j))], width=16>

Operator groups

“Groups” represent one of the fundamental properties of an operator for eigensolving. A group consists of all terms that share the same off-diagonal structure, i.e. terms that have off-diagonal operators, e.g. X, Y, +, or -, on the same qubit indices. We can see the number of groups easily:

H.num_groups

25

In this case we have a single group that contains all the diagonal components, i.e. the off-diagonal structure is trivial since they are diagonal, and then 24 more groups, each one corresponding to XX and YY on each of the 24 edges in the graph. We can verify this using the split operator pieces from before:

diag_H.num_groups

1

offdiag_H.num_groups

24

The importance of groups of terms come from the fact that, for a given row, each group represents a matrix-element in the Hamiltonian. We see that the diagonal Hamiltonian has a single group, indicating the obvious fact that a diagonal must have a single matrix element per row. The off-diagonal Hamiltonian has 24 groups, and thus we can expect at most 24 off-diagonal elements per row. The “at most” here represents the fact that, in a subspace, not all matrix-elements need be found within the subspace. That there are 48 total terms in the off-diagonal operator, but only 24 groups is due to the fact that each group is comprised of the XX and YY terms on each edge of the graph.

We can slice the Hamiltonian by computing the group pointers, indicating the start and stop indices of terms within each group:

group_ptrs = offdiag_H.group_ptrs()
group_ptrs

array([ 0,  2,  4,  6,  8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32,
       34, 36, 38, 40, 42, 44, 46, 48], dtype=uint64)

We see that, for example, the 0th group starts with term 0 and goes to term 2, ie. terms 0 and 1 are in the group. We can slice the operator to see this:

offdiag_H[0:2]

<QubitOperator[('X:0 X:1', (-0.5+0j)), ('Y:0 Y:1', (-0.5+0j))], width=16>

showing that indeed, the group is made up of XX and YY terms on the [0, 1] edge of the graph. Other groups are similar.

Removing constants from an operator

Physics says that it is always possible to add a constant energy to the Hamiltonian. This has the effect of shifting all of the operator eigenenergies, but does not change the energy splittings. It is best to remove this constant before eigensolving. Let us add a constant to our lattice Hamiltonian

Hc = fq.QubitOperator.from_constant(H.width, 100)
Hc += H

We not only want to remove the constant term, but would also like to return the energy value as well, e.g. to add back later:

H0, constant_energy = Hc.remove_constant_terms()
constant_energy

100.0

We can verify that H0 has no constant energy remaining:

H0.constant_energy()

0.0