Fixed bond dimension operations

ITensor Julia Questions
,
1

Hi,

I wanted to ask for help as I am trying to reproduce the technique of using Chebysev polynomials of MPOs of a certain Hamiltonian to estimate the DOS from arXiv:1909.01398v1.

The Chebysev polynomials are computed using the recursive relation:

T_{n+1}(H) = 2H T_n(H) - T_{n-1}(H), \quad \text{with } T_0(H) = I \text{ and } T_1(H) = H.

where H is a MPO of the Hamiltonian. As the authors say in the paper, you need to fix the bond dimension if you want to sum these MPOs as they must be the same kind of tensor, so I guess one should try to “extend” the bond dimension of the identity MPO and first order polynomials.

I am having trouble to find a way to get I and H to have the desired fixed tensor structure to perform those operations in the documentation or by myself.

If someone could help me I would really appreciate it.

(btw: Sorry if the question is too naive, I am still getting used to ITensor and TNs by myself)

2

Hi Jorge,
The ITensor toolkit ought to be a good fit for what you’re trying to do.

To make H as a tensor network, you can use the “OpSum system” which converts data about products of operators (somewhat similar to Latex input) into MPO tensor networks automatically. You can see an example of the use of this system at this link:
https://itensor.github.io/ITensors.jl/dev/tutorials/DMRG.html
(Just the part leading up to and including the line that reads H = MPO(os,sites).)

Then for the identity operator in the same Hilbert space, you could make this as:

os_I = OpSum()
os_I += "Id",1
I = MPO(os_I,sites)

where adding just the "Id",1 operator is enough because our system understands it to mean, implicitly, the identity operator on site 1 times the identity operator on site 2, etc. times identity on all the other sites. (This is generally implied in all physics operator notation, that any sites not carrying a non-trivial operator are actually just a product with the identity operator on that site.)

Then once you have made H and I as MPO’s, you can perform operations with them. Note that for the usual Chebyshev MPO technique, if I recall correctly, you do not actually want to compute T_n(H) itself, but the action of T_n(H) on a reference state (represented as an MPS) which is much more efficient.

1 Like
3

I’ll try to implement it that way, thank you very much!

4

Thanks again for the suggestion. I developed a function that reproduces the Chebysev polynomials as functions:

function Chebyshev_polynomial_mpo(sites, H::MPO, BondDimension::Int, MinBondDimension::Int, precision::Number, order::Int)
    T_list = Vector{Function}(undef, order + 1)  # Create an empty list of functions

    os_I = OpSum()
    os_I += "Id", 1
    I = MPO(os_I, sites)

    T_list[1] = x -> truncate!(I; maxdim=BondDimension, mindim=MinBondDimension, cutoff=precision)  # T_0(x) = I
    T_list[2] = x -> truncate!(H; maxdim=BondDimension, mindim=MinBondDimension, cutoff=precision)  # T_1(x) = H

    for k in 2:order
        T_list[k + 1] = x -> truncate!(2 * apply(T_list[k](x), H; maxdim=BondDimension, cutoff=precision) - T_list[k - 1](x); maxdim=BondDimension, mindim= MinBondDimension, cutoff=precision)
    end

    return T_list
end

and it works without giving error if applied to H, (an MPO represntation of the LR Ising chain). If I print the function applied to H I obtain the following:


julia> Chebyshev_polynomial_mpo(sites, H, 10, 10, 1e-6, 5)[6](H)
MPO
[1] ((dim=2|id=745|"S=1/2,Site,n=1")', (dim=2|id=745|"S=1/2,Site,n=1"), (dim=4|id=245|"Link,l=1"))
[2] ((dim=2|id=454|"S=1/2,Site,n=2")', (dim=2|id=454|"S=1/2,Site,n=2"), (dim=4|id=17|"Link,l=2"), (dim=4|id=245|"Link,l=1"))
[3] ((dim=2|id=110|"S=1/2,Site,n=3")', (dim=2|id=110|"S=1/2,Site,n=3"), (dim=4|id=17|"Link,l=2"))


julia> Chebyshev_polynomial_mpo(sites, H, 10, 10, 1e-6, 5)[5](H)
MPO
[1] ((dim=2|id=745|"S=1/2,Site,n=1")', (dim=2|id=745|"S=1/2,Site,n=1"), (dim=4|id=631|"Link,l=1"))
[2] ((dim=2|id=454|"S=1/2,Site,n=2")', (dim=2|id=454|"S=1/2,Site,n=2"), (dim=4|id=94|"Link,l=2"), (dim=4|id=631|"Link,l=1"))
[3] ((dim=2|id=110|"S=1/2,Site,n=3")', (dim=2|id=110|"S=1/2,Site,n=3"), (dim=4|id=94|"Link,l=2"))


julia> Chebyshev_polynomial_mpo(sites, H, 10, 10, 1e-6, 5)[4](H)
MPO
[1] ((dim=2|id=745|"S=1/2,Site,n=1")', (dim=2|id=745|"S=1/2,Site,n=1"), (dim=4|id=310|"Link,l=1"))
[2] ((dim=2|id=454|"S=1/2,Site,n=2")', (dim=2|id=454|"S=1/2,Site,n=2"), (dim=4|id=451|"Link,l=2"), (dim=4|id=310|"Link,l=1"))
[3] ((dim=2|id=110|"S=1/2,Site,n=3")', (dim=2|id=110|"S=1/2,Site,n=3"), (dim=4|id=451|"Link,l=2"))


julia> Chebyshev_polynomial_mpo(sites, H, 10, 10, 1e-6, 5)[3](H)
MPO
[1] ((dim=2|id=745|"S=1/2,Site,n=1")', (dim=2|id=745|"S=1/2,Site,n=1"), (dim=4|id=501|"Link,l=1"))
[2] ((dim=2|id=454|"S=1/2,Site,n=2")', (dim=2|id=454|"S=1/2,Site,n=2"), (dim=3|id=320|"Link,l=2"), (dim=4|id=501|"Link,l=1"))
[3] ((dim=2|id=110|"S=1/2,Site,n=3")', (dim=2|id=110|"S=1/2,Site,n=3"), (dim=3|id=320|"Link,l=2"))


julia> Chebyshev_polynomial_mpo(sites, H, 10, 10, 1e-6, 5)[2](H)
MPO
[1] ((dim=2|id=745|"S=1/2,Site,n=1")', (dim=2|id=745|"S=1/2,Site,n=1"), (dim=4|id=44|"Link,l=1"))
[2] ((dim=2|id=454|"S=1/2,Site,n=2")', (dim=2|id=454|"S=1/2,Site,n=2"), (dim=4|id=846|"Link,l=2"), (dim=4|id=44|"Link,l=1"))
[3] ((dim=2|id=110|"S=1/2,Site,n=3")', (dim=2|id=110|"S=1/2,Site,n=3"), (dim=4|id=846|"Link,l=2"))


julia> Chebyshev_polynomial_mpo(sites, H, 10, 10, 1e-6, 5)[1](H)
MPO
[1] ((dim=2|id=745|"S=1/2,Site,n=1")', (dim=2|id=745|"S=1/2,Site,n=1"), (dim=2|id=464|"Link,l=1"))
[2] ((dim=2|id=454|"S=1/2,Site,n=2")', (dim=2|id=454|"S=1/2,Site,n=2"), (dim=2|id=747|"Link,l=2"), (dim=2|id=464|"Link,l=1"))
[3] ((dim=2|id=110|"S=1/2,Site,n=3")', (dim=2|id=110|"S=1/2,Site,n=3"), (dim=2|id=747|"Link,l=2"))

Here comes my question: If the link indices are of different dimension depending on the order, how it is possible that the sum operation performs correctly? What is the code doing that I am missing?

Thanks again for the help, I really appreciate it!

5

You don’t need the link dimensions of two MPOs to be the same in order to sum them. There are multiple algorithms to sum MPOs, the simplest one conceptually is that you direct sum each MPO tensor over the link dimensions.

6

mm I see. I was think more of a trying to do a sum over tensors of the same kind(with the bond dimension fixed) from what I understood of [1909.01398] Probing thermalization through spectral analysis with matrix product operators page 2, last paragraph. Maybe it was a misunderstanding of mine.

Therefore if you have two equal MPOs when performing the operation “+” you sum each element with its corresponding one, but if those MPOs are not exactly the same “+” means direct tensor sum?

7

If the two MPOs happen to be the same or share bases, after you direct sum them you can then compress them back down to a smaller bond dimension (if they are exactly the same, you can compress them back down the original bond dimension of each of the MPOs you are summing).

Compressing MPOs that represent Hamiltonians (or more generally, sums of quasi-local operators) is a little bit subtle, see for example [1909.06341] Local Matrix Product Operators: Canonical Form, Compression, & Control Theory, but if the MPOs don’t have too many sites you can use ITensorMPS.truncate!. That is the algorithm that is implemented when you call +(a, b; alg="directsum", cutoff=1e-12) for MPOs a and b, see the docs. You can also try the alg="densitymatrix" backend which is similar to the algorithm described here.

8

Also, when summing MPOs one is never summing each element with its corresponding one (ie one is never doing an element wise sum). That kind of sum would not correspond to forming the sum of two linear operators, which is what is wanted.

9

I’ll give it a look, thank you very much!

10

Thanks for the clarification, I appreciate it!

11

This topic was automatically closed 10 days after the last reply. New replies are no longer allowed.