Monte Carlo Sampling Expectation Values of Obsevables

download script

The example below demonstrates how to use the *_basis_general methods Op_bra_ket(), representative() and get_amp(), which do not require computing all basis states, to do Monte-Carlo sampling in a symmetry-reduced subspace.

Depending on the problem of interest, sampling in the symmetry-reduced Hilbert space, as opposed to sampling in the full Hilbert space, may or may not be advantages (in terms of speed and efficiency).

Physics Setup

The expectation value of an operator \(H\) in a state \(\psi\) can be written as

\[\langle\psi|H|\psi\rangle = \sum_{s,s'} \psi_s^\ast H_{ss'}\psi_{s'} = \sum_s |\psi_s|^2 E_{s},\qquad E_s = \frac{1}{\psi_s}\sum_{s'}H_{ss'}\psi_{s'},\]

where \(\{|s\rangle\}_s\) can be any basis, in particular the Fock (\(z\)-) basis used in QuSpin.

The above expression suggests that one can use sampling methods, such as Monte Carlo, to estimate the expectation value of \(\langle\psi|H|\psi\rangle\) using the quantity \(E_s\) (sometimes referred to as the local energy) evaluated in samples drawn from the probability distribution \(p_s=|\psi_s|^2\). If we have: (i) a function \(s\mapsto\psi_s\) which computes the amplitudes for every spin configuration \(s\), and (ii) the matrix elements \(H_{ss'}\), then

\[\langle\psi|H|\psi\rangle \approx \frac{1}{N_\mathrm{MC}}\sum_{s\in\mathcal{S}} E_s,\]

where \(\mathcal{S}\) contains \(N_\mathrm{MC}\) spin configurations sampled from \(p_s=|\psi_s|^2\).

Since this procedure does not require the the state \(|\psi\rangle\) to be normalized, ideas along these lines allow to look for variational approximations to the wavefunction \(\psi_s\) in large system sizes, for instance with the help of restricted Boltzmann machines [arXiv:1606.02318].

The code below performs Monte-Carlo (MC) sampling in the symmetry-reduced Hilbert space. New states are proposed by swapping the spin values on random lattice sites; this update can take a representative state outside the symmetry-reduced basis, and we apply the basis.representative() function to project it back. To satisfy Detailed Balance in the MC acceptance/rejection condition, we use the function basis.get_amp() to compute the probability ratio in the full (symmetry-free) basis without actually constructing it.

In the example below, we assume that we already have a quantum state \(\psi_s\) in the Fock basis, and we sample the expectation value of an operator H using the *_basis_general methods Op_bra_ket(), representative(), and get_amp(). These methods do not require to compute the full basis, and thus allow to develop techniques for reaching system sizes beyond exact diagonalization (cf. the basis optional argument make_basis, as well as block_order which can deliver extra speed).

Script

download script

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
from __future__ import print_function, division
import sys,os
# line 4 and line 5 below are for development purposes and can be removed
qspin_path = os.path.join(os.getcwd(),"../../")
sys.path.insert(0,qspin_path)
###########################################################################
#                            example 11                                   #
#  In this script we demonstrate how to use QuSpin's methods of           #
#  the general_basis class which do not require explicit calculation      #
#  of the basis itself. Using the J1-J2 model on a square lattice, we     #
#  show how  to estimate the energy of a state using Monte-Carlo sampling.#
###########################################################################
from quspin.operators import hamiltonian
from quspin.basis import spin_basis_general
from quspin.operators._make_hamiltonian import _consolidate_static
import numpy as np
from scipy.special import comb
np.random.seed(1) #fixes seed of rng
from time import time # timing package
#
ti = time() # start timer
###### define model parameters ######
J1=1.0 # nn interaction
J2=0.5 # nnn interaction
Lx, Ly = 4, 4 # linear dimension of 2d lattice
N_2d = Lx*Ly # number of sites
#
###### setting up user-defined symmetry transformations for 2d lattice ######
sites = np.arange(N_2d) # site labels [0,1,2,....]
x = sites%Lx # x positions for sites
y = sites//Lx # y positions for sites
#
T_x = (x+1)%Lx + Lx*y # translation along x-direction
T_y = x +Lx*((y+1)%Ly) # translation along y-direction
#
T_a = (x+1)%Lx + Lx*((y+1)%Ly) # translation along anti-diagonal
T_d = (x-1)%Lx + Lx*((y+1)%Ly) # translation along diagonal
#
P_x = x + Lx*(Ly-y-1) # reflection about x-axis
P_y = (Lx-x-1) + Lx*y # reflection about y-axis
P_d = y + Lx*x # reflection about diagonal
#
Z   = -(sites+1) # spin inversion
#
###### setting up operator string for Hamiltonian matrix elements H_{ss'} ######
# setting up site-coupling lists for the J1-J2 model on a 2d square lattice
J1_list=[[J1,i,T_x[i]] for i in range(N_2d)] + [[J1,i,T_y[i]] for i in range(N_2d)]
J2_list=[[J2,i,T_d[i]] for i in range(N_2d)] + [[J2,i,T_a[i]] for i in range(N_2d)]
# setting up opstr list
static=[["xx",J1_list],["yy",J1_list],["zz",J1_list],  ["xx",J2_list],["yy",J2_list],["zz",J2_list]]
# convert static list to format which is easy to use with the basis_general.Op and basis_general.Op_bra_ket methods. 
static_formatted = _consolidate_static(static)
#
###### setting up basis object without computing the basis (make=False) ######
basis = spin_basis_general(N_2d, pauli=0, make_basis=False, 
			Nup=N_2d//2, 
			kxblock=(T_x,0), kyblock=(T_y,0),
			pxblock=(P_x,0), pyblock=(P_y,0), pdblock=(P_d,0),
			zblock=(Z,0),
			block_order=['zblock','pdblock','pyblock','pxblock','kyblock','kxblock'] # momentum symmetry comes last for speed
		)
print(basis) # examine basis: contains a single element because it is not calculated due to make_basis=False argument above.
print('basis is empty [note argument make_basis=False]')
#
###### define quantum state to compute the energy of using Monte-Carlo sampling ######
#
# auxiliary basis, only needed for probability_amplitude(); not needed in a proper variational ansatz.
aux_basis = spin_basis_general(N_2d, pauli=0, make_basis=True, 
			Nup=N_2d//2, 
			kxblock=(T_x,0), kyblock=(T_y,0),
			pxblock=(P_x,0), pyblock=(P_y,0), pdblock=(P_d,0),
			zblock=(Z,0),
			block_order=['zblock','pdblock','pyblock','pxblock','kyblock','kxblock'] # momentum symmetry comes last for speed
		)
# set quantum state to samplee from to be GS of H 
H = hamiltonian(static,[],basis=aux_basis,dtype=np.float64)
E, V = H.eigsh(k=2, which='SA') # need NOT be (but can be) normalized
psi=(V[:,0] + V[:,1])/np.sqrt(2)
#
##### define proposal function #####
#
def swap_bits(s,i,j):
	""" Swap bits i, j in integer s.

	Parameters
	-----------
	s: int
		spin configuration stored in bit representation.
	i: int
		lattice site position to be swapped with the corresponding one in j.
	j: int
		lattice site position to be swapped with the corresponding one in i.

	"""
	x = ( (s>>i)^(s>>j) ) & 1
	return s^( (x<<i)|(x<<j) )
#
##### define function to compute the amplitude `psi_s` for every spin configuration `s` #####
basis_state_inds_dict=dict()
for s in aux_basis.states:
	basis_state_inds_dict[s]=np.where(aux_basis.states==s)[0][0]
def probability_amplitude(s,psi):
	''' Computes probability amplitude `psi_s` of quantum state `psi` in z-basis state `s`.

	Parameters
	----------
	s: array_like(int)
		array of spin configurations [stored in their bit representation] to compute their local energies `E_s`.
	psi_s: array
		(unnormalized) probability amplitude values, corresponding to the states `s`. 

	'''
	return psi[[basis_state_inds_dict[ss] for ss in s]]
#
##### define function to compute local energy `E_s` #####
#
def compute_local_energy(s,psi_s,psi):
	"""Computes local energy E_s for a spin configuration s.

	Parameters
	----------
	s: array_like(int)
		array of spin configurations [stored in their bit representation] to compute their local energies `E_s`.
	psi_s: array
		(unnormalized) probability amplitude values, corresponding to the states `s` in the symmetry-reduced basis. 
	psi: array
		(unnormalized) array which encodes the mapping $s \\to \\psi_s$ (e.g. quantum state vector) in the symmetry-reduced basis.
		
	"""
	# preallocate variable
	E_s=np.zeros(s.shape,dtype=np.float64)
	#
	# to compute local energy `E_s` we need matrix elements `H_{ss'}` for the operator `H`.
	# These can be computed by looping overthe static list without constructing the operator matrix. 
	for opstr,indx,J in static_formatted:
		# for every state `s`, compute the state it connects to `s'`, and the corresponding matrix element `ME`
		ME,bras,kets = basis.Op_bra_ket(opstr,indx,J,np.float64,s,reduce_output=False)
		# performs sum over `s'`
		E_s+=ME * probability_amplitude(bras,psi)
	# normalize by `psi_s`
	E_s/=psi_s
	return E_s 
#
##### perform Monte Carlo sampling from `|psi_s|^2` ##### 
#
# draw random spin configuratio
s=[0 for _ in range(N_2d//2)] + [1 for _ in range(N_2d//2)]
np.random.shuffle(s)
s = np.array(s)
s=''.join([str(j) for j in s])
print('random initial state in bit representation:', s)
# transform state in bit representation
s=int(s,2)
print('same random initial state in integer representation:', s)
# compute representative of state `s` under basis symmetries (here only Z-symmetry)
s=basis.representative(s)
print('representative of random initial state in integer representation:', s)
# compute amplitude in state s
psi_s=probability_amplitude(s,psi)
#
psi_s_full_basis=np.copy(psi_s)
# overwrite the symmetry-reduced space psi_s_full_basis with its amplitude in the full basis
basis.get_amp(s,amps=psi_s_full_basis,mode='representative') # has the advantage that basis need not be made.
#
# define MC sampling parameters
equilibration_time=200
autocorrelation_time=N_2d
# number of MC sampling points
N_MC_points = 1000
#
##### run Markov chain MC #####
#
# compute all distinct site pairs to swap
#
# preallocate variables
E_s=np.zeros(N_MC_points,dtype=np.float64)
MC_sample=np.zeros(N_MC_points,dtype=basis.dtype)
#
j=0 # set MC chain counter
k=0 # set MC sample counter 
while k<N_MC_points:
	# propose new state t by swapping two random bits
	t=s
	while t==s: # repeat until a different state is reached
		# draw two random sites
		site_i=np.random.randint(0,N_2d)
		site_j=np.random.randint(0,N_2d)
		# swap bits in spin configurations
		t = swap_bits(s,site_i,site_j) # new state t corresponds to integer with swapped bites i and j
	#
	# compute representatives or proposed configurations to bring them back to symmetry sector
	t=basis.representative(t)
	psi_t=probability_amplitude(t, psi)
	# CAUTION: symmetries break detailed balance which we need to restore by using the amplitudes in the full basis.
	psi_t_full_basis=np.copy(psi_t)
	# overwrite the symmetry-reduced space psi_t_full_basis with its amplitude in the full basis
	basis.get_amp(t,amps=psi_t_full_basis,mode='representative') # has the advantage that basis need not be made.
	#
	### accept/reject new state
	#
	# use amplitudes psi_t_full_basis and psi_s_full_basis to restore detailed balance
	eps=np.random.uniform(0,1)
	if eps * np.abs(psi_s_full_basis)**2 <= np.abs(psi_t_full_basis)**2: 
		s=t
		psi_s=psi_t
		psi_s_full_basis=psi_t_full_basis
	#
	# wait for MC chain to quilibrate and collect uncorrelated samples
	if (j>equilibration_time) and (j%autocorrelation_time==0):
		# compute local energy
		print('computing local energy E_s for MC sample {0:d}'.format(k))
		E_s[k] = compute_local_energy(s,psi_s,psi)
		# update sample
		MC_sample[k]=s
		# update MC samples counter
		k+=1
	#
	j+=1 # update MC chain counter
#
##### compute MC-sampled average energy #####
# compute energy expectation and MC variance
E_mean=np.mean(E_s)
E_var_MC=np.std(E_s)/np.sqrt(N_MC_points)
#
# compute exact expectation value
E_exact=H.expt_value(psi/np.linalg.norm(psi))
#####   compute full basis   #####
# so far the functions representative(), get_amp(), and Op_bra_ket() did not require to compute the full basis
basis.make(Ns_block_est=16000)
print(basis) # after the basis is made, printing the basis returns the states
# compare results
print('mean energy: {0:.4f}, MC variance: {1:.4f}, exact energy {2:.4f}'.format(E_mean, E_var_MC, E_exact) )
print("simulation took {0:.4f} sec".format(time()-ti))