This notebook is for generating the figures for my thesis
Sam Greydanus. May 2017. MIT License
Code inspired by https://github.com/simple-dmrg/simple-dmrg/
% matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import copy
import time
from scipy.sparse import kron, identity
from scipy.sparse.linalg import eigsh # Lanczos routine from ARPACK
np.random.seed(seed=123) # for reproducibility
class FreeSite():
def __init__(self):
self.length = 1 # length
self.ops = ops = {}
# build operator dictionary
ops["H"] = np.zeros((2,2)) # local Hamiltonian np.random.randn(2,2)
ops["Sz"] = np.array([[0.5, 0], [0, -0.5]]) # z spin (S^z) operator
ops["Sp"] = np.array([[0.0, 1.0], [0.0, 0.0]]) # raising (S^+) operator
def get_dim(self):
return list(self.ops.values())[0].shape[0] # all ops should have same dimensionality
def enlarge(self, site):
'''Enlarge block by a single site'''
D1, H1, Sz1, Sp1 = self.get_dim(), self.ops['H'], self.ops['Sz'], self.ops['Sp'] # this block
D2, H2, Sz2, Sp2 = site.get_dim(), site.ops['H'], site.ops['Sz'], site.ops['Sp'] # another block (ie free site)
enlarged = copy.deepcopy(self)
enlarged.length += site.length
ops = enlarged.ops
ops['H'] = kron(H1, identity(D2)) + kron(identity(D1), H2) + self.interaction_H(site)
ops['Sz'] = kron(identity(D1), Sz2)
ops['Sp'] = kron(identity(D1), Sp2)
return enlarged
def interaction_H(self, site):
'''Given another block, returns two-site term in the
Hamiltonain that joins the two sites.'''
Sz1, Sp1 = self.ops["Sz"], self.ops["Sp"] # this block
Sz2, Sp2 = site.ops["Sz"], site.ops["Sp"] # another block
J = 1.; Jz = 1.
# J = np.random.rand(); Jz = np.random.rand()
join_Sp = (J/2)*(kron(Sp1, Sp2.conjugate().transpose()) + kron(Sp1.conjugate().transpose(), Sp2))
join_Sz = Jz*kron(Sz1, Sz2)
return (join_Sp + join_Sz)
def rotate_ops(self, transformation_matrix):
# rotate and truncate each operator.
new_ops = {}
for name, op in self.ops.items():
new_ops[name] = self.rotate_and_truncate(op, transformation_matrix)
self.ops = new_ops
@staticmethod
def rotate_and_truncate(S, O):
'''Transforms the operator to a new (possibly truncated) basis'''
return O.conjugate().transpose().dot(S.dot(O)) # eqn 7 in arXiv:cond-mat/0603842v2
def step(sys, env, free_site, m):
'''One DMRG step keeping maximum of m states in new basis'''
# enlarge blocks
sys_enl = sys.enlarge(free_site)
env_enl = sys_enl if env is sys else env.enlarge(free_site)
sys_ops, env_ops = sys_enl.ops, env_enl.ops
sys_enl_m, env_enl_m = sys_enl.get_dim(), env_enl.get_dim()
# find superblock hamiltonian
H_E1 = kron(sys_ops["H"], identity(env_enl_m))
H_E2 = kron(identity(sys_enl_m), env_ops["H"])
H_SS = sys_enl.interaction_H(env_enl)
H_super = H_E1 + H_E2 + H_SS
# find the superblock ground state. "SA"= smallest amplitude e'value
(energy,), psi0 = eigsh(H_super, k=1, which="SA")
# construct reduced density matrix. we want (sys, env) -> (row, column)
psi0 = psi0.reshape([sys_enl_m, -1])
rho = np.dot(psi0, psi0.conjugate().transpose())
# diagonalize reduced density matrix, sort e'vecs by e'vals
evals, evecs = np.linalg.eigh(rho)
possible_eigenstates = []
for eval, evec in zip(evals, evecs.transpose()):
possible_eigenstates.append((eval, evec))
possible_eigenstates.sort(reverse=True, key=lambda x: x[0]) # largest e'values first
# build transformation matrix from m overall most significant e'vecs
m = min(len(possible_eigenstates), m)
transformation_matrix = np.zeros((sys_enl_m, m))
for i, (eval, evec) in enumerate(possible_eigenstates[:m]):
transformation_matrix[:, i] = evec
error = 1 - sum([x[0] for x in possible_eigenstates[:m]])
sys_enl.rotate_ops(transformation_matrix)
return sys_enl, energy, error
def graphic(sys, env, sys_label='l'):
assert sys_label in ('l', 'r'), 'Label must be "l" or "r"'
graphic = ('=' * sys.length) + '**' + ('-' * env.length)
return graphic if sys_label is "l" else graphic[::-1]
# global plot settings
fs = [8,5]
cs = 4
L=16
m_default = 16
cl = [-.5,.5]
m_vals = [4,6,8,10,12]
assert L % 2 == 0, "L must be even"
infin_results = []
for m in m_vals:
# build system using infinite system
free_site = FreeSite()
block = copy.deepcopy(free_site) # block starts out as one free site
infin_blocksize = []
infin_energies = []
infin_errors = []
while 2 * block.length < L:
# one DMRG step, save block to mem
# print(graphic(block, block))
block, energy, error = step(block, block, free_site, m=m)
# print("E/L =", energy / (block.length * 2))
infin_blocksize.append(block.length * 2)
infin_energies.append(energy / (block.length * 2))
infin_errors.append(error)
infin_results.append(((m, (infin_blocksize, infin_energies, infin_errors, block.ops['H']))))
print('Finished infinite algorithm for for m={}'.format(m))
Finished infinite algorithm for for m=4 Finished infinite algorithm for for m=6 Finished infinite algorithm for for m=8 Finished infinite algorithm for for m=10 Finished infinite algorithm for for m=12
free_site = FreeSite()
block = free_site
for i in range(2,17):
block = block.enlarge(free_site)
(actual_e,), psi0 = eigsh(block.ops['H'], k=1, which="SA")
actual_e /= 16
actual_e_error = 0
f1 = plt.figure(figsize=fs)
plt.title("Infinite DMRG $E_{0}$ estimates")
plt.xlabel("Block size ($L$)")
plt.ylabel("Ground state energy ($E_{0}$)")
k = 0
ee = 1
labels = []
for i in range(len(m_vals)):
x = infin_results[i][1][0][k:]
y = infin_results[i][1][1][k:]
yerr = infin_results[i][1][2][k:]
yerr = np.cumsum(yerr)
plt.errorbar(x, y, yerr=yerr, errorevery=ee, capsize=cs)
labels.append(r'$m = %i$' % (m_vals[i]))
y = [actual_e]*len(x)
plt.errorbar(x, y, color='k')
labels.append('exact')
plt.legend(labels, ncol=2, loc='upper right')
plt.show() ; f1.savefig('./figures/infin-full-e0.pdf', bbox_inches='tight')
f2 = plt.figure(figsize=fs)
plt.title("Infinite DMRG $E_{0}$ estimates (near convergence)")
plt.xlabel("Block size ($L$)")
plt.ylabel("Ground state energy ($E_{0}$)")
# colormap = plt.cm.gist_ncar
# plt.gca().set_color_cycle([colormap(i) for i in np.linspace(0, 0.9, len(m_vals)-2)])
k = int((L-10)/2)
labels = []
for i in range(len(m_vals)):
x = infin_results[i][1][0][k:]
y = infin_results[i][1][1][k:]
yerr = infin_results[i][1][2][k:]
yerr = np.cumsum(yerr)
plt.errorbar(x, y, yerr=yerr, errorevery=ee, capsize=cs)
labels.append(r'$m = %i$' % (m_vals[i]))
y = [actual_e]*len(x)
plt.errorbar(x, y, color='k')
labels.append('exact')
plt.legend(labels, ncol=2, loc='upper right')
plt.show() ; f2.savefig("./figures/infin-convergence-e0.pdf", bbox_inches='tight')
m_warmup=12
m_sweep_list=[12, 14, 16]
block_mem = {} # memory for Block objects
assert L % 2 == 0, "L must be even"
# build system using infinite system
free_site = FreeSite()
block = copy.deepcopy(free_site) # block starts out as one free site
block_mem["l", block.length] = block
block_mem["r", block.length] = block
print("building finite system to size L")
while 2 * block.length < L:
# one DMRG step, save block to mem
block, energy, error = step(block, block, free_site, m=m_warmup)
block_mem["l", block.length] = block
block_mem["r", block.length] = block
fin_results = []
# now that we've built the system, sweep back and forth
sys_label, env_label = "l", "r"
sys = block; del block # rename block
print("sweeping over system with finite algorithm")
i = 0
for m in m_sweep_list:
print("\t using m={}".format(m))
fin_step = [] ; fin_energies = [] ; fin_errors = []
while True:
i += 1
env = block_mem[env_label, L - sys.length - 2] # load env block from memory
if env.length == 1: # reverse direction
sys, env = env, sys
sys_label, env_label = env_label, sys_label
# one DMRG step
# print(graphic(sys, env, sys_label))
sys, energy, error = step(sys, env, free_site, m=m)
fin_step.append(i)
fin_energies.append(energy/L)
fin_errors.append(error/L)
# print("E/L: {} \t error: {}".format(energy/L, error))
# save new block
block_mem[sys_label, sys.length] = sys
# full sweep?
if sys_label == "l" and 2 * sys.length == L: break
fin_results.append(((m, (fin_step.copy(), fin_energies.copy(), fin_errors.copy(), sys.ops['H'].copy()))))
building finite system to size L sweeping over system with finite algorithm using m=12 using m=14 using m=16
f3 = plt.figure(figsize=fs)
plt.title("Finite DMRG $E_{0}$ estimates")
plt.xlabel("step #")
plt.ylabel("Ground state energy ($E_{0}$)")
k = 0
ee = 3
labels = []
yerr_prev = 0
for i in range(len(m_sweep_list)):
x = fin_results[i][1][0][k:]
y = fin_results[i][1][1][k:]
yerr = fin_results[i][1][2][k:]
yerr = np.cumsum(yerr) + yerr_prev
plt.errorbar(x, y, yerr=yerr, errorevery=ee, capsize=cs)
labels.append(r'$m = %i$' % (m_sweep_list[i]))
yerr_prev = yerr[-1]
x = range(fin_results[i][1][0][-1])
y = [actual_e]*len(x)
plt.errorbar(x, y, color='k')
labels.append('exact')
plt.legend(labels, ncol=1, loc='upper right')
plt.show() ; f3.savefig("./figures/fin-e0.pdf", bbox_inches='tight')
m = m_default
block_mem = {} # memory for Block objects
assert L % 2 == 0, "L must be even"
# build system using infinite system
free_site = FreeSite()
block = copy.deepcopy(free_site) # block starts out as one free site
block_mem["l", block.length] = block
block_mem["r", block.length] = block
while 2 * block.length < L:
# one DMRG step, save block to mem
block, energy, error = step(block, block, free_site, m=m)
block_mem["l", block.length] = block
block_mem["r", block.length] = block
fin_results = []
infin_block = block
# now that we've built the system, sweep back and forth
m_sweep_list=[m, m, m]
sys_label, env_label = "l", "r"
sys = block; del block # rename block
print("sweeping over system with finite algorithm")
for m in m_sweep_list:
print("\t using m={}".format(m))
while True:
env = block_mem[env_label, L - sys.length - 2] # load env block from memory
if env.length == 1: # reverse direction
sys, env = env, sys
sys_label, env_label = env_label, sys_label
# one DMRG step
sys, energy, error = step(sys, env, free_site, m=m)
# save new block
block_mem[sys_label, sys.length] = sys
# full sweep?
if sys_label == "l" and 2 * sys.length == L: break
fin_block = sys
sweeping over system with finite algorithm using m=16 using m=16 using m=16
f4 = plt.figure(figsize=[10,7])
plt.subplot(231) ; plt.title("$H$") ; plt.ylabel("Infinite DMRG") ; plt.imshow(infin_block.ops['H'])
plt.clim(*cl)
plt.subplot(232) ; plt.title("$S^z$") ; plt.imshow(infin_block.ops['Sz'])
plt.clim(*cl)
plt.subplot(233) ; plt.title("$S^+$") ; plt.imshow(infin_block.ops['Sp'])
plt.clim(*cl)
plt.subplot(234) ; plt.title("$H$") ; plt.ylabel("Finite DMRG") ; plt.imshow(fin_block.ops['H'])
plt.clim(*cl)
plt.subplot(235) ; plt.title("$S^z$") ; plt.imshow(fin_block.ops['Sz'])
plt.clim(*cl)
plt.subplot(236) ; plt.title("$S^+$") ; im = plt.imshow(fin_block.ops['Sp'])
plt.clim(*cl)
cax = f4.add_axes([0.95, 0.15, 0.03, 0.7])
cb = f4.colorbar(im, cax=cax, orientation='vertical')
cb.set_label('color scale')
plt.show() ; f4.savefig("./figures/compare-ops.pdf", bbox_inches='tight')
f5 = plt.figure(figsize=fs)
plt.subplot(121) ; plt.title("Infinite DMRG Hamiltonian") ; plt.imshow(infin_block.ops['H'])
plt.clim(*cl)
plt.subplot(122) ; plt.title("Finite DMRG Hamiltonian") ; im = plt.imshow(fin_block.ops['H'])
plt.clim(*cl)
cax = f5.add_axes([0.95, 0.15, 0.03, 0.7])
cb = f5.colorbar(im, cax=cax, orientation='vertical')
cb.set_label('color scale')
plt.show() ; f5.savefig("./figures/compare-Hs.pdf", bbox_inches='tight')
runs = []
for i in range(5):
runtimes = [] ; x = []
free_site = FreeSite()
block = free_site
for j in range(2,21):
block = block.enlarge(free_site)
start = time.time()
(energy,), psi0 = eigsh(block.ops['H'], k=1, which="SA")
runtimes.append(time.time()-start) ; x.append(j)
runs.append(runtimes)
print("{}x".format(i+1))
1x 2x 3x 4x 5x
f6 = plt.figure(figsize=[6,5])
plt.title("Scaling of $\mathbf{eigsh()}$ computation time", fontsize=16)
plt.xlabel("Number of sites ($N$)", fontsize=14)
plt.ylabel("Compute time (seconds)", fontsize=14)
s = 10
runtimes = np.vstack(runs)
ymeans = np.mean(runtimes, axis=0)[s:]
ystds = np.std(runtimes, axis=0)[s:]
x_ = x[s:]
plt.errorbar(x_, ymeans, yerr=ystds, capsize=cs)
plt.show() ; f6.savefig('./figures/eigsh-runtimes.pdf', bbox_inches='tight')
m = 20
runs = []
for i in range(5):
runtimes = []
L_vals = [10,20,30,40,50,60,70,80,90,100]
for L in L_vals:
# build system using infinite system
start = time.time()
free_site = FreeSite()
block = copy.deepcopy(free_site) # block starts out as one free site
while 2 * block.length < L:
block, energy, error = step(block, block, free_site, m=m)
runtimes.append(time.time()-start)
print('\tfinished infinite algorithm m={}, L={}'.format(m, L))
runs.append(runtimes)
print("{}x".format(i+1))
finished infinite algorithm m=20, L=10 finished infinite algorithm m=20, L=20 finished infinite algorithm m=20, L=30 finished infinite algorithm m=20, L=40 finished infinite algorithm m=20, L=50 finished infinite algorithm m=20, L=60 finished infinite algorithm m=20, L=70 finished infinite algorithm m=20, L=80 finished infinite algorithm m=20, L=90 finished infinite algorithm m=20, L=100 1x finished infinite algorithm m=20, L=10 finished infinite algorithm m=20, L=20 finished infinite algorithm m=20, L=30 finished infinite algorithm m=20, L=40 finished infinite algorithm m=20, L=50 finished infinite algorithm m=20, L=60 finished infinite algorithm m=20, L=70 finished infinite algorithm m=20, L=80 finished infinite algorithm m=20, L=90 finished infinite algorithm m=20, L=100 2x finished infinite algorithm m=20, L=10 finished infinite algorithm m=20, L=20 finished infinite algorithm m=20, L=30 finished infinite algorithm m=20, L=40 finished infinite algorithm m=20, L=50 finished infinite algorithm m=20, L=60 finished infinite algorithm m=20, L=70 finished infinite algorithm m=20, L=80 finished infinite algorithm m=20, L=90 finished infinite algorithm m=20, L=100 3x finished infinite algorithm m=20, L=10 finished infinite algorithm m=20, L=20 finished infinite algorithm m=20, L=30 finished infinite algorithm m=20, L=40 finished infinite algorithm m=20, L=50 finished infinite algorithm m=20, L=60 finished infinite algorithm m=20, L=70 finished infinite algorithm m=20, L=80 finished infinite algorithm m=20, L=90 finished infinite algorithm m=20, L=100 4x finished infinite algorithm m=20, L=10 finished infinite algorithm m=20, L=20 finished infinite algorithm m=20, L=30 finished infinite algorithm m=20, L=40 finished infinite algorithm m=20, L=50 finished infinite algorithm m=20, L=60 finished infinite algorithm m=20, L=70 finished infinite algorithm m=20, L=80 finished infinite algorithm m=20, L=90 finished infinite algorithm m=20, L=100 5x
f7 = plt.figure(figsize=[6,5])
plt.title("Scaling of DMRG computation time ($m=20$)", fontsize=16)
plt.xlabel("Number of sites ($N$)", fontsize=14)
plt.ylabel("Compute time (seconds)", fontsize=14)
s = 0
runtimes = np.vstack(runs)
ymeans = np.mean(runtimes, axis=0)[s:]
ystds = np.std(runtimes, axis=0)[s:]
x = L_vals[s:]
plt.errorbar(x, ymeans, yerr=ystds, capsize=cs)
plt.show() ; f7.savefig('./figures/dmrg-runtimes.pdf', bbox_inches='tight')