More efficiently compute the expectation value of an MPO that scale as O(N^4) (Binder cumulant)?

DMRG and Numerical Methods
,
1

Hello,

I am curious what is the most efficient way to compute the expectation value of an MPO which contains N^4 terms (like the Binder cumulant).

I am asking this becasuse I am recently going through this paper (Matrix product state techniques for two-dimensional systems at finite temperature), where the Binder cumulant is computed to detect the phase boundaries. And I am trying to do exactly the same thing on another spin-1/2 Heisenberg type model on the square lattice. However, I found it became very hard to efficiently compute the fourth order cumulants.

I am using Julia version ITensor. The bond dimension is D, and the N_y's I choose are 4,6,8 (for D=2, N=16, 36, 64). And I tried to compute the fourth order cumulants in two different ways.

(1) Construct an MPO that contains all the “necessary” terms. By saying necessary, I mean that for all the repeated terms, I only count once and multiply the term with a corresponding prefactor. And then use
‘’’
inner(psi, MPO, psi)
‘’’

(2) Make for nested loops and compute each term individually by doing something similar to what “correlation_matrix” function does, but with four sites (not using inner() function, but only scalar() function).

What I get from the two methods are basically as follows.

(1) works great for N_y=4. But N_y=6 then becomes a disaster. And the scaling of the computational cost seems to be greated than O(N^4).

(2) is less efficient than (1) for N_y=4. But I guess the scaling of its computational complexity is under control, which is O(N^4)=O((N_y^2)^4).

Do you have any advice on accelerating the computation of such observables? And if I may ask, I am also curious about how you computed the Binder cumulant in the above paper I mentioned.

Best,
Yi

2

Just found another paper which is relevant. Efficient evaluation of high-order moments and cumulants in tensor network states There, they basically computed the generating function exp(a*O) by applying Trotter gates.

3

Thanks for your patience with my slow reply here. It’s a really interesting question about how to do this efficiently, because it does seem like there should be a nice solution. The good news is that I think there is, if I’m understanding correctly what you want to do.

First of all, that does seem like a very useful and relevant paper – glad you found it.

But more to the point, if you are looking to construct operators M^2 and M^4 efficiently, where for concreteness let’s say that

M = \sum_j Z_j

where Z_j means the Pauli matrix \sigma^z on site j, then the good news is that there is a very efficient, exact MPO constructor for these operators.

The MPO for M^n in this case has a bond dimension of n+1 and basically implementing a rule like “put the first operator on site i1, then the second on some later site i2, etc. until we have put n operators then just put identity operators after that”. So that only takes a bond dimension of n+1 because the internal “states” of the MPO are (1) having put zero operators, (2) having put one operator, … (n) having put (n-1) operators, (n+1) having put all operators. You can use ITensor’s OpSum system to check this statement for small lattice sizes.

The actual construction of the MPO is a bit tedious but also straightforward.

If we take the example of M^2, we can write this as:

M^2 = \sum_{i,j} Z_i Z_j = 2 \sum_{i < j} Z_i Z_j + \sum_{j} Z^2_j

Then if you know how to write operators like the above as MPOs consisting of “operator valued matrices”, then here is an MPO pattern that makes the above M^2 operator:

\begin{bmatrix} I_j & 0 & 0 \\ Z_j & I_j & 0\\ Z^2_j & 2 Z_j & I_j \end{bmatrix}

with boundary conditions being the vectors [0, 0, 1] on the left and [1, 0, 0]^\dagger on the right.

A similar, but more complicated construction holds for M^4.

Please let me know if you want to pursue this approach further and have more questions about it.

1 Like
4

By the way, here’s a short demo code to verify that there is an n+1 dimensional MPO for the operator M^n. But as you mentioned, this approach of using OpSum doesn’t scale well to really large systems, so that’s why I mentioned the “direct” construction approach instead.

using ITensors

let
  N = 20
  s = siteinds("S=1/2",N)

  os2 = OpSum()
  for i=1:N, j=1:N
    os2 += "Z",i,"Z",j
  end
  M2 = MPO(os2,s)

  os4 = OpSum()
  for i=1:N, j=1:N, k=1:N, l=1:N
    os4 += "Z",i,"Z",j,"Z",k,"Z",l
  end
  M4 = MPO(os4,s)

  @show linkdims(M2)
  @show linkdims(M4)

  return
end

Output:

linkdims(M2) = [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3]
linkdims(M4) = [5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5]
5

Hi Miles,

Thank you for the detailed reply! I see what you mean. But I am thinking in a slightly different way. For M^n, at each site, one can choose to put 0, 1, \cdots, n operators. And this gives n+1 possibilities, which means there are n+1 “channels” from this site to other sites, i.e., the bond dimension is n+1. I am not sure whether this is the correct way to think of it.

Regarding the “operator valued matrices”, I have seen something similar in your iDMRG code, but I never know how to construct the matrices and make use of them to compute the MPOs. It would be wonderful if you could point me to some literatures or perhaps related answers.

Thanks again!

Best,
Yi

6

By the way, you can take powers of OpSum to easily create these MPOs:

using ITensors

function sum_z(n)
  os = OpSum()
  for j in 1:n
    os += "Z", j
  end
  return os
end

n = 10
s = siteinds("S=1/2", n)
h = sum_z(n)
H = @time MPO(h, s)
@show maxlinkdim(H)
H² = @time MPO(Ops.expand(h^2), s)
@show maxlinkdim(H²)
H³ = @time MPO(Ops.expand(h^3), s)
@show maxlinkdim(H³)
H⁴ = @time MPO(Ops.expand(h^4), s)
@show maxlinkdim(H⁴)
H⁵ = @time MPO(Ops.expand(h^5), s)
@show maxlinkdim(H⁵)

This outputs:

  0.000758 seconds (3.12 k allocations: 826.562 KiB)
maxlinkdim(H) = 3
  0.004852 seconds (27.59 k allocations: 5.191 MiB)
maxlinkdim(H²) = 3
  0.038790 seconds (214.14 k allocations: 27.082 MiB, 19.01% gc time)
maxlinkdim(H³) = 4
  0.294138 seconds (1.79 M allocations: 151.750 MiB, 8.14% gc time)
maxlinkdim(H⁴) = 5
  3.404573 seconds (17.71 M allocations: 1.162 GiB, 13.94% gc time, 0.11% compilation time)
maxlinkdim(H⁵) = 6
3 Likes
7

Thanks for pointing it out! And the scaling looks very nice. I will try it out.

best,

8

Nice to see the Ops.expand trick above – I knew Matt had upgraded the OpSum system recently but didn’t know about that feature :^)

@slowlight, to answer your question about what operator-valued matrices mean and how they relate to MPOs, here is a small primer showing the example of an MPO that would make the sum \sum_j Z_j:

IMG_CC88C1188D69-1
IMG_CC88C1188D69-11000×1216 92.7 KB

This way of thinking about MPOs (and MPS) was introduced by various people, and is related to “finite state automata” ideas. Here are some references:

2 Likes
9

Hi Miles and Matt,

Thank you so much for the messages! And apology for the late reply. It took me some time to digest all the information. These really help me understand the MPOs stuff!

best,

10

Its an old thread, but there is also a relevant article here: [1906.03833] Cumulants and Scaling Functions of Infinite Matrix Product States

3 Likes
11

Thanks Ian,

This is still very helpful as I have recently started to try to play with iDMRG as well.

Best,
Yi

1 Like
12

Thanks, Ian