矩阵转置

前几天同事问我一个问题：矩阵的在位转置或者是原地转置(transpose inplace)，意思就是尽量不分配额外的内存空间，把一个矩阵转置。
矩阵的存储结构是用一块线性内存来存储数据，然后记录宽度、高度，大致如下

class Matrix { double* data; int width, height; };

【问题分析】
因为没有额外空间可用，所以改变数据只能通过“交换”，也就是把矩阵中的两个元素交换位置，这个操作可以不使用任何额外存储空间。那么通过交换操作，可以衍生出来的更复杂的操作就是循环移位，大致意思就是一组数据向左或者向右移动一个位置，比如原来是
{ x a1 a2 a3 ... an y }，那么左移位后就是 { a1 a2 ... an y x }，显然通过连续的数据交换是可以实现这个操作的。那么我们的算法中可以使用的数据移动操作就只有这两种：swap(i, j)完成数据交换，move_data(src, dst)完成{src, dst}之间的数据循环移位。（要说明的是，显然数据循环移位操作，不仅仅局限于连续的数据）

矩阵转置，实质上是对存储索引的重新排列(permutation)，和排序有点类似，但不一样的是，排序是值排序，算法决定一个数据元素的最终位置的时候只要知道这个数据的数值就可以了，但是转置是对索引位置的排序，我们不关心数据值的大小，只关心它在源序列中的位置，而如果不用额外存储空间，这个信息是很难保存的。

比如对于一个3*4的矩阵

1  2  3  4
5  6  7  8
9 10 11 12转置后为：

1 5  9
2 6 10
3 7 11
4 8 12观察在内存中的存储形式，源矩阵是
S = { 1 2 3 4 5 6 7 8 9 10 11 12 }

目标矩阵是
D = { 1 5 9 2 6 10 3 7 11 4 8 12 }

转置问题就是对于S，求出D即可。S是先验可知的，只是和矩阵的宽度高度有关，D是S的一个置换。

【正推】

比较容易想到的思路就是对S中的每个元素S[i]，求出它在D中的位置j1，我们把S[j1]先保存一下，然后再看S[j1]应该在D中的位置，比如是j2，那么再看S[j2]应该被放到的位置j3，这样继续下去，就得到了一个子序列 { i, j1, j2, ... i }，最后一个一定是i，也就是说这个子序列最后必然要回到i，很容易想到啊，最后总有一个数据要填充到第一个数据的位置上，因为就这么些个数据，不可能出来一个洞洞啊，呵呵。我们把这个子序列称为T(i)

那么对这样一个子序列完成一个move_data操作之后，这些位置上的数据就安放好了，然后我们再寻找下一个子序列，从i+1开始，循环下去，直到所有的数据都放好了，就收工了。

这个算法的一个问题就是当我完成一个子序列T(i)后，如果i+1包含在T(i)中，那么T(i+1)和T(i)是一样的，必须要跳过i+1，看i+2，如果i+2也在T(i)中，那就继续跳过，所以我需要一个标志数组来标志整个序列中那些位置的数据是已经放好的了，那些数据是还没有放好的，我只要对那些没有放好的数据做move_data。那么这个额外的标志数组的大小和矩阵是一样大的，呵呵额外存储空间啊，不行，所以这个方法失败。

不用标志数组也可以，我可以事先先把整个序列划分成若干个互不重叠的子序列，然后再进行操作，但是划分结果保存在那里呢？呵呵，所以，这个方法还是不行。

【反推】

换个思路，也许就好了。

如果我们从目标矩阵开始呢，我们考虑D中的一个位置i，这个i应该由S中的S[j]来填充，那么我们直接把S[i..j]之间做一个循环移位，这样i位置就被正确的设置了，同时[i+1,j]之间的数据仍然是按原来在S中的顺序。如果我们考虑把i从第一个位置开始依次循环递加，那么每次操作之后，我都可以保证[1, i]之间是放好的，[i, n]之间是没有放好的，呵呵看出来好处了么？我们不再需要任何额外的存储空间来保存哪些数据是已经放好的了。

虽然原理是很容易，但是写出程序还是要费一点心思，关键是在于对于[i, n]之间的这些待放置数据，它们之间的相对顺序虽然没有乱，但是位置都变动了，所以你要能够一直跟踪这些位置变动信息。这个比看上去的要复杂一点。

最后的程序如下

void TransposeInplace() { int from, to, offset, step; int i, j; for (i = to = offset = 0, step = width; i < width; ++i, --step, offset += height-1) { for (j = 0, from = i + offset; j < height; ++j, ++to, from += step) { _move_data(from, to); } } swap(width, height); }

显然，如果以move_data为单元操作，这个算法的复杂度是O(n)的，考虑到move_data本身也是O(n)的，那么这个算法的复杂度是O(n^2)

【再快一点】

能不能更快一些呢？呵呵

算法上的加速不是很容易了，但是技巧上的加速还是可以的。

move_data是O(n)的，但是swap就是O(1)的了，上面的算法中对于每一个数据都要通过move_data来进行放置，所以我想到的一个加速就是尽可能的用swap来取代move_data

假设m = min(width, height)，矩阵左上角的，大小为m×m 的方阵显然是可以通过交换来完成转置的，如果我们先用swap对这个小的子方阵做好转置，再对剩下的数据放置，也许就可以快一些了。

比如，对前面用过的那个3*4的矩阵

1  2  3  4
5  6  7  8
9 10 11 12首先变化成

1 5  9  4
2 6 10  8
3 7 11 12也就是左边的3*3的小方阵转置一下，右边的部分不动，这样变换后，子序列{ 1, 5, 9, 2, 6, 10, 3, 7, 11 }有序，相对位置不需要变动了，剩下的工作只是把子序列 {4, 8, 12} 移到矩阵末尾去。呵呵，这个要快很多吧。

再仔细推敲这个算法，可以写出程序来，这个程序的编写仍然是以跟踪数据索引的位置变动为难点。比第一个要稍微困难一点，但是也不致于写不出来。

程序如下

void TransposeInplaceQuick() { int m = width < height ? width : height; int i, j; for (i = 0; i < m; ++i) { for (j = i+1; j < m; ++j) { int i_old = _index(j, i); int i_new = _index(i, j); swap(data[i_old], data[i_new]); } } if (m == width) { for (i = m; i < height; ++i) { for (j = 0; j < width; ++j) { _move_data(_index(j, i), j * (width + i - m + 1) + i); } } } else { int last = _index(width-1, height-1); int n = 0; for (j = width - 1; j >= m; --j) { int step = n; for (i = height - 1; i >= 0; --i, --last) { _move_data(_index(j, i) - step, last); step -= (width - 1 - j); } n += height - 1; } } swap(width, height); }

这个算法的时间复杂度仍然是O(n^2)的，但是要比第一个快一些。平均看来时间应该是前一个算法的一半。