How to make MPO out of an MPS (outer product/projector)?

Question

Hello,

Is there an easy way to form a projector |psi><psi| from an MPS |psi>?
I have been looking through the documentation and some of files manually but could not find such function.

Thanks

Answer 1

Hi,

What would you like to use it for? I ask because in most use cases, you would not need to explicitly form the MPO (which would have a bond dimension of m^2 if formed from an MPS with bond dimension m). Instead, you can work directly with the tensors of the MPS. If you let us know your use case we can help you determine how to do that.

Cheers,
Matt

Comments

Hi Matt,

So here is what I want to achieve:
I have an MPS that suppose to represent (approximately) the wavefunction @@|\psi> = |g>\otimes|\phi_{1}> + |b>\otimes|\phi_2>@@, where @@|g>@@ and @@|b>@@ are orthogonal. I want to extract @@|\phi_1>@@.
What I mean by the tensor product, is that @@|g>@@ and @@|\phi_1>@@ describe different set of sites.
I have access to an MPS of @@|g>@@. To be more precise, what I really have is an MPS of the same number of sites as @@|\psi>@@, but @@|g>@@ describes a localized state, so I should be able to just throw away the irrelevant parts (although I dont know how to that either) .

To give some physical motivation in case it is not very clear, @@|\psi>@@ is the final state after a scattering experiment. @@|g>@@ or @@|b>@@ represent the (localized) bound states of the scatterer and @@|\phi_{1,2}>@@ are the output fields.

Would taking the overlap of a portion of @@|g>@@ with the corresponding portion of @@|\psi>@@ be the right way of doing this? In practice, @@|g>@@ is supported on 10-20 sites, and the whole @@|\psi>@@ is 300-400 sites.

This is probably not very difficult but I am new to both MPS and c++.

Hope it makes sense,
Thanks!

---

Thanks for the detailed explanation! If I understand correctly, it seems like you have two steps that you want to perform:

1. It sounds like you have some state @@|\psi'> = |g>\otimes|\phi'>@@, and you want to extract @@|g>@@. In that case, to figure out where the division between @@|g>@@ and @@|\phi'>@@ are, you can loop through the sites of your MPS and calculate the entanglement entropy as shown here (http://itensor.org/docs.cgi?page=formulas/entanglement_mps ). You can take the dividing line to be where the entanglement is minimal. For example, say that link 15 of your MPS @@|\psi'>@@ has minimal entanglement. In that case, you can take the first 15 tensors of @@|\psi'>@@ to be your state @@|g>@@ (after truncating the bond dimension of link 15 according to the singular values).

2. Once you save the tensors of @@|g>@@ from step 1. (perhaps saved as a `vector<ITensor>`), you can contract them over the corresponding sites of the MPS @@|\psi>@@, and the non-contracted MPS tensors over the rest of the sites in the system will be (approximately) the state @@|\phi_1>@@.

I am giving a fairly general picture for what you need to do, but I think it should give you an idea for how to get started. Please don't hesitate to ask for more specific guidance.

Cheers,
Matt

---

Hello Matt,

I am afraid I will need a bit more specific guidance, regarding the contraction part.

As far as I understand, to take an overlap I need to perform the same operations described here http://itensor.org/docs.cgi?page=tutorials/correlations except without sandwiching the additional operators. I am not sure if this is the correct approach for the "partial overlap" which I need.
Here is something I tried:

    auto ir = commonIndex(gs.A(190),gs.A(191),Link); //index to right of bra
    psi.Aref(190) *= dag(prime(gs.A(190),ir));
    Print(psi.A(190));
    for (int j = 191;j <= 209;++j){
        psi.Aref(190) *= psi.A(j);
        psi.Aref(190) *= dag(prime(gs.A(j),Link));
    }
    auto il = commonIndex(gs.A(209),gs.A(210),Link); //index to left of bra
    psi.Aref(190) *= dag(prime(gs.A(210),il));

where gs is the @@|psi'>@@ state. Now, I realized that the primes are unnecessary, because the Link indexes on psi and gs are different.
This also means that I have dangling indexes on the leftmost and rightmost tensors, which dont get contracted.  I dont know how to make ITensor contract them.

Supposing that I manage to contract these indexes, psi.A(190) will be a scalar, and I am not sure whether all this makes sense.
Sorry for the elementary questions!

---

Sorry for the late response, I didn't notice your comment until recently. I think you have the right idea. Here is an example code that may be useful for you to work from to get the details of what to do with all of the extra dangling bonds:

    // N is the total system size
    int N = 20;
    
    // Here are the sites that |psi'> and |psi> live on
    auto sites = SpinOne(N);
    
    // Here is |psi> = |g>|phi1> + |b>|phi2>
    // (This is a product state but yours will not be)
    auto psi = MPS(sites);
    
    // Here is |psi'> = |g>|phi'>
    // (This is a product state but yours will not be)
    auto psip = MPS(sites);
    
    // nmin is a stand in for where the minimal entanglement
    // cut of |psi'> is. You will need to calculate the entanglement
    // entropy at each cut to figure out where this bond is for your
    // example
    auto nmin = 10;
    
    // Move the orthogonality center to the minimal entanglement bond
    psip.position(nmin);
    
    // We truncate the bond nmin to 1 to divide |psi'> into
    // |g> and |phi'>
    psip.svdBond(nmin,psip.A(nmin)*psip.A(nmin+1),Fromright,{"Maxm=",1});
    
    // Get the new bond, so that we can get rid of it from the new
    // psip.A(nmin)
    auto lmin = commonIndex(psip.A(nmin),psip.A(nmin+1));
    psip.Aref(nmin) *= setElt(lmin(1));
    
    // Now that the bond is truncated, the site 1...nmin represent
    // the approximation for |g>.
    // Now we can contract <g|psi>:
    auto gpsi = dag(psip.A(1))*psi.A(1);
    for(int i = 2; i<=nmin; i++)
      gpsi *= dag(psip.A(i))*psi.A(i);
    
    // Now incorporate the overlap of <g|psi> into the rest of |psi>
    // The sites nmin+1...N of psi now represent an approximation
    // to |phi1>
    psi.Aref(nmin+1) *= gpsi;

Let me know if you have any questions about any of the steps!

---

Hi Matt,

Thanks a lot for the detailed solution! I have come up with my own, just before you posted the comment.

To eliminate the extra Link index I used the delta tensor, like this:
auto spec = svd(tbond,U,S,V,{"Maxm",1});//,"ShowEigs",true});
edgeR = U*S;
edgeR *= delta(commonIndex(S,V));

I suppose that should be equivalent to what you used
    psip.Aref(nmin) *= setElt(lmin(1));

Are those equivalent?

Here is my full function to cut an MPS by half and store the two halves as new MPS:

    void cutMPS(MPS& psi,MPS& psi_cut_L,MPS& psi_cut_R, int bond){
    //right cut begins with bond+1, left cut includes up to bond
    psi.position(bond);
    auto tbond = psi.A(bond)*psi.A(bond+1);
    ITensor U = psi.A(bond);
    ITensor S,V;
    auto spec = svd(tbond,U,S,V,{"Maxm",1});//,"ShowEigs",true});
    int N_L = psi_cut_L.N();
    int N_R = psi_cut_R.N();
    ITensor edgeR,edgeL,temp;
    if ( bond != N_L ) {
        throw std::invalid_argument( "bond and psi_cut.L dont match" );
    }
    edgeR = U*S;
    edgeR *= delta(commonIndex(S,V)); //remove the unneccesery Link index of size one
    psi_cut_L.setA(N_L,edgeR*delta(findtype(edgeR,Site),findtype(psi_cut_L.A(N_L),Site)));
    for (int j = 1;j < N_L; ++j){
        temp = psi.A(j);
        psi_cut_L.setA(j,temp*delta(findtype(temp,Site),findtype(psi_cut_L.A(j),Site)));
    }

    edgeL = V;
    edgeL *= delta(commonIndex(V,S)); //remove the unneccesery Link index of size one
    psi_cut_R.setA(1,edgeL*delta(findtype(edgeL,Site),findtype(psi_cut_R.A(1),Site)));
    for (int j = 2;j <= N_R; ++j){
        temp = psi.A(j+bond);
        psi_cut_R.setA(j,temp*delta(findtype(temp,Site),findtype(psi_cut_R.A(j),Site)));
    }
    
    psi_cut_L.position(1);
    normalize(psi_cut_L);
    psi_cut_R.position(1);
    normalize(psi_cut_R);
    }