Artificially increase bond dimension of MPS with conserved QNs

ITensor Julia Questions
, ,
1

Hi everyone,
I’m trying to artificially increase the bond dimensions of a product-state MPS to a given bond dimension N (since I want to use this state in a TDVP evolution algorithm with fixed bond dimensions). I start from

sites = siteinds("Qubit", 10; conserve_parity=true)
initstate = MPS(sites, "0")

Normally, without QNs, I would use the method described here to “artificially” enlarge the bond dimension to the desired size.
Taking quantum numbers into account, I defined the following function:

function enlarge!(v::MPS, new_d)
    N = length(v)
    if hasqns(v[1])
        lflux = sum(flux, v[1:end-1])
        enlarged_links = Vector{ITensors.QNIndex}(undef, N - 1)
        for bond in (N - 1):-1:1
            enlarged_links[bond] = dag(Index(lflux => new_d; tags=tags(commonind(v[bond], v[bond + 1]))))
            lflux -= flux(v[bond])
        end
    else
        enlarged_links = [
            Index(new_d; tags=tags(commonind(v[bond], v[bond + 1]))) for bond in 1:(N - 1)
        ]
    end
    @show first(linkinds(v), 5)
    @show first(enlarged_links, 5)

    for bond in 1:(N - 1)
        bond_index = commonind(v[bond], v[bond + 1])
        v[bond] = v[bond] * delta(dag(bond_index),enlarged_links[bond])
        v[bond + 1] = v[bond + 1] * delta(bond_index, dag(enlarged_links[bond]))
    end
    return v
end

I took inspiration from the MPS constructor method

MPS(eltype::Type{<:Number}, sites::Vector{<:Index}, states_)

in ITensors/…/src/lib/ITensorMPS/src/mps.jl to define the new link indices, since this constructor is what gets called when I create the product state above; I figured I would get the same index structure.
Unfortunately this does not happen: the @show statements print

julia> enlarge!(initstate, 10)

first(linkinds(v), 5) = Index{Vector{Pair{QN, Int64}}}[(dim=1|id=315|"Link,l=1") <In>
 1: QN("Parity",1,2) => 1, (dim=1|id=773|"Link,l=2") <In>
 1: QN("Parity",0,2) => 1, (dim=1|id=639|"Link,l=3") <In>
 1: QN("Parity",1,2) => 1, (dim=1|id=593|"Link,l=4") <In>
 1: QN("Parity",0,2) => 1, (dim=1|id=474|"Link,l=5") <In>
 1: QN("Parity",1,2) => 1]
first(enlarged_links, 5) = Index{Vector{Pair{QN, Int64}}}[(dim=10|id=583|"Link,l=1") <In>
 1: QN("Parity",0,2) => 10, (dim=10|id=654|"Link,l=2") <In>
 1: QN("Parity",0,2) => 10, (dim=10|id=118|"Link,l=3") <In>
 1: QN("Parity",0,2) => 10, (dim=10|id=72|"Link,l=4") <In>
 1: QN("Parity",0,2) => 10, (dim=10|id=973|"Link,l=5") <In>
 1: QN("Parity",0,2) => 10]

and as you can see the alternating parity in the link indices is lost in the enlarged MPS.
So I ask you: what isn’t working in my method? Or is there another way I can enlarge the bond dimensions of an MPS?

2

There are many ways of expanding the bond dimension of an MPS with tradeoffs in terms of what properties you want the new basis to have, computational cost of performing the expansion, etc. and it is even an active area of research to find new ways.

You could try the ITensorMPS.expand function mentioned here: GitHub - ITensor/ITensorMPS.jl: MPS and MPO methods based on ITensor (ITensors.jl), execute julia> ?expand at the REPL to see instructions on how to use it.

3

I’ve already tried expand in the following way:

N=10
s = siteinds("Qubit", N; conserve_parity=true)
v = MPS(s, "Dn")
h = OpSum()
for n in 1:N
       h += "S+ * S-", n
end
for n in 1:N-1
       h += "S+", n, "S-", n+1
       h -= "S-", n, "S+", n+1
end
H = MPO(h, s)
state_t = expand(v, H; alg="global_krylov")

but I always get this error:

ERROR: BoundsError: attempt to access 0-element Vector{Pair{ITensors.QuantumNumbers.QN, Int64}} at index [1]
Stacktrace:
  [1] throw_boundserror(A::Vector{Pair{ITensors.QuantumNumbers.QN, Int64}}, I::Tuple{Int64})
    @ Base ./essentials.jl:14
  [2] getindex
    @ ./essentials.jl:915 [inlined]
  [3] combineblocks(qns::Vector{Pair{ITensors.QuantumNumbers.QN, Int64}})
    @ ITensors /opt/julia/packages/ITensors/AbUYR/src/qn/qnindex.jl:417
  [4] combineblocks
    @ /opt/julia/packages/ITensors/AbUYR/src/qn/qnindex.jl:479 [inlined]
  [5] combiner(inds::Tuple{ITensors.Index{…}, ITensors.Index{…}}; dir::ITensors.QuantumNumbers.Arrow, tags::String)
    @ ITensors /opt/julia/packages/ITensors/AbUYR/src/qn/qnitensor.jl:449
  [6] combiner
    @ /opt/julia/packages/ITensors/AbUYR/src/qn/qnitensor.jl:445 [inlined]
  [7] #combiner#322
    @ /opt/julia/packages/ITensors/AbUYR/src/tensor_operations/itensor_combiner.jl:7 [inlined]
  [8] eigen(A::ITensors.ITensor, Linds::Vector{…}, Rinds::Vector{…}; mindim::Int64, maxdim::Int64, cutoff::Float64, use_absolute_cutoff::Nothing, use_relative_cutoff::Nothing, ishermitian::Bool, tags::ITensors.TagSets.GenericTagSet{…}, lefttags::Nothing, righttags::Nothing, plev::Nothing, leftplev::Nothing, rightplev::Nothing)
    @ ITensors /opt/julia/packages/ITensors/AbUYR/src/tensor_operations/matrix_decomposition.jl:365
  [9] contract(::NDTensors.BackendSelection.Algorithm{…}, A::MPO, ψ::MPS; cutoff::Float64, maxdim::Int64, mindim::Int64, normalize::Bool, kwargs::@Kwargs{})
    @ ITensors.ITensorMPS /opt/julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/mpo.jl:762
 [10] contract
    @ /opt/julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/mpo.jl:691 [inlined]
 [11] product(alg::NDTensors.BackendSelection.Algorithm{…}, A::MPO, ψ::MPS; kwargs::@Kwargs{…})
    @ ITensors.ITensorMPS /opt/julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/mpo.jl:604
 [12] product
    @ /opt/julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/mpo.jl:603 [inlined]
 [13] #apply#453
    @ /opt/julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/mpo.jl:600 [inlined]
 [14] product
    @ /opt/julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/mpo.jl:599 [inlined]
 [15] expand(::NDTensors.BackendSelection.Algorithm{…}, state::MPS, operator::MPO; krylovdim::Int64, cutoff::Float64, apply_kwargs::@NamedTuple{…})
    @ ITensorTDVP /opt/julia/packages/ITensorTDVP/O4qn7/src/expand.jl:159
 [16] expand(::NDTensors.BackendSelection.Algorithm{:global_krylov, @NamedTuple{}}, state::MPS, operator::MPO)
    @ ITensorTDVP /opt/julia/packages/ITensorTDVP/O4qn7/src/expand.jl:147
 [17] expand(state::MPS, reference::MPO; alg::String, kwargs::@Kwargs{})
    @ ITensorTDVP /opt/julia/packages/ITensorTDVP/O4qn7/src/expand.jl:31
 [18] top-level scope
    @ REPL[24]:1
Some type information was truncated. Use `show(err)` to see complete types.

and now that I took more time to investigate this error I realize that it comes from just doing H * v. I didn’t expect that… I’ll try to understand what’s happening here.

4

Starting with an all down state:

v = MPS(s, "Dn")

likely isn’t what you want, the state will always remain in that state because of the QN conservation. You should start in a state in the symmetry sector you are interested in, i.e. start in a Néel state if you want to study the 0 spin sector.

Probably the error you are getting is because that state is being projected to the norm zero state by the Hamiltonian, ideally we would have a better error message (or output a norm zero state).

5

You are absolutely right. I chose the first v and H that came to mind without thinking about their action. I can get something from expand if I use a proper v and H, and the bond dimension does increase a bit on some sites.
However I don’t think this method is for me since I can’t control the bond dimension of the output. I’d like to be able to get an MPS of a predetermined bond dimension.

In the end I realized the mistake in my original method: I was computing the flux of the MPS elements v[n], but the MPS constructor I took the code from computes it from the state ITensors the MPS is built from. The states do not yet have the link indices attached to them, while the v[n] elements do, and I’m guessing this is the cause of the different QNs I’m getting.
I can recover the original states (without the links) by doing

states = Vector{ITensor}(undef, N)                                                   
states[1] = v[1] * dag(onehot(linkind(v,1) => 1))                                    
for n in 2:N-1                                                                       
    states[n] = v[n] * onehot(linkind(v,n-1) => 1) * onehot(dag(linkind(v, n)) => 1) 
end                                                                                  
states[N] = v[N] * onehot(linkind(v, N-1) => 1)

and compute flux on states instead of v. This seems to get me what I want in the end, even if it’s probably a bit convoluted.