In-place matrix transposition

On a computer, one can often avoid explicitly transposing a matrix in memory by simply accessing the same data in a different order. For example, software libraries for linear algebra, such as BLAS, typically provide options to specify that certain matrices are to be interpreted in transposed order to avoid data movement.

However, there remain a number of circumstances in which it is necessary or desirable to physically reorder a matrix in memory to its transposed ordering. For example, with a matrix stored in row-major order, the rows of the matrix are contiguous in memory and the columns are discontiguous. If repeated operations need to be performed on the columns, for example in a fast Fourier transform algorithm (e.g. Frigo & Johnson, 2005), transposing the matrix in memory (to make the columns contiguous) may improve performance by increasing memory locality. Since these situations normally coincide with the case of very large matrices (which exceed the cache size), performing the transposition in-place with minimal additional storage becomes desirable.

Also, as a purely mathematical problem, in-place transposition involves a number of interesting number theory puzzles that have been worked out over the course of several decades.

For example, consider the 2×4 matrix:

${\displaystyle {\begin{bmatrix}11&12&13&14\\21&22&23&24\end{bmatrix}}.}$

In row-major format, this would be stored in computer memory as the sequence (11, 12, 13, 14, 21, 22, 23, 24), i.e. the two rows stored consecutively. If we transpose this, we obtain the 4×2 matrix:

${\displaystyle {\begin{bmatrix}11&21\\12&22\\13&23\\14&24\end{bmatrix}}}$

which is stored in computer memory as the sequence (11, 21, 12, 22, 13, 23, 14, 24).

If we number the storage locations 0 to 7, from left to right, then this permutation consists of four cycles:

(0), (1 2 4), (3 6 5), (7)

That is, the value in position 0 goes to position 0 (a cycle of length 1, no data motion). Next, the value in position 1 (in the original storage: 11, 12, 13, 14, 21, 22, 23, 24) goes to position 2 (in the transposed storage 11, 21, 12, 22, 13, 23, 14, 24), while the value in position 2 (11, 12, 13, 14, 21, 22, 23, 24) goes to position 4 (11, 21, 12, 22, 13, 23, 14, 24), and position 4 (11, 12, 13, 14, 21, 22, 23, 24) goes back to position 1 (11, 21, 12, 22, 13, 23, 14, 24). Similarly for the values in position 7 and positions (3 6 5).

In the following, we assume that the N×M matrix is stored in row-major order with zero-based indices. This means that the (n,m) element, for n = 0,...,N−1 and m = 0,...,M−1, is stored at an address a = Mn + m (plus some offset in memory, which we ignore). In the transposed M×N matrix, the corresponding (m,n) element is stored at the address a' = Nm + n, again in row-major order. We define the transposition permutation to be the function a' = P(a) such that:

${\displaystyle Nm+n=P(Mn+m)\,}$ for all ${\displaystyle (n,m)\in [0,N-1]\times [0,M-1]\,.}$

This defines a permutation on the numbers ${\displaystyle a=0,\ldots ,MN-1}$.

It turns out that one can define simple formulas for P and its inverse (Cate & Twigg, 1977). First:

${\displaystyle P(a)={\begin{cases}MN-1&{\text{if }}a=MN-1,\\Na{\bmod {(}}MN-1)&{\text{otherwise}},\end{cases}}}$

where "mod" is the modulo operation.

Second, the inverse permutation is given by:

${\displaystyle P^{-1}(a')={\begin{cases}MN-1&{\text{if }}a'=MN-1,\\Ma'{\bmod {(}}MN-1)&{\text{otherwise}}.\end{cases}}}$

(This is just a consequence of the fact that the inverse of an N×M transpose is an M×N transpose, although it is also easy to show explicitly that P⁻¹ composed with P gives the identity.)

As proved by Cate & Twigg (1977), the number of fixed points (cycles of length 1) of the permutation is precisely 1 + gcd(N−1,M−1), where gcd is the greatest common divisor. For example, with N = M the number of fixed points is simply N (the diagonal of the matrix). If N − 1 and M − 1 are coprime, on the other hand, the only two fixed points are the upper-left and lower-right corners of the matrix.

The number of cycles of any length k>1 is given by (Cate & Twigg, 1977):

${\displaystyle {\frac {1}{k}}\sum _{d|k}\mu (k/d)\gcd(N^{d}-1,MN-1),}$

where μ is the Möbius function and the sum is over the divisors d of k.

Furthermore, the cycle containing a=1 (i.e. the second element of the first row of the matrix) is always a cycle of maximum length L, and the lengths k of all other cycles must be divisors of L (Cate & Twigg, 1977).

For a given cycle C, every element ${\displaystyle x\in C}$ has the same greatest common divisor ${\displaystyle d=\gcd(x,MN-1)}$.

This theorem is useful in searching for cycles of the permutation, since an efficient search can look only at multiples of divisors of MN−1 (Brenner, 1973).

Laflin & Brebner (1970) pointed out that the cycles often come in pairs, which is exploited by several algorithms that permute pairs of cycles at a time. In particular, let s be the smallest element of some cycle C of length k. It follows that MN−1−s is also an element of a cycle of length k (possibly the same cycle).

The following briefly summarizes the published algorithms to perform in-place matrix transposition. Source code implementing some of these algorithms can be found in the references, below.

for n = 0 to N - 1
    for m = n + 1 to N
        swap A(n,m) with A(m,n)

for each length>1 cycle C of the permutation
    pick a starting address s in C
    let D = data at s
    let x = predecessor of s in the cycle
    while x ≠ s
        move data from x to successor of x
        let x = predecessor of x
    move data from D to successor of s