Hi Kevin,
I’m the creator of ITensorMPOConstruction the independent library you alluded to. I spent a lot of time thinking about MPO construction and I’ve never seen that paper, so thank you for bringing it to my attention. It does seem to be primarily focused on MPO compression for infinite systems, which is not something I am familiar with.
I have only had a chance to skim it, but my basis understanding is that they treat the starting and ending identity strings differently from the “bulk” terms and truncate the bulk after bringing it into canonical form. ITensorMPOConstruction does not treat the identity strings differently whereas ITensorMPS does treat the starting and ending identity string differently and only does the SVD on the bulk. Both algorithms, should construct the MPO in left cannonical form if the local operator basis is orthonormal.
So in that sense the algorithm in ITensorMPS seems to match the algorithm from that paper, but it does the “compression” online as the MPO is built, instead of building a naive MPO first and then compressing it.
To that end, I constructed the Hamiltonian from Equation 34 the bond dimensions of which they plot in Figure 2
H = \sum_{k, n, m} \frac{1}{(k - n)^2 (n - m)^2} Z_k Z_n Z_m + \sum_{n, m} \frac{1}{(n - m)^4} Z_n Z_m \ .
They claim that a naive construction has a bond dimension of 1000, and they obtained a bond dimension of roughly 15 using their truncation scheme with a cutoff of 10^{-4} (note: their truncation cutoff likely does not directly translate to the cutoff in ITensors). ITensorMPS produces an MPO of bond dimension 18 and ITensorMPOConstruction an MPO of bond dimension 34. This situation seems very similar to the situation I encountered for the Haldane-Shastry Hamiltonian, which I documented extensively.
Code: click to expand
using ITensors, ITensorMPS, ITensorMPOConstruction
function get_OpSum(N)
os = OpSum{Float64}()
J(n, m) = 1.0 / (n - m)^2
for k in 1:N
for n in 1:N
for m in 1:N
k == n && continue
n == m && continue
os .+= J(k, n) * J(n, m), "Z", k, "Z", n, "Z", m
end
end
end
for n in 1:N
for m in 1:N
n == m && continue
os .+= J(n, m)^2, "Z", n, "Z", m
end
end
basis = [["I", "Z"] for _ in 1:N]
return os, siteinds("Qubit", N; conserve_qns=true), basis
end
let
os, sites, basis = get_OpSum(30)
@time H = MPO(os, sites; cutoff=1e-15)
@show maxlinkdim(H)
@time Hnew = MPO_new(os, sites; splitblocks=true, basis_op_cache_vec=basis)
@show maxlinkdim(Hnew)
end
return nothing
As to your question about reducing the MPO bond dimension for complicated systems, ITensorMPOConstruction will construct an exact MPO of the optimal bond dimension up to numeric precision. ITensorMPS has a similar guarantee. For the systems I have studied, primarily the momentum space Fermi-Hubbard Hamiltonian and electronic structure Hamiltonian, further numeric compression of the MPO after it has already been constructed is counterproductive. A large truncation error is required to meaningfully decrease the bond dimension, and the resulting energies obtained with DMRG are poor. For example, see Figure 9 of my paper Scaling up the transcorrelated density matrix renormalization group, where a truncation cutoff of 10^{-6} is required to reduce the bond dimension by a factor of two.
Aside on truncation
I claim, as I say at the top of the README, that “ITensorMPOConstruction is not designed to construct approximate compressed MPOs”. But this really depends on what you mean by compressed. In my mind there are three different starting places for a Hamiltonian
H = \sum_{i = 1}^M c_i O_i
- As a truly naive MPO of bond dimension M, one for each term.
- Using the MPO produced by the minimum vertex cover algorithm. This is the minimum MPO bond dimension required to represent the Hamiltonian for any coefficients c_i.
- Using the MPO produced by the rank decomposition algorithm used in ITensorMPS and by default in ITensorMPOConstruction. This is the the minimum bond dimension required for those specific coefficients c_i.
I think using option 2 or 3 as your starting point for compression is the only fair choice.
Update on the Hamiltonian from Equation 34
I tried using the vertex cover algorithm, which does no numerical compression of any kind, and obtained an MPO of bond dimension 32. Interestingly this is less than the MPO from the QR decomposition, which in exact arithmetic should not be possible. Upon further investigation it appears to be due to the numerical noise present when the terms in the Hamiltonian span six orders of magnitude. It appears from the trends that the vertex cover produces an MPO of bond dimension N + 2 for N qubits, so this is actually a fairly compact MPO. Certainly this scaling could be greatly reduced if you were to tolerate some approximation error.