1. PyTeNet Basics

This tutorial shows how to construct a matrix product state (MPS) with random entries, a Hamiltonian in matrix product operator (MPO) form, and compute the corresponding average energy.

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

# seed random number generator for reproducibility
rng = np.random.default_rng(42)

1.1. Matrix product states

Create a matrix product state with random tensor entries:

# physical dimension (of each lattice site);
# this could be spin 1, 0, -1 for example, or 0, 1, 2 particles per lattice site
d = 3

# virtual bond dimensions (note the leading and trailing 1)
b = [1, 4, 15, 13, 7, 1]

# set all physical and virtual bond quantum numbers to zero
# (effectively disabling quantum numbers)
qsite = np.zeros(d, dtype=int)
qbonds = [np.zeros(bi, dtype=int) for bi in b]

# now create the matrix product state
mps = ptn.MPS(qsite, qbonds, fill="random", rng=rng)

The actual tensors of the MPS are stored in the mps.a instance variable (list of NumPy arrays).

type(mps.a)

list

type(mps.a[0])

numpy.ndarray

# axis ordering convention of MPS tensors is
# ("left virtual bond", "physical axis", "right virtual bond")
mps.a[0].shape

(1, 3, 4)

# show random entries of leftmost tensor as illustration
mps.a[0]

array([[[ 0.06220011+0.01347846j, -0.21228587+0.23009715j,
          0.15318521+0.09542994j,  0.19199197-0.17540234j],
        [-0.39825339+0.07527094j, -0.26580628-0.19573109j,
          0.02609531+0.17931292j, -0.06455275-0.01019108j],
        [-0.00342952-0.03773487j, -0.17412686-0.13899416j,
          0.17950636+0.24955021j,  0.15876611-0.0315432j ]]])

# the number of lattice sites is (by definition) the length of mps.a;
# by construction same as length of virtual bond dimension list - 1
nsites = mps.nsites
(nsites, len(mps.a), len(b)-1)

(5, 5, 5)

Here are our virtual bond dimensions (same as b above):

mps.bond_dims

[1, 4, 15, 13, 7, 1]

The to_vector() function computes the vector representation of a matrix product state:

x = mps.to_vector()
# the length of x is equal to the physical Hilbert space dimension,
# i.e., d^nsites with nsites the number of lattice sites:
(x.shape, d**nsites)

((243,), 243)

# in general, the MPS is not normalized
np.linalg.norm(x)

np.float64(0.008359386283800499)

Call mps.orthonormalize() to left- or right-orthonormalize the mps.a tensors (via QR decompositions). This function returns the hitherto norm.

# returns previous norm
mps.orthonormalize(mode="left")

np.float64(0.008359386283800499)

# now the norm is 1, as expected
x = mps.to_vector()
np.linalg.norm(x)

np.float64(1.0000000000000002)

1.2. Hamiltonian as matrix product operator

PyTeNet supports some common quantum Hamiltonians in MPO form, see hamiltonian.py. The construction method there is very general and straightforward to adapt (see the hamiltonian_mpo.ipynb tutorial). Here we use the Bose-Hubbard model with the following Hamiltonian as an illustration.

\[H = -t \sum_{j=1}^{L-1} \left(b^{\dagger}_j b_{j+1} + \text{h.c.}\right) + \frac{1}{2} u \sum_{j=1}^L \hat{n}_j \left(\hat{n}_j - 1\right) - \mu \sum_{j=1}^L \hat{n}_j\]

# Hamiltonian parameters
t  =  1.0
u  =  4.0
mu = -0.5
# construct MPO (allowed local occupancies are 0, 1, ..., d - 1)
h_mpo = ptn.bose_hubbard_1d_mpo(nsites, d, t, u, mu)

As for matrix product states, the actual tensors of the MPO are stored in the h_mpo.a instance variable (list of NumPy arrays):

type(h_mpo.a)

list

type(h_mpo.a[0])

numpy.ndarray

# axis ordering convention of MPO tensors is
# ("left virtual bond", "physical output axis", "physical input axis", "right virtual bond")
(h_mpo.a[0].shape, h_mpo.a[1].shape, h_mpo.a[2].shape)

((1, 3, 3, 4), (4, 3, 3, 4), (4, 3, 3, 4))

# tensors are quite sparse
h_mpo.a[0]

array([[[[0.        , 1.        , 0.        , 0.        ],
         [0.        , 0.        , 1.        , 0.        ],
         [0.        , 0.        , 0.        , 0.        ]],

        [[1.        , 0.        , 0.        , 0.        ],
         [0.        , 1.        , 0.        , 0.5       ],
         [0.        , 0.        , 1.41421356, 0.        ]],

        [[0.        , 0.        , 0.        , 0.        ],
         [1.41421356, 0.        , 0.        , 0.        ],
         [0.        , 1.        , 0.        , 5.        ]]]])

The to_matrix() function computes the matrix representation of a matrix product operator:

h_mat = h_mpo.to_matrix()

type(h_mat)

numpy.ndarray

h_mat.shape

(243, 243)

1.3. Expectation values

Now let’s compute the average energy (expectation value) \(\langle\psi \vert H \vert \psi\rangle\) via the MPS formalism:

ptn.mpo_average(mps, h_mpo)

np.complex128(9.188269028617416+1.3183898417423734e-16j)

# should agree with matrix-vector representation
np.vdot(x, np.dot(h_mat, x))

np.complex128(9.188269028617416+0j)