Permutation average of an array

General Usage
1

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