2. Quantum Hamiltonians as matrix product operators

In this tutorial, we will first introduce so-called 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 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 as_matrix() function forms the outer product, i.e., the matrix representation of the operator chain (without the leading and trailing identities):

opchain.as_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 a MPO version of a sum of operator chains? A useful mental picture are 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 besides the tracks, each displaying a local operator (like the Pauli matrices in the example above). The operator chain (train) collects the markers which 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: its markers solely display the local identity operation. Operator trains share this identity track before and after the sites which they act on.

To convert the operator chains to an MPO, we use an intermediate OpGraph data structure, which 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)
L = 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], L, 0)
# internal consistency check
print(opgraph.is_consistent(verbose=True))

True

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

# construct MPO representation of the operator graph
mpo = ptn.MPO.from_opgraph(qd, 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.as_matrix(opmap)),
    np.identity(2**(L - (opc.istart + opc.length))))
        for opc in [opchain, opchain2, opchain3])
# difference should be zero
np.linalg.norm(mpo.as_matrix() - opchains_mat)

0.0

2.3. Constructing quantum Hamiltonians

The module hamiltonian.py uses the mechanism to construct quantum Hamiltonians as MPOs. As 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(L, 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., perform 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_mpo(L, 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, like preserving the total spin or particle number of a physical system, for example. 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 qd.) Likewise, each virtual bond is also associated with a quantum number (qD in PyTeNet). Now every MPO tensor obeys a sparsity pattern dictated by the quantum numbers: The sum of first physical and left virtual bond quantum number of each non-zero tensor entry is equal to the sum of second physical and right virtual bond quantum number. (Here “first” and “second” refers to the row and column dimension 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 appears naturally when constructing the MPO representation of a Hamiltonian, and quantum numbers provide the means for actually recognizing and understanding it. In practical terms, quantum numbers allow to optimize common operations like 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.) By the same functionality as above, the quantum numbers are copied from the operator chains to the graph and finally into the virtual bond quantum numbers of the MPO, endowing 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_mpo(L, J, D, h)

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

array([ 1, -1])

# virtual bond dimensions
mpo_xxz.bond_dims

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

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

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

# sparsity pattern consistency check
ptn.is_qsparse(mpo_xxz.A[3], [mpo_xxz.qd, -mpo_xxz.qd, mpo_xxz.qD[3], -mpo_xxz.qD[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.