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.
[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)
[ ]: