# Tullio Order of Operations
**URL:** https://discourse.julialang.org/t/tullio-order-of-operations/85538
**Category:** General Usage
**Tags:** tullio
**Created:** [August 9, 2022, 2:31pm UTC](https://discourse.julialang.org/t/tullio-order-of-operations/85538 "2022-08-09T14:31:27Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@andrewkhardy](https://discourse.julialang.org/u/andrewkhardy)
#### Post date: [August 9, 2022, 2:31pm UTC](https://discourse.julialang.org/t/tullio-order-of-operations/85538/1 "2022-08-09T14:31:27Z")
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.
```julia
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,:,:])
```
---
### Author:  [@tbeason](https://discourse.julialang.org/u/tbeason)
#### Post date: [August 9, 2022, 4:07pm UTC](https://discourse.julialang.org/t/tullio-order-of-operations/85538/2 "2022-08-09T16:07:35Z")
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`.
---
### Author:  [@mcabbott](https://discourse.julialang.org/u/mcabbott)
#### Post date: [August 9, 2022, 6:27pm UTC](https://discourse.julialang.org/t/tullio-order-of-operations/85538/3 "2022-08-09T18:27:14Z")
> [@tbeason](#):
>
> Probably you don’t need to construct those
Yes, you can use the indices directly in the expression:
```julia
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`.
> [@andrewkhardy](#):
>
> I need to take the sum of inverses and instead am taking the inverse of sums.
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
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
```
---
### Author:  [@andrewkhardy](https://discourse.julialang.org/u/andrewkhardy)
#### Post date: [August 9, 2022, 7:15pm UTC](https://discourse.julialang.org/t/tullio-order-of-operations/85538/4 "2022-08-09T19:15:17Z")
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.
---
### Author:  [@mcabbott](https://discourse.julialang.org/u/mcabbott)
#### Post date: [August 9, 2022, 7:30pm UTC](https://discourse.julialang.org/t/tullio-order-of-operations/85538/5 "2022-08-09T19:30:47Z")
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
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}"
```
---
### Author:  [@andrewkhardy](https://discourse.julialang.org/u/andrewkhardy)
#### Post date: [August 10, 2022, 1:38pm UTC](https://discourse.julialang.org/t/tullio-order-of-operations/85538/6 "2022-08-10T13:38:06Z")
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
```julia
@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)
```julia
98.897 ms (16 allocations: 34.40 MiB)
273.807 ms (1500052 allocations: 267.12 MiB)
```
Can I do any better than that?