Tullio broadcasting issue

New to Julia
, ,
1 
@tensor tmp1[m,s,n,p] := -λ[l,m,s]*ω[l,k]*λ[k,n,p] + λ[l,m,s]*g[l,n]*ones(2)[p])
@tullio tmp2[m,s,n,p] := -λ[l,m,s]*ω[l,k]*λ[k,n,p] + λ[l,m,s]*g[l,n]

tmp1 turns out to be different from tmp2. It seems that @tullio performs no broadcasting in the p-index, or is there something I missed?

2

Hi @swanchristmas, can you provide a minimal reproducable example? i.e. Please read: make it easier to help you

For what it’s worth, if I try to construct such an example, I see no difference here:

julia> using Tullio, TensorOperations

julia> let 
           N, P = 2,3
           A = rand(N)
           B = rand(P)
           C = rand(N)
           @tullio out1[n, p] := A[n] * B[p] + C[n]
           @tullio out2[n, p] := A[n] * B[p] + C[n] * ones(P)[p]
           @tensor out3[n, p] := A[n] * B[p] + C[n] * ones(P)[p]
          
           @info "" out1 out2 out3
       end
┌ Info: 
│   out1 =
│    2×3 Matrix{Float64}:
│     0.670981  0.719456  0.856811
│     0.577745  0.641696  0.822899
│   out2 =
│    2×3 Matrix{Float64}:
│     0.670981  0.719456  0.856811
│     0.577745  0.641696  0.822899
│   out3 =
│    2×3 Matrix{Float64}:
│     0.670981  0.719456  0.856811
└     0.577745  0.641696  0.822899

Is it possible there’s an error or typo somewhere in your code that’s causing the discrepancy you’re seeing?

3

Ah okay, nevermind, I can reproduce the discrepancy you’re seeing by introducing a contracted index:

julia> let 
           N, M, P = 2,4, 3
           A = rand(N, M)
           B = rand(M, P)
           C = rand(N)
           @tullio out1[n, p] := A[n, m] * B[m, p] + C[n]
           @tullio out2[n, p] := A[n, m] * B[m, p] + C[n] * ones(P)[p]
           @tensor out3[n, p] := A[n, m] * B[m, p] + C[n] * ones(P)[p]
          
           @info "" out1 out2 out3
       end
┌ Info: 
│   out1 =
│    2×3 Matrix{Float64}:
│     2.68322  2.8603   2.62903
│     4.80266  4.25346  4.36135
│   out2 =
│    2×3 Matrix{Float64}:
│     2.68322  2.8603   2.62903
│     4.80266  4.25346  4.36135
│   out3 =
│    2×3 Matrix{Float64}:
│     1.31961  1.49669  1.26541
└     2.12742  1.57822  1.68612

It seems like what’s happening here is that instead of calculating

map(Iterators.product(1:N, 1:P)) do (n, p)
    sum(1:M) do m
        A[n,m] * B[m,p]
    end + C[n]
end

Tullio is instead calculating

map(Iterators.product(1:N, 1:P)) do (n, p)
    sum(1:M) do m
        A[n,m] * B[m,p] + C[n]
    end
end

That’s rather unexpected, I didn’t know it did this. cc @mcabbott

I guess this behaviour makes some sense, since Tullio is meant to represent general reductions, but there should maybe be some warnings about this in the README.

1 Like
4

This is the correct diagnosis, Tullio always sums the entire RHS. It does not think that + is special. This issue has a longer description.

2 Likes
5

I’ve opened a PR trying to document this: Add note about reduction scope by MasonProtter · Pull Request #183 · mcabbott/Tullio.jl · GitHub