Is it possible to use MetropolisHamiltonian sampler in dissipative models?

[GillianGrayson]

Hi!

As far as I understood MetropolisHamiltonian can be used, when Hilbert space has symmetries (docs).
When I tried to use MetropolisHamiltonian with DissipativeIsing1d model I got the error:

Traceback (most recent call last):
  File "/home/user/Work/dl/ndm/ising1d/4_samplers_test.py", line 62, in <module>
    vs = nk.vqs.MCMixedState(sa, ndm, n_samples=n_samples, n_samples_diag=n_samples_diag)
  File "/home/user/anaconda3/envs/py38/lib/python3.8/site-packages/netket/vqs/mc_mixed_state.py", line 156, in __init__
    self._diagonal = MCState(
  File "/home/user/anaconda3/envs/py38/lib/python3.8/site-packages/netket/vqs/mc_state.py", line 223, in __init__
    self.sampler = sampler
  File "/home/user/anaconda3/envs/py38/lib/python3.8/site-packages/netket/vqs/mc_state.py", line 267, in sampler
    self.sampler_state = self.sampler.init_state(
  File "/home/user/anaconda3/envs/py38/lib/python3.8/site-packages/netket/sampler/base.py", line 213, in init_state
    return sampler._init_state(wrap_afun(machine), parameters, key)
  File "/home/user/anaconda3/envs/py38/lib/python3.8/site-packages/netket/sampler/metropolis.py", line 274, in _init_state
    rule_state = sampler.rule.init_state(sampler, machine, params, key_rule)
  File "/home/user/anaconda3/envs/py38/lib/python3.8/site-packages/netket/sampler/rules/hamiltonian.py", line 50, in init_state
    raise ValueError(
ValueError: 
            The hilbert space of the sampler (Spin(s=1/2, N=6)) and the hilbert space
            of the operator (DoubledHilbert(Spin(s=1/2, N=6))) for HamiltonianRule must be the same.

The script with the error:

import netket as nk


# Model params
L = 6
gp = 0.3
Vp = 2.0

# Ansatz params
beta = 2
alpha = 2
n_samples = 5000
n_samples_diag = 2000
n_iter = 1000

# Create graph
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 = []

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))

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

# Neural quantum state model: Positive-Definite Neural Density Matrix using the ansatz from Torlai and Melko
ndm = nk.models.NDM(
    alpha=alpha,
    beta=beta,
)

# Metropolis Local Sampling
sa = nk.sampler.MetropolisHamiltonian(lind.hilbert, hamiltonian=lind)

# Variational state
vs = nk.vqs.MCMixedState(sa, ndm, n_samples=n_samples, n_samples_diag=n_samples_diag)
vs.init_parameters(nk.nn.initializers.normal(stddev=0.01))

The problem rises during initializing of sampler_diag with non-DoubledHibert space.

Motivation: in issue #850 I tried to consider a dissipative system with Many-Body Localization. This system has symmetry: half of the spins are in the ground state, half are in the excited state. I assumed that using the MetropolisHamiltonian, I could indicate the correct exchange graph in the doubled Choi space that corresponds to symmetry (by passing the corresponding matrix into the MetropolisHamiltonian).

[PhilipVinc]

Yes your intuition was correct (though you will have extremely bad performance, because MetropolisHamiltonian is slow and becomes even slower the more complex the operator is... and Lindlab operators are quite slow. They could use some love.
I think a better take would be to create the graph for MetropolisExchange manually)

Anyhow, that's an interesting error...
This line

  File "/home/user/anaconda3/envs/py38/lib/python3.8/site-packages/netket/vqs/mc_mixed_state.py", line 156, in __init__
    self._diagonal = MCState(

Is telling you that it is erroring during the initialisation of the internal MCState used to sample the diagonal of the density matrix.

For ease of use, when you don't specify the kwarg sampler_diag=... when constructing MCMixedState, NetKet tries to adapt the sampler used for the full density matrix to the diagonal by doing

        if sampler_diag is None:
            sampler_diag = sampler.replace(hilbert=hilbert_physical)

This works for simple samplers such as MetropolisLocal because they have no additional data, but utterly fails with MetropolisHamiltonian because it de-sync those two states.

The fix for you would be to also specify a sampler defined on the physical hilbert space for the diagonal such as MetropolisLocal or MetropolisHamiltonian(hamiltonian).

--

Instead, to avoid such errors in the future, we might try to default to MetropolisLocal.