Tullio Order of Operations

General Usage
1

I’m attempting to use Tullio to perform a sum over matrices and running into difficulty. I need to take the sum of inverses and instead am taking the inverse of sums. Is there a way to properly write order of operations, or do I need to use a for-loop? Here is a minimal working example.

using Tullio, LoopVectorization, LinearAlgebra
GMatrix = Array{Float64}(undef,10,3,3)
x = fill(1,10)
Id = Matrix{Float64}(I, 3, 3)
H = Array{Float64}(undef,5,3,3)
filldata = collect(LinRange(-pi, pi, 5))
H[:,1,2] = filldata
H[:,1,3] = filldata
H[:,2,1] = filldata
H[:,3,1] = filldata

@tullio GMatrix[i,j,k] = inv(x[i]*Id[j,k] - H[l,j,k])
print(GMatrix[1,:,:])
2

I think when you use a function inside a Tullio expression it is applied to each element individually. You may need to consult the Tullio readme where it talks about “larger expressions” if you are trying to sum the inverse of matrices.

EDIT: Also if your real code involves these same x and Id variables, it seems like you could much more efficient write this because x[i]*Id[j,k] is always either 1 or 0. Probably you don’t need to construct those and can further simplify your Tullio expression and construction of GMatrix.

3

Yes, you can use the indices directly in the expression:

julia> H .= rand.();

julia> @tullio GMatrix[i,j,k] = inv(x[i]*Id[j,k] - H[l,j,k]);

julia> GMatrix ≈ @tullio alt[i,j,k] := 1/(x[i]*(j==k) - H[l,j,k])
true

Notice BTW that you do not initialise all entries in H, some are garbage from the undef.

This expression is the sum of inverses. The entire right hand side is being summed over l (since this does not appear on the left). The inv is acting on one number, not a matrix, before the sum.

julia> @tullio alt[i,j,k] := inv(x[i]*(j==k) - H[l,j,k]) verbose=true;
...
┌ Info: left index ranges
│   i = Base.OneTo(10)
│   j = Base.OneTo(3)
└   k = Base.OneTo(3)
┌ Info: reduction index ranges
└   l = Base.OneTo(5)

julia> @tullio s1 := inv(H[1,1,i])
4.997426578816242

julia> sum(inv.(H[1,1,:]))
4.997426578816242

julia> inv(sum(H[1,1,:]))
0.5423160388877815
1 Like
4

I see my mistake.

I would actually like to take the MATRIX inverse of the j,k indices. Is this possible in Tullio? I was looking at initializing static SMatrix, but as a novice, was struggling to figure out how to create filled arrays of them.

5

OK. No, then Tullio is not really the thing. It does not, except perhaps accidentally, think about slices of arrays.

To use broadcasting, you want something like inv.(x' .* Ref(I) .- eachslice(H, dims=1)), which makes a matrix of matrices, which you then glue together & sum on the right axis. And index-notation way is this:

julia> using TensorCast, LinearAlgebra

julia> @reduce mat[i,j,k] := sum(l) inv(x[i]*Ref(I) - H[l,:,:])[j,k];  # that you need Ref(I) is probably a bug, sorry

julia> summary(mat)
"10×3×3 Array{Float64, 3}"
3 Likes
6

Super generous of you to write TensorCast and then also answer questions. If I can ask a performance question, I saw some advice about using static arrays to increase inversion time. I’m now comparing

@cast B[k,l]{i,j} :=  x[k]*Id[i,j] - H[l,i,j]
@btime @reduce GMatrix[k] := sum(l) inv(B[k,l])
@btime @reduce mat[i,j,k] := sum(l) inv(x[i]*Id[:,:] - H[l,:,:])[j,k]

and it seems that casting each 3x3 matrix as a static array gives me about a 3x speed-up ( for a 1000 x 500 array)

98.897 ms (16 allocations: 34.40 MiB)
273.807 ms (1500052 allocations: 267.12 MiB)

Can I do any better than that?