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,:,:])
```

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`

.

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

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.

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

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?