Using the sum of exponentials to fit long-range interactions

[Ang]

This question comes from a GitHub discussion. According to the answer there, the OpSum function from the Julia version will automatically compress the MPO. Is that correct?

The original question is as follows.

Several years ago, the idea of fitting long-range interactions with sums of exponentials was discussed on the ITensor forum[ Constructing long-range hamiltonians with MPO's? - ITensor Support Q&A ]. I am opening this issue because I would like to learn more about the technical details of how this approach is implemented in ITensorMPS. My question may be a bit roundabout to explain, but I’ll try to make it as clear as possible below.

Example

Let me use an example to illustrate my question. If I want to construct the MPO of a simple Hamiltonian H=\sum_{i<j} \frac{1}{|j-i|^\alpha} \sigma_i^z \sigma_j^z, I can approximate it by the sum of several Hamiltonians

H = \sum_{m=1}^M x_m H_m = \sum_{m=1}^M x_m \left[ \sum_{i<j} \lambda_k^{|i-j|} \sigma_{i}^{z} \sigma_{j}^{z} \right] ，

where the coefficients {x_m} can be determined by an optimization procedure. The question is: how should one construct the Hamiltonian with exponentially decaying interactions using ITensor?

Direct coding approach

One straightforward way is to explicitly write code such as:

H_m = OpSum()

for i in 1:L
    for j in 1:i-1  
        xi = (i - 1) % L
        xj = (j - 1) % L
        Vij = lambda^abs(xi-xj)
        H_m += Vij, "Sz", i, "Sz", j
    end
end

H_MPO_m = MPO(Float64, H_m, sites)

Then, one could combine all {H_{\text{MPO},m}} with the coefficients {x_m} to obtain the MPO for H.

My question

However, according to Pirvu, Murg, Cirac, and Verstraete, New J. Phys. 12, 025012 (2010), the MPO representation of each H_m has a compact analytical form with small bond dimension.

• Does the MPO constructor in ITensor automatically produce such a compact form, or will it generate a more redundant representation?
• Is there a specialized function in ITensor (similar to fit_with_sum_of_exp in TeNPy) that can directly construct MPOs with exponentially decaying interactions?
• If not, does one need to manually construct the MPO by explicitly writing out the B matrices (as shown in the figure I attached below)?

[corbett5]

I am not aware of anything like fit_with_sum_of_exp. However, the algorithm ITensors uses in theory will produce an exact MPO of optimal bond dimension (up to floating point error). But this optimality is not guaranteed when performing approximate MPO construction. That said, in practice I would imagine that in this situation the approximation (by tuning the cutoff parameter of MPO(Float64, H_m, sites; cutoff)), would work very well. Or alternatively by calling truncate! on the fully “exactly” constructed MPO.

[miles]

Ben is correct that the algorithm used by the ITensor OpSum to MPO constructor / converter will use an approach that accomplishes the same thing as the method of fitting to exponentials. The exponential fitting can be viewed as a kind of special case of the approach we use, and as Ben said very often our code will produce the optimal bond dimension which is the main thing you want.