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.