Inconsistent estimation from steady state sampling

[qetuo2]

Hi,

Recently I have experimented the Examples/DissipativeIsing/ising1d.py code with netket.3.3.2.post2. I try to verify the sampling for LdaggerL and other observables in the following code:

import netket as nk
import numpy as np

# 1D Lattice
L = 4
gp = 0.3
Vp = 2.0

g = nk.graph.Hypercube(length=L, n_dim=1, pbc=False)

# Hilbert space of spins on the graph
hi = nk.hilbert.Spin(s=1 / 2, N=g.n_nodes)

# The hamiltonian
ha = nk.operator.LocalOperator(hi)

# List of dissipative jump operators
j_ops = []

# Observables
obs_sx = nk.operator.LocalOperator(hi)
obs_sy = nk.operator.LocalOperator(hi, dtype=complex)
obs_sz = nk.operator.LocalOperator(hi)

for i in range(L):
    ha += (gp / 2.0) * nk.operator.spin.sigmax(hi, i)
    ha += (
        (Vp / 4.0)
        * nk.operator.spin.sigmaz(hi, i)
        * nk.operator.spin.sigmaz(hi, (i + 1) % L)
    )
    # sigma_{-} dissipation on every site
    j_ops.append(nk.operator.spin.sigmam(hi, i))
    obs_sx += nk.operator.spin.sigmax(hi, i)
    obs_sy += nk.operator.spin.sigmay(hi, i)
    obs_sz += nk.operator.spin.sigmaz(hi, i)


#  Create the liouvillian
lind = nk.operator.LocalLiouvillian(ha, j_ops)

ma = nk.models.NDM(
    beta=3,
)

# Metropolis Local Sampling
sa = nk.sampler.MetropolisLocal(lind.hilbert)

# Optimizer
op = nk.optimizer.Sgd(0.01)
sr = nk.optimizer.SR(diag_shift=0.01)

vs = nk.vqs.MCMixedState(sa, ma, n_samples=2000, n_samples_diag=500)
vs.init_parameters(nk.nn.initializers.normal(stddev=0.01))

ss = nk.SteadyState(lind, op, variational_state=vs, preconditioner=sr)

obs = {"Sx": obs_sx, "Sy": obs_sy, "Sz": obs_sz}

out = ss.run(n_iter=200, out="test", obs=obs)

rho = vs.to_matrix()
ll_exact = float(np.sum(np.abs(lind.to_sparse()@rho.reshape(-1))**2))
sz_exact = np.trace(obs_sz.to_sparse()@rho)
print('ll_exact', ll_exact)
print('sz_exact', sz_exact)
n_batch = 100
vs.n_samples=10000
vs.n_samples_diag=10000

def compute_mean(vs, n_batch, ops):
  ll_list = []
  for _ in range(n_batch):
    vs.sample()
    tmp = vs.expect(ops)
    ll_list.append(float(tmp.mean))
  ll_list = np.array(ll_list)
  ll_avg = np.mean(ll_list)
  ll_std = np.std(ll_list) / np.sqrt(n_batch)
  return ll_avg, ll_std, ll_list

def compute_mean_diag(vs, n_batch, ops):
  ll_list = []
  for _ in range(n_batch):
    vs.diagonal.sample()
    tmp = vs.expect(ops)
    ll_list.append(complex(tmp.mean))
  ll_list = np.array(ll_list)
  ll_avg = np.mean(ll_list)
  ll_std = np.std(ll_list) / np.sqrt(n_batch)
  return ll_avg, ll_std, ll_list

ll_avg, ll_std, ll_list = compute_mean(vs, n_batch, lind.H@lind)
sz_avg, sz_std, sz_list = compute_mean_diag(vs, n_batch, obs_sz)
print('ll_sampling', ll_avg, ll_std)
print('sz_sampling', sz_avg, sz_std)

I have the following outcomes:

ll_exact 0.1279338892236631
sz_exact (-3.0008300177896667+0j)
ll_sampling 0.21384986400238362 0.0012018498004683997
sz_sampling (-3.002326+0j) 0.0022461819249562118

While the sz observable estimation is correct, the LdaggerL estimation seems to be quite off. May I know what is the correct way to estimate LdaggerL in this case? Thanks!

[qetuo2]

In addition, the 200 steps optimization shows the following:
200/200 [00:38<00:00, 5.25it/s, LdagL=0.195 ± 0.014 [σ²=0.205]
which is close to the resampling ll_sampling value, while is different from ll_exact.

[PhilipVinc]

Are you sure you are normalising correctly the estimation of LdagL?
The estimator is defined in my PRL paper and it is normalised by the L2 norm of the density matrix, not the trace squared as you are doing in ll_exact = float(np.sum(np.abs(lind.to_sparse()@rho.reshape(-1))**2))

I also suggest you to use current master (which will soon be released as v3.4) because I changed the algorithm used to compute LdagL to one that is a bit more stable and faster.