Permute Multidimensional Array without allocation

General Usage
1

Hi guys,

I am trying to permute/shuffle columns of a multidimensional array, an example is here:

using BenchmarkTools

const arra = zeros(Int, 4, 3, 2)
arra[:,:, 1] = 1:12
arra[:,:, 2] = 13:24

shuffle_col = [3, 1, 1]

@benchmark arra[:, shuffle_col, :] #solution 1 - works but allocates
@benchmark @view arra[:, shuffle_col, :] #solution 2 - much faster but arra is now a subarray

Both solutions work but the first one allocates, and the second one forces me to continue to work with a subarray (which I would like to avoid). I tried to put everything in a function in a loop, but my solution takes longer than the two above for bigger arrays (and needs a copy to work to begin with):

function permute_col_array!(array::AbstractArray, permutation::Vector{Int}) #, dimension::Int
        temp = copy(array)
        for (i, permute) in enumerate(permutation)
            array[:,i,:] = temp[:,permute,:]
        end
end

@benchmark permute_col_array!(arra, shuffle_col)

Does anyone here know if I can permute multidimensional arrays without allocation? I would gain significant improvements in functions if I could do that. There is a base function permute! in Julia but it only works on vector elements, not on whole columns.

Best regards,

1 Like
2

Did you test your function? I think it must overwrite some elements before they are used.

3

Thanks it should work now! But my try is now even worse than before as I needed the copy to not overwrite it in a loop.

1 Like
4

I’ve been stuck trying to solve this problem too! I have a use case for this kind of mechanism that I haven’t been able to crack: Complex Numbers Wrappers

5

AFAIK permuting arrays in place in general is not a trivial problem. If allocation (and not copying) is the costly part per se, perhaps you should just preallocate for output.

6

Preallocation is unfortunately not possible as I need to shuffle the array in a loop and the object of interest is the array itself ( a filtering problem).

7

Could you elaborate on this ? In place permutations are really difficult to optimize. In 3D you have to compute some permutation cycles… In general out of place (with pre-allocation + ptr swap) permutations are much easier to implement and perform better. If it is simple permutedims you can tile.

1 Like
8

Here would be a dummy code of what I have to do:

#Input
observations = Matrix(T,P)

#Algorithm
function myproblem(observations, k, ....)
    #initialize trajectory
    trajectories = Matrix(T,M,N) 
    for i in 1:T
        #get column order
        weights = somefunction(observations[1:i,:], trajectories[:,(i-k):i,:] )
        #Reshuffle trajectories
        trajectories = trajectories[:,weights,:]
    end
    return trajectories
end

I can do preallocation for trajectories, but so far I have not found a way for trajectories[:,weights,:] to not allocate a new array at each iteration unfortunately. I also cannot simply make a output container as I only know the output after looping through all columns in this case.

Also, looping column by column performs better here than row by row, as I fill the trajectories via columns in this case.

9

(Realize that array[:,i,:] = temp[:,permute,:] allocates an a extra copy of the temp[:,permute,:] slice unless you use @views.)

There actually is an undocumented Base.permutecols!!(A, p) function that permutes the whole columns of a matrix A by p (also overwriting the permutation array p). It’s not very long, and you could base something for 3d arrays off it pretty easily: https://github.com/JuliaLang/julia/blob/99d1d6723d72d43074b5d2765e63d08439989e8e/base/combinatorics.jl#L66-L95

4 Likes
10

Thank you a lot! I tried to get it to run in Julia 1.2, but failed so far. I will think about how I could extend it to more dimension without copying at least p.