Monte Carlo Sampling Expectation Values of Obsevables
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
1###########################################################################
2# example 11 #
3# In this script we demonstrate how to use QuSpin's methods of #
4# the general_basis class which do not require explicit calculation #
5# of the basis itself. Using the J1-J2 model on a square lattice, we #
6# show how to estimate the energy of a state using Monte-Carlo sampling.#
7###########################################################################
8from quspin.operators import hamiltonian
9from quspin.basis import spin_basis_general
10from quspin.operators._make_hamiltonian import _consolidate_static
11import numpy as np
12from scipy.special import comb
13
14np.random.seed(1) # fixes seed of rng
15from time import time # timing package
16
17#
18ti = time() # start timer
19###### define model parameters ######
20J1 = 1.0 # nn interaction
21J2 = 0.5 # nnn interaction
22Lx, Ly = 4, 4 # linear dimension of 2d lattice
23N_2d = Lx * Ly # number of sites
24#
25###### setting up user-defined symmetry transformations for 2d lattice ######
26sites = np.arange(N_2d) # site labels [0,1,2,....]
27x = sites % Lx # x positions for sites
28y = sites // Lx # y positions for sites
29#
30T_x = (x + 1) % Lx + Lx * y # translation along x-direction
31T_y = x + Lx * ((y + 1) % Ly) # translation along y-direction
32#
33T_a = (x + 1) % Lx + Lx * ((y + 1) % Ly) # translation along anti-diagonal
34T_d = (x - 1) % Lx + Lx * ((y + 1) % Ly) # translation along diagonal
35#
36P_x = x + Lx * (Ly - y - 1) # reflection about x-axis
37P_y = (Lx - x - 1) + Lx * y # reflection about y-axis
38P_d = y + Lx * x # reflection about diagonal
39#
40Z = -(sites + 1) # spin inversion
41#
42###### setting up operator string for Hamiltonian matrix elements H_{ss'} ######
43# setting up site-coupling lists for the J1-J2 model on a 2d square lattice
44J1_list = [[J1, i, T_x[i]] for i in range(N_2d)] + [
45 [J1, i, T_y[i]] for i in range(N_2d)
46]
47J2_list = [[J2, i, T_d[i]] for i in range(N_2d)] + [
48 [J2, i, T_a[i]] for i in range(N_2d)
49]
50# setting up opstr list
51static = [
52 ["xx", J1_list],
53 ["yy", J1_list],
54 ["zz", J1_list],
55 ["xx", J2_list],
56 ["yy", J2_list],
57 ["zz", J2_list],
58]
59# convert static list to format which is easy to use with the basis_general.Op and basis_general.Op_bra_ket methods.
60static_formatted = _consolidate_static(static)
61#
62###### setting up basis object without computing the basis (make=False) ######
63basis = spin_basis_general(
64 N_2d,
65 pauli=0,
66 make_basis=False,
67 Nup=N_2d // 2,
68 kxblock=(T_x, 0),
69 kyblock=(T_y, 0),
70 pxblock=(P_x, 0),
71 pyblock=(P_y, 0),
72 pdblock=(P_d, 0),
73 zblock=(Z, 0),
74 block_order=[
75 "zblock",
76 "pdblock",
77 "pyblock",
78 "pxblock",
79 "kyblock",
80 "kxblock",
81 ], # momentum symmetry comes last for speed
82)
83print(
84 basis
85) # examine basis: contains a single element because it is not calculated due to make_basis=False argument above.
86print("basis is empty [note argument make_basis=False]")
87#
88###### define quantum state to compute the energy of using Monte-Carlo sampling ######
89#
90# auxiliary basis, only needed for probability_amplitude(); not needed in a proper variational ansatz.
91aux_basis = spin_basis_general(
92 N_2d,
93 pauli=0,
94 make_basis=True,
95 Nup=N_2d // 2,
96 kxblock=(T_x, 0),
97 kyblock=(T_y, 0),
98 pxblock=(P_x, 0),
99 pyblock=(P_y, 0),
100 pdblock=(P_d, 0),
101 zblock=(Z, 0),
102 block_order=[
103 "zblock",
104 "pdblock",
105 "pyblock",
106 "pxblock",
107 "kyblock",
108 "kxblock",
109 ], # momentum symmetry comes last for speed
110)
111# set quantum state to samplee from to be GS of H
112H = hamiltonian(static, [], basis=aux_basis, dtype=np.float64)
113E, V = H.eigsh(k=2, which="SA") # need NOT be (but can be) normalized
114psi = (V[:, 0] + V[:, 1]) / np.sqrt(2)
115
116
117#
118##### define proposal function #####
119#
120def swap_bits(s, i, j):
121 """Swap bits i, j in integer s.
122
123 Parameters
124 -----------
125 s: int
126 spin configuration stored in bit representation.
127 i: int
128 lattice site position to be swapped with the corresponding one in j.
129 j: int
130 lattice site position to be swapped with the corresponding one in i.
131
132 """
133 x = ((s >> i) ^ (s >> j)) & 1
134 return s ^ ((x << i) | (x << j))
135
136
137#
138##### define function to compute the amplitude `psi_s` for every spin configuration `s` #####
139basis_state_inds_dict = dict()
140for s in aux_basis.states:
141 basis_state_inds_dict[s] = np.where(aux_basis.states == s)[0][0]
142
143
144def probability_amplitude(s, psi):
145 """Computes probability amplitude `psi_s` of quantum state `psi` in z-basis state `s`.
146
147 Parameters
148 ----------
149 s: array_like(int)
150 array of spin configurations [stored in their bit representation] to compute their local energies `E_s`.
151 psi_s: array
152 (unnormalized) probability amplitude values, corresponding to the states `s`.
153
154 """
155 return psi[[basis_state_inds_dict[ss] for ss in s]]
156
157
158#
159##### define function to compute local energy `E_s` #####
160#
161def compute_local_energy(s, psi_s, psi):
162 """Computes local energy E_s for a spin configuration s.
163
164 Parameters
165 ----------
166 s: array_like(int)
167 array of spin configurations [stored in their bit representation] to compute their local energies `E_s`.
168 psi_s: array
169 (unnormalized) probability amplitude values, corresponding to the states `s` in the symmetry-reduced basis.
170 psi: array
171 (unnormalized) array which encodes the mapping $s \\to \\psi_s$ (e.g. quantum state vector) in the symmetry-reduced basis.
172
173 """
174 # preallocate variable
175 E_s = np.zeros(s.shape, dtype=np.float64)
176 #
177 # to compute local energy `E_s` we need matrix elements `H_{ss'}` for the operator `H`.
178 # These can be computed by looping overthe static list without constructing the operator matrix.
179 for opstr, indx, J in static_formatted:
180 # for every state `s`, compute the state it connects to `s'`, and the corresponding matrix element `ME`
181 ME, bras, kets = basis.Op_bra_ket(
182 opstr, indx, J, np.float64, s, reduce_output=False
183 )
184 # performs sum over `s'`
185 E_s += ME * probability_amplitude(bras, psi)
186 # normalize by `psi_s`
187 E_s /= psi_s[0]
188 return E_s
189#
190##### perform Monte Carlo sampling from `|psi_s|^2` #####
191#
192# draw random spin configuratio
193s = [0 for _ in range(N_2d // 2)] + [1 for _ in range(N_2d // 2)]
194np.random.shuffle(s)
195s = np.array(s)
196s = "".join([str(j) for j in s])
197print("random initial state in bit representation:", s)
198# transform state in bit representation
199s = int(s, 2)
200print("same random initial state in integer representation:", s)
201# compute representative of state `s` under basis symmetries (here only Z-symmetry)
202s = basis.representative(s)
203print("representative of random initial state in integer representation:", s)
204# compute amplitude in state s
205psi_s = probability_amplitude(s, psi)
206#
207psi_s_full_basis = np.copy(psi_s)
208# overwrite the symmetry-reduced space psi_s_full_basis with its amplitude in the full basis
209basis.get_amp(
210 s, amps=psi_s_full_basis, mode="representative"
211) # has the advantage that basis need not be made.
212#
213# define MC sampling parameters
214equilibration_time = 200
215autocorrelation_time = N_2d
216# number of MC sampling points
217N_MC_points = 1000
218#
219##### run Markov chain MC #####
220#
221# compute all distinct site pairs to swap
222#
223# preallocate variables
224E_s = np.zeros(N_MC_points, dtype=np.float64)
225MC_sample = np.zeros(N_MC_points, dtype=basis.dtype)
226#
227j = 0 # set MC chain counter
228k = 0 # set MC sample counter
229while k < N_MC_points:
230 # propose new state t by swapping two random bits
231 t = s
232 while t == s: # repeat until a different state is reached
233 # draw two random sites
234 site_i = np.random.randint(0, N_2d)
235 site_j = np.random.randint(0, N_2d)
236 # swap bits in spin configurations
237 t = swap_bits(
238 s, site_i, site_j
239 ) # new state t corresponds to integer with swapped bites i and j
240 #
241 # compute representatives or proposed configurations to bring them back to symmetry sector
242 t = basis.representative(t)
243 psi_t = probability_amplitude(t, psi)
244 # CAUTION: symmetries break detailed balance which we need to restore by using the amplitudes in the full basis.
245 psi_t_full_basis = np.copy(psi_t)
246 # overwrite the symmetry-reduced space psi_t_full_basis with its amplitude in the full basis
247 basis.get_amp(
248 t, amps=psi_t_full_basis, mode="representative"
249 ) # has the advantage that basis need not be made.
250 #
251 ### accept/reject new state
252 #
253 # use amplitudes psi_t_full_basis and psi_s_full_basis to restore detailed balance
254 eps = np.random.uniform(0, 1)
255 if eps * np.abs(psi_s_full_basis) ** 2 <= np.abs(psi_t_full_basis) ** 2:
256 s = t
257 psi_s = psi_t
258 psi_s_full_basis = psi_t_full_basis
259 #
260 # wait for MC chain to quilibrate and collect uncorrelated samples
261 if (j > equilibration_time) and (j % autocorrelation_time == 0):
262 # compute local energy
263 print("computing local energy E_s for MC sample {0:d}".format(k))
264 E_s[k] = compute_local_energy(s, psi_s, psi)[0]
265 # update sample
266 MC_sample[k] = s[0]
267 # update MC samples counter
268 k += 1
269 #
270 j += 1 # update MC chain counter
271#
272##### compute MC-sampled average energy #####
273# compute energy expectation and MC variance
274E_mean = np.mean(E_s)
275E_var_MC = np.std(E_s) / np.sqrt(N_MC_points)
276#
277# compute exact expectation value
278E_exact = H.expt_value(psi / np.linalg.norm(psi))
279##### compute full basis #####
280# so far the functions representative(), get_amp(), and Op_bra_ket() did not require to compute the full basis
281basis.make(Ns_block_est=16000)
282print(basis) # after the basis is made, printing the basis returns the states
283# compare results
284print(
285 "mean energy: {0:.4f}, MC variance: {1:.4f}, exact energy {2:.4f}".format(
286 E_mean, E_var_MC, E_exact
287 )
288)
289print("simulation took {0:.4f} sec".format(time() - ti))