I am looking for SVD for non-contiguous bipartitions of an MPS. However, even before that if I wish to compute entropy for 2:rest bipartition (i.e., 2:1,3,4,5…N), one way to do that is by applying a swap gate that exchanges the sites (1,2)–> (2,1) and then uses the usual SVD along 1:rest. First I constructed the swapped MPS phi from psi. And then used the SVD for phi along 1:rest. When I am trying to do so, it is showing some problem.
INPUT:
psi.position(1);
auto phi=MPS(N);
phi.ref(1) = psi(2);
phi.ref(2) = psi(1);

auto l1 = leftLinkIndex(phi,1);
auto s1 = siteIndex(phi,1);

auto [UU,S,V] = svd(phi(1),{l1,s1});

OUT PUT:
From line 568, file indexset.cc

In findIndex: more than one Index found, consider using findInds instead

In findIndex: more than one Index found, consider using findInds instead
Aborted (core dumped)

Any input will be highly appreciated.

Also, do you think finding the Schmidt coefficients for non-contiguous bipartitions SWAP operation is the only way out? I tried first constructing the reduced density matrix and then diagonalizing it but even for moderately large system size, I am having memory problems when I am using diagHermitian.

Hi, thanks for reposting this, actually because I had meant to get back to you before but it kind of fell off of my radar screen.

To start with, let’s talk about part (2) of your question. As long as the region “A” you are wanting the entanglement of is a small group of neighboring sites, then makind the reduced density matrix for these sites is the best choice, and doesn’t involve any swap gates etc. But this approach will scale exponentially in the number of sites making up region A so would only work for up to maybe 10 sites or so. How big of a region A do you want the entanglement for?

Regarding your code and the bug you’re getting, it’s hard to know just from reading the code what the problem is. Did you identify which line of the code is throwing the error?

My aim is to compute Schmidt values for a bipartitionA= [1, 3, 5,7…], B=[2, 4, 6,…] and the size of both the blocks are 10 (Totals system size is say, N=20). When I am trying to compute the reduced density matrix it is showing out of memory. That’s why I tried the swapping option.

But how efficient would be that one?

For Swapping I started with an empty MPS, say phi=MPS(N), and then did the following two things,
a) Swapped the elements:
psi.position(1);
auto phi=MPS(N);
phi.ref(1) = psi(2);
phi.ref(2) = psi(1);

b) Then also changed the site ordering for the swapped MPS
phi.ref(1) = phi(1) * delta(dag(sites(2)),sites(1));
phi.ref(2) = phi(2) * delta(dag(sites(1)),sites(2));

for the case of making the rdm directly by contracting MPS tensors (and tracing non-contiguous sites out that don’t appear in the rdm you are making) indeed you can run into very high costs. This approach scales as d^(2M) where d is the dimension of your site indices and M is the size of the region whose rdm you are trying to compute. So in practice it is limited to a small number of sites (M < 10).

swapping could indeed work for certain states, but for an MPS you have to remember that there are virtual or bond indices connecting the tensors, so you cannot just switch the two tensors and expect that to work. A better approach is to act with a swap operator on neighboring sites (similar to how you would act with a Trotter gate in a TEDB simulation). After merging and swapping the two sites you then SVD the tensor back apart to restore the MPS form.