I have the following problem: I’m given a square matrix `M in M_{N,N}`

and a left/right group action on `M_{N,N}`

given by explicit permutations of the row/column spaces. I need to find a fixed point of this action, i.e. an average of `M`

under this action. Code:

```
begin
N = 500
perms = [randperm(N) for _ in 1:3*N];
M = Matrix(sprand(N,N, 0.1));
end
julia> function perm_avg_loop(result, P, preps)
lp = size(P,1)
for p in preps
@inbounds for j in 1:lp
k = p[j]
for i in 1:lp
result[i,j] += P[p[i], k]
end
end
end
return result
end
perm_avg_loop (generic function with 1 method)
julia> function perm_avg_view(result, P, preps)
lp = size(P,1)
for p in preps
result .+= view(P, p, p)
end
return result
end
perm_avg_view (generic function with 1 method)
julia> @time pl = perm_avg_loop(zeros(size(M)), M, perms);
0.564013 seconds (23.23 k allocations: 3.038 MiB)
julia> @time pl = perm_avg_loop(zeros(size(M)), M, perms);
0.538265 seconds (7 allocations: 1.908 MiB)
julia> @time pv = perm_avg_view(zeros(size(M)), M, perms);
0.776470 seconds (45.63 k allocations: 3.796 MiB)
julia> @time pv = perm_avg_view(zeros(size(M)), M, perms);
0.725294 seconds (4.51 k allocations: 2.068 MiB)
julia> all(pl .≈ pv)
true
```

Is there a better way of doing this for larger (in practice `N~4000`

or so)?