[julia] Out of memory problem

Question

Dear all,

I have obtained the ground state MPS wavefunction, Now the problem is: when calculating the correlators, I encounter the "out of memory error" and the program get killed.

Note this problem is similar to this one, but in my case the correct code for calculating the correlator has been implemented, but the problem remains...

I guess there may be some way to save the intermediate data into disk? Can someone give me more details on how to do it, or give me some links about this method.

Comments

Hi JW, could you provide some more information? It sounds like your DMRG calculation succeeded without running out of memory, correct? But then later you were calling the function correlation_matrix and then you got an out of memory error at that point? Is that right?

One thing to please check is: what version of ITensors.jl are you using?

Best,
Miles

---

Hi miles, thank you again for the help! I am using Julia version 1.7.1, and ITensors v0.2.12

Yes, the DMRG calculation is finished. And then I use my own function to do some correlator measurement. Here I do not use the built-in correlation_matrix because my site type is customized, and correlation_matrix can not work

In case you want to check my measurement function. I attach my code below:



   
                sites = [
                    if mod(n,3) == 2
                        Index([QN()=>2]; tags="S=1/2, Site, n=$n")
                    else
                        siteind("Fermion", n; conserve_qns=true)
                    end
                    for n = 1:3*ns
                        ]
                state = [
                    if mod(n,3) == 2
                        "Up"
                    else
                        (mod(n,fl*3) == 1 || mod(n,fl*3) == 0) ? "Occ" : "Emp"
                    end
                    for n = 1:3*ns
                        ]
            
                psii = randomMPS(sites, state)
            
                # Sz-Sz
                function szz(psi::MPS)
                    szz = zeros(ns,ns)
                    ts = ITensor(1.)
                    i1 = 0
                    for i in 1:3*ns
                        if mod(i,3) == 2
                            i1 += 1
                            j1 = i1
                            orthogonalize!(psi, i)
                            szz[i1,i1] = scalar(ts*psi[i]*op(szo,sites,i)*prime(op(szo,sites,i))*dag(prime(prime(psi[i],"Site"),"Site")))
                            lind = linkind(psi,i)
                            tsi = ts*psi[i]*op(szo,sites,i)*dag(prime(prime(psi[i],"Site"),lind))
                            for j = i+1:3*ns
                                tsi *= psi[j]
                                if mod(j,3) == 2 # spin site only
                                    j1 += 1
                                    lind = linkind(psi,j-1)
                                    szz[i1,j1] = scalar(tsi*op(szo,sites,j)*dag(prime(prime(psi[j],"Site"),lind)))
                                    szz[j1,i1] = conj(szz[i1,j1])
                                end
                                tsi *= dag(prime(psi[j],"Link")) # no JW string
                            end
                        end
                    end
                    return szz
                end
                
                # Sz
                function sz(psi::MPS)
                    sz = zeros(ns)
                    ts = ITensor(1.)
                    i1 = 0
                    for i in 1:3*ns
                        if mod(i,3) == 2
                            i1 += 1
                            orthogonalize!(psi,i)
                            sz[i1] = scalar(ts*psi[i]*op(szo,sites,i)*dag(prime(psi[i],"Site")))
                        end
                    end
                    return sz
                end

                # Nf
                function n(psi::MPS)
                    n = zeros(2*ns)
                    ts = ITensor(1.)
                    i1 = 0
                    for i in 1:3*ns
                        if mod(i,3) != 2
                            i1 += 1
                            orthogonalize!(psi, i)
                            n[i1] = scalar(ts*psi[i]*op("N",sites,i)*dag(prime(psi[i],"Site")))
                        end
                    end
                    return n
                end

                # Nf-Nf
                function nn(psi::MPS)
                    nn = zeros(2*ns,2*ns)
                    ts = ITensor(1.)
                    i1 = 0
                    for i in 1:3*ns
                        if mod(i,3) != 2
                            i1 += 1
                            j1 = i1
                            orthogonalize!(psi,i)
                            nn[i1,i1] = scalar(ts*psi[i]*op("N",sites,i)*prime(op("N",sites,i))*dag(prime(prime(psi[i],"Site"),"Site")))
                            lind = linkind(psi,i)
                            tsi = ts*psi[i]*op("N",sites,i)*dag(prime(prime(psi[i],"Site"),lind))
                            for j = i+1:3*ns
                                tsi *= psi[j]
                                if mod(j,3) != 2
                                    j1 += 1
                                    lind = linkind(psi,j-1)
                                    nn[i1,j1] = scalar(tsi*op("N",sites,j)*dag(prime(prime(psi[j],"Site"),lind)))
                                    nn[j1,i1] = conj(nn[i1,j1])
                                    tsi *= dag(prime(psi[j],"Link"))
                                else
                                    tsi *= dag(prime(psi[j],"Link"))
                                end
                            end
                        end
                    end
                    return nn
                end

                # Cdag-C
                function ff(psi::MPS)
                    ff = zeros(2*ns,2*ns)
                    ts = ITensor(1.)
                    i1 = 0
                    for i in 1:3*ns
                        if mod(i,3) != 2
                            i1 += 1
                            j1 = i1
                            orthogonalize!(psi,i)
                            ff[i1,i1] = scalar(ts*psi[i]*op("N",sites,i)*dag(prime(psi[i],"Site")))
                            lind = linkind(psi,i)
                            tsi = ts*psi[i]*op("Cdag",sites,i)*dag(prime(prime(psi[i],"Site"),lind))
                            for j in i+1:3*ns
                                tsi *= psi[j]
                                if mod(j,3) != 2
                                    j1 += 1
                                    lind = linkind(psi, j-1)
                                    ff[i1,j1] = scalar(tsi*op("C",sites,j)*dag(prime(prime(psi[j],"Site"),lind)))
                                    ff[j1,i1] = conj(ff[i1,j1])
                                    # with JW string
                                    tsi *= op("F",sites,j)*dag(prime(psi[j]))
                                else
                                    # no JW string
                                    tsi *= dag(prime(psi[j],"Link"))
                                end
                            end
                        end
                    end
                    return ff
                end


                # Fz
                # fzi = Nli - Nri / 2
                function fz(psi::MPS)
                    fz = zeros(ns)
                    ts = ITensor(1.)
                    i1 = 0
                    for i in 1:3*ns
                        if mod(i,3) == 1
                            i1 += 1
                            orthogonalize!(psi,i)
                            nup = scalar(ts*psi[i]*op("N",sites,i)*dag(prime(psi[i],"Site")))
                            orthogonalize!(psi,i+2)
                            ndn = scalar(ts*psi[i+2]*op("N",sites,i+2)*dag(prime(psi[i+2],"Site")))
                            fz[i1] = (nup-ndn)/2
                        end
                    end
                    return fz
                end
                
                # Fx & Fy
                # sxi = fli fri + fri fli / 2
                # syi = -i fli fri + i fri fli / 2
                function fxfy(psi::MPS)
                    fx = zeros(ns)
                    fy = zeros(ns)
                    ts = ITensor(1.)
                    i1 = 0
                    for i in 1:3*ns
                        if mod(i,3) == 1
                            i1 += 1
                            orthogonalize!(psi,i)
                            lind1 = linkind(psi,i)
                            lind2 = linkind(psi,i+1)
                            tsi = ts*psi[i]*op("Cdag",sites,i)*dag(prime(prime(psi[i],"Site"),lind1))
                            tsi *= psi[i+1]*dag(prime(psi[i+1],"Link"))
                            tsi *= psi[i+2]*op("C",sites,i+2)*dag(prime(prime(psi[i+2],"Site"),lind2))
                            lri = scalar(tsi)
                            tsi = ts*psi[i]*op("C",sites,i)*dag(prime(prime(psi[i],"Site"),lind1))
                            tsi *= psi[i+1]*dag(prime(psi[i+1],"Link"))
                            tsi *= psi[i+2]*op("Cdag",sites,i+2)*dag(prime(prime(psi[i+2],"Site"),lind2))
                            rli = scalar(tsi)
                            fy[i1] = (-lri+rli)/2
                            fx[i1] = (lri+rli)/2
                        end
                    end
                    return fx,fy
                end

                # Cdga-Sz-C
                function rungTunneling(psi::MPS)
                    rungTunneling = zeros(ns)
                    ts = ITensor(1.)
                    i1 = 0
                    for i in 1:3*ns
                        if mod(i,3) == 1
                            i1 += 1
                            orthogonalize!(psi,i)
                            lind1 = linkind(psi,i)
                            lind2 = linkind(psi,i+1)
                            tsi = ts*psi[i]*op("Cdag",sites,i)*dag(prime(prime(psi[i],"Site"),lind1))
                            tsi *= psi[i+1]*op(szo,sites,i+1)*dag(prime(psi[i+1]))
                            tsi *= psi[i+2]*op("C",sites,i+2)*dag(prime(prime(psi[i+2],"Site"),lind2))
                            rungTunneling[i1] = scalar(tsi)
                        end
                    end
                    return rungTunneling
                end

---

Hi JW, I see. Yes if you are using custom code to do this sort of thing then I think I can guess what might be causing the memory issue. The answer in this case is almost always a tensor contraction or set of tensor contractions that is not being done in the best or right order. So for example, it could be that two MPS tensors are being contracted together in a way that 4 virtual indices stick out afterward, which would lead to the algorithm having chi^4 memory scaling and cost scaling which can be really bad.

One culprit here is that the order of operations is not always so obvious in Julia. (In C++ it's always simply left to right.)

So what I recommend is explicitly putting parenthesis in chains of tensor contractions such as this one from your code: scalar(ts*psi[i]*op("N",sites,i)*dag(prime(psi[i],"Site")))
(I just picked that one at random, it could be another one or multiple that are leading to the poor scaling.)

Or instead of parenthesis you can break the contractions across multiple lines like:
C = A * B
E = C * D
...

Hope that helps,
Miles

---

Also, I wasn't sure what you meant about your site type being customized and thus correlation_matrix not working. If you define overloads of op to tell ITensor about the custom operator names that go with your sites, then correlation_matrix should automatically work and recognize these custom op names that you defined.

Miles

---

Thank you for your information, I actually assume that correlation_matrix won't work if the site site is not of the standard ones. And my impression on correlation_matrix is that it is still not reliable, but it turns out sometimes my own function is even more unreliable...

---

After the bug fix I made following the forum discussion you linked, I think correlation_matrix is now pretty reliable. You could always switch away from it if it doesn't work well. It's definitely reliable in terms of giving correct results.

Also yes please do try it with your custom operator(s). If they are correctly defined then correlation_matrix will automatically "know" about them through the `op` system.

---

Hi Miles, I just realized that in some cases the Jordan-Wigner string may be added when measuring correlators. And for some custom site type (say, spin-fermion-spin sites), this string may even be discontinuous. And my naive guess is that the built-in correlation_matrix can not deal with some special cases like this.

I probably will check the source code of correlation_matrix later. But anyway, I have adopt your suggestion and find that the problem is solved!

Best,
Junsen

---

Hi Junsen,
Yes you're right to wonder or be skeptical about the behavior of correlation_matrix for fermions. I think actually it would work even for the case you described because it will put the "F" operator on all sites, but then if "F" is defined to just be the identity for spin sites (say) then all would be ok and still work. But it's complicated I agree and one would have to do some tests on cases where the answer is known to check everything.

But because of these kinds of complicated matters, I am still working on a new fermion system for ITensor that will completely avoid having to use things like "F" operators and so on. The new system is working well but there are some remaining technical challenges with things like the SVD I'm working through.

Still, I'd suggest just trying correlation_matrix. On the other hand, feel welcome to use your own code of course and I would suggest checking all the contraction orderings as I mentioned as the most likely root cause of the memory issue you had.

Miles

---

Oh I had missed your comment saying that your problem is now solved. That is great! Hope the rest of your project continues to go smoothly.

Answer 1

(Please see discussion above.)