# Permutation average of an array

**URL:** <https://discourse.julialang.org/t/permutation-average-of-an-array/17753>\
**Category:** General Usage\
**Created:** [November 19, 2018, 8:34pm UTC](https://discourse.julialang.org/t/permutation-average-of-an-array/17753 "2018-11-19T20:34:20Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [November 19, 2018, 8:34pm UTC](https://discourse.julialang.org/t/permutation-average-of-an-array/17753/1 "2018-11-19T20:34:20Z")

</div>

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:

```julia
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)?