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.

[1]:

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

[2]:

# 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:

[3]:

# 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)
D = [1, 4, 15, 13, 7, 1]

# set all physical and virtual bond quantum numbers to zero
# (effectively disabling quantum numbers)
qd = np.zeros(d, dtype=int)
qD = [np.zeros(Di, dtype=int) for Di in D]

# now create the matrix product state
mps = ptn.MPS(qd, qD, fill='random', rng=rng)

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

[4]:

type(mps.A)

[4]:

list

[5]:

type(mps.A[0])

[5]:

numpy.ndarray

[6]:

# the first dimension is the physical dimension,
# and the second and third dimension the left and right virtual bond dimensions
mps.A[0].shape

[6]:

(3, 1, 4)

[7]:

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

[7]:

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 ]]])

[8]:

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

[8]:

(5, 5, 5)

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

[9]:

mps.bond_dims

[9]:

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

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

[10]:

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

[10]:

((243,), 243)

[11]:

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

[11]:

0.00873358512801628

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

[12]:

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

[12]:

0.008733585128016278

[13]:

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

[13]:

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 illustration.

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

[14]:

# Hamiltonian parameters
t  =  1.0
U  =  4.0
mu = -0.5
# construct MPO (allowed local occupancies are 0, 1, ..., d - 1)
BH = ptn.bose_hubbard_mpo(d, L, t, U, mu)

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

[15]:

type(BH.A)

[15]:

list

[16]:

type(BH.A[0])

[16]:

numpy.ndarray

[17]:

# the first and second dimensions are the physical dimensions,
# and the third and fourth dimension the left and right virtual bond dimensions
(BH.A[0].shape, BH.A[1].shape, BH.A[2].shape)

[17]:

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

[18]:

# tensors are quite sparse
BH.A[0]

[18]:

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

        [[ 0.        +0.j,  0.        +0.j, -1.        +0.j,
           0.        +0.j]],

        [[ 0.        +0.j,  0.        +0.j,  0.        +0.j,
           0.        +0.j]]],


       [[[-1.        +0.j,  0.        +0.j,  0.        +0.j,
           0.        +0.j]],

        [[ 0.        +0.j,  1.        +0.j,  0.        +0.j,
           0.5       +0.j]],

        [[ 0.        +0.j,  0.        +0.j, -1.41421356+0.j,
           0.        +0.j]]],


       [[[ 0.        +0.j,  0.        +0.j,  0.        +0.j,
           0.        +0.j]],

        [[-1.41421356+0.j,  0.        +0.j,  0.        +0.j,
           0.        +0.j]],

        [[ 0.        +0.j,  1.        +0.j,  0.        +0.j,
           5.        +0.j]]]])

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

[19]:

BH_mat = BH.as_matrix()

[20]:

type(BH_mat)

[20]:

numpy.ndarray

[21]:

BH_mat.shape

[21]:

(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:

[22]:

ptn.operator_average(mps, BH)

[22]:

(9.249829654930819-4.0592529337857286e-16j)

[23]:

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

[23]:

(9.24982965493082+1.3877787807814457e-16j)

[ ]: