Confusing/bugged? index convention for braiding tensors.

[Gertian]

In the following code, which contracts two tensors in various ways depicted below, I'd expect all diagrams to be equivalent.

Nevertheless I get a spacemismatch for d-type diagram for reasons that totally evade me. I'm not sure if this is intended behavior or a bug but either way I think it deserves clarification.

using TensorKit, TensorOperations, Test

sp = ℂ^2
T = TensorMap(rand, ComplexF64, sp, sp)

#the normal diagram
a = @planar opt=true T[a;b]*T[b;a]

#the diagram with one winding, both conventions for the labeling of the braiding
b = @planar opt=true T[d;a]*T[b;c]*τ[c b;a d]
c = @planar opt=true T[d;a]*T[b;c]*τ[a c;d b]
@test a ≈ b ≈ c

#the diagrams with one pulling trough --> two taus --> both conventions lead to 4 diagrams
d = @planar opt=true T[f;a]*T[c;d]*τ[d b;c e]*τ[e b;a f]   #this one gives a space mismatch for reasons that evade me !
e = @planar opt=true T[f;a]*T[c;d]*τ[d b;c e]*τ[a e;f b]
f = @planar opt=true T[f;a]*T[c;d]*τ[c d;e b]*τ[e b;a f]
g = @planar opt=true T[f;a]*T[c;d]*τ[c d;e b]*τ[a e;f b]

Resulting error for the d-type diagram :

ERROR: SpaceMismatch("(ℂ^2 ⊗ ℂ^2) ← ProductSpace{ComplexSpace, 0}() ≠ permute((ℂ^2 ⊗ ℂ^2) ← (ℂ^2 ⊗ ℂ^2)[(1, 2), (3, 4)] * ℂ^2 ← ℂ^2[(2, 1), ()], ((1, 2), ())")
Stacktrace:
 [1] planarcontract!(::TensorMap{ComplexSpace, 2, 0, Trivial, Matrix{ComplexF64}, Nothing, Nothing}, ::TensorKit.BraidingTensor{ComplexSpace, Matrix{Float64}}, ::Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}, ::TrivialTensorMap{ComplexSpace, 1, 1, Matrix{ComplexF64}}, ::Tuple{Tuple{Int64, Int64}, Tuple{}}, ::Tuple{Tuple{Int64, Int64}, Tuple{}}, ::One, ::Zero)
   @ TensorKit ~/.julia/dev/TensorKit/src/tensors/braidingtensor.jl:232
 [2] _planarcontract!(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any)
   @ TensorKit ~/.julia/dev/TensorKit/src/planar/postprocessors.jl:57
 [3] top-level scope
   @ REPL[652]:1

[Gertian]

Some further info : I used @expand from MacroTool.jl to expand the macro for the failing contraction. The result (after some cleaning) is :

(sealion, nightingale) = (numout(T), numin(T))
    numout(T) == 1 && (TensorKit.numin)(T) == 1 || throw(TensorOperations.IndexError("incorrect number of input-output indices: (1, 1) instead of " * string((sealion, nightingale)) * " for T."))
    (dugong, fish) = (numout(B), numin(B))
    numout(B) == 1 && (TensorKit.numin)(B) == 1 || throw(TensorOperations.IndexError("incorrect number of input-output indices: (1, 1) instead of " * string((dugong, fish)) * " for B."))

    weasel = TensorKit.BraidingTensor(space(B, 1), (space(B, 2))')  #this is tau2 ##To get this to work correctly we need to replace space(B, 1) with space(B, 1)' !
    vicuña = TensorKit.BraidingTensor(space(B, 1), (space(B, 2))')  #this is tau1

    #define the scalartype of the contracted thing
    porpoise = TensorOperations.promote_contract(TensorOperations.scalartype(weasel), TensorOperations.scalartype(T))
    #contract tau2=weasel with T 
    heron = TensorOperations.tensoralloc_contract(porpoise, ((1, 2), ()), weasel, ((1, 2), (3, 4)), :N, T, ((2, 1), ()), :N, false)
    heron = TensorKit._planarcontract!(heron, ((1, 2), ()), weasel, ((1, 2), (3, 4)), T, ((2, 1), ()), One(), Zero())

    #again find the correct scalartype
    okapi = TensorOperations.promote_contract(TensorOperations.scalartype(B), TensorOperations.scalartype(vicuña))
    #contract tau1=vicuña with B
    antelope = TensorOperations.tensoralloc_contract(okapi, ((), (2, 1)), B, ((), (1, 2)), :N, vicuña, ((3, 1), (2, 4)), :N, false)
    antelope = TensorKit._planarcontract!(antelope, ((), (2, 1)), B, ((), (1, 2)), vicuña, ((3, 1), (2, 4)), One(), Zero())

    #again find the correct scalartype
    alpaca = TensorOperations.promote_contract(TensorOperations.scalartype(antelope), TensorOperations.scalartype(heron))
    #contract T*tau and B*tau
    monkey = TensorOperations.tensoralloc_contract(alpaca, ((), ()), antelope, ((), (1, 2)), :N, heron, ((1, 2), ()), :N, false)
    monkey = TensorKit._planarcontract!(monkey, ((), ()), antelope, ((), (1, 2)), heron, ((1, 2), ()), One(), Zero())

    #we don't need T*tau and B*tau anymore, free them up from memory !
    TensorOperations.tensorfree!(antelope)
    TensorOperations.tensorfree!(heron)

    #free and typeset all other things. 
    rabbit = TensorOperations.scalartype(monkey)
    cattle = TensorOperations.tensoralloc_add(rabbit, ((), ()), monkey, :N, true)
    cattle = TensorKit._planaradd!(cattle, ((), ()), monkey, One(), Zero())
    TensorOperations.tensorfree!(monkey)
    lion = TensorOperations.tensorscalar(cattle)
    TensorOperations.tensorfree!(cattle)
    d = lion

Which leads to a spacemismatch because weasel=tau2 has incorrectly inferred spaces ! Indeed, we can make the generated code run fie by replacing weasel = TensorKit.BraidingTensor(space(B, 1), (space(B, 2))') with weasel = TensorKit.BraidingTensor(space(B, 1)', (space(B, 2))') .

I think that the problem lies within the _construct_braidingtensors(ex::Expr) function. It seems to correctly infers which spaces are associated to the various tau tensors but then still constructs a braiding tensor with the wrong spaces. I can't figure out how to patch this though.

[lkdvos]

hmm, that's interesting. I think it should have pulled the spaces for weasel from T, and somehow this didn't happen. I'll look into it tomorrow!

[Gertian]

@lkdvos, I went trought the code a bit and I think it correctly recognizes where to pull the spaces from but it doesn't seem to know if they should be conjugate or not.

Thanks for looking into it !

[lkdvos]

This should be fixed and tagged as of v0.12.1

[Gertian]

For the following contraction : @diffset @planar opt=true edges[WEST, row, cop][L1 E2 e4;L4] := envs.edges[WEST, row, col][L2 W4 w4; L3] peps_above[row, col][P1; N1 E1 S1 W1] conj(peps_below[row, col][P3; n1 e1 s1 w1]) P[L3 N2 n4; L4] Q[L1;L2 S3 s3] τ[w3 W3;W4 w4]τ[n3 N1;N2 n4]τ[n2 W2;W3 n3]τ[n1 w2;w3 n2]τ[e2 S2;S3 e1] τ[e3 s2;s3 e2]τ[e4 E1;E2 e3]τ[S1 s1;s2 S2]τ[P2 W1;W2 P1]τ[P3 w1;w2 P2] that you can see below : I still get a non-planar error when I remove the opt=true flag.

To be fair this is not really problematic but is there any change the opt=true interferes with automatic differentiation as this would pose a problem later down the line...

[lkdvos]

The opt=true flag tells TensorOperations to optimize the contraction order. If you don't specify this, @planar will just multiply pairs in the order you provide them, which in your case would lead to a non-planar contraction. In particular, planarity is not invariant under order changes :wink: .

This should not be relevant for automatic differentiation, but is definitely highly relevant for performance. In particular, at a later stage you might even want to consider specifying which indices are large, and which are not

[Gertian]

Thanks for your swift response !

How would this lead to a non planar diagram ? All crossings have been eliminated with tau tensors so I don't see how this is affected by the order of contraction. I could see that the order is ill defined due to my usage of names like w1, w2, W1, W2 etc but this doesn't seem to bother other contractions in in my code.

[lkdvos]

Consider the following diagram, and compare contracting orders: (b, a, b' c) versus (a, ...).

@planar D[-1; -2] := A[-1; a b c] * B[b; b'] * C[a b' c; -2]

In particular, any order starting with contracting index a ot c does not work, because you can't connect A and C while B is still squeezed in between.

In this case however, because you're not using NCON index notation, but rather just the letters, @planar decides to just multiply in the order you supplied:

@planar D[-1; -2] := (A[-1; a b c] * B[b; b']) * C[a b' c; -2]

which would work, however this would not:

@planar D[-1; -2] := (A[-1; a b c] * C[a b' c; -2]) * B[b; b']) 

[lkdvos]

As an aside, I am not sure if the "optimal" contraction order is guaranteed to be planar (I honestly don't think it is), but at least intuitively it feels like the example above would be less likely to occur in an optimal contraction path.

[Gertian]

In particular, any order starting with contracting index a ot c does not work, because you can't connect A and C while B is still squeezed in between.

Meaning that it has to permute A and C to be next to each other and that it introduces crossings by doing so ?

[lkdvos]

It's a bit hard to show you an example, because I want to say that it is not possible :p but basically, when you contract two tensors, you need to permute them such that all contracted indices are on one side, and all uncontracted on the other, so you cannot achieve this when contracting A and C, because index b is in the way.

[Gertian]

Ok, that clarifies it ! Thank you !