2. Quantum Hamiltonians as matrix product operators

In this tutorial, we first introduce the concept of operator chains and demonstrate how to convert them to MPO form. This mechanism facilitates a concise construction of quantum Hamiltonians as MPOs.

# import NumPy and the PyTeNet package
import numpy as np
import pytenet as ptn

2.1. Operator chains

An “operator chain” defined in opchain.py consists of an outer product of local operators, i.e.,

\[c \cdot \mathrm{op}_i \otimes \mathrm{op}_{i+1} \otimes \cdots \otimes \mathrm{op}_{i+n-1},\]

with \(c\) a coefficient and \(\mathrm{op}_i\) acting on lattice site \(i\) (and identity operations on the remaining lattice sites). Here is a basic illustration:

# Pauli matrices
sigma_x = np.array([[0.,  1.], [1.,  0.]])
sigma_y = np.array([[0., -1j], [1j,  0.]])
sigma_z = np.array([[1.,  0.], [0., -1.]])

# store local operators in a dictionary
opmap = {
    0: np.identity(2),
    1: sigma_x,
    2: sigma_y,
    3: sigma_z
}

# create the symbolic operator chain 0.7 * I I X X Z I I ...
opchain = ptn.OpChain([1, 1, 3], [0, 0, 0, 0], 0.7, 2)
# The second parameter contains quantum numbers to be interleaved with the operators,
# see the section below for an explanation.
# The last argument is the first site which the operator chain acts on other than identities.

Our opchain simply stores the arguments used for constructing it as instance variables:

opchain.oids

[1, 1, 3]

opchain.qnums

[0, 0, 0, 0]

opchain.coeff

0.7

opchain.istart

2

The to_matrix() function forms the outer product, i.e., the matrix representation of the operator chain (without the leading and trailing identities):

opchain.to_matrix(opmap)

array([[ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0.7,  0. ],
       [ 0. , -0. ,  0. , -0. ,  0. , -0. ,  0. , -0.7],
       [ 0. ,  0. ,  0. ,  0. ,  0.7,  0. ,  0. ,  0. ],
       [ 0. , -0. ,  0. , -0. ,  0. , -0.7,  0. , -0. ],
       [ 0. ,  0. ,  0.7,  0. ,  0. ,  0. ,  0. ,  0. ],
       [ 0. , -0. ,  0. , -0.7,  0. , -0. ,  0. , -0. ],
       [ 0.7,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ],
       [ 0. , -0.7,  0. , -0. ,  0. , -0. ,  0. , -0. ]])

2.2. Converting operator chains to MPOs

How can we construct an MPO version of a sum of operator chains? A useful mental picture is railway tracks running from east to west (i.e., from the rightmost to the leftmost lattice site), possibly interconnected by switches. Each operator chain corresponds to a train running once from east to west along one route through the tracks. (There is no risk of collision since trains operate at different hours.) Markers are placed at uniform intervals along the tracks, each displaying a local operator (such as the Pauli matrices in the example above). The operator chain (train) collects the markers that it encounters during its ride (preserving order).

We can exploit the locality of typical operator chains (e.g., acting non-trivially only on two neighboring lattice sites) by creating a special “identity” track whose markers display only the local identity operation. Operator trains share this identity track before and after the sites on which they act.

To convert the operator chains to an MPO, we use an intermediate OpGraph data structure that represents the overall track layout as a graph. We can then assemble an MPO from it. The number of tracks running in parallel (at a given longitude) is precisely the virtual bond dimension.

# first create two additional operator chains:
# 0.5 * I Y Z I I ...
opchain2 = ptn.OpChain([2, 3], [0, 0, 0], 0.5, 1)
# -1.3 * I I I X I I ...
opchain3 = ptn.OpChain([1], [0, 0], -1.3, 4)

# overall system size (number of lattice sites)
nsites = 8

# construct the graph;
# last argument is the operator index (ID) in the above `opmap` of the identity operation
opgraph = ptn.OpGraph.from_opchains([opchain, opchain2, opchain3], nsites, 0)
# internal consistency check
print(opgraph.is_consistent(verbose=True))

True

# physical quantum numbers (see below)
qsite = [0, 0]

# construct MPO representation of the operator graph
mpo = ptn.MPO.from_opgraph(qsite, opgraph, opmap)

# need at most 3 virtual bonds
mpo.bond_dims

[1, 1, 2, 3, 3, 1, 1, 1, 1]

# consistency check
# include leading and trailing identities to construct reference matrix represenation
opchains_mat = sum(np.kron(np.kron(
    np.identity(2**opc.istart),
    opc.to_matrix(opmap)),
    np.identity(2**(nsites - (opc.istart + opc.length))))
        for opc in [opchain, opchain2, opchain3])
# difference should be zero
np.linalg.norm(mpo.to_matrix() - opchains_mat)

np.float64(0.0)

2.3. Constructing quantum Hamiltonians

The module hamiltonian.py uses the mechanism to construct quantum Hamiltonians as MPOs. As an illustration, for the Ising Hamiltonian

\[H = \sum_{j=1}^{L-1} J\,\sigma^z_j \sigma^z_{j+1} + \sum_{j=1}^L \left( h\,\sigma^z_j + g\,\sigma^x_j \right)\]

(where we have omitted the \(\otimes\) symbol between \(\sigma^z_j\) and \(\sigma^z_{j+1}\) for brevity), the function ising_MPO(nsites, J, h, g) in hamiltonian.py creates an operator chain for \(J\,\sigma^z_0 \sigma^z_1\) and one for \(h\,\sigma^z_0 + g\,\sigma^x_0\), and then shifts them along the lattice (i.e., performs the sum over \(j\)). The same procedure works for other quantum Hamiltonians and is not restricted to nearest-neighbor interactions.

# Hamiltonian parameters
J =  1.0
h = -0.4
g =  0.7
# construct Ising Hamiltonian as MPO
mpo_ising = ptn.ising_1d_mpo(nsites, J, h, g)

# virtual bond dimensions
mpo_ising.bond_dims

[1, 3, 3, 3, 3, 3, 3, 3, 1]

2.4. Quantum numbers and conservation laws

Let’s discuss how (additive) quantum numbers enter the story so far. In many cases, a Hamiltonian respects conservation laws, such as the conservation of total spin or particle number in a physical system. How to exploit this within the matrix product operator formalism is not entirely obvious, but here is how it works: we first identify the physical quantum numbers at each lattice site, like \(\pm \frac{1}{2}\) for spin-up or spin-down. (PyTeNet stores them in variables denoted qsite.) Likewise, each virtual bond is also associated with a quantum number (qbonds in PyTeNet). Now every MPO tensor obeys a sparsity pattern dictated by the quantum numbers: The sum of the first physical and left virtual bond quantum number of each non-zero tensor entry is equal to the sum of the second physical and right virtual bond quantum number. (Here, “first” and “second” refer to the row and column dimensions of a site-local operator.) PyTeNet provides the is_qsparse utility function to probe such a sparsity pattern. Note that we do not have to manually enforce it; instead, the sparsity pattern naturally arises when constructing the MPO representation of a Hamiltonian, and quantum numbers provide the means to recognize and understand it. In practical terms, quantum numbers enable the optimization of common operations, such as singular value decompositions, via partitioning into non-zero blocks.

How can we obtain the virtual bond quantum numbers in the first place? Let’s illustrate this via the XXZ Heisenberg model, with Hamiltonian represented as

\[H = \sum_{j=1}^{L-1} \left( \tfrac{1}{2} J\,S^{+}_j S^{-}_{j+1} + \tfrac{1}{2} J\,S^{-}_j S^{+}_{j+1} + \Delta\,S^z_j S^z_{j+1} \right) - \sum_{j=1}^L h\,S^z_j\]

where \(S^{\pm}_j\) are the spin raising and lowering operators at site \(j\), respectively, and \(S^z_j = \frac{1}{2} \sigma^z_j\). As one might guess, the raising and lowering operators change the spin quantum number by \(1\) (such that the overall net effect of both \(S^{+}_j S^{-}_{j+1}\) and \(S^{-}_j S^{+}_{j+1}\) is zero). We can translate this into code by sandwiching the virtual bond quantum number \(1\) between \(S^{+}_j\) and \(S^{-}_{j+1}\), and likewise \(-1\) between \(S^{-}_j\) and \(S^{+}_{j+1}\). (To avoid any numerical rounding issues, quantum numbers are generally stored as integers, and thus we multiply all spin quantum numbers by \(2\) in the code.) Using the same functionality as above, the quantum numbers are copied from the operator chains to the graph and, finally, to the virtual bond quantum numbers of the MPO, endowing the MPO representation of the Hamiltonian with quantum numbers. See hamiltonian.py for the concrete implementation.

# Hamiltonian parameters
J =  1.0
D =  0.8
h = -0.1
# construct XXZ Heisenberg Hamiltonian as MPO
mpo_xxz = ptn.heisenberg_xxz_1d_mpo(nsites, J, D, h)

# physical quantum numbers (multiplied by 2)
mpo_xxz.qsite

array([ 1, -1])

# virtual bond dimensions
mpo_xxz.bond_dims

[1, 4, 5, 5, 5, 5, 5, 4, 1]

# virtual bond quantum numbers of third bond (multiplied by 2)
mpo_xxz.qbonds[3]

array([ 2,  0, -2,  0,  0])

# sparsity pattern consistency check
ptn.is_qsparse(mpo_xxz.a[3], (mpo_xxz.qbonds[3], mpo_xxz.qsite, -mpo_xxz.qsite, -mpo_xxz.qbonds[4]))

True

Note that we can effectively disable quantum numbers by setting them all to zero, since the above sparsity pattern rule is then always trivially satisfied.