Proliferating matrix elements in LocalOperators

[attila-i-szabo]

I have been trying to compute the expectation values of LocalOperators of the form Sx @ Sz, where Sx and Sz are linear combinations of Pauli X and Z operators acting on all N lattice sites. Sx connects each element of the Hilbert space to N others, and Sz is diagonal, both of these work out with the output of get_conn. However, their product ends up having N^2 non-unique connections:

N = 32
hi = nk.hilbert.Spin(0.5, N)
LC = np.ones(N)
Sx = nk.operator.LocalOperator(
    hi, LC[:, None, None] * np.array([[0, 1], [1, 0]]), [[i] for i in range(N)]
)
Sz = nk.operator.LocalOperator(
    hi, LC[:, None, None] * np.array([[1, 0], [0, -1]]), [[i] for i in range(N)]
)
x = np.sign(np.random.normal(size=N)) # a random vector of +/- 1
a, b = Sx.get_conn(x)
a.shape
>>> (32, 32)
a, b = Sz.get_conn(x)
a.shape
>>> (1,32)
a,b = (Sz@Sx).get_conn(x)
a.shape
>>> (1024, 32)
np.unique(a, axis=0).shape
>>> (32, 32)

and the same for Sx@Sz.

Is there a way to force LocalOperator to aggregate these matrix elements? I imagine something like this must be in place for the purely diagonal operators already...

[PhilipVinc]

Hmm.. is this a bug?

In [34]: (Sz@Sx).operators[1]
Out[34]:
array([[ 0.,  1.,  1.,  0.],
       [ 1.,  0.,  0., -1.],
       [ 1.,  0.,  0., -1.],
       [ 0., -1., -1., -0.]])

LocalOperator is built under the assumption of Local (or few body) terms.
PauliStringshould be better at those 'row sparse' but not local operators.

As for conflating matrix terms... What matrix terms would you conflate?

[attila-i-szabo]

I've had a quick look at LocalOperator.get_conn_flattened and I see that it explicitly treats diagonal matrix elements separately. One could rewrite it roughly as follows:

• One outer loop over the quantum numbers
  • find connections and matrix elements naively (don't separate diagonal and offdiagonal)
  • perform np.unique on the connections and use the indices returned by it to aggregate the matrix elements
  • find the length of the section at the end
• Keep the connections and the matrix elements in a list of arrays until the end, then assemble them using np.concatenate at the very end

Would this be slower than the current code under numba?

[attila-i-szabo]

LocalOperator is built under the assumption of Local (or few body) terms. PauliStringshould be better at those 'row sparse' but not local operators.

These operators are 2-body, and there is many of them, so I don't think PauliString would help.
Okay, PauliString could be a workaround. However, it would be a hassle for a more complex operator of this kind (e.g. Sx @ Sz @ Sx) to work out all the coefficients, and I guess it doesn't work with ones that contain S+ or S-?

As for conflating matrix terms... What matrix terms would you conflate?

Sx_i @ Sz_j produce the same connections for every value of j, because every Sz is diagonal. These could be aggregated

Hmm.. is this a bug?

In [34]: (Sz@Sx).operators[1]
Out[34]:
array([[ 0.,  1.,  1.,  0.],
       [ 1.,  0.,  0., -1.],
       [ 1.,  0.,  0., -1.],
       [ 0., -1., -1., -0.]])

What do you mean?

[PhilipVinc]

What do you mean?

I would expect that Sx@Sy yields mostly 4x4 matrices with only 1 nonzero entry per row, not two...

However, it would be a hassle for a more complex operator of this kind (e.g. Sx @ Sz @ Sx)

PauliString supports products of PauliStrings, so if you write Sx and Sz as PauliStrings (simply pass it in ["XIII", "IXIII", "IIXII", ...]) you can compose them as you want, like you are dong now for LocalOperators... Or is there something I',m missing?

and I guess it doesn't work with ones that contain S+ or S-?

You can still write those operators as "Sx +1j Sy" but indeed, it won't be efficient.
I would like to add those kind of operators to a "GeneralisedPauliOperator" but this was never implemented.

I see that it explicitly treats diagonal matrix elements

This code is quite old, but the reason we do that is to avoid duplicating the diagonal term on every term of the sum as we do for the off diagonal terms.
I'm not sure that getting rid of this would be a good idea because, at least for the diagonal, it can be precomputed easily (and we also have the diagonal offset to consider)

perform np.unique on the connections and use the indices returned by it to aggregate the matrix elements

That's a good idea actually!
I don't know how it would perform. But in general I think that we should strive to reduce number of connected entries. In any case this code has lower runtime complexity than evaluating the NN ansatz.

I would be curious to see this implemented.

[attila-i-szabo]

I would expect that Sx@Sy yields mostly 4x4 matrices with only 1 nonzero entry per row, not two...

This is fine, it just combines Sx_i @ Sz_j and Sx_j @ Sz_i into the same matrix.

PauliString supports products of PauliStrings, so if you write Sx and Sz as PauliStrings (simply pass it in ["XIII", "IXIII", "IIXII", ...]) you can compose them as you want, like you are dong now for LocalOperators... Or is there something I',m missing?

Yep, that's true.

You can still write those operators as "Sx +1j Sy" but indeed, it won't be efficient.
I would like to add those kind of operators to a "GeneralisedPauliOperator" but this was never implemented.

Actually, it seems that the zero matrix elements of "Sx + 1j Sy" get eliminated, which is nice. I think this solves my immediate issue.

perform np.unique on the connections and use the indices returned by it to aggregate the matrix elements

That's a good idea actually! I don't know how it would perform. But in general I think that we should strive to reduce number of connected entries. In any case this code has lower runtime complexity than evaluating the NN ansatz.

I would be curious to see this implemented.

I could give this a go when I have time. I imagine constructing the connected matrix elements is the cheapest part of calculating expectation values, so even if this is a bit slower, it might be better for code readability.

[gcarleo]

Is this a bug then?

[PhilipVinc]

No

[gcarleo]

It is true though that PauliString is much more efficient at removing duplicate connected elements, I wonder if for LocalOperators that can be converted to PauliString it makes sense to do the conversion by default?