Contracting planar tensors with changeable legs
Hi Jutho and Lukas,
I am trying to contract two planar tensors, but the corresponding contracted legs are some variables that can be changed as an input parameter of a function. For example, I want to contract
@planar LT[(-1, -2, -3); -4] := L[-3; 1] * T[1; (-2, -1, -4)]
but in the following sense
ind1 = (-1, -2, -3)
ind2 = (-2, -1, -4)
@planar LT[ind1; -4] := L[-3; 1] * T[1; ind2]
Is there a way to achieve this?
Edit:
Found a planarcontract!
Sorry I am still confused with the conventions used in planarcontract!. According to the documentation TensorOperations.tensorcontract!, in my naive thinking, if I want to bend leg 2 and 3 in LT after contraction, it seems that I should write something like
V = Vect[IsingAnyon](:I => 1, :σ => 1) # error
# V = ℂ^2 # no error
T = rand(V ← V ⊗ V ⊗ V)
L = rand(V ← V)
ind1 = (-1, -2, -3)
ind2 = (-2, -1, -4)
LT = transpose(T, ((3,2,1), (4,)))
LT = TensorKit.planarcontract!(LT, L, ((1,), (2,)), T, ((1,), (2,3,4)), ((3,2,1), (4,)), 1, 0)
But this is not correct for sector being IsingAnyon. It seems that for non-bosonic sector, the planarcontract! will call reorder_indices, which restricts the length(pAB[1]) == length(pA[1]).
Find the tensorexpr in MPSKit.
It is true that the whole planar contractions are a bit less expressive than the corresponding general tensor contractions (even when the restrictions of being planar are fulfilled), in the sense that planarcontract! is more restrictive and there is no ncon equivalent for dynamic indices.
But let us focus on your use case: Is it only the ind1 and ind2 that are variable in the expression @planar LT[ind1; -4] := L[-3; 1] * T[1; ind2] ? Then definitely you can just write this as transpose(L * T, some_cyclic_permutation), which might be the most efficient approach.
Thanks for your explanation, here the example is only a special case. I want that the contracted leg of T to vary also. Now I am trying to use the tensorexpr like
l = MPSKit.tensorexpr(:L, (-1,), (:x,))
t = MPSKit.tensorexpr(:T, (:x,), (-4, -3, -2))
lt = MPSKit.tensorexpr(:LT, (-1,-2,-3), (-4,))
LT = eval(macroexpand(@__MODULE__, :(return @planar $lt := $l * $t)))
As Jutho mentioned, we are still lacking a bit of extra features for the planar contraction schemes. At some point, I tried to improve on this, but got sidetracked and didn't immediately need these features so that got stalled (https://github.com/Jutho/TensorKit.jl/pull/124).
Depending on the specifics of what you need, there are some ways to work around that.
If you really only need binary contractions, the easiest and cleanest solution would be to manually rewrite the contraction in terms of planarcontract! etc.
As you discovered, our implementation works somewhat interestingly, and for planarcontract!(C, A, pA, B, pB, pAB, alpha, beta), this should be interpreted as:
C = beta * C + alpha * transpose(transpose(A, pA) * transpose(B, pB), pAB).
Typically though, you want to specify the indices using labels instead of permutations, so you have to go from labels to permutations. To achieve this, we want something similar to the implementation from tensorcontract! here:
https://github.com/Jutho/TensorOperations.jl/blob/4d0d7f7abebc24d858031e1bb4caa0d7b4d2a29e/src/implementation/functions.jl#L205-L209
Ie, for labels IA, IB, IC you get pA, pB, pAB = contract_indices(IA, IB, IC).
However, note that there are multiple permutations that lead to the same indices.
For example, permuting the codomain indices of pA and relabeling the pAB indices accordingly should not have any effect on the result.
However, for planarcontract! we have to choose the specific option that is planar, which is what I attempted to implement here: https://github.com/Jutho/TensorKit.jl/blob/236e3f1b0ae90d70dd3a3f751374faf199151e34/src/planar/indices.jl#L86-L179
In fact, once that works, it is relatively easy to support ncon-like syntax, so this could then be used to implement larger networks.
Alternatively, when you need larger networks with sums etc, but not too many different ones, it might be more convenient to simply "hand-craft" the required expression and give that to @planar, similar to how you started above.
Because @planar is a macro, this gets evaluated at parse-time. Therefore, we need to build the expression before that. To achieve this, there are two options:
- We dynamically build an expression at runtime, and then use
@evalto evaluate it. - We use a generated function to build an expression that depends only on static information, and use that.
From what I can tell, your suggested implementation looks correct, and should evaluate to the correct answer. The use of eval however means that building the expression and finding the functions etc has to be called every time the function is invoked, which might become a performance issue at some point.
This second approach is something you might find examples of in MPSKit and PEPSKit, which has the benefit of being type-stable and doesn't have to re-evaluate the expressions. The drawback being that it is quite heavy on the compiler, so if you have many different contractions this also might not be ideal.
Thanks for your long reply! I'll read it carefully.