Average over the columns of the matrix in Julia

user_19104784 0
Viewed 97

I have a big matrix with float entries of the form

[ a b c d 
  e f g h
  i j k l
  m n o p ]

some of the values are outliers, so I wanted to average each of the entries with values with its recent k entries in the corresponding column and preserve the shape. In other words to have something like this for k = 3:

[        a                  b                   c                  d 
      (e + a)/2          (f + b)/2          (g + c)/2          (h + d)/2
    (e + a + i)/3      (f + b + j)/3      (g + c + k)/3      (h + d + l)/3
    (e + i + m)/3      (f + j + n)/3      (g + k + o)/3      (h + l + p)/3   ] 


 etc.
julia matrix
2 Answers

You can do this with RollingFunctions and mapslices:

julia> a = reshape(1:16, 4, 4)
4×4 reshape(::UnitRange{Int64}, 4, 4) with eltype Int64:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> using RollingFunctions

julia> mapslices(x -> runmean(x, 3), a, dims = 1)
4×4 Matrix{Float64}:
 1.0  5.0   9.0  13.0
 1.5  5.5   9.5  13.5
 2.0  6.0  10.0  14.0
 3.0  7.0  11.0  15.0
user_24998920

I didn't know about RollingFunctions, but a regular loop is 4X faster. I'm not sure if it's some kind of type instability caused by mapslices?.

function runmean(a,W)
    A = similar(a)
    for j in axes(A,2), i in axes(A,1)
        l = max(1, i-W+1)
        A[i,j] = mean(a[k,j] for k=l:i)
    end
    A 
end

Testing yields:

using RollingFunctions 
@btime mapslices(x -> runmean(x, 3), A, dims = 1) setup=(A = rand(0.0:9,1000,1000))
@btime runmean(A,3) setup=(A = rand(0.0:9,1000,1000))

  15.326 ms (10498 allocations: 23.45 MiB)
   4.410 ms (2 allocations: 7.63 MiB)
user_39431700
Related