How to convert contractions from TensorOperations to ITensor ?

ITensor Julia Questions
1

I am trying to write a code for an iterative solver of a set of non linear equations. Doing this needs computation of a lot of Tensor Contractions which I had been doing using [TensorOperations.jl] . My code is running fine but I also want to try to do the same using other Tensor frameworks. [ITensor.jl] is what I want to try out now to accomplish the same jobs. Before a small function to calculate a contraction and add it to a given workspace tensor looked like below

using ITensors,TensorOperations
using ITensors: permute!
include("ccd-helper.jl")
nv,nocc = 2,4;
function addres(R2u::Array{Float64,4},nv,nocc,T2::Array{Float64,4})
    println("TensorOperations is being used")
    g_voov = deserialize("g_voov.jlbin")
    @tensor Trm2[a_1,a_2,i_1,i_2] := - g_voov[a_1,i_3,i_1,a_3] * T2[a_2,a_3,i_3,i_2]
    @tensor R2u[a_1,a_2,i_1,i_2] += Trm2[a_1,a_2,i_1,i_2]
    @tensor Trm5[a_1,a_2,i_1,i_2] := 2*g_voov[a_1,i_3,i_1,a_3] * T2[a_2,a_3,i_2,i_3]
    @tensor R2u[a_1,a_2,i_1,i_2] += Trm5[a_1,a_2,i_1,i_2]
    return R2u
end

Now, I cannot understand how to accomplish the same using [ITensor.jl] as elegantly as possible. The problem seems to stem from the fact that “Indices” are treated differently here and are just not dummy variables as they were when using [TensorOperations.jl]. I tried to implement the same in [ITensors.jl] as shown below

function addres(R2u::ITensor,nv,nocc,T2::ITensor)
    println("ITensor is being used")
    a_1 = Index(nv,"a_1")
    a_2 = Index(nv,"a_2")
    i_1 = Index(nocc,"i_1")
    i_2 = Index(nocc,"i_2")
    i_3 = Index(nocc,"i_3")
    a_3 = Index(nv,"a_3")
    g_voov = ITensor(deserialize("g_voov.jlbin"),a_1,i_3,i_1,a_3)
    T2 = replaceinds(T2, inds(T2) => (a_2, a_3, i_3, i_2))
    R2u = replaceinds(R2u, inds(R2u) => (a_1, a_2, i_1, i_2))
    Trm2 = ITensor(a_1,a_2,i_1,i_2)
    Trm5 = ITensor(a_1,a_2,i_1,i_2)
    Trm2 = g_voov * T2
    Trm2 = permute(Trm2,a_1,a_2,i_1,i_2)
    R2u += Trm2
    Trm5 = 2*g_voov * T2
    Trm5 = permute(Trm5,a_1,a_2,i_1,i_2)
    R2u += Trm5
    return R2u
    println("All Good")
end

But when trying to compare the outputs of the functions using the code below

initialize_cc_variables()
T2 = initialize_t2_only();
a = Index(nv,"a");
b = Index(nv,"b");
i = Index(nocc,"i");
j = Index(nocc,"j");
iT2 = ITensor(initialize_t2_only(),a,b,i,j);
R2u = zeros(Float64,nv,nv,nocc,nocc);
iR2u = ITensor(R2u,a,b,i,j);
R2u = addres(R2u,nv,nocc,T2)
iR2u = addres(iR2u,nv,nocc,iT2)
@show iR2u
display(R2u)

They seem to be giving different results. Is there any way to accomplish this code conversion without a significant amount of work ? (If you want to run the code exactly as shown above you can use the attached file which has been included at the start of the code, that file contains all the helper functions. The file t1.jl contains the current faulty code.)
ccd-helper.jl (26.6 KB)
t1.jl (1.7 KB)

2

I think that is basically the right strategy, I’m not sure if we would be able to do a much better job of debugging your code than you would, I don’t see any obvious mistakes but I only glanced through it.

Part of the struggle you are having may just stem from the fact that designing an algorithm and code around temporary index labels like in TensorOperations.jl vs. permanent index labels like those of ITensors.jl is fundamentally quite different, and requires you to think through the design of your algorithm in a different way. However, we believe that once you “buy into” permanent index labels it can make it easier to develop more general tensor network algorithms with less code.

1 Like
3

Also, note that the calls to permute shouldn’t be necessary since ITensor automatically permutes the dimensions based on the Index labels.

4

There are other simplifications to your code that I can see, for example initializing Trm2 = ITensor(a_1,a_2,i_1,i_2) isn’t necessary since it is just overwritten by Trm2 = g_voov * T2. So I think you should try to go through the ITensor code and simplify things down, that could make it easier to debug.

5

Also, lines like T2 = replaceinds(T2, inds(T2) => (a_2, a_3, i_3, i_2)) make me pretty suspicious, perhaps that can be achieved in a different way. You could pretty easily be setting the indices incorrectly, since if T2 was input with a different index ordering that would make that code incorrect. So instead of just trying to translate from TensorOperations to ITensor, maybe you should step back a little bit and think more critically about what you are trying to achieve with the code (i.e. think about redesigning it from scratch with ITensor’s design philosophy around labelled indices in mind).

1 Like
6

Two other suggestions that may help:

  1. building on Matt’s point about rethinking / redesigning, in ITensor you should think about how each ITensor maps from an input space (defined by a collection of indices) to an output space. If the input and output space are ‘physically’ or the same, then a good strategy can be to have the output space consist of the same indices but with the ‘primelevel’ of those indices raised, say to the value 1. You can make the prime of an index by writing i' or higher levels by i'', i''', or prime(i, level). The nice thing about primes is that it is easy to map back to the original space by calling functions like noprime on a resulting tensor. If the spaces are physically different, you should have the output space consist of a distinct set of indices and not use primes of the input indices. Only in more complicated situations, though hopefully yours will be the two cases above, it may be helpful to replace indices or map them to other indices using delta tensors.

  2. You may be doing this already, but drawing tensor diagrams is enormously helpful when using ITensor and setting up calculations. The whole design of ITensor is based on the idea that code should contain nearly the same information that tensor diagrams contain, no more and no less.

1 Like
7

Turns out that is exactly where the problem was ! I should have permuted the indices before the second tensor contraction. I also took into account the advice about the lack of need to manually permute the Indices to arrive at the following code

function addres_gvoov(R2u::Array{Float64,4},nv,nocc,T2::Array{Float64,4})
    g_voov = deserialize("g_voov.jlbin")
    @tensor Trm2[a_1,a_2,i_1,i_2] := - g_voov[a_1,i_3,i_1,a_3] * T2[a_2,a_3,i_3,i_2]
    @tensor R2u[a_1,a_2,i_1,i_2] += Trm2[a_1,a_2,i_1,i_2]
    @tensor Trm5[a_1,a_2,i_1,i_2] := 2*g_voov[a_1,i_3,i_1,a_3] * T2[a_2,a_3,i_2,i_3]
    @tensor R2u[a_1,a_2,i_1,i_2] += Trm5[a_1,a_2,i_1,i_2]
    return R2u
end

function addres_gvoov(R2u::ITensor,nv,nocc,T2::ITensor)
    i_3 = Index(nocc,"i_3")
    a_3 = Index(nv,"a_3")
    g_voov = ITensor(deserialize("g_voov.jlbin"),a_1,i_3,i_1,a_3)
    Trm2 = -g_voov * replaceinds(T2,inds(T2)=>(a_2,a_3,i_3,i_2))
    R2u += Trm2
    Trm5 = 2*g_voov * replaceinds(T2,inds(T2)=>(a_2,a_3,i_2,i_3))
    R2u += Trm5
    return R2u
end

These are the two versions of the same function that deal with TensorOperations and ITensor respectively. And now both give the expected correct result. However, as rightly pointed out, it is probably a good idea to rethink the structure of my algorithm and how to modify it to be compatible with ITensors.jl instead of just blindly trying to translate the current version.

1 Like
8

Thanks for the update. Note that as a reference for you and others reading this, I would recommend avoiding calls like:

replaceinds(T2,inds(T2)=>(a_2,a_3,i_3,i_2))

where you replace based on inds(T2) since that can lead to bugs where the new indices don’t properly correspond to the memory layout of the underlying data if you input an ITensor T2 with the same indices but different index ordering. Anyway, glad to see you were able to both simplify your ITensor code and catch the bug you were seeing.

9

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