Inplace Dimension Permutation in O(1) space

But why? Why not! it was fun. Okay, with that out of the way...

I wanted to ramble a bit about the permutation of a multidimensional array in $O(1)$ space. That is, without extra auxiliary buffers^[1], using only the actual memory of the array itself.

It's a very well known and studied problem, so I'm not claiming any sort of novelty. It started just as a simple curiosity, an interesting problem to try and play around with (without looking at the solution beforehand). What took me by surprise was the depth of the rabbit hole. It was more fun than anticipated and touched on some pretty interesting math.

^[1]"swept under the rug" would be more precise...

This is the classic case of doing a simple thing in a complicated way. Because boredom.

Also, small note, this is not going to be super formal for exposition's sake. Ok, let's go!

^[2] In particular all code here was run on Julia 1.12.6

^[3] honeypot em-dashes for all those pesky AI detection bots.

Matrices and Memory

Let's assume we have a matrix $A=(a_{i,j})$, with dimensions say $3\times5$. The transpose is obtained by considering the elements with the indices swapped, so $A^T=(a_{j,i})$.

But our computers reason monodimensionally and cartesian indices, from the point of view of a computer, are just a silly abstraction humans need. In fact our matrix $A$ is represented in memory as a sequence (let's assume contiguous) of bits representing 15 numbers. In particular if we have:

$$A= \begin{bmatrix} a & b & c & d & e \cr f & g & h & i & j \cr k & l & m & n & o \end{bmatrix}$$

Then we will have an in-memory linear representation of:

$$A = [a\quad f\quad k\quad b\quad g\quad l\quad c\quad h\quad m\quad d\quad i\quad n\quad e\quad j\quad o]$$

No, I'm not having a stroke, julia is a column-major language, which means that arrays are stored column by column instead of row by row like in the case of a row-major language.

By comparison, the transpose of our matrix is denoted with $A^T$ and is a $5\times 3$ matrix:

$$A^T=\begin{bmatrix} a & f & k \cr b & g & l \cr c & h & m \cr d & i & n \cr e & j & o \cr \end{bmatrix}$$

with a linear in-memory representation:

$$A^T = [a\quad b\quad c\quad d\quad e\quad f\quad g\quad h\quad i\quad j\quad k\quad l\quad m\quad n\quad o]$$

The elements of the transposed matrix are a permutation of the elements of our original matrix. At least linearly. To obtain the actual transposed matrix we would also need to interpret this strip of memory with the correct binning of elements. Where with binning we mean grouping together which elements make up a column.

In Code:

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

3×5 Matrix{Symbol}:
 :a  :b  :c  :d  :e
 :f  :g  :h  :i  :j
 :k  :l  :m  :n  :o

@show vec(A) # linear representation

vec(A) = [:a, :f, :k, :b, :g, :l, :c, :h, :m, :d, :i, :n, :e, :j, :o]

permutedims(A) # the transpose

5×3 Matrix{Symbol}:
 :a  :f  :k
 :b  :g  :l
 :c  :h  :m
 :d  :i  :n
 :e  :j  :o

@show vec(permutedims(A)) # linear repr of the transpose

vec(permutedims(A)) = [:a, :b, :c, :d, :e, :f, :g, :h, :i, :j, :k, :l, :m, :n, :o]

# reinterpreted linear rapresentation into 3 "buckets" (columns) of 5 elements (rows)
reshape(vec(permutedims(A)), (5,3))

5×3 Matrix{Symbol}:
 :a  :f  :k
 :b  :g  :l
 :c  :h  :m
 :d  :i  :n
 :e  :j  :o

Permutations

As a small refresher, a permutation is a function $σ: S \rightarrow S$ from a set $S$ to the same set $S$ such that every element is mapped to another element without overlaps or leftover elements. In technical terms, a bijection. A permutation $\sigma$ is equivalent to a rearrangement of the elements of $S$ such that the $i$-th element is replaced by the corresponding element given by $\sigma(i)$.

A usual notation to represent permutations is the cyclic notation:

$$\sigma =(g\quad \sigma(g)\quad \sigma(\sigma(g))\quad \ldots)$$

But we will implicitly use also the "one line" notation:

$$\sigma = \sigma(1) \sigma(2) \sigma(3) \,\ldots\, \sigma(N)$$

.

In particular, the "one line" notation can be more practical since a permutation can then be represented simply as an array or tuple $T$ where

$$T[i] = \sigma(i)$$

.

Another straight-forward observation is that we don't really care about the elements themselves for this problem, but simply about their position in memory.

So, in our simple case of a transposition we would like to have a permutation that given the position in memory of the $(i,j)$-th element, it will give us the position in memory of the $(j, i)$-th element.

To cut this presentation short, we'll skip ahead and tackle the more general presentation of the problem. But I encourage the curious reader to stop here and take some time to try and figure out what would be the formula of such permutation.

Multidimensional Arrays and Arbitrary Permutations

Let's go ahead and define the more generic version of the problem:

Let $A$ be a $d$-dimensional array with dimensions $D = (N_i)_{i\in[d]}$ and column-major order, and $\tau$ a permutation of the dimensions. We indicate with $A^\tau$ the $d$-dimensional array with dimensions given by $\tau(D) = (N_{\tau(i)})_{i\in[d]}$.

Let $L=\prod N_i$ be the number of elements in the array. We want to find a permutation $\sigma$ such that $\sigma([L])$ has the same arrangement of elements of $A^{τ}$. With the notation $[X]$ we indicate the set of number from $1$ to $X$ included, so $[X] = \{1,2,3,\ldots,X\}$.

We begin by noticing that the linear ordering of the elements of $A$ is the co-lexicographical order (first index varies fastsest) of the cartesian indices. That is, with $i_k \in \{1,\ldots,N_k\}$, the sequence runs

$$\begin{aligned} &(1, 1, \ldots)(2, 1, \ldots)\ldots(N_1, 1, \ldots)(1, 2,\ldots)\ldots\\ &(N_1,N_2,1,\ldots)(1,1,2,\ldots)\ldots(N_1,N_2,\ldots,N_d)\\ \end{aligned}$$

Given a linear index $l$ and a shape $D=(N_1,\ldots,N_d)$ we want to recover the corresponding cartesian index $I=(i_1,\ldots,i_d)$. We define the function

$$M_D^{-1}:\left[\prod N_i\right] \longrightarrow \prod [N_i]$$ $$\qquad l \mapsto M_D^{-1}(l) = \bigl(m_{D,i}^{-1}(l)\bigr)_{i\in[d]}$$

For example, it will map $1$ to $(1,1,1,...)$ and $2$ to $(2,1,1,...)$. The reason why we define it as the inverse of another function is mainly because when I was exploring this problem I picked the wrong convention and now I am too lazy to change everything. Sorry for the confusion.

The key quantity is the stride of each dimension: $S_{D,k}$. The $k$-th stride is the distance between two subsequent elements along the $k$ dimension, or alternatively is the number of linear positions gained by incrementing $i_k$ by one. In column-major order:

$$S_{D, k} \;=\; \prod_{i=1}^{k-1}N_i\quad\forall k\in[d] \qquad S_{D,1} = 1 \text{ (empty product).}$$

The actual function is given by

$$m^{-1}_{D,k}(l) = 1 + \left\lfloor\frac{l-1}{S_{D,k}}\right\rfloor \quad \mathrm{mod}_1 \quad N_k$$

introducing the two 1-base helpers:

$$\mathrm{fld}_1(a,b) \;=\; 1 + \left\lfloor \tfrac{a-1}{b}\right\rfloor, \qquad (x \;\mathrm{mod}_1\; n) \in \{1,\ldots,n\}.$$

where the $\mathrm{mod}_1$ is the modulo that maps into the range $(0, N_k]$, which is more suited for base 1 computations. Comparing to the usual $mod$ that maps into the range $[0, N_k)$. i.e.

$$K\quad \mathrm{mod}_1\quad K = K$$ $$K\quad \mathrm{mod}\quad K = 0$$

and $\mathrm{fld}_1$ is the 1-based integer floored division.

the function becomes:

$$m^{-1}_{D,k}(l) \;=\; \mathrm{fld}_1(l,\, S_{D,k}) \;\;\mathrm{mod}_1\;\; N_k,$$

in other words: "divide by the stride, then wrap around the dimension size."

The empty product for the first stride is a bit of an annoyance, but one way around it is by introducing a fake dimension $N_0 = 1$ which will allow us to start the index $i$ from 0 instead of 1 without changing the result but allowing for a more generic and branchless implementation.

Ok, let's break it up a little bit and write some code to help sediment the ideas. Below is a simple implementation of our $M_{D}^{-1}$ function.

function Minv(dims::NTuple{N, Int}) where {N}
	# pad with a leading 1 (N_0) so that the prod(dims[1:k-1])
	# falls out of prod(dims0[1:k]) without branching on k = 1.

	# or prod(dims[0:k-2]) and prod(dims0[0:k-1])
	# if you're foaming by the mouth already
	dims0 = (1, dims...)

	# we return the function m_D^{-1} (lowercase m)
	l -> ntuple(N) do k
		Sₖ = prod(dims0[1:k]) 		# our k-th stride (1-based)
		mod1(fld1(l, Sₖ), dims0[k+1])	# +1 offsets the 0th element
	end
end

We've laid down all the boring bits, let's continue. The inverse operation (pick your conventions correctly people) consists of, given a cartesian index, obtaining the linear index given a tuple of dimensions. In particular we can define:

$$M_D:\prod [N_i]\longrightarrow\left[\prod N_i\right]$$ $$\qquad(i_1,i_2,\ldots) \to m_D(i_1, i_2, \ldots)$$

with

$$m_D(i_1, i_2, \ldots) = 1 + \sum_{k=1}^{d}(i_k - 1)S_{D,k}$$

If you are particularly bored you can check that indeed $(M_D \circ M_D^{-1}) = id$.

Again, lets write some code to help concretize

function M(dims::NTuple{N}) where {N}
	dims0 = (1, dims...)
	# returns a function that takes a tuple
	# of cartesian indices and spits out a linear index
	(I) -> begin
		1 + sum(1:N) do k
			Sₖ = prod(dims0[1:k])
			(I[k]-1)*Sₖ
		end
	end
end

Let's check if the identity holds:

# let's check if the identity holds!
D = size(A)
m, minv = M(D), Minv(D)

@assert all(l -> (m ∘ minv)(l) == l, 1:length(A))

Hurray! (doesn't print anything, but the page wouldn't compile if it didn't work)

Now, let's recall our final objective: Given our permuted array $A^\tau$ by the permutation $\tau$, we want to know for each of its linear indices, what is the corresponding linear index in the non permuted array $A$.

To approach this problem the general strategy will be:

• given a linear index $l'$, obtain the cartesian index associated to the permuted dimensions,

• apply the inverse permutation to that cartesian index

• map back to the linear index.

In a more schematic way we have:

In particular, we're interested in the inverse permutation $\sigma^{-1}$, which according to the diagram is defined by:

$$\sigma^{-1} = (M_D \circ \tau^{-1} \circ M_{\tau(D)}^{-1})$$

We have all the intermediate pieces to build our permutation, In code this would look like this:

#Small helper function
applyperm(x::NTuple{N}, perm::NTuple{N, Int}) where {N} =
    ntuple(i -> x[perm[i]], N)

function σinv_naive(D::NTuple{N, Int}, τ::NTuple{N, Int}) where {N}
    permD = applyperm(D, τ)
    invτ = invperm(τ)

    # we construct the function directly by concatenating other functions!
    return M(D) ∘ ((I) -> applyperm(I, invτ)) ∘ Minv(permD)
end

Great! We can now show the permutation in its entirety!

$$\sigma^{-1} = \begin{pmatrix} 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14 & 15 \\ 1 & 4 & 7 & 10 & 13 & 2 & 5 & 8 & 11 & 14 & 3 & 6 & 9 & 12 & 15 \end{pmatrix}$$

Let's see how long it takes to compute the whole permutation:

# First benchmark - naive approach
f(A) = begin
    s = σinv_naive(size(A), (2, 1))
    @inbounds for l in eachindex(A)
        Core.donotdelete(s(l)) # prevent the compiler from being too smart
    end
end
show(stdout, "text/plain", @benchmark f(A))

BenchmarkTools.Trial: 10000 samples with 9 evaluations per sample.
 Range (min … max):  2.926 μs …  14.324 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.991 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   3.039 μs ± 229.390 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▂▆████▆▅▅▃▃▃▂▁▂▁            ▁ ▂▂▃▃▂▁▁                      ▂
  ▇██████████████████▇▇▇▆▆█▇▆█▇█▇████████▇▇▇▇▇▆▆▆▃▅▆▅▄▅▅▄▄▃▅▅ █
  2.93 μs      Histogram: log(frequency) by time      3.56 μs <

 Memory estimate: 1.41 KiB, allocs estimate: 60.

We can definitely do better! Function composition is cute, but when done like this, it tends to be heavy.

Expanding the permutation

If you've had enough math for now, look away, skip to the next section for more code. Otherwise carry on.

The composition for $\sigma^{-1}$ is simple enough that it can be expanded analytically. Given a position $l'$ we have:

$$\sigma^{-1}(l') = (m_D ∘ τ^{-1})(m_{\tau(D)}^{-1}(l')) =\,\ldots$$

with

$$m_{\tau(D)}^{-1}(l') = \left( \mathrm{fld}_1(l', S_{\tau(D),k})\quad \mathrm{mod}_1\quad N_{\tau(k)}\right)_{k\in[d]}$$

let's call this tuple $(a_k)_{k∈[d]}$ for ease of notation.

The denominators of the floored divisions are the strides of the permuted dimensions. For example, if an array has dimensions $(3, 4, 2)$ then its strides will be: $(1, 3, 12)$ (remember we start at $0$ and go up until $k-1$). If we were to shuffle the order of the dimensions, say $(2, 4, 3)$, then the strides will be different: $(1, 2, 8)$

The first step gave us a tuple of numbers $(a_k)_{k∈[d]}$, which we now want to reshuffle back using the inverse of our permutation $\tau$:

$$...\,= m_D\left( \tau^{-1}( (a_k)_{k\in[d]})\right) = m_D((a_{\tau^{-1}(k)})_{k\in[d]}) =$$ $$=1+\sum_{k=1}^{d} (a_{\tau^{-1}(k)}-1)S_{D,k}=$$

and if we expand back the $a_k$'s we get:

$$=1 + \sum_{k=1}^{d} \left( -1 + \mathrm{fld}_1\left(l', S_{\tau(D),\tau^{-1}(k)}\right)\quad \mathrm{mod}_1\quad N_{\tau^{-1}(\tau(k))}\right)S_{D,k}$$

Let's now implement this expanded version of $\sigma$. A keen eye will notice that the strides depend only on the dimensions and not on the actual element $l'$, which means that we can precompute them.

function σinv(dims::NTuple{N}, perm::NTuple{N}) where {N}
    permdims = applyperm(dims, perm)
    iperm = invperm(perm)
    # we can precompute all the strides
    strides = accumulate(*, (1, dims...))
    permStrides = accumulate(*, (1, permdims...))

    # 1 + Σ_{k=1}^{N} ( fld1( gen,  Π_{i=1}^{τ^-1(k)-1} Nτ(i)) mod1 Nτ(k) ) - 1)*Π_{i=1}^{k-1} Ni
    return it -> begin
        1 + sum(Base.OneTo(N)) do k
            aₖ = mod1(fld1(it, permStrides[iperm[k]]), dims[k])
            (aₖ - 1) * strides[k]
        end
    end
end

# Second benchmark - optimized approach
f2(A) = begin
	s = σinv(size(A), (2, 1))
	@inbounds for l in eachindex(A)
		Core.donotdelete(s(l)) # prevents the compiler from being too smart
	end
end
show(stdout, "text/plain", @benchmark f2(A))

BenchmarkTools.Trial: 10000 samples with 965 evaluations per sample.
 Range (min … max):  83.031 ns … 197.409 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     83.204 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   84.243 ns ±   4.326 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▄     ▁▃                                                    ▁
  ████▇▇████▆▆▆▇▇▅▆▆▄▄▅▄▃▄▄▃▅▅▄▄▄▄▄▄▄▂▅▄▅▄▅▄▆▅▄▄▃▄▅▅▄▄▅▃▄▄▄▄▄▄ █
  83 ns         Histogram: log(frequency) by time       105 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

Much better! Almost 40 times faster, but more importantly: no allocations. That means, no GC involved.

Swapping the elements?

Ok great, now what? We know where elements should go, but we need to actually put them there.

It's easy! Just swap the elements according to our $\sigma$!

# define our sigma and its inverse
siginv = σinv(size(A), (2, 1))
iS = siginv.(eachindex(A))
S = invperm(iS)

B = deepcopy(A)
for l in eachindex(B)
    B[siginv(l)] = B[l]
end
# we swapped the elements, but we need to tell the program
# that now this array is a 5x3 and no longer a 3x5!
Bt = reshape(B, (5, 3))

5×3 Matrix{Symbol}:
 :a  :n  :m
 :l  :k  :f
 :m  :h  :g
 :f  :g  :f
 :k  :f  :o

Err... No. That doesn't look right. Repeated simbols? What's going on?

The issue comes from the fact that we're operating inplace. If you had another copy of your array, then, this approach would work. But since we don't, what happens is that a swap might put the wrong element in a slot that will need to be swapped later.

Let's recall the full definition of the memory permutation for our running example, as well as the inverse:

$$\sigma^{-1} = \begin{pmatrix} 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14 & 15 \\ 1 & 4 & 7 & 10 & 13 & 2 & 5 & 8 & 11 & 14 & 3 & 6 & 9 & 12 & 15 \end{pmatrix}$$

$$\sigma = \begin{pmatrix} 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14 & 15 \\ 1 & 6 & 11 & 2 & 7 & 12 & 3 & 8 & 13 & 4 & 9 & 14 & 5 & 10 & 15 \end{pmatrix}$$

Our $\sigma^{-1}$ is telling us that at position $l$ in memory of our transposte, we need to put the element at position $\sigma^{-1}(l)$ of the current array.

So for example with l = 3

$$A^T[3] \leftarrow A[7]$$

Then, once we reach

$$\sigma(3) = 11$$

we will need to do the following swap:

$$A^T[11] \leftarrow A[3]$$

But $A^t$ and $A$ share the same memory, so after the first swap

$$A[3] = A^T[3] = A[7]$$

Leading to:

$$A^T[11] \leftarrow A[7]$$

instead! Which is wrong! Since we've been going in order and we swapped already

$$A^T[7] \leftarrow A[5]$$

We see that what's happening is that the order of our swaps can't be the naive from 1 to N but should be given instead by the permutation itself:

$$A[l] \leftarrow A[\sigma^{-1}(l)] \leftarrow A[\sigma^{-1}(\sigma^{-1}(l))] \leftarrow \ldots A[(\sigma^{-1})^k(l)]$$

Swapping the elements: The right way

A permutation can be decomposed into disjoint cycles. Meaning that each cycle will only touch a subset of all the elements and partition the whole interval.

In particular, these cycles are the orbits of the cyclic group $\langle\sigma\rangle = \{1, \sigma, \sigma^2, \ldots\}$ that acts on our set of "positions" $[L]$. In other words we can always write $\sigma$ as:

$$\sigma = (g_1,\,\sigma(g_1),\, \sigma(\sigma(g_1)),\,\ldots)\,(g_2,\, \sigma(g_2),\, \sigma(\sigma(g_2)),\,\ldots)\,\ldots$$

for some $g_i$'s that act as the generators of the cycles.

Computing these cycles for our permutation can be done in a brute-forcy way with a technique called cycle chasing:

1. Pick a random element, for example $2$ (this is actually the optimal choice^[4]),

2. apply the permutation until you obtain $2$ again. Marking all elements you visit along the way.

3. Then pick the next available element not yet visited, go to 2. If no other elements are available end.

^[4] It is proved (Cate & Twigg, 1977) that the cycle generated by the second element is always the cycle of maximum length $G$, and any other cycle must have length that divides $G$.

In code this looks like this:

function σ_cycles(dims::NTuple{N, Int}, perm::NTuple{N, Int}) where {N}
    available = trues(prod(dims)) # a bool vector with prod(dims) trues

    g = 2

    # first and last element are always fixed (source: trust me)
    available[1] = available[end] = false
    cycles = Vector{Int}[[1], [prod(dims)]]

    f = σinv(dims, perm)
    while any(available)
        available[g] = false
        cycle = Int[g]
        next = f(g)

        while next != g
            available[next] = false
            push!(cycle, next)
            next = f(next)
        end

        g = findfirst(available)
        push!(cycles, cycle)
    end

    return cycles
end

In our case, we have that our permutation can be split into the following disjoint cycles:

$$\sigma^{-1} = (2\hspace{0.5em}4\hspace{0.5em}10\hspace{0.5em}14\hspace{0.5em}12\hspace{0.5em}6)(3\hspace{0.5em}7\hspace{0.5em}5\hspace{0.5em}13\hspace{0.5em}9\hspace{0.5em}11)$$

(we're almost done, I swear! Please bear with me!)

We need one extra bit, that is the fact that any cycle can be represented as a sequence of transpositions! A transposition is simply a cycle of length $2$. i.e.: $(1\,\,\,\,\,2)$ is a transposition.

In particular, we have:

$$(a_1\,\,\,\,\,a_2\,\,\,\,\,a_3\,\,\,\,\,\ldots\, a_n) = (a_1\,\,a_n)(a_1\,\,a_{n-1})\,\ldots\,(a_1\,\,a_3)(a_1\,\,a_2)$$

This is the last ingredient!

We generate the cycles, and then split them up into a sequence of transpositions (swaps) to get the correct behaviour!

function overcomplicated_permutedims!(array, perm)
    cycles = σ_cycles(size(array), perm)
    for c in cycles
        length(c) == 1 && continue

        g, rest... = c
        for i in reverse(rest)
            array[g], array[i] = array[i], array[g]
        end
    end
    return reshape(array, applyperm(size(array), perm))
end

let
    B1 = deepcopy(A)
    B2 = deepcopy(A)

    B_t = overcomplicated_permutedims!(B1, (2, 1))
    B_t_true = permutedims(B2) # compare with the sane way

    @assert all(B_t .== B_t_true)
end

Huzzah!

We successfully swapped all elements to obtain the correct transpose. And without using extra arrays. Or did we? cue V-Sauce music

The astute among you might have noticed something fishy is going on. I promised $O(1)$ space, but I am making use of multiple auxiliary arrays to generate and store the cycles. So yeah, I lied about the last ingredient being, well..., the last.

The last ingredient (really this time)

To overcome this last hurdle, we'll take advantage of the great julia's metaprogramming capabilities. In particular, we'll make use of @generated functions.

Instead of badly parroting a half-assed definition I'll just quote the opening paragraph from the Julia docs:

A very special macro is @generated, which allows you to define so-called generated functions. These have the capability to generate specialized code depending on the types of their arguments with more flexibility and/or less code than what can be achieved with multiple dispatch. While macros work with expressions at parse time and cannot access the types of their inputs, a generated function gets expanded at a time when the types of the arguments are known, but the function is not yet compiled.

Together with the use of the Val type, we can hoist the dimensions and permutation into the typesystem and use that information to emit the body of a function with all the swaps hardcoded.

# here we define it without the @generated to show the generated AST.
function _truly_overcomplicated_permutedims!(
    ::Type{<:AbstractArray{T, N}}, dims::NTuple{N, Int}, perm::NTuple{N, Int}
) where {T, N}

    isperm(perm) || throw(ArgumentError("input is not a permutation"))

    C = σ_cycles(dims, perm)

    # in generated functions we don't ouput a value,
    # we output the AST of a function that will be compiled and run.
    exs = Expr[]
    for c in C
        length(c) == 1 && continue # skip fixed points

        g, rest... = c
        for i in reverse(rest)
            # interpolate the indices directly in the expressions
            ex = :((array[$g], array[$i]) = (array[$i], array[$g]))
            push!(exs, ex)
        end
    end

    push!(exs, :(reshape(array, $(applyperm(dims, perm)))))
    return Expr(:block, exs...)
end


# the arguments of a generated function are types, not values.
# we pass them along, generate the AST which will be compiled and run.
@inline @generated truly_overcomplicated_permutedims!(
    array::AbstractArray, ::Val{dims}, ::Val{perm}
) where {dims, perm} = _truly_overcomplicated_permutedims!(array, dims, perm)

# some helper functions to help hoist information in the type system.
truly_overcomplicated_permutedims!(arr::AbstractArray, perm) =
    truly_overcomplicated_permutedims!(arr, Val(size(arr)), Val(perm))

truly_overcomplicated_permutedims!(arr::AbstractArray, perm::Val) =
    truly_overcomplicated_permutedims!(arr, Val(size(arr)), perm)

The AST generated looks something like this:

begin
    (array[2], array[6]) = (array[6], array[2])
    (array[2], array[12]) = (array[12], array[2])
    (array[2], array[14]) = (array[14], array[2])
    (array[2], array[10]) = (array[10], array[2])
    (array[2], array[4]) = (array[4], array[2])
    (array[3], array[11]) = (array[11], array[3])
    (array[3], array[9]) = (array[9], array[3])
    (array[3], array[13]) = (array[13], array[3])
    (array[3], array[5]) = (array[5], array[3])
    (array[3], array[7]) = (array[7], array[3])
    reshape(array, (5, 3))
end

@assert truly_overcomplicated_permutedims!(deepcopy(A), (2,1)) == permutedims(A)

Ole'!

BenchmarkTools.Trial: 10000 samples with 383 evaluations per sample.
 Range (min … max):  248.475 ns …  48.187 μs  ┊ GC (min … max): 0.00% … 99.25%
 Time  (median):     258.376 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   268.159 ns ± 500.443 ns  ┊ GC (mean ± σ):  2.31% ±  1.39%

   ▃█▄      ▄▅                                                   
  ▃███▆▅▅▄▆███▇▅▄▅▅▇▆▆▅▄▄▄▃▃▃▃▄▅▃▂▂▂▃▃▃▃▂▂▂▂▂▂▂▂▂▁▂▁▂▁▁▁▁▁▁▁▁▁▁ ▃
  248 ns           Histogram: frequency by time          297 ns <

 Memory estimate: 112 bytes, allocs estimate: 3.

Where are those allocations coming from? Well, we're playing a bit with fire here, the compiler is good, but it's not magic. Those allocations come from us trying to hoist a runtime known value into the type system as if it were known at compile time. To really get rid of every last bit of allocations we need to work with quantities that are actually known at compile time.

Let's try to do that:

using StaticArrays
SA = SizedArray{Tuple{size(A)...}}(A)

const P = (2, 1)
const Pval = Val(P)

show(stdout, "text/plain", @benchmark truly_overcomplicated_permutedims!($SA, $Pval))

BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
 Range (min … max):  9.500 ns … 45.208 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     9.625 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   9.661 ns ±  0.672 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

     ▄   █   █   ▅   ▃   ▃   ▂  ▁   ▂   ▂                    ▂
  ▆▁▁█▁▁▁█▁▁▁█▁▁▁█▁▁▁█▁▁▁█▁▁▁█▁▁█▁▁▁█▁▁▁█▁▁▁▇▁▁▁▅▁▁▁▆▁▁▁▅▁▁▄ █
  9.5 ns       Histogram: log(frequency) by time     10.1 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

Compared to the simple, efficient and sane way:

BenchmarkTools.Trial: 10000 samples with 999 evaluations per sample.
 Range (min … max):  10.468 ns … 36.203 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     10.553 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   10.661 ns ±  1.110 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

            █     ▂                                            
  ▂▁▁▁▁▃▁▁▁▁█▁▁▁▁▁█▁▁▁▁▃▁▁▁▁▃▁▁▁▁▁▃▁▁▁▁▂▁▁▁▁▂▁▁▁▁▁▂▁▁▁▁▂▁▁▁▁▂ ▂
  10.5 ns         Histogram: frequency by time        10.9 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

And there you go, a completely in-place transpose without extra allocations. And apparently for small enough arrays it's on par with the inplace one with an extra destination array.

The Catch

Admit it, you were waiting for this one.

Shh... Listen. Do you hear those muffled sounds? It's the compiler choking to death after you tried to run this piece of code with a 50-dimensional array and 10000 elements. I hope you like watching paint dry.

The catch here is that the @generated function body contains one swap (2 expressions) per cycle element. And if we consider all cycles we have $O(N)$ expressions where $N$ is the number of elements in the array. So far so good, if it weren't for the fact that the compiler roughly scales quadratically in the number of expressions [citation needed].

This means that for large arrays, if it doesn't die midway, when it finishes compiling, you'll be noticeably older. And god forbid you change dimensions or permutations. For each new pair, a new function is generated and compiled.

But not all is lost: A neat trick is to tell the generating function to emit hardcoded swaps only up until there are no more to emit or if you reach a maximum number of expressions. If any swaps are left over, emit a for loop that cycles over the remaining elements (basically a partial manual loop unrolling).

This of course violates a bit the O(1) promise, for large $N$'s, but allows for a safe fallback that doesn't completely hang the system.

@generated function _even_more_overcomplicated_permutedims!(
	array::AbstractArray{T,N},
	::Val{dims},
	::Val{perm}
) where {T, N, dims, perm}
	isperm(perm) || throw(ArgumentError("input is not a permutation"))

	pdims = ntuple(i -> dims[perm[i]], N)

	C = σ_cycles(dims, perm)

	max_exprs = 1000
	nexpr = 0

	breaknext = false
	rest_i = 1

	exs = Expr[]
	for (i,cycle) in enumerate(C)
		length(cycle) == 1 && continue

		leader = cycle[1]
		for j in length(cycle):-1:2

			push!(exs, quote
				@inbounds array[$leader], array[$(cycle[j])] = array[$(cycle[j])], array[$leader]
			end)
			nexpr += 1
			if nexpr >= max_exprs
				breaknext=true
			end
		end

		rest_i = i+1
		breaknext && break
	end

	if rest_i <= length(C)
		push!(exs, quote
			cycles = $(C[rest_i:end])
			@inbounds for cycle in cycles
				leader = first(cycle)
				for i in length(cycle):-1:2
					arr[leader], 	arr[cycle[i]] = arr[cycle[i]], arr[leader]
				end
			end
		end)
	end

	push!(exs, :(reshape(array, $pdims))
	Expr(:block, exs...)
end

Other permutations and higher dimensions?

So far, we've only played with a toy example of a 2 dimensional array and very precise dimensions. But this approach works for any dimensions and sizes. In principle.

Let's quickly test that all this also works for other dimensions and permutations! Some kind of fuzzy testing at home.

# multi-dim array benchmark example

begin
    while true
        global ndims = rand(2:5)
        global dims = ntuple(i -> rand(1:10), ndims)
        global arr = rand(dims...)
        global arrstatic = SizedArray{Tuple{dims...}}(arr)
        global perm = Tuple(Random.shuffle(1:ndims))

        # So I can publish within today
        prod(dims) <= 2000 && break
    end

    println(
        """Case:
        	ndims = $ndims
        	dims = $dims
        """
    )

    @assert truly_overcomplicated_permutedims!(deepcopy(arr), Val(perm)) ==
        permutedims(arr, perm)


    arrcopy = similar(arr, applyperm(dims, perm)...)

    perm_inplace = @benchmark permutedims!($arrcopy, $arr, $perm)

    complicated_perm_inplace = (@benchmark
        truly_overcomplicated_permutedims!(arr, $(Val(perm)))
        setup = (arr = deepcopy($arrstatic))
    )

    show(stdout, "text/plain", perm_inplace)
    println()
    show(stdout, "text/plain", complicated_perm_inplace)

end

Case:
	ndims = 4
	dims = (6, 9, 2, 6)

BenchmarkTools.Trial: 10000 samples with 303 evaluations per sample.
 Range (min … max):  277.363 ns … 376.238 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     277.776 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   279.204 ns ±   6.089 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▅▃             ▂▃                                            ▁
  ██████▇▇▆▄▆▆▇▇▇▇██▇▄▅▆▆▅▅▅▆▆▆▆▅▁▃▁▆▆▅▅▄▃▁▄▄▃▄▃▁▄▃▁▄▄▃▁▄▄▄▄▃▅▅ █
  277 ns        Histogram: log(frequency) by time        312 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.
BenchmarkTools.Trial: 10000 samples with 839 evaluations per sample.
 Range (min … max):  146.454 ns … 272.895 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     147.298 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   147.998 ns ±   3.277 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▇▇██▅▁     ▂▃▃▄▂                                              ▂
  ██████▅▇▆▇▇█████▇▅▆▅▆▇▆▆▇▇▇▆▆▄▄▅▅▄▄▅▃▃▃▄▅▄▅▅▅▅▄▄▄▂▂▄▂▅▄▅▃▅▅▅▅ █
  146 ns        Histogram: log(frequency) by time        165 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.