netket: Errors when using 'exact' representation of density matrices

A worrisome email I received:

Since the 3.0 update some observables in my calculations switched signs (the MC estimates).

When calculating expectation values of an operator that is not symmetric in a steady state solution I get a different sign when calculating
rho = nk.exact.steady_state(lind)
tr(rho@op.to_dense())
0.15
ss.estimate(op)
-0.15
When investigating this, I didn’t find any fault in the MC nk code, but one possible explanation could be that operator.to_dense() returns a transposed matrix or is this your convention that operators are applied to the left?

He also showed me a case that I further reduced a bit. here is my explanation:

In nette, for reasons unbeknown to me, we have the convention of labelling the Hilbert space from down to up, from low to high. So the local Hilbert space of Pauli matrices is [down, up], as testified by

>>> nk.hilbert.Spin(0.5).all_states()
array([[-1.],
       [ 1.]])

However, if I now create the sigma plus operator I get:

>>> nk.operator.spin.sigmap(nk.hilbert.Spin(0.5), 0).operators
[array([[0., 1.],
       [0., 0.]])]

And clearly in this conventiont this is the sigma_minus operator. The same issue happens with non symmetric matrices like sigma_y.

The error is on my part, as I had taken that code from some older code of mine which I guess assumed opposite ordering. Interestingly enough, boson operators are correct…

This is consistent with the fact that I long have seen that sigma_y had the opposite sign in some calculations compared to my old article, but after investigating I thought the error was on my old article…

Still, while the fix is clear, I have to double check the Impact on the ordering of the lindblad jump operators because I think I might have implemented a transposition there otherwise I don’t understand how it can work…

Note: THIS WILL BREAK CODES

About this issue

Original URL
State: closed
Created 3 years ago
Comments: 16 (11 by maintainers)

Most upvoted comments

@jmellak I’m confident I solved the issue, but could you double check that PR #617 fixes the issue?

PhilipVinc on Apr 12, 2021

Upon more debugging it seems to me that:

nk.exact.steady_state returns the transpose of the steady_state. I say that by comparing the steady state computed by it with the steady state given by qutip. (I’m sure of that)
nk.MCMixedState.to_matrix() returns the transpose of the state (I’m 95% sure of this)
nk.MCMixedState.estimate(obs) works correctly

The fix would be fixing the two methods above…

@jmellak I think your ‘exact’ code is wrong. I double checked my previous results in this paper, which is the same model you have in your example, and sigma_y must be positive, so what is wrong here is the exact part, not the monte-carlo…

This would explain the inconsistency you see.

PhilipVinc on Apr 12, 2021

But again, I stress that this is not a bug, one can well call the basis states [-3435,+231212], define the usual Pauli matrices, and still have consistent results. (modulo knowing that the quantum numbers used are not the eigenstates of Sigma_z). The issue I think here arises only in the density matrix, to be investigated further …

gcarleo on Apr 10, 2021

Unfortunate naming on my part - op is an operator representing an observable (an asymmetric one, e.g. Sigma_y).

I attached a minimal demo script based on the dissipative ising1d example, showing what exactly I mean. testOperators.py.zip The output I get:
testOperators_output.txt Can you reproduce this? Takes a minute to run.

Here

print(f'RBM result for the observable: {ss.estimate(obs)}')
print(f'Exact result from ED:          {np.trace(rho_0@obs.to_dense())}')

Gives the output:

… RBM result for the observable: 0.063+0.021j ± 0.052 [σ²=1.360, R̂=0.9978] Exact result from ED: (-0.08141487086313647+1.5265566588595902e-16j) …

In the second half of the script I recreated the netket mechanism for calculating rho.expect(operator) to better understand what is going on. Using this code and removing the Monte-Carlo part I get the exact same result as the ED result rho_0 but with a different sign.

… Reproduction of expectation value using rho_0: (0.08141487086313648-1.1102230246251565e-16j) Result from exact diagonalization rho_0 : (-0.08141487086313647+1.5265566588595902e-16j) …

But if I transpose all operator matrices, which would change the sign of the observable, the lindblad operator is not correct. I hope you can reproduce my results with this and maybe tell me where I messed up?

jmellak on Apr 9, 2021

Mmmm… no, it seems that I am exaggerating. The historical definition of the Pauli matrices is correct, and the only bug is me transposing them in operator.spin.

It’s simply that we are using an opposite convention for the states, where [-1] corresponds to the up state which is… weird and confusing

PhilipVinc on Apr 9, 2021