Unexpected scaling of OpSum

DMRG and Numerical Methods
1

Hello!
I am trying to build a MPO which includes 4-sites non-local interaction. The first thing I do is define a OpSum and add all the terms I want. Given that it’s a 4 site interaction, I expect the time to construct this to scale like N^4, where N is the number of sites in the system. However, what I see is significantly worse. See the following example

using ITensorMPS
using Combinatorics

function H(N)
    ampo = OpSum()
  for (k1,k2,k3,k4) in collect(combinations(1:N, 4))
    ampo += 2*(N-k4),"Z",k1,"Z",k2,"Z",k3,"Z",k4
  end
  return ampo
end


Nvec = [10,20,30,40,50,60]

for n in Nvec
    time4 = @elapsed xx = H(n)
    println("Time for N = $(n): $(time4)")
end

The timing appears to grow much faster than N^4: for example, N=30 gives me ~1.5 seconds, and N=60 gives ~550 seconds. This is clearly a problem for me, because going to larger system sizes becomes prohibitive.
Is my expectation of a N^4 scaling wrong? I expected the MPO function to potentially take a while, but that’s already happening for OpSum

2

Hopefully others can chime in but check out

which may tackle the problem you have

3

Hi Ryan,

thanks, it does look like OpIDSum would be what I need. Unfortunately I cannot import the ITensorMPOConstruction package, I believe it wants ITensorMPS in version 0.2, while I have the latest version, and I get a “Unsatisfiable requirements detected for package ITensorMPS” error

4

Hi Bernardo,
I’d suggest you chat with the maintainer (Ben Corbett) of the ITensorMPOConstruction package about lifting the restriction on the ITensorMPS version. It’s probably not an intentional restriction, and just an accidental one.

About the existing OpSum to MPO constructor: it just wasn’t designed to handle complicated chemistry Hamiltonians. We are currently working on a rewrite of it where one goal of the rewrite is to improve the scaling.

5

Let me add one thing: you can separate an N^4 Hamiltonian into the sum of N different N^3 terms–say just by pulling one of the index sums outside. Each N^3 MPO construction may behave much better than the N^4 one. In particular, the N^3 coefficients matches the number of variables in an NxN MPO of length N–and, in fact, it is possible to construct this MPO directly without using a compression algorithm. ITensor can use an array of MPOs just as easily as a single MPO. See THE JOURNAL OF CHEMICAL PHYSICS 145, 014102 (2016)

2 Likes
6

First off simply change ampo += ... to ampo .+=, this way you are directly modifying ampo each time instead of constructing a copy and then reassigning it. With this one character modification here are my timings

julia> include("foo.jl"); # ampo += ...
Time for N = 10: 0.018505916
Time for N = 30: 3.156324
Time for N = 40: 32.901548875

julia> include("foo.jl"); # ampo .+= ...
Time for N = 10: 0.014833458
Time for N = 30: 0.033778666
Time for N = 40: 0.095046791
Time for N = 50: 0.271597917
Time for N = 60: 0.841291958

That said, I have encountered situations where this change simply isn’t fast enough, or when memory is limited, that’s why I created OpIDSum.

I just merged a PR to ITensorMPOConstruction that should get rid of the version conflicts (for now at least). Below is the example code modified to also use ITensorMPOConstruction.MPO_new both from an OpSum and OpIDSum. TLDR is that for this case OpSum is probably fast enough, and that MPO_new can construct the MPO for N = 110 in the time it takes ITensorMPS.MPO to construct for N = 50.

Code 
using ITensorMPS
using Combinatorics
using ITensorMPOConstruction

function H_OpSum(N)::OpSum{Float64}
  os = OpSum{Float64}()
  for (k1,k2,k3,k4) in collect(combinations(1:N, 4))
    os .+= 2*(N-k4),"Z",k1,"Z",k2,"Z",k3,"Z",k4
  end

  return os
end

function H_OpIDSum(N)::Tuple{OpIDSum{Float64}, Vector{Vector{String}}}
  operatorNames = [["I", "Z"] for _ in 1:N]
  opZ(n) = OpID(2, n)

  os = OpIDSum{Float64}()
  for (k1,k2,k3,k4) in collect(combinations(1:N, 4))
    push!(os, Float64(2*(N-k4)), opZ(k1), opZ(k2), opZ(k3), opZ(k4))
  end

  return os, operatorNames
end

for N in [10:10:120]
  sites = siteinds("Qubit", N)

  println("N = $N")
  if N <= 50
    @time "\tOpSum construction  " os = H_OpSum(N)
    @time "\tMPO construction    " mpo = MPO(os, sites)
    @time "\tMPO_new construction" mpo = MPO_new(os, sites)
  end

  @time "\tOpIDSum construction" os, operatorNames = H_OpIDSum(N)
  @time "\tMPO_new construction" mpo = MPO_new(os, sites, operatorNames; basis_op_cache_vec=operatorNames)
  println()
end
Output 
N = 10
    OpSum construction  : 0.000198 seconds (8.41 k allocations: 554.844 KiB)
    MPO construction    : 0.011609 seconds (149.35 k allocations: 16.557 MiB)
    MPO_new construction: 0.001477 seconds (11.27 k allocations: 1.067 MiB)
    OpIDSum construction: 0.000036 seconds (452 allocations: 55.156 KiB)
    MPO_new construction: 0.001240 seconds (10.66 k allocations: 1.308 MiB)

N = 20
    OpSum construction  : 0.003953 seconds (193.81 k allocations: 12.399 MiB)
    MPO construction    : 0.406268 seconds (4.61 M allocations: 438.247 MiB)
    MPO_new construction: 0.017774 seconds (155.02 k allocations: 17.642 MiB)
    OpIDSum construction: 0.000392 seconds (9.74 k allocations: 1.029 MiB)
    MPO_new construction: 0.016200 seconds (113.98 k allocations: 19.220 MiB)

N = 30
    OpSum construction  : 0.031350 seconds (1.10 M allocations: 70.375 MiB)
    MPO construction    : 3.718616 seconds (33.19 M allocations: 2.870 GiB, 9.76% gc time)
    MPO_new construction: 0.100924 seconds (834.97 k allocations: 118.611 MiB)
    OpIDSum construction: 0.001746 seconds (54.88 k allocations: 6.059 MiB)
    MPO_new construction: 0.088726 seconds (592.16 k allocations: 124.450 MiB)

N = 40
    OpSum construction  : 0.092522 seconds (3.66 M allocations: 233.050 MiB)
    MPO construction    : 17.120962 seconds (134.23 M allocations: 11.427 GiB, 6.97% gc time)
    MPO_new construction: 0.419496 seconds (2.78 M allocations: 493.040 MiB, 0.35% compilation time)
    OpIDSum construction: 0.006401 seconds (182.86 k allocations: 18.643 MiB)
    MPO_new construction: 0.364457 seconds (1.94 M allocations: 510.816 MiB, 9.73% gc time)

N = 50
    OpSum construction  : 0.265503 seconds (9.21 M allocations: 583.985 MiB, 21.69% gc time)
    MPO construction    : 70.386649 seconds (400.07 M allocations: 34.228 GiB, 6.95% gc time)
    MPO_new construction: 1.170543 seconds (6.95 M allocations: 1.326 GiB, 2.96% gc time, 0.12% compilation time)
    OpIDSum construction: 0.015412 seconds (460.69 k allocations: 41.196 MiB)
    MPO_new construction: 1.002872 seconds (4.85 M allocations: 1.372 GiB, 3.58% gc time)

N = 60
    OpIDSum construction: 0.086254 seconds (975.38 k allocations: 83.676 MiB, 63.58% gc time)
    MPO_new construction: 2.551916 seconds (10.29 M allocations: 3.235 GiB, 3.50% gc time, 0.07% compilation time)

N = 70
    OpIDSum construction: 0.092260 seconds (1.83 M allocations: 157.595 MiB, 23.97% gc time)
    MPO_new construction: 6.453138 seconds (19.29 M allocations: 6.803 GiB, 10.16% gc time, 0.03% compilation time)

N = 80
    OpIDSum construction: 0.142721 seconds (3.16 M allocations: 280.324 MiB, 15.90% gc time)
    MPO_new construction: 12.330305 seconds (33.23 M allocations: 12.989 GiB, 6.08% gc time, 0.02% compilation time)

N = 90
    OpIDSum construction: 0.243893 seconds (5.11 M allocations: 479.688 MiB, 15.81% gc time)
    MPO_new construction: 22.095196 seconds (53.65 M allocations: 22.364 GiB, 6.99% gc time, 0.01% compilation time)

N = 100
    OpIDSum construction: 0.347497 seconds (7.84 M allocations: 656.863 MiB, 16.94% gc time)
    MPO_new construction: 40.768082 seconds (82.30 M allocations: 37.863 GiB, 7.64% gc time, 0.00% compilation time)

N = 110
    OpIDSum construction: 0.499944 seconds (11.55 M allocations: 1011.284 MiB, 15.51% gc time)
    MPO_new construction: 67.607090 seconds (121.15 M allocations: 59.888 GiB, 11.14% gc time, 0.00% compilation time)

N = 120
    OpIDSum construction: 0.793480 seconds (16.43 M allocations: 1.342 GiB, 19.97% gc time)
    MPO_new construction: 142.217085 seconds (172.36 M allocations: 86.625 GiB, 30.62% gc time, 0.00% compilation time)
1 Like
7

I don’t think this is always a good idea. For example, the interaction term in the momentum space Fermi-Hubbard Hamiltonian

\sum_{p, q, k} \hat{c}^\dagger_{p - k, \uparrow} \hat{c}^\dagger_{q + k, \downarrow} \hat{c}_{q, \downarrow} \hat{c}_{p, \uparrow} \ ,

has O(N^3) terms and can be represented as an MPO of bond dimension O(N). By factoring out either the sum over p, q or k you still wind up with N MPO’s of bond dimension O(N). This behavior is essentially unchanged for the transcorrelated momentum space Fermi-Hubbard Hamiltonian, which has O(N^5) terms (https://doi.org/10.1063/5.0028608). If there’s a smarter way of grouping terms than simply factoring out a sum I would be very interested.

Output 
N = 10
  All terms bond dim : 96
  Max, median bond dim after factoring out p: 32, 29.0
  Max, median bond dim after factoring out q: 32, 29.0
  Max, median bond dim after factoring out k: 100, 50.0
N = 20
  All terms bond dim : 196
  Max, median bond dim after factoring out p: 62, 62.0
  Max, median bond dim after factoring out q: 62, 62.0
  Max, median bond dim after factoring out k: 400, 144.0
N = 40
  All terms bond dim : 396
  Max, median bond dim after factoring out p: 122, 122.0
  Max, median bond dim after factoring out q: 122, 122.0
  Max, median bond dim after factoring out k: 1600, 484.0
N = 80
  All terms bond dim : 796
  Max, median bond dim after factoring out p: 242, 242.0
  Max, median bond dim after factoring out q: 242, 242.0
Code 
using ITensorMPS
using ITensorMPOConstruction
using Statistics

function get_mpo_bond_dim(N; pRange=1:N, qRange=1:N, kRange=1:N)::Int
  sites = siteinds("Electron", N)
  operatorNames = [
    [
      "I",
      "Cdn",
      "Cup",
      "Cdagdn",
      "Cdagup",
      "Ndn",
      "Nup",
      "Cup * Cdn",
      "Cup * Cdagdn",
      "Cup * Ndn",
      "Cdagup * Cdn",
      "Cdagup * Cdagdn",
      "Cdagup * Ndn",
      "Nup * Cdn",
      "Nup * Cdagdn",
      "Nup * Ndn",
    ] for _ in 1:N
  ]

  os = OpSum{Float64}()
  for p in pRange
    for q in qRange
      for k in kRange
        pmk = mod1(p - k, N)
        qpk = mod1(q + k, N)
        os .+= 1.0, "Cdagup", pmk, "Cdagdn", qpk, "Cdn", q, "Cup", p
      end
    end
  end

  return maxlinkdim(MPO_new(os, sites; basis_op_cache_vec=operatorNames))
end

for N in [10, 20, 40, 80]
  println("N = $N")
  println("  All terms bond dim : ", get_mpo_bond_dim(N))
  
  pdims = [get_mpo_bond_dim(N; pRange=i:i) for i in 1:N]
  println("  Max, median bond dim after factoring out p: $(maximum(pdims)), $(median(pdims))")

  qdims = [get_mpo_bond_dim(N; qRange=i:i) for i in 1:N]
  println("  Max, median bond dim after factoring out q: $(maximum(qdims)), $(median(qdims))")

  kdims = [get_mpo_bond_dim(N; kRange=i:i) for i in 1:N]
  println("  Max, median bond dim after factoring out k: $(maximum(kdims)), $(median(kdims))")
end
8

Thanks! I ended up being able to use the ITensorMPOConstruction package before the change, I got around the version compatibility issues by doing a clean installation of julia. With this, I am able to compute the operators I need in a reasonable amount of time

1 Like