Ordering of exact diagonalization basis

[bart]

I am trying to match the results of TeNPy with a few other codes and am trying to figure out the TeNPy convention for the ordering of the basis in ExactDiag. For example, if we have SpinHalfFermionSites with 3 sites and (2, 2) electrons, then in total the Hilbert space has dimension 9 (N and Sz are conserved). Usually, I would expect that the Hilbert space is ordered in ascending binary numbers. For example, I would expect the statevector in this case, to be returned in the following order (using SpinHalfFermionSite notation):
[3, 3, 0]
[3, 1, 2]
[1, 3, 2]
[3, 2, 1]
[3, 0, 3]
[1, 2, 3]
[2, 3, 1]
[2, 1, 3]
[0, 3, 3]
However, it is actually returned in this order (which I can deduce by converting all the product states one-by-one):
[2, 3, 1]
[2, 1, 3]
[0, 3, 3]
[3, 2, 1]
[3, 0, 3]
[3, 3, 0]
[3, 1, 2]
[1, 2, 3]
[1, 3, 2]
Does anyone know what is the TeNPy convention for the ordering? Is it sorted in some binary representation?

Sample code to recover the basis ordering by brute force:

Python: Select all

basis_ordering = []
for i, state in enumerate(product_states):
    prod_mps = MPS.from_product_state(ED.model.lat.mps_sites(), state)     
    prod_statevector = list(ED.mps_to_full(prod_mps).to_ndarray())
    idx = prod_statevector.index(1)
    basis_ordering.append(idx)

[Johannes]

Hi Bart,
it's essentially defined by the following line in ExactDiag.__init__()

Python: Select all

self._pipe = npc.LegPipe(legs, qconj=1, sort=(not sparse), bunch=(not sparse))

Unless you use the `sort=False` keyword, the pipe internally sorts by charge values, more precises such that it's `qmap` is lexiographically sorted - that's why you get that strange order in this case.
The task of the Pipe is essentially to 1) store the orginal legs it was contracted from and 2) the permutations done, such that the combine_legs can be un-done with the `split_legs`. Which order the pipe chooses is considered an implementation detail (it's a somewhat arbitrary convention how we do that), but we guarantee that the `split_legs` can undo the permutation done.

When you call psi_full.to_ndarray() where psi_full is e.g. the ground state returend by exact_diag.groundstate(), you will get a full single leg with d^L entries. The straightforward solution is to call split_legs() first, which undoes whatever permutation was done for the combine_legs, such that
psi_full.split_legs().to_ndarray() is simply a shape [d]*L array instead - then the mapping of indices should be very clear.
This is a bit inefficient (since it fills on all the zeros that are guaranteed by to be zero by charge conservation) - if you need to avoid that as well, you instead need to directly get the data block out of the np_conserved array, and undo the permutation yourself.
Technically, the permutation data is stored in LegPipe._perm and ._qmap, but it's probably hard to read that out and disentangle the code, if you don't know the structure. The more bruteforce but simple way would be to just call LegPipe.map_incoming_inds() for all combinations of indices you need, and get the permutation from that. (Given that we are currently rewriting the np_conserved for version 2.0, I would recommend to not invest too much time into understanding the LegPipe's internal details)

[bart]

Ah I see, it is lexiographically sorted by `qmap`. Thank you, I will try out these suggestions to restore the conventional ordering of the statevectors. As long as one of these methods is faster than explicitly checking/converting each product state MPS, then it should be an improvement.

[Jakob]

Hi bart, just FYI that this will be resolved in the next version of tenpy, developed at https://github.com/tenpy/cyten.
There, we consider this permutation in ``Tensor.to_numpy()`` such that it behaves as expected, i.e. such that ``.combine_legs(...).to_numpy()`` is equivalent to ``.to_numpy().reshape(...)``, *without* any permutations

[bart]

I see, that's great to hear. In that case, I will not dig too deep into this. Thank you!

[bart]

Hi Johannes and Jakob,

Thank you again for your answers and suggestions to this question. I have had a chance to think about this now and I think that these suggestions do not solve the problem. I will try to explain more clearly with a fully self-contained minimal example. Here is my example script...

Python: Select all

from tenpy.models.hubbard import FermiHubbardChain
from tenpy.algorithms.exact_diag import ExactDiag
from tenpy.networks.mps import MPS

M = FermiHubbardChain({"L": 2})

product_states = [[3, 0], [1, 2], [2, 1], [0, 3]]
basis_map = []
for state in product_states:
    prod_mps = MPS.from_product_state(M.lat.mps_sites(), state)
    ED = ExactDiag(M, charge_sector=[2, 0])
    prod_statevector = list(ED.mps_to_full(prod_mps).to_ndarray())
    print(f"{state} => {prod_statevector}")
    basis_map.append(prod_statevector.index(1))
print(f"basis map = {basis_map}")

...and here is its output...

Code: Select all

[3, 0] => [0.0, 0.0, 1.0, 0.0]
[1, 2] => [0.0, 0.0, 0.0, 1.0]
[2, 1] => [1.0, 0.0, 0.0, 0.0]
[0, 3] => [0.0, 1.0, 0.0, 0.0]
basis map = [2, 3, 0, 1]

In the above example, I consider a L=2 FermiHubbardChain and I fix the symmetry sector to (N, Sz)=(2, 0). There are 4 possible product states in this symmetry sector [[3, 0], [1, 2], [2, 1], [0, 3]]. I would now like to deduce the basis ordering for a statevector given from the mps_to_full function. To this end, I can convert the individual MPS product states into statevectors (total Hilbert space dimension of 4). I can do this one-by-one and see which entry each product state corresponds to, which finally gives me a basis map. What I would like to do is find a way to deduce basis_map a priori, without having to explicitly convert each product state MPS into a statevector.

Previously suggested methods:

• I do not notice any difference applying .split_legs() to the statevector from ED.mps_to_full(prod_mps).
• As far as I understand, ED._pipe._perm will give me the permutation of the complete Hilbert space (size 16) to the one which is sorted by symmetry sectors. Indeed, if I print ED._mask, I can tell which of these entries are in the (2, 0) sector. However, this does not help me find basis_map, as far as I can tell.
• As far as I understand, ED._pipe.map_incoming_flat(state) is simply another way to derive ED._pipe._perm.

Do you know how I can find the basis_map in my example script, without explicitly converting every possible product state? In other words, what is the default ordering of the exact diagonalization basis in TeNPy, within a symmetry sector?