randomMPS does not produce random state with relabeled quantum numbers.

ITensor Julia Questions
,
1

I’m working on finding the ground state the Fermi-Hubbard Hamiltonian, so the electron site type is the natural choice. Specifically I am interested in the ground state with a given particle number and total spin, so I am using the two native quantum numbers associated with the electron site.

The DMRG accuracy is improved by using non-standard site orderings (mapping from a physical DOF to the position in the MPS). However the use of the electron site type limits the possible re-orderings since spin up and spin down sites are always bundled together. I recently tried using fermion sites instead as this should allow me to mix spin up sites and spin down sites as I see fit. I seem to have gotten everything working except for the initial state that I use as input to DMRG. I am using the quantum number version of randomMPS, and I found in the electron case that providing a large initial link dimension significantly improved the DMRG results. However when using fermion sites randomMPS simply spits out the same input product state.

This certainly comes down to my use of quantum numbers and my ordering of the sites. What seemed natural to me is to rename the “total number of fermions” quantum number into Ndn(Nup) for the spin down(up) sites as done in the documentation. If I don’t rename the quantum numbers then randomMPS produces a non-product state. If I order the sites differently so that all the spin up sites come first then randomMPS also produces a non-product state.

Much appreciated!

using ITensors

N = 4
nUp = 1
nDn = 1

sites = [isodd(i) ? siteind("Fermion"; conserve_qns=true, qnname_nf="Ndn") :
                    siteind("Fermion"; conserve_qns=true, qnname_nf="Nup") for i = 1:N]

# Construct the initial state
initialState = ["Emp" for i = 1:N]

for i in 1:nDn
  initialState[2 * (i - 1) + 1] = "Occ"
end

for i in 1:nUp
  initialState[2 * (i - 1) + 2] = "Occ"
end

@show initialState

psi = randomMPS(sites, initialState; linkdims=100)
@show(psi)

@show norm(psi - MPS(sites, initialState))

The output is

initialState = ["Occ", "Occ", "Emp", "Emp"]
┌ Warning: MPS center bond dimension is less than requested (you requested 100, but in practice it is 1. This is likely due to technicalities of truncating quantum number sectors.
└ @ ITensors ~/.julia/packages/ITensors/HjjU3/src/mps/mps.jl:149
psi = MPS
[1] ((dim=2|id=524|"Fermion,Site") <Out>
 1: QN("Ndn",0,-1) => 1
 2: QN("Ndn",1,-1) => 1, (dim=1|id=983|"Link,l=1") <Out>
 1: QN(("Ndn",0,-1),("Nup",1,-1)) => 1)
[2] ((dim=2|id=884|"Fermion,Site") <Out>
 1: QN("Nup",0,-1) => 1
 2: QN("Nup",1,-1) => 1, (dim=1|id=400|"Link,l=2") <Out>
 1: QN(("Ndn",0,-1),("Nup",0,-1)) => 1, (dim=1|id=983|"Link,l=1") <In>
 1: QN(("Ndn",0,-1),("Nup",1,-1)) => 1)
[3] ((dim=2|id=562|"Fermion,Site") <Out>
 1: QN("Ndn",0,-1) => 1
 2: QN("Ndn",1,-1) => 1, (dim=1|id=819|"Link,l=3") <Out>
 1: QN("Nup",0,-1) => 1, (dim=1|id=400|"Link,l=2") <In>
 1: QN(("Ndn",0,-1),("Nup",0,-1)) => 1)
[4] ((dim=2|id=572|"Fermion,Site") <Out>
 1: QN("Nup",0,-1) => 1
 2: QN("Nup",1,-1) => 1, (dim=1|id=819|"Link,l=3") <In>
 1: QN("Nup",0,-1) => 1)

norm(psi - MPS(sites, initialState)) = 0.0

My version information

julia> versioninfo()
Julia Version 1.9.0
Commit 8e630552924 (2023-05-07 11:25 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin22.4.0)
  CPU: 16 × Intel(R) Core(TM) i9-9880H CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, skylake)
  Threads: 1 on 16 virtual cores

julia> using Pkg; Pkg.status("ITensors")
Status `~/.julia/environments/v1.9/Project.toml`
  [9136182c] ITensors v0.3.41
2

My guess is that it is something about how we are creating the random unitaries that generate the random state, and having two alternating fermion flavors is constraining the unitaries to be trivial.

A workaround is to create the states for the different fermion flavors separately, and then interleave those states into a single state:

N = 10
chi = 10
s1 = siteinds("Fermion", N ÷ 2; conserve_qns=true, qnname_nf="Ndn")
s2 = siteinds("Fermion", N ÷ 2; conserve_qns=true, qnname_nf="Nup")
psi1 = randomMPS(s1, j -> isodd(j) ? "Emp" : "Occ"; linkdims=chi)
psi2 = randomMPS(s2, j -> isodd(j) ? "Emp" : "Occ", linkdims=chi)
@show linkdims(psi1)
@show linkdims(psi2)

psi_paired = MPS([psi1[j] * psi2[j] for j in 1:(N ÷ 2)])
psi_paired = orthogonalize(psi_paired, 1)
psi = MPS(mapreduce(vcat, 1:(N ÷ 2)) do j
  psi_j, psi_jp1 = factorize(psi_paired[j], [linkinds(psi_paired, j - 1); s1[j]])
  return [psi_j, psi_jp1]
end)
@show linkdims(psi)

It’s a bit unfortunate that we need this workaround, ideally we could automatically detect this kind of situation and fix it in randomMPS but I’m not sure how feasible that is.

Something else you can do is emulate how randomMPS works and create your own random circuit and apply it to a product state, in which case you can make that circuit however you want to. In that way you can make sure to manually make a circuit which accounts for having two alternating fermion flavors.

1 Like
3

This works great, thank you! A quick test run shows that I am getting both better energies and more rapid convergence with this new construction. I don’t fully understand what is going on, but I’ll play around with it and see if I can modify it to support arbitrary mappings.

The only immediate concern I have is that the resulting MPS does not have the requested bond dimension. When I run the example code you posted I have the following output

linkdims(psi1) = [2, 4, 4, 2]
linkdims(psi2) = [2, 4, 4, 2]
linkdims(psi) = [2, 4, 8, 16, 22, 16, 8, 4, 2]

As you can see the bond dimension of the final MPS is 22 and not the requested 10. For a problem with 16 spin up and 16 spin down sites with a requested bond dimension of 100 the maximum bond dimension winds up being above 14,000. Is there a way to construct the final MPS with the requested bond dimension without having to construct a much larger MPS and then truncate it?

Thanks again!

4

Yeah generically this approach would lead to an MPS with bond dimensions that are a product of the two MPS bond dimensions. You may just have to set of bond dimensions of the separate MPS to square root of the desired bond dimension. Otherwise you could just truncate the MPS with:

psi = truncate(psi, chi)

after the fact, or pass a maxdim to factorize in my code above, which won’t give an “optimal” truncation but it may not matter too much since they are random MPS anyway.

@miles any thoughts on whether this would be easy to fix inside randomMPS?

5

@mtfishman your solution works just fine when the sites alternate between spin down and spin up, but doesn’t work when the sites are arranged in a more arbitrary order.It also looks to me like extending it to support an arbitrary order would be computationally prohibitive.

It turns out that the random MPS wasn’t as important as I thought, if I slowly grow the bond dimension and perform many sweeps that performs just as well. However I do feel like this is a limitation of randomMPS that could be documented with say a warning message. Similar to how you get a warning if you request a much larger bond dimension than necessary.

6

Good question @mtfishman. I don’t actually know if this issue would be straightforward to fix in the current code or not. But for some time I’ve wanted to overhaul the QN-conserving randomMPS anyway, since it has other drawbacks (e.g. mainly that it produces states with strictly finite-ranged correlations, instead of having some exponential decay of correlations as it should). So in the process of changing the algorithm it might very well be possible to fix this issue too.

I’ve updated our issue about improving the QN case of randomMPS to mention this post. It’s an issue that is assigned to me but a blocker for it was that we did not have a QN-conserving QR factorization. Now that we have that, it would be a good time for me to rewrite that code (and could probably be written generically to use nearly the same, better, code as for the dense case).